I recently rewrote my R code so that it can be used with new dataset. In doing so, I changed the nomenclature of my sample Groups (used in Design). For some reason, the "Control vs. Treated" comparison is now being made, instead of the intended "Treated vs. Control".
Before -> Now
Auxin_00min -> Control
Auxin_20min -> Auxin
Note that alphabetically, the prior nomenclature naturally ordered Auxin_00min before Auxin_20min. However, I had always included a relevel step, just in case, and always got the expected comparison (Auxin_20min vs. Auxin_00min).
control_name = "Auxin_0min"
dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~Treatment)
relevel(dds$Treatment, ref = control_name)
## [1] Auxin_0min Auxin_0min Auxin_0min Auxin_20min Auxin_20min Auxin_20min
## Levels: Auxin_0min Auxin_20min
dds <- DESeq(dds)
results <- results(dds, alpha = 0.05)
resOrdered <- results[order(results$padj, -results$log2FoldChange), ] #Order by increasing padj, decreasing LFC
resOrdered
## log2 fold change (MLE): Treatment Auxin 20min vs Auxin 0min
## Wald test p-value: Treatment Auxin 20min vs Auxin 0min
## DataFrame with 6530 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue
## <numeric> <numeric> <numeric> <numeric> <numeric>
## YLR044C 1794.59 -3.39792 0.0829905 -40.9434 0.00000e+00
## YHR174W 1641.12 -2.90740 0.0798101 -36.4289 1.48365e-290
## YGR254W 1776.95 -3.07650 0.0852199 -36.1007 2.21288e-285
Alphabetically, the new nomenclature orders Auxin before Control, but the relevel step should solve that issue. I receive no errors, but the results reveal that the incorrect "Control vs Auxin" comparison has been made.
UPDATE: changing groups with relevel() worked in no capacity - not even when trying to set ref="Auxin". However, changing 'Auxin' to 'Treated' so that it once again fell alphabetically after 'Control' caused the expected 'Treatment vs. Control' comparison to be made again. So I'm fairly certain there is some issue with relevel.
control_name = "Control"
dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~Treatment)
relevel(dds$Treatment, ref = control_name)
## [1] Control Control Control Auxin Auxin Auxin
## Levels: Control Auxin
dds <- DESeq(dds)
results <- results(dds, alpha = 0.05)
resOrdered <- results[order(results$padj, -results$log2FoldChange), ] #Order by increasing padj, decreasing LFC
resOrdered
## log2 fold change (MLE): Treatment Control vs Auxin
## Wald test p-value: Treatment Control vs Auxin
## DataFrame with 6530 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue
## <numeric> <numeric> <numeric> <numeric> <numeric>
## YLR044C 1794.59 3.39792 0.0829905 40.9435 0.00000e+00
## YHR174W 1641.12 2.90740 0.0798101 36.4289 1.48361e-290
## YGR254W 1776.95 3.07650 0.0852199 36.1007 2.21281e-285
Hi Michael,
Those were copy-paste errors. I'm running this within Rmarkdown and was copying and pasting from the completed report. The first step in the Rmarkdown code clear the environment, so there's no carry over from prior work. The report completes successfully and I get all of my exports.
Since changing the word 'Auxin' to 'Treatment' reestablishes the original alphabetic order and simultaneously reverts the comparison to the expected 'Treatment vs. Control', I'm fairly certain that the issue is with relevel.
Here's a summary of different Group namings and the comparisons they yielded:
All of these were run with relevel set to the value listed under (Control).
I just noticed you are not storing your relevel command. That function (nor others in R) does not modify arguments in place.
E.g. if you have
x, then:Typically will not modify
x(unless the developers are doing something very unexpected vis a vis the language).lol. Did I mention I’m a beginner? 4 hours of troubleshooting to find out I needed to save as a variable. I believe that is the weekend calling me. :P