Question

How to setup a two factor design in DESeq2

0

Entering edit mode

dequattro.concetta ▴ 10

@dequattroconcetta-21510

Last seen 2.2 years ago

Italy

Hi! I have to perform a differential expression analysis for a RIP-seq experiment. In my design I have two variables treatment (miR and control) and Immunoprecipitation (Ago, IgG and TL). The sample table looks likes as follow:

samples treatment   IP  condition
sample1 miR Ago Ago.miR
sample2 miR Ago Ago.miR
sample3 miR Ago Ago.miR
sample4 miR IgG IgG.miR
sample5 miR IgG IgG.miR
sample6 miR IgG IgG.miR
sample7 miR TL  TL.miR
sample8 miR TL  TL.miR
sample9 miR TL  TL.miR
sample10    Control Ago Ago.Control
sample11    Control Ago Ago.Control
sample12    Control Ago Ago.Control
sample13    Control IgG IgG.Control
sample14    Control IgG IgG.Control
sample15    Control IgG IgG.Control
sample16    Control TL  TL.Control
sample17    Control TL  TL.Control
sample18    Control TL  TL.Control

I want to create a two-factorial design considering the variables treatment (miR and control) and Immunoprecipitation. In paticular I would to compare: Ago miR vs Ago Control adjusting for the IgG factor. The comparisons should looks like as follow: (Ago.miR – IgG.miR) vs (Ago. control– IgG. control).

I was wondering if the right way to get this result is to create the dds model and extract the result as follow:

enter code here
ddsHTSeq <- DESeqDataSetFromHTSeqCount(sampleTable = samples, directory = directory, design= ~ treatment * IP)
res <- results(dds, contrast=list(c("treatment_miR_vs_Control", "IP_Ago_vs_IgG")))

Thank you,

Concetta

DESeq2 two-factorialdesign • 2.9k views

ADD COMMENT • link updated 3.2 years ago by Michael Love 41k • written 3.2 years ago by dequattro.concetta ▴ 10

score 0 · Answer 1 · 2021-02-15

0

Entering edit mode

Michael Love 41k

@mikelove

Last seen 1 day ago

United States

The interaction term will represent the ratio of ratios, see the example in the vignette under Interactions.

So you only pull out the interaction term with name in results().

ADD COMMENT • link 3.2 years ago Michael Love 41k

0

Entering edit mode

Thank you for your reply. Unfortunately I did not understand very well how to extract the result. In my experiment I have RNA-seq experiments corresponding to TL and RIP-seq experiments corresponding to IgG and Ago. I would like to identify microRNAs targets comparing Ago vs Control. However in my comparisons I want to correct first for the aspecific binding to the beads considering the IgG factor and then I want to to correct for the transcriptome changes considering the TL factor. I have read this older post about RIP-seq (https://support.bioconductor.org/p/61509/) and I changed my command into:

ddsHTSeq <- DESeqDataSetFromHTSeqCount(sampleTable = samples, directory = directory, design= ~ IP + treatment + IP : treatment)
dds_1 <- DESeq(ddsHTSeq)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship 
final dispersion estimates
fitting model and testing
resultsNames(dds_1)
[1] "Intercept"                    "IP_IgG_vs_Ago"                "IP_TL_vs_Ago"                
[4] "treatment_Control_vs_miR" "IPIgG.treatmentControl"      "IPTL.treatmentControl”

I was wondering how I can extract the result of Ago miR vs Ago Control corrected for IgG and the results of Ago miR vs Ago Control corrected for TL. Thank you for your help.

ADD REPLY • link 3.2 years ago dequattro.concetta ▴ 10

0

Entering edit mode

results() can pull out the coefficents by name.

results(dds, name="...")

When you have interaction terms (the 5th and 6th term above), these represent ratios of ratios as explained in the vignette diagram.

I might recommend you work with a statistician to help interpret the linear model terms.

Also see the vignette on how to set factor levels. Here you have miR as the reference level, instead of Control.

ADD REPLY • link 3.2 years ago Michael Love 41k