Hi everyone, I am reanalyzing a published RNA-seq dataset from GEO as a learning exercise to better understand the DESeq2 workflow. I downloaded the original FASTQ files, performed quality control, alignment, read counting, and differential expression analysis using the standard DESeq2 pipeline. My experimental design was: Factor 1 = Control (3 biological replicates) Factor 2 = Simvastatin treatment (3 biological replicates) The analysis completed successfully, and the adjusted p-values and significant genes appear reasonable. However, I noticed that the direction of the log2FoldChange is opposite to what is reported in the published paper. For example, genes that the authors describe as upregulated in the treatment group appear as negative log2FoldChange in my DESeq2 results, while genes reported as downregulated appear positive. Interestingly, when I reverse the comparison (using the treatment as the reference level), my results become almost identical to those reported in the publication. This makes me wonder: Could the authors have used the treatment group as the reference level in DESeq2 without explicitly stating it? Is there any way to determine the reference level used in a published RNA-seq analysis if the paper does not mention it? Has anyone encountered a similar situation when reproducing published DESeq2 analyses? I would appreciate any advice on how to verify this issue before concluding that the difference is only due to the reference level rather than an error in my analysis. Thank you!
•
written
by
ashraffouad1998 •
1 answer
You don't have to guess what the authors did, just check what you did. resultsNames(dds) prints the contrast explicitly (something like condition_Simvastatin_vs_Control), which tells you what a positive log2FC means in your own table. And from now on pass it explicitly, results(dds, contrast = c("condition","Simvastatin","Control")), then the default level ordering stops mattering at all.
For working out whether you or the paper has the sign right, use the biology instead of the reference level. Simvastatin inhibits HMG-CoA reductase, the cells sense the sterol drop and activate SREBP-2, so HMGCR, LDLR, HMGCS1, INSIG1 and the rest of that feedback set should come out UP in the treated samples. Look at the sign on those genes in your results -- if they're positive, your direction is correct no matter what the figure legend claims. That's a much firmer check than trying to reverse-engineer someone's factor levels.
Also worth knowing this isn't rare: published figures do sometimes have the direction flipped, or the text describes it one way while the supplementary table is the other way round. If your sterol genes behave sensibly, trust your own run.
Log in to answer this question.
-
DESeq2
written by t.ru •
-
fishpond RNA-seq analysis: Two genotype and four treatments
written by scsc185 •
-
Single Factor Design and Dispersion and BCV
written by plewis8 •
-
RNA-Seq differential gene expression without biological and technical replicates.
written by utsafar
-
How many genes are expressed in each condition: control / treatment group
written by ek699 •
-
DEXSeq output from Galaxy
written by anand.gururajan •
-
Differential Expression Help
written by mlai2567 •
-
Question about analysis design in DeSeq2
written by jillianrjaycox •
-
Differential expressed genes in DESeq
written by Lila M
-
Question about the DESeq2: log2foldchange contradictory
written by tiago211287
If you just eyeball the normalized counts, what does that look like?
A/B vs B/A mistakes can easily happen - from simple oversight to subtly devious bugs in R code when strings are turned into factors.