4
$\begingroup$

I'm looking at Broad Institute's orignal GSEA-P algorithm R script which I downloaded here: http://software.broadinstitute.org/gsea/downloads.jsp.

I'm trying to adapt their GSEA.1.0.R script to process datasets that have 1 gene expression profile as for class 1, and N gene expression profiles for class 2. The bottleneck seems to be the GSEA.MetaboliteRanking() method. I looked and it appears they are using "signal to noise ratio" (SNR) metric to rank genes in patient profiles.

I was wondering how the signal to noise is calculated from gene expression profiles?

Note: I am using/modifying their script because it's the only implementation of the GSEA algorithm that I could find where they expose their code.

i.e., what is the mathematical reasoning behind the following lines of code? Note, there are little to no comments in the code, but the following variables, from what I can tell, are:

A: Matrix of gene expression values (rows are genes, columns are samples) 
class.labels: Phenotype of class distinction of interest. A vector of 
              binary labels having first the 1s and then the 0s 



P = class.labels1 * subset.mask
n1 = sum(P[,1])         
M1 = A %*% P
M1 = M1/n1      
A2 = A*A        
S1 = A2 %*% P   
S1 = S1/n1 - M1*M1    
S1 = sqrt(abs((n1/(n1-1)) * S1))

P = class.labels2 * subset.mask
n2 = sum(P[,1])           
M2 = A %*% P           
M2 = M2/n2          
A2 = A*A           
S2 = A2 %*% P      
S2 = S2/n2 - M2*M2 
S2 = sqrt(abs((n2/(n2-1)) * S2))
$\endgroup$
1
$\begingroup$

The noise is comming from the subset.mask, which is created above in a loop with the number of permutations.

for (r in 1:nperm) { #L90
....
        subset.mask[, r] <- as.numeric(c(subset.class1, subset.class2)) # L107

So by multiplying the random selection of subsets by the expression we get the "noise". Later on line 251 we get the ratio of noise/signal by dividing the signal from the contrast and the signal by the noise.

Those transformation on the middle seem to come from normalizing the signal in order to be able to compare the resampling, but I don't fully understand why they do this. Also note that the R scripts are outdated and no longer maintained, so they can contain errors.

|improve this answer|||||
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.