I am performing an analysis on RNA-Seq data, where only genes relating to specific pathways are of interest to the researcher.
1) One option is normalizing, estimating the dispersion and performing DE with DESeq2 as usual (since DESeq's assumption that most genes are not differentially expressed pertains to the whole set of genes, not to the subset).
Following that, it would be possible to manually select only the relevant subset of genes, and apply FDR only to this specific subset, ( based on the p-values calculated when taking all genes into account).
This is somewhat analogous imho to what independent-filtering does. Independent-filtering after calculating p-values for all genes, subsets the list for only those with mean higher than a certain cuttoff, maximizing that cutoff. The explicitly stated goal is to increase the number of significantly DE genes, the rationale being that genes with genes with low expression are not interesting in the first place.
Here, what defines which genes are interesting is not the mean level, but the inclusion in a specific set.
Would the described process be suitable, and is there another one if not?
Hi! Just wanted to know the update or any suggestions about your queries above.
As I am also doing DEGs using RNA-seq on a subset of protein family.
The power of common RNA-seq tools comes from borrowing information from all genes to most accurately estimate dispersions along the baseMean gradient. Later subsetting genes to a few in order to get smallest possible FDRs for them after using all genes to achieve this "power" sounds odd to me. I would do a standard analysis (after prefiltering for counts, see edgeR and DESeq2 manuals) and then take the stats as they're returned. It's cherrypicking if you do lots of custom subsetting.
Sorry to revive this ancient thread but I'm currently discussing the same question with some people in my lab and would love to somebody pointing out the error in this line of reasoning, because so far I actually don't see it.
As far as I understand, DESeq2's normalization/dispersion correction shrinks variance & logFCs by using information from the whole dataset. This results in more realistic raw p-values and effect size estimates, thereby also reducing type I errors, as explicitly stated in the seminal paper from Love et al. When running thousands of hypothesis tests, FDR is corrected for to keep type I error at 5% (or however i set my FDR threshold).
If I am now interested in a set of target genes (defined a priori due to a certain biological questions, not after looking at the data), and I just happen to have a whole genome dataset, isn't it simply the most robust approach to run DESeq2 to get variance/effect size shrinkage (again, thereby also reducing type I error, which is kind of the opposite of p-hacking?) before I run my statistics? FDR correcting over 1000s of genes would then however also be incorrect because i am only interested in 30 hypothesis. So to not inflate type II error, I extract the raw p-values from DESeq2s output and run benjamini hochberg over those.
The alternative would be to just use the raw data without normalization, run Wilcoxon Rank Sum tests over those and have a higher risk of false positives and negatives given high dispersion rates that stem from my typically small sample size. Subsequent multiple comparison correction on these test would also only use the number of genes in my set, so I don't 'gain' anything in terms of statistical rigor when using this approach. In both variants FDR-correction is based on the number of genes that I'm interested in. I only lose the normalization/dispersion shrinkage that should increase the robustness of my inference.
I would really appreciate some insight on this question :D
Best regards