Hi,

First of all, thanks for this great GSVA package.

When I ran GSVA (mx.diff=T) on healthy and diseased samples together, the enrichment scores are quite different from when I ran GSVA on healthy and diseased samples separately.

For example, when I plotted the distribution of enrichment scores for diseased samples when the samples were ran using GSVA separately from the healthy samples, i get a bimodal (a peak at somewhere ES > 0 and one at somewhere ES < 0 ).

When the diseased samples were ran using GSVA together with the healthy samples, the distribution of ES for diseased is mostly unimodal, for example ES > 0 (for genes in gene set that should be upregulated in disease vs. in healthy).

I looked at Fhat GSVA paper equation (1) and see that it is dependent on the distribution of the gene in all samples. So i tried to map out the intermediate steps to understand what each output is-- basically I see that when a gene that is expressed highly when compared to other genes within a diseased sample itself, but lowly expressed as compared to healthy samples and similarly expressed to other disease samples, the ranking of the gene dropped within the sample. However, when no healthy samples were provided, the rank of the gene is higher within the sample.

So, when only diseased samples are included, does it "exaggerate" the difference of the gene expression present in the gene set between diseased samples, hence resulting in the bimodal distribution of ES?

When diseased samples vs. healthy are ran together, technically healthy gene expression in the gene set should be quite different from diseased, hence it gives some weight to the genes in the diseased /healthy that results in the ES separating out disease from healthy.

Do i have the right picture or am I misunderstanding something?

Thanks,

Jenn

hi Jenn,

Indeed, the output of GSVA depends on the number of input samples because the first step of the algorithm calculates an expression statistic that attempts to bring distinct gene expression profiles to a common scale. As explained in the paper, the motivation for this step is that sample-by-gene specific effects, caused for instance by differences in GC content or gene length, may result in identical expression values meaning slightly different things in different genes. While this is not relevant in a classical differential expression analysis, we thought this is important when you want to put together those expression values in a gene set. Because this first step is calculated for every gene across all samples, it will lead to different results depending on how many samples are available, which by the way, we recommend in the paper to be at least 10, but in general the larger the sample size, the more robust will be this first step to differences in the type and number of samples. The algorithm runs in an "unsupervised" manner simply meaning that it does not take into account any phenotype, it does not know whether a sample belongs to a healthy or to a diseased individual.

Regarding the distribution of GSVA enrichment scores, to give an informed opinion I would have to look at your actual data or see at least the plots of the distributions you're talking about. Do we talk about the distribution of scores across genes? across samples? have lowly-expressed genes been filtered out? how many samples are we talking about? what kind of samples are they? has data been normalized? etc.. ==========================Edited answer after discussion below==========

Since I don't have your data, I only can try to simulate it, you can run this by yourself:

library(GSVA)

p <- 10000 ## number of genes
n <- 100   ## number of samples
nGrp1 <- 50 ## number of samples in group 1
nGrp2 <- n - nGrp1 ## number of samples in group 2

## consider three disjoint gene sets
geneSets <- list(set1=sample(paste("g", 1:p, sep=""), size=200, replace=FALSE))

## sample data from a normal distribution with mean 0 and st.dev. 1
y <- matrix(rnorm(n*p), nrow=p, ncol=n,
            dimnames=list(paste("g", 1:p, sep="") , paste("s", 1:n, sep="")))

## consider samples divided in two groups and defined group2 as follows
grp2 <- (nGrp1+1):n

## genes in set1 are expressed at higher levels in the
## second group of samples
y[geneSets$set1, grp2] <- y[geneSets$set1, grp2] + 2
dim(y)
## [1] 10000 100

## run GSVA in all samples
es <- gsva(y, geneSets)

## subset data and the run GSVA
y2 <- y[, grp2]
dim(y2)
## [1] 10000   50
es1 <- gsva(y2, geneSets)
dim(es1)
## [1]  1 50

## run GSVA and then subset data
es2 <- gsva(y, geneSets)[, grp2, drop=FALSE]
dim(es2)
## [1]  1 50

## plot the results
par(mfrow=c(3, 1), mar=c(4, 5, 3, 2))
plot(density(es[1, ]), main="All data")
plot(density(es1[1, ]), main="Subset data")
plot(density(es2[1, ]), main="Subset results")

As the plots show, only when using all the samples you get the bimodal distribution for a gene set which is differentially expressed in half of the samples. In the other two scenarios, where either you subset one half of the data and calculate GSVA enrichment scores, or you calculate the scores on all the data and look at the distribution of the corresponding half of the scores, the distribution is unimodal.

cheers,

robert.