The first thing that comes to mind is to use the Standard Deviation (SD). This will give you a measure of the variability in the spread of the data from the averages. Using something like the Coefficient of Variation (CV), which is the SD divided by the average, you can express the same value relative to the average. I think the latter is more appropriate in this case, since it is easier to compare across features/genes.

Hi Mahmoud,

Would it be valid to use the "gene read count" file from here?

https://gtexportal.org/home/datasets

It seems that genes are found in there listed by their ensemble gene IDs, with the gene read counts for each sample within gtex. Are these the values to enter into the analysis you described in your initial answer?

Many thanks,

Dave

Yes, David Green

In case of the RNA-Seq data, the read counts are what you can use to represent the gene expression.

Depending on the dataset was processed, you might want to do some quality control before you apply any analysis though. I think the ones to consider are filtering, normalization and transformation.

Also, because you are trying to estimate and compare the variance of the genes, you might want to consider the type of replicates in this dataset; biological or technical.

Let me warn you against using unnormalized, unfiltered read counts for a few reasons.

1) Even though differences in read depth between samples will be accounted for by dividing the gene SD by the gene mean, this does not adjust for differences in read length, which can greatly impact mean read counts. Using an expression metric that has been normalized by gene length (e.g. TPM, FPKM or RPKM), will take care of this issue. TPM is the preferred normalization these days, and if you are using GTEx data, they have TPMs as one of the data types you can download.

2) CV obviously depends very strongly on mean expression such that lowly expressed genes are going to tend to have much higher CV. Therefore, we cannot trust the CV of very lowly expressed genes, and so these should be filtered out of your analysis. If you want a well tested expression filtering scheme, you could use GTEx's scheme (Genes were selected based on expression thresholds of >0.1 TPM in at least 20% of samples and ≥6 reads in at least 20% of samples).

3) There could also be a concern with the heteroscedasticity of the expression data. This is the phenomenon of variance depending strongly on mean expression. Now one may hope that this will be fully handled in the CV calculation since you are normalizing the variance to the mean expression, but that is only true if the association between variance and mean is linear. You can reduce heteroscedasticity concerns by either log2 or VST transforming your normalized expression values. For a more detailed discussion on this, check out this blog post: https://seqqc.wordpress.com/2015/02/16/should-you-transform-rna-seq-data-log-vst-voom/

I've attached an image with some scatterplots from the GTEx Brain-Caudate RNA-seq dataset comparing CV to mean expression using counts, log2(counts), TPMs or log2(TPMs). They show the difference that normalizing and log2 transforming your data can make in your CV results, as well as the crazy CV values given by lowly expressed genes. I'd most trust the results represented in the log2(TPM) analysis (after also filtering out lowly expressed genes).

Best of luck!