Requiring (At Least Some) Methods Testing for Every Project

It may currently be a little hard to find, but I wanted to point out a couple links relevant to showing the value in testing different RNA-Seq methods for every project:

• SourceForge repository for public data analysis
• I still have a ways to go before being able to start working on a paper, but you can see how I am progressing here
• I think the Target_Recovery_Status.xlsx file (for checking recovery of the known genetic perturbation in an experiment) is the most relevant for showing that you could not choose 1 method out of edgeR, DESeq2, and limma-voom to maximally recover the known gene knock-down or over-expression
• I am also experimenting with having a completely public log for notes and analysis
• Acknowledgement for GitHub RNA-Seq gene expression template
• Includes some papers with modified methods
• While the newer analysis tends to have smaller samples sizes, you can see noticable differences between methods in a much larger cohort in this post

On Biostars (which you can see in a variety of responses, including but not limited to this one), I would generally give the following recommendations:

• If possible, test calculating p-values with edgeR, DESeq2, and limma-voom
• I would recommend having an independently calculated expression method (like FPKM, Fragment Per Kilobase per Million), in order to help assess method selection
• For example, you might see an extremely obvious change in expression for a gene (such as the one that you altered), but it might not have a significant p-value (or have a missing p-value) for one of the methods.
• While the optimal strategy for discovery may not necessarily be the one that most stringently recovers previous results, you may be able to tell some strategies clearly don't work well on your data.
• I would also recommend using this gene expression measurement to create heatmaps to compare clustering of replicates
• I would typically use this instead of exporting normalized counts from the method to calculate the p-value, but testing clustering of replicates (without defining the groups in the normalization) is another possible way to compare strategies.
• Sometimes this can be a bit qualitative.  However, if you define your gene lists / enrichment as a "hypothesis," then I think this is made up for my having independent validation for your claim.
• I do realize this treads the line between p-hacking and needing to test methods due to limits in precision (which I mention a little bit in this comment and this Twitter discussion).  However, as scientists, I think this is part of why it is extremely important to be transparent and admit errors as soon as we discover them (in the interests of training ourselves to be as objective as possible).
• Robustness of identifying a result with different methods may also give you some extra confidence in the results (unless the methods are not really independent, for example)
• If you test alternative normalization, make sure you have a visualization before and after applying that normalization (to try and assess the likelihood of over-fitting in your adjustment)
• I also think it is important that these are open-source, freely available programs (so that you can have the ability to determine what works best for your individual project)

In general, these posts may also be relevant to the discussion of limits to precision in the genomics methods:

• What is the difference between a "Pipeline" and a "Template"?
• Emphasizing "Hypothesis Generation" in Genomics
• For benefits for re-analysis of mutation calling and/or variant annotation, you may also be interested in a some other posts (like this or this)

Again, it is going to be a while, but I do hope to eventually have a preprint to cover the above points (as well as some other observations that I have had from working on a variety of projects for RNA-Seq gene expression analysis).

Change Log:

11/6/2019 - public post
6/3/2020 - add link for earlier (larger) RNA-Seq benchmark

6/7/2022 - minor formatting change

Variance Stabilization and Pseudocounts / Rounding Factors

To be fair, I first want to point out that most of this content was previously posted in an answer Biostars thread.

While I could tell that thread got a fair number of views overall, there wasn't much of a reaction (at least not shortly after posting something on Twitter).  While patience is a virtue, I discovered something when I started working with my own code that I didn't initially notice with the DESeq2 example code.  Plus, I thought some formatting differences may make things more clear, and it might help to separate this point from the slightly different topic of the original Biostars thread.

So, here are those (and additional) results presented in a slightly different way.

First, I thought it was interesting that the plot for rlog values looked a lot like log-transformed values with a pseudocount of 100 (from the meanSdPlot() function in the vsn package, as used in the DESeq2 vignette):

In other words, I am saying the middle plot above (log-transformation with a higher pseudocount) looks a lot more like the rlog plot (on the far right, above).  As a minor point, it may also be worth mentioning that I reformatted the plot (to look more like the other plots that I create below), if you don't immediately recognize the plots above from the Biostars thread.  In the original thread, I thought this was worth mentioning because the log-transformed expression values are independent for each sample, but the rlog values for a given sample will vary depending upon what other samples are processed (so, all other things being equal, I thought the independent per-sample normalization would be better).

However, what I didn't initially notice is what happens if I plot the regular mean versus SD:

In both cases, the higher pseudocount causes there to be less of a bump in the SD values (for lower count genes).  However, with this alternative plot, you can notice more of a difference between the log-transformed values and the rlog values (above, middle-plot versus right-plot).

I previously discussed fold-change values the context of "rounding factors" (in this paper, and this blog post).  While the mean-versus-SD trend is similar, there are some differences from changing the order of operations (CPM with rounding factor in 2nd and 4th columns, with the qualitatively lowest increase in SD with the larger rounding factor in the 4th column):

I frequently find having independently calculated expression values (such as log2(FPKM+X) values) to be helpful in assessing results from differential expression programs (for each project).  However, to be fair, there is in fact some noticeable differences between the log-transformed expression with the higher psuedocount (or rounding factor) and the rlog values from DESeq2.

The code to create the plots above can be downloaded here.  As mentioned in the Biostars thread, I also owe a big thanks to Mike Love (because the 1st part of that code is copied from the DESeq2 vignette).