Question: Block x Treatment interaction test with DeSeq
0
7.5 years ago by
Miguel Gallach10 wrote:
Hi list, I am interested in testing Genetic x Environment interaction for RNA- Seq data. However, if I understood correctly, DeSeq do not test for interaction, right? My experimental design consists in comparing expression of two different populations, two replicates per population (Pop A1, Pop A2 vs. Pop B1 and Pop B2), at different temperatures (T1 and T2). The next table represents my data frame and my procedure: myData: A1.T1 A2.T1 B1.T1 B2.T1 A1.T2 A2.T2 B1.T2 B2.T2 Gene1 count count count count count count count count Gene2 count count count count count count count count Gene n count count count count count count count count Desing = data.frame (treatment = ("T1"," T1"," T1"," T1"," T2"," T2"," T2"," T2"), block = c("A1","A2","B1","B2","A1","A2","B1","B2")) cds = newCountDataSet(myData, desing) cds = estimateSizeFactors(cds) cds = estimateVarianceFunctions(cds, method = "pooled") fit0 = nbinomFitGLM(cds, count ~ block) fit1 = nbinomFitGLM(cds, count ~ block + treatment) fit2 = nbinomFitGLM(cds, count ~ block + treatment + block:treatment) pvals2 = nbinomGLMTest (fit2, fit1) padj2 = p.adjust (pvals2, method = "BH") pvals = nbinomGLMTest (fit1, fit0) padj = p.adjust (pvals, method = "BH") myData = data.frame (myData, pvals, padj) May I use the same cds to test fit0, fit1 and fit2? If understand ok, this is not correct since I need to estimate different SizeFactor and VarianceFuntions for each GLM, wright? If I am wright, is then there anyway to test for Block, treatments and Block:Treatment with DeSeq? Thank you very much for your help. Miguel [[alternative HTML version deleted]]
deseq • 857 views
modified 7.5 years ago by Simon Anders3.5k • written 7.5 years ago by Miguel Gallach10
Answer: Block x Treatment interaction test with DeSeq
0
7.5 years ago by
Simon Anders3.5k
Zentrum für Molekularbiologie, Universität Heidelberg
Simon Anders3.5k wrote:
Hi Miguel On 11/28/2011 10:16 AM, Miguel Gallach wrote: > I am interested in testing Genetic x Environment interaction for RNA-Seq > data. However, if I understood correctly, DeSeq do not test for > interaction, right? Of course, DESeq can test for interaction. > My experimental design consists in comparing expression of two different > populations, two replicates per population (Pop A1, Pop A2 vs. Pop B1 and > Pop B2), at different temperatures (T1 and T2). The next table represents > my data frame and my procedure: You table got a bis messed up, so I rather write down how I understood your design data to look like. library population temperature A1.T1 A T1 A2.T1 A T1 B1.T1 B T1 B2.T1 B T1 A1.T2 A T2 A2.T2 A T2 B1.T2 B T2 B2.T2 B T2 This is now assuming that by different "population", you mean something like different genotype or ecotype or colony, and that of either of these, you have four aliquots/cultures/etc, two of which you let grow at temperature T1 and the other two at T2. If you give DESeq the columns 'population' and 'temperature' as a 'design' data frame and the call 'estimateDispersions' with 'method="pooled"', it will consider those sample pairs which have the same population and temperature values as replicates and estimate dispersions accordingly. Then, you can fit models like fit0 <- nbinomFitGLM(cds, count ~ population) fit1 <- nbinomFitGLM(cds, count ~ population + treatment ) fit2 <- nbinomFitGLM(cds, count ~ population + treatment + population:treatment ) The last one is equivalent to fit2 <- nbinomFitGLM(cds, count ~ population * treatment ) Now depending on what you want to test, you may for example do pvals <- nbinomGLMTest( fit1, fit0 ) if you consider 'population' as a blocking factor (i.e., a nuisance covariate whose influence you want to get rid of), and pvals <- nbinomGLMTest( fit2, fit1 ) if you consider both factors as biologically interesting and want to see their interactions. Simon

Dear Simon - I know this is a very late follow up - but there is a reason I need to relate directly to this post.

I am interested in replicating this (i.e. ~ population + treatment + population:treatment, in my case ~ gender + disease + gender:disease, 22 samples in total (each one from a different subject)) with DESeq2. DESeq2 appears to no longer use nbinomFitGLM - does the function have a new name in DESeq2, or is the approach now different?

Could you take me through it step by step from the cds object ( I make it using DESeqDataSetFromHTSeqCount)

cds <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design = ~gender+disease)

cds <- DESeq(cds)

My disease factor has 2 levels: ASD and LTC, and I'm interested in how the disease effects may be different in male and female samples. I have 11 in each of the diesase groups: ASD = 9x male, 2x female, LTC = 4x male, 7x female)

I tried to follow the steps in the manual for standard and expanded model matrices, but got quickly lost.

Thanks,

Matt