Search
Question: DESeq2 for Ribosomal profiling analysis [two-factor designs]
0
gravatar for xiaozhuo
3 months ago by
xiaozhuo0
xiaozhuo0 wrote:

Hi Mike,

I am using DESeq2 to analyze differential translation efficiency (TE) from Ribosomal profiling (Riboseq) and RNAseq datasets. Translation efficiency = normalized Riboseq read density / normalized RNASeq read density. I have four groups of samples with duplicates each group for both Riboseq and RNAseq: 1) wild type in physiological condition, 2) wild type in stress condition, 3) mutant in physiological condition, and 4) mutant in stress condition. We aim to identify gene lists where those genes have significantly changed translation efficiency (TE) between different conditions. There are helpful posts giving me examples for just comparing TE between single factor, such as treatment vs control: DESeq2 for Ribosomal profiling analysis (design = ~SeqType + Genotype + SeqType:Genotype). 

However, when there are two factors, such as "WT and mutant" plus "Treatment and control", Can we use DESeq2 for TE analysis? Shall I divide them into two single-factor-analyses, or I can combine them into one two-factor analysis using 

> design(dds) ~SeqType + Genotype + Stress + seqType:Genotype + seqType:Stress

Exp. Desgin

SeqType

Genotype

Stress

Ribo_WT_CTL

Riboseq

WT

CTL

Ribo_WT_Stress

Riboseq

WT

Stress

Ribo_mutant_CTL

Riboseq

mutant

CTL

Ribo_ mutant _Stress

Riboseq

mutant

Stress

RNA_WT_CTL

RNAseq

WT

CTL

RNA_WT_Stress

RNAseq

WT

Stress

RNA_ mutant _CTL

RNAseq

mutant

CTL

RNA_ mutant _Stress

RNAseq

mutant

Stress

 

If using design = ~ SeqType + Genotype + Stress + SeqType:Genotype + SeqType:Stress , how can I get specific set of TE comparison between the following groups?

1) TE_WT_Stress vs TE_WT_CTL:

a1 <- results(dds, contrast=list("SeqType_Riboseq_vs_Rnaseq", "Stress_Stress_vs_CTL", SeqType.Genotype")))

 

2) TE_Mutant_Stress vs TE_Mutant_CTL:

a2 <- results(dds, contrast=list("SeqType_Riboseq_vs_Rnaseq", "Stress_Stress_vs_CTL", SeqType.Genotype_Mutant")))

 

3) TE_Mutant_CTL vs TE_WT_CTL:

a3 <- results(dds, contrast=list("SeqType_Riboseq_vs_Rnaseq", "Genotype_mutant_vs_WT", SeqType.Stress")))

 

4) TE_Mutant_Stress vs TE_WT_Stress:

a4 <- results(dds, contrast=list("SeqType_Riboseq_vs_Rnaseq", "Genotype_mutant_vs_WT", SeqType.Stress_Stress")))

 

Apologies if this is full of mistake and confusing. 

Thanks for the help!

Best,

Xiaozhuo

ADD COMMENTlink modified 3 months ago by Michael Love19k • written 3 months ago by xiaozhuo0
2
gravatar for Michael Love
3 months ago by
Michael Love19k
United States
Michael Love19k wrote:

hi Xiaozhuo,

I think the easiest would be to combine genotype and stress into a single factor, group using factor(paste0(...)), and use a design of ~group + group:type. Then you will have a baseline difference in counts for each group, and then 4 terms in resultsNames(dds) that gives the Ribo/RNA ratio for each group. You can contrast these to find differences in the Ribo/RNA ratio across groups:

results(dds, contrast=list("groupMutStress.typeRibo",
                           "groupWTStress.typeRibo"))

Although it will just say, e.g. typeRibo in the interaction term, it is the Ribo/RNA ratio.

ADD COMMENTlink modified 3 months ago • written 3 months ago by Michael Love19k

hi Mike, 

Great! Thanks for your help as always! For "grouping" function, I notice you had another helpful post before: https://support.bioconductor.org/p/67600/#67612 (DESEq2 comparison with mulitple cell types under 2 conditions). Following this post and reply, here is the code for "grouping" and getting the Ribo/RNA ratio across groups, 

1) Could you tell me whether the code is correct for the Ribo/RNA ratio?

Note: The code coldata$group <- factor(paste0(coldata$genotype, coldata$Stress)) didn't work for my case, so I used dds$group <- factor(paste0(dds$Genotype, dds$Stress)) instead.

dds <- DESeqDataSetFromMatrix(countData = countMat, colData = designMat, design = ~SeqType+Genotype+Stress+SeqType:Genotype+Seq:Stress)
dds$SeqType <- factor(dds$SeqType, levels=c("RNAseq","Riboseq"))
dds$Genotype <- factor(dds$Genotype, levels=c("WT","Mutant")) 
dds$Stress <- factor(dds$Stress, levels=c("CTL","Stress")) 
dds$group <- factor(paste0(dds$Genotype, dds$Stress))
design(dds) <- ~ group + group:Seqtype
dds <- DESeq(dds)
resultsNames(dds)

And here is the resultsNames(dds) output lists given by this code:

"Intercept", "group_MutantStress_vs_MutantCTL", "group_WTCTL_vs_MutantCTL", "group_WTStress_vs_MutantCTL", "groupMutantCTL.SeqTypeRiboseq", "groupMutantStress.SeqTypeRiboseq", "groupWTCTL.SeqTypeRiboseq", "groupWTStress.SeqTypeRiboseq"  

# e.g. TE_WT_Stress vs TE_WT_CTL:

results(dds, contrast=list("groupWTStress.SeqTypeRiboseq", "groupWTCTL.SeqTypeRiboseq"))

# e.g. TE_Mutant_Stress vs TE_Mutant_CTL:

results(dds, contrast=list("groupMutantStress.SeqTypeRiboseq", "groupMutantCTL.SeqTypeRiboseq"))

# e.g. TE_Mutant_CTL vs TE_WT_CTL:

results(dds, contrast=list("groupMutantCTL.SeqTypeRiboseq", "groupWTCTL.SeqTypeRiboseq")) 

# e.g. TE_Mutant_Stress vs TE_WT_Stress:

results(dds, contrast=list("groupMutantStress.SeqTypeRiboseq", "groupWTStress.SeqTypeRiboseq")) 

 

Again, thanks for your time and support!

Best,

Xiaozhuo

ADD REPLYlink modified 3 months ago • written 3 months ago by xiaozhuo0
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 133 users visited in the last hour