--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SMART/DiffExpAnal/DESeqTools/pairwiseSERE.R	Tue Apr 30 14:33:21 2013 -0400
@@ -0,0 +1,41 @@
+# pairwiseSERE
+# compute pairwise SERE statistics
+
+# input : counts
+# output : matrix of SERE values
+
+# created october 19th, 2012
+# Marie-Agnes Dillies
+
+
+pairwiseSERE <- function( counts ){
+  
+  sere <- matrix( NA, ncol=ncol(counts), nrow=ncol(counts) )
+  for (i in 1:ncol(counts)){
+    for (j in 1:ncol(counts)){
+      sere[i,j] <- sigfun_Pearson( counts[,c(i,j)] )
+    }
+  }
+  colnames(sere) <- rownames(sere) <- colnames(counts)
+  return( formatC(sere, format="f", digits=2) )
+}
+
+sigfun_Pearson <- function(observed) {
+  #calculate lambda and expected values
+  laneTotals<- colSums(observed);
+  total <- sum(laneTotals)
+  fullObserved <- observed[rowSums(observed)>0,];
+  fullLambda <- rowSums(fullObserved)/total;
+  fullLhat <- fullLambda > 0;
+  fullExpected<- outer(fullLambda, laneTotals);
+
+  #keep values
+  fullKeep <- which(fullExpected > 0);
+  
+  #calculate degrees of freedom (nrow*(ncol -1) >> number of parameters - calculated (just lamda is calculated >> thats why minus 1)
+  #calculate pearson and deviance for all values
+  oeFull <- (fullObserved[fullKeep] - fullExpected[fullKeep])^2/ fullExpected[fullKeep] # pearson chisq test
+  dfFull <- length(fullKeep) - sum(fullLhat!=0);
+  
+  return(c(sqrt(sum(oeFull)/dfFull)));
+}