[BioC] DESeq2

[Michael Love]

hi Aroa,

On Mon, May 12, 2014 at 5:22 AM, AROA SUÁREZ VEGA <asuav at unileon.es> wrote:
> Hi Michael,
>
>
>
> Now that the new version is available I have some doubts about how to do the
> contrast.
>
>
> I run the following commands:
>
> ##Loading DESeq2
>
> source("http://bioconductor.org/biocLite.R")
>
> biocLite("DESeq2")
>
> library("DESeq2")
>
>
>
> setwd("~/DESeq2")
>
> milk<- read.csv("global_count_genes_milk.csv", header= TRUE, row.names=1,
> sep=" ")
>
> milk<-subset(milk,select=c("A1_D10","A2_D10","A3_D10","A4_D10","C1_D10","C2_D10","C3_D10","C4_D10","A1_D50","A2_D50","A3_D50","A4_D50","C1_D50","C2_D50","C3_D50","C4_D50","A2_D120","A3_D120","A4_D120","C1_D120","C3_D120","C4_D120","A1_D150","A2_D150","A3_D150","A4_D150","C1_D150","C2_D150","C3_D150","C4_D150"))
>
> milk_design=data.frame(row.names = colnames(milk), breed = c("A", "A",
> "A","A","C","C","C","C","A","A","A","A","C","C","C",
> "C","A","A","A","C","C","C","A","A","A","A","C","C","C","C"), day =
> c("D10","D10","D10","D10","D10","D10","D10","D10","D50","D50","D50","D50","D50","D50","D50","D50","D120","D120","D120","D120","D120","D120","D150","D150","D150","D150","D150","D150","D150","D150"))
>
>
>
> ##Creating DESeqDataSetFromMatrix
>
> dds<- DESeqDataSetFromMatrix(countData= milk, colData= milk_design, design=
> ~ breed + day)
>
> dds$breed<- factor(dds$breed, levels=c("A","C"))
>
> dds$day<- factor(dds$day, levels=c("D10","D120","D150","D50"))
>
>

The following code is unnecessary, from here...

>
> ##Normalization procedure
>
> dds<-estimateSizeFactors(dds)
>
> dds<-estimateDispersions(dds)
>
> dds<-nbinomWaldTest(dds)
>
>

...to here. These lines perform the same operation as the DESeq() call
below, except using the wrong design formula. Furthermore, the results
from these calls above are erased by calling DESeq(), except for the
estimateSizeFactors(dds) call, which produces the same result as below
because it doesn't depend on the design.

>
> ##Setting model and performing analysis
>
> ddsX<-dds
>
> design (ddsX)<-~breed+day+breed:day
>
> ddsX<-DESeq(ddsX)
>
>
>
> Running the model with interactions I obtain the following results names:
>
> resultsNames(ddsX)
>
> [1] "Intercept"           "breedA"          "breedC"
>
>  [4] "dayD10"              "dayD120"             "dayD150"
>
>  [7] "dayD50"              "breedA.dayD10"   "breedC.dayD10"
>
> [10] "breedA.dayD120"  "breedC.dayD120" "breedA.dayD150"
>
> [13] "breedC.dayD150" "breedA.dayD50"   "breedC.dayD50"
>
>
>
> I know that if I run the following contrast:
>
>
> res_CvsA<-results(ddsX, contrast=list("breedC","breedA"))
>
>
> I obtain the differentially expressed genes between C and A for all the
> points.
>

Yes. And this is equivalent to using contrast=c("breed","C","A").

>
> And if I run:
>
>
> resD10_CvsA<-results(ddsX, contrast=list("breedC.dayD10","breedA.dayD10"))
>
>
> I obtain the interaction effect in addition to the main effect of breed C
> over breed A.

Yes.

> But, if I want to obtain only the differentially expressed
> genes at that point without taken into account the baseline.

This is then the main effect and interaction effect added together:

results(ddsX, contrast=list(c("breedC","breedC.dayD10"),c("breedA","breedA.dayD10")))

> In my
> experiment I want to know the differentially expressed genes between the two
> breeds but it would be also interesting to obtain the differentially
> expressed genes between the two breeds at each time point. How do I have to
> proceed?
>
>
> We also want to obtain the genes differentially expressed for each breed
> comparing the different time points.
>

For example, the day 120 vs 10 effect for breed A is:

results(ddsX, contrast=list(c("dayD120","breedA.dayD120"),c("dayD10","breedA.dayD10")))

Mike

>
> Thank you very much in advance.
>
>
> Aroa
>
>
> Aroa Suárez Vega
> PhD student
> Dpto. Producción Animal
> Facultad de Veterinaria
> Campus de Vegazana s/n
> 24071 Leon
> Telef. 987291000 Ext. 5296
>
>