Entering edit mode
I am having trouble getting specific coefficient names after running DESeq2. The coefficients I get when I run resultsNames(dds) are "Treat1", "Treat2", etc instead of "Treat_B_vs_A", "Treat_C_vs_A", etc. Has anyone else had this issue? The code that I was running and the output are below:
samples$Treat <- factor(samples$Treat, levels = c("A","B","C","D","E","F"))
samples$Treat
# A A A B B B C C C D D D E E E F F F
# Levels: A B C D E F
dds <- DESeqDataSetFromTximport(txi,
colData = samples,
design = ~ Subject + Treat)
keep <- rowSums(counts(dds) >= 10) >= 3
dds <- dds[keep,]
dds$Treat <- relevel(dds$Treat, ref = "A")
dds <- DESeq(dds)
resultsNames(dds)
# "Intercept" "Subject1" "Subject2" "Treat1" "Treat2" "Treat3" "Treat4" "Treat5"
sessionInfo( )
# R version 3.6.0 (2019-04-26)
# Platform: x86_64-apple-darwin15.6.0 (64-bit)
# Running under: macOS 10.15.7
Thank you!
Hi Michael,
Yes, Subject is also a simple factor.
samples$Subjectreturns:12 14 15 12 14 15 12 14 15 12 14 15 12 14 15 12 14 15 12 14 15 12 14 15 Levels: 12 14 15
class(samples)returns"data.frame"
The strange thing is that I had it working properly last week, but it stopped giving specific coefficients when I ran it this week. I am not sure what I did differently - nothing I don't think.
Check if your installation is still ok:
BiocManager::valid()I get:
Bioconductor version '3.9'
Versions of Bioconductor packages used for the analysis: tximport_1.12.3 tximportData_1.12.0 DESeq2_1.24.0 apeglm_1.6.0 IHW_1.12.0 rhdf5_2.28.1
I did the package update recommended by Bioconductor and re-ran DESeq2 but am still getting the weird coefficient names.
I'm not sure, could you try with the same data on a different machine?
Also what do you get with:
I have 8 conditions with 3 biol reps each. I releveled the treatment factor so that every sample is compared to A (control).
As a (slight) aside, how would I specify the comparison B vs C if Treat2 is B vs A and Treat 2 is C vs A? Is this correct? contrast = c("Treat", "Treat3", "Treat2")
A_rep1 1 -1 -1 -1 -1 -1 -1 -1
A_rep2 1 -1 -1 -1 -1 -1 -1 -1
A_rep3 1 -1 -1 -1 -1 -1 -1 -1
B_rep1 1 1 -1 -1 -1 -1 -1 -1
B_rep2 1 1 -1 -1 -1 -1 -1 -1
B_rep3 1 1 -1 -1 -1 -1 -1 -1
C_rep1 1 0 2 -1 -1 -1 -1 -1
C_rep2 1 0 2 -1 -1 -1 -1 -1
C_rep3 1 0 2 -1 -1 -1 -1 -1
D_rep1 1 0 0 3 -1 -1 -1 -1
D_rep2 1 0 0 3 -1 -1 -1 -1
D_rep3 1 0 0 3 -1 -1 -1 -1
E_rep1 1 0 0 0 4 -1 -1 -1
E_rep2 1 0 0 0 4 -1 -1 -1
E_rep3 1 0 0 0 4 -1 -1 -1
F_rep1 1 0 0 0 0 5 -1 -1
F_rep2 1 0 0 0 0 5 -1 -1
F_rep3 1 0 0 0 0 5 -1 -1
G_rep1 1 0 0 0 0 0 6 -1
G_rep2 1 0 0 0 0 0 6 -1
G_rep3 1 0 0 0 0 0 6 -1
H_rep1 1 0 0 0 0 0 0 7
H_rep2 1 0 0 0 0 0 0 7
H_rep3 1 0 0 0 0 0 0 7
attr(,"assign")
[1] 0 1 1 1 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$Treat
[1] "contr.helmert"
I see -- the issue is that you've set a non-default setting for contrasts in your
.Rprofile.I'd recommend to take that out if you want default contrasts from DESeq2.
Ahh that's right, I used helmert contrasts for a different project and didn't realize it changed it in the .Rprofile. I set it back to default and now it works. Mystery solved! Thanks so much for your prompt help and for writing these tools!