|
18
|
1 # pairwiseSERE
|
|
|
2 # compute pairwise SERE statistics
|
|
|
3
|
|
|
4 # input : counts
|
|
|
5 # output : matrix of SERE values
|
|
|
6
|
|
|
7 # created october 19th, 2012
|
|
|
8 # Marie-Agnes Dillies
|
|
|
9
|
|
|
10
|
|
|
11 pairwiseSERE <- function( counts ){
|
|
|
12
|
|
|
13 sere <- matrix( NA, ncol=ncol(counts), nrow=ncol(counts) )
|
|
|
14 for (i in 1:ncol(counts)){
|
|
|
15 for (j in 1:ncol(counts)){
|
|
|
16 sere[i,j] <- sigfun_Pearson( counts[,c(i,j)] )
|
|
|
17 }
|
|
|
18 }
|
|
|
19 colnames(sere) <- rownames(sere) <- colnames(counts)
|
|
|
20 return( formatC(sere, format="f", digits=2) )
|
|
|
21 }
|
|
|
22
|
|
|
23 sigfun_Pearson <- function(observed) {
|
|
|
24 #calculate lambda and expected values
|
|
|
25 laneTotals<- colSums(observed);
|
|
|
26 total <- sum(laneTotals)
|
|
|
27 fullObserved <- observed[rowSums(observed)>0,];
|
|
|
28 fullLambda <- rowSums(fullObserved)/total;
|
|
|
29 fullLhat <- fullLambda > 0;
|
|
|
30 fullExpected<- outer(fullLambda, laneTotals);
|
|
|
31
|
|
|
32 #keep values
|
|
|
33 fullKeep <- which(fullExpected > 0);
|
|
|
34
|
|
|
35 #calculate degrees of freedom (nrow*(ncol -1) >> number of parameters - calculated (just lamda is calculated >> thats why minus 1)
|
|
|
36 #calculate pearson and deviance for all values
|
|
|
37 oeFull <- (fullObserved[fullKeep] - fullExpected[fullKeep])^2/ fullExpected[fullKeep] # pearson chisq test
|
|
|
38 dfFull <- length(fullKeep) - sum(fullLhat!=0);
|
|
|
39
|
|
|
40 return(c(sqrt(sum(oeFull)/dfFull)));
|
|
|
41 }
|
