Mercurial > repos > yufei-luo > s_mart

comparison SMART/DiffExpAnal/DESeqTools/pairwiseSERE.R @ 31:0ab839023fe4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author m-zytnicki
date Tue, 30 Apr 2013 14:33:21 -0400
parents 94ab73e8a190
children
comparison
equal deleted inserted replaced
30:5677346472b5 31:0ab839023fe4
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 }