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))