Compare more than one factor
1
0
Entering edit mode
Yuya Liang ▴ 20
@yuya-liang-14202
Last seen 5.7 years ago
United States

 

Hello everyone,

I have questions about how to write design for DESeqDataSetFromMatrix and using "name" and "contrast" argument in results().

Here's how my data looks like:

> sample_info
                   X Genotype     Trt
1     Flood_4610_D14     4610 Control
2     Flood_4610_D14     4610 Control
3     Flood_4610_D14     4610 Control
4     Flood_4610_D14     4610 Control
5    Flood_Rondo_D14    Rondo Control
6    Flood_Rondo_D14    Rondo Control
7    Flood_Rondo_D14    Rondo Control
8    Flood_Rondo_D14    Rondo Control
9  Drought_Rondo_D14    Rondo Drought
10 Drought_Rondo_D14    Rondo Drought
11 Drought_Rondo_D14    Rondo Drought
12 Drought_Rondo_D14    Rondo Drought
13  Drought_4610_D14     4610 Drought
14  Drought_4610_D14     4610 Drought
15  Drought_4610_D14     4610 Drought
16  Drought_4610_D14     4610 Drought​

 

The aim of my study is to identify what genes express differently in 4610 compare to Rondo under drought treatment. I wrote "design = ~ Genotype + Trt + Genotype:Trt" in DESeqDataSetFromMatrix()

 

My question is: How to see the comparison for "Drought_4610" vs "Drought_Rondo" in results()? (I can only do either Drought vs Control or 4610 vs Rondo. How to consider both factor at the same time?)

 

Thank you,

 
deseq2 • 929 views
ADD COMMENT •
0
Entering edit mode

I followed A: DESEq2 comparison with mulitple cell types under 2 conditions to create group.

DESeq.ds$group <- factor(paste0(DESeq.ds$Genotype, DESeq.ds$Trt))
design(DESeq.ds) <- ~ group
dds <- DESeq(DESeq.ds)
resultsNames(dds)​

But my resultsNames(dds)​ only has:

> resultsNames(dds)

[1] "Intercept"                         "group_4610Drought_vs_4610Control"  "group_RondoControl_vs_4610Control"

[4] "group_RondoDrought_vs_4610Control"

 

Doesn't have "group_RondoDrought_vs_4610Drought".

ADD REPLY •
0
Entering edit mode
@mikelove
Last seen 28 minutes ago
United States

You can use the contrast argument. Take a look at ?results
ADD COMMENT •
0
Entering edit mode

Thank you! I just realized I can simply use `contrast = c("group", "4610Drought", "RondoDrought"` to compare. I have another question. The DEGs I got here, did they compare with their own control first? ((i.e. do 4610Drought vs 4610Control and RondoDrought vs RondoControl to see what gene different under drought in each genotype. After that, compare two genotypes.) 

Or just compare "4610Drought" and "RondoDrought" without considering "4610Contorl" and "RondoControl"?
 

 > res <- results(dds, contrast = c("group", "4610Drought", "RondoDrought"))
> head(res)
log2 fold change (MLE): group 4610Drought vs RondoDrought 
Wald test p-value: group 4610Drought vs RondoDrought 
DataFrame with 6 rows and 6 columns
                       baseMean log2FoldChange     lfcSE       stat     pvalue      padj
                      <numeric>      <numeric> <numeric>  <numeric>  <numeric> <numeric>
OS07G0656900       6.897109e+02     -0.4559368 0.3011922 -1.5137737 0.13008324 0.6777180
OS10G0188900       1.714817e+03     -0.2148277 0.1513032 -1.4198492 0.15565159 0.7143898
EPlOSAG00000042159 5.019592e-01      1.6213232 2.4343327  0.6660237 0.50539597        NA
OS02G0192150       3.408006e+00      0.6896324 0.9069645  0.7603742 0.44703097 0.9054431
EPlOSAG00000042157 6.907324e-02      0.2376521 4.3280446  0.0549098 0.95621032        NA
EPlOSAG00000042155 1.591523e+01     -0.8045365 0.4864350 -1.6539445 0.09813879 0.6225005

ADD REPLY •
0
Entering edit mode

Thank you, Michael! 
But can I set two different ref level for two genotypes respectively? 


DESeq$group <- relevel(DESeq$group, ref = c("4610Control", "RondoControl"))
Error in DESeq$group : object of type 'closure' is not subsettable

ADD REPLY •
1
Entering edit mode

After re-reading your above comment, I realize that I think you want to test the difference of differences here. The easiest way to do this is to first replace the design with ~genotype*treatment,

design(dds) <- ~genotype * treatment

and then use the following:

dds <- DESeq(dds, test="LRT", reduced=~genotype + treatment)
res <- results(dds)

This design is equivalent to ~group, but while the first makes it easier to perform pairwise comparisons, the second will give you a test of the difference of differences (this is called an interaction).

ADD REPLY •
0
Entering edit mode

Thank you, Michael!!
I was reading your previous post DESEQ2 results and contrast (and a relevel question). That's really clearly explaining. After thinking about design matrix, I suddenly understand.
   
Before reducing= ~ genotype + treatment, should I still set the reference level?

ADD REPLY •
0
Entering edit mode

You should generally always set the reference level.

ADD REPLY •
0
Entering edit mode

Dear Michael,

I have another problem. I was trying to visualize DEGs by heatmap but somehow the graph looks weird. Doesn't have the different expression.

df <- as.data.frame(colData(DESeq.ds)[,c("Genotype","Trt")])
mat <- assay(DESeq.rlog) head(rownames(mat))
idx <- rownames(res)[ which(res$padj < 0.1) ]
pheatmap(mat[idx, ], cluster_rows=FALSE, show_rownames=TRUE, cluster_cols=FALSE, annotation_col=df)

Thank you, Yuya

ADD REPLY •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 810 users visited in the last hour
Content Search
Users
Tags
Badges
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6