Hi Michael, I followed the vignette for the section on htseq-count and your suggestion as best as I could and got some results. But I am not sure if I did everything right, so I would really appreciate if you could verify if the following steps is the correct way to analyze this data set: There are 4 samples for two patients: Patient A: abc.htseq.count (Tumor), pqr.htseq.count (Normal) Patient B: def.htseq.count (Tumor), xyz.htseq.count (Normal) > sampleFiles=list.files(directory) > sampleFiles [1] "abc.htseq.count" "def.htseq.count" [3] "pqr.htseq.count" "xyz.htseq.count" > > sampleCondition=c("Tumor","Tumor","Normal","Normal") > samplePatient=c("A","B","A","B") > sampleTable=data.frame(sampleName=sampleFiles, fileName=sampleFiles, condition=sampleCondition, patient=samplePatient) > dds=DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory=directory, design= ~ patient + condition) > > colData(dds) DataFrame with 4 rows and 2 columns condition patient <factor> <factor> abc.htseq.count Tumor A def.htseq.count Tumor B pqr.htseq.count Normal A xyz.htseq.countt Normal B > > dds$condition [1] Tumor Tumor Normal Normal Levels: Normal Tumor > > dds$patient [1] A B A B Levels: A B > > dds = DESeq(dds) estimating size factors estimating dispersions gene-wise dispersion estimates mean-dispersion relationship final dispersion estimates fitting model and testing > > result=results(dds) > > resOrdered2=result[order(result$padj),] > > head(resOrdered2,20) DataFrame with 20 rows and 6 columns baseMean log2FoldChange lfcSE stat pvalue padj <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> GAL1 873531.57 -11.645190 0.7175710 -16.22863 3.164692e-59 5.627771e-55 ... ... I did not get the normalization factors in the output, whereas it does show up in the results in the vignette, so I am not sure if I missed some step. Also, would you expect to get meaningful results with just two matched pairs of samples with this analysis? I am wondering if it is better to just report the FC rather than the p-values. Thank you, - Pankaj ________________________________ From: Michael Love [michaelisaiahlove@gmail.com] Sent: Thursday, April 10, 2014 4:08 PM To: Pankaj Agarwal Cc: bioconductor@r-project.org Subject: Re: [BioC] DESeq2: Paired Test hi Pankaj, You should follow the multifactor analysis section of the vignette. You would include the patient effect in the design like so: design(dds) <- ~ patient + condition Where 'patient' and 'condition' are factors which are columns in colData(dds), and condition has levels: normal, tumor (where normal is the base level). Then results(dds) will build a result table for tumor vs normal, controlling for the patient effect. Mike On Thu, Apr 10, 2014 at 3:52 PM, Pankaj Agarwal <p.agarwal@duke.edu<mailto:p.agarwal@duke.edu>> wrote: Hi, This query is in reference to the following post https://stat.ethz.ch/pipermail/bioconductor/2010-April/032869.html where it was mentioned that "DESeq does not support paired tests (yet)". I have a need for performing rna-seq data analysis with matched tumor/normal pairs of data and was wondering if this capability is now available in the latest version of DESeq2. I searched the list site and also the latest vignette, but could not find anything indicating this capability is now available in DEseq. I would be grateful if, in case this feature is not yet available, someone could advice me how to do the analysis. I have samples from two patients, and from each patient there is a "normal" and a "tumor" rna-seq sample data. Thanks, - Pankaj [[alternative HTML version deleted]] _______________________________________________ Bioconductor mailing list Bioconductor@r-project.org<mailto:bioconductor@r-project.org> https://stat.ethz.ch/mailman/listinfo/bioconductor Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor [[alternative HTML version deleted]]

hi Pankaj, The output of summarizeOverlaps() is a SummarizedExperiment, which can be used in object construction with: dds <- DESeqDataSet(se, ~ condition) for a SummarizedExperiment, se. See the first part of either of the vignettes, and the manual page for ?DESeqDataSet. Or if you can read in any arbitrary counts from a file (using read.table, read.csv, etc.), then check the manual page for ?DESeqDataSetFromMatrix Mike On Wed, Apr 30, 2014 at 4:37 PM, Pankaj Agarwal <p.agarwal at="" duke.edu=""> wrote: > Hi, > > I am trying to use DESeq2 for a set of counts obtained via > summerizedOverlaps and have the counts for 4 samples in one file. Earlier I > had used DESeq2 with counts from htseq-count, which reports the counts for > each samples in a separate file. Now that I have these in one file, is > there a method analogous to "DESeqDataSetFromHTSeqCount" for constructing > the dds object. With htseq-count I had used the following set of commands: > > "directory" contains the 4 count files. > >> sampleFiles=list.files(directory) >> sampleFiles >> [1] "abc.htseq.count" "def.htseq.count" >> [3] "pqr.htseq.count" "xyz.htseq.count" >> >> sampleCondition=c("Tumor","Tumor","Normal","Normal") >> samplePatient=c("A","B","A","B") >> sampleTable=data.frame(sampleName=sampleFiles, fileName=sampleFiles, >> condition=sampleCondition, patient=samplePatient) >> dds=DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, >> directory=directory, design= ~ patient + condition) > > Thanks, > > - Pankaj > ________________________________ > From: Michael Love [michaelisaiahlove at gmail.com] > > Sent: Thursday, April 17, 2014 9:18 AM > To: Pankaj Agarwal > Subject: RE: [BioC] DESeq2: Paired Test > > Hi Pankaj > > For future questions to Bioconductor maintainers, always cc the mailing > list. This helps avoids duplicate questions, which we get a lot of. Answer > below, > > On Apr 17, 2014 8:49 AM, "Pankaj Agarwal" <p.agarwal at="" duke.edu=""> wrote: >> >> Thanks for your help and response. >> >> >> I did not get the normalization factors in the output, whereas it does >> >> show up in the results in the vignette, so I am not sure if I missed >> >> some step. >> >> >I'm confused, which section of the vignette are you referencing? >> >> I was referring to Section 1.5 which report "sizeFactor" as the third >> column in the table, but when I execute the same function I don?t see that >> column: > > You have to run DESeq or estimateSizeFactor first, see ?DESeq > >> >> colData(dds) >> DataFrame with 4 rows and 2 columns >> condition patient >> <factor> <factor> >> abc.htseq.count Tumor A >> def.htseq.count Tumor B >> pqr.htseq.count Normal A >> xyz.htseq.countt Normal B >> >> Is the third column in Section 1.5 the output of estimateSizeFactors( dds >> )? > > Yes > >> >> Thanks, >> >> - Pankaj >> >> -----Original Message----- >> From: Michael Love [mailto:michaelisaiahlove at gmail.com] >> Sent: Wednesday, April 16, 2014 4:37 PM >> To: Pankaj Agarwal >> Subject: Re: [BioC] DESeq2: Paired Test >> >> hi Pankaj, >> >> On Wed, Apr 16, 2014 at 3:59 PM, Pankaj Agarwal <p.agarwal at="" duke.edu=""> >> wrote: >> > Hi Michael, >> > >> > I followed the vignette for the section on htseq-count and your >> > suggestion as best as I could and got some results. >> > But I am not sure if I did everything right, so I would really >> > appreciate if you could verify if the following steps is the correct >> > way to analyze this data set: >> > >> > There are 4 samples for two patients: >> > Patient A: abc.htseq.count (Tumor), pqr.htseq.count (Normal) Patient >> > B: def.htseq.count (Tumor), xyz.htseq.count (Normal) >> > >> >> sampleFiles=list.files(directory) >> >> sampleFiles >> > [1] "abc.htseq.count" "def.htseq.count" >> > [3] "pqr.htseq.count" "xyz.htseq.count" >> >> >> >> sampleCondition=c("Tumor","Tumor","Normal","Normal") >> >> samplePatient=c("A","B","A","B") >> >> sampleTable=data.frame(sampleName=sampleFiles, fileName=sampleFiles, >> >> condition=sampleCondition, patient=samplePatient) >> >> dds=DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, >> >> directory=directory, design= ~ patient + condition) >> >> >> >> colData(dds) >> > DataFrame with 4 rows and 2 columns >> > condition patient >> > <factor> <factor> >> > abc.htseq.count Tumor A >> > def.htseq.count Tumor B >> > pqr.htseq.count Normal A >> > xyz.htseq.countt Normal B >> >> >> >> dds$condition >> > [1] Tumor Tumor Normal Normal >> > Levels: Normal Tumor >> >> >> >> dds$patient >> > [1] A B A B >> > Levels: A B >> >> >> >> dds = DESeq(dds) >> > estimating size factors >> > estimating dispersions >> > gene-wise dispersion estimates >> > mean-dispersion relationship >> > final dispersion estimates >> > fitting model and testing >> >> >> >> result=results(dds) >> >> >> >> resOrdered2=result[order(result$padj),] >> >> >> >> head(resOrdered2,20) >> > DataFrame with 20 rows and 6 columns >> > baseMean log2FoldChange lfcSE stat pvalue >> > padj >> > <numeric> <numeric> <numeric> <numeric> <numeric> >> > <numeric> >> > GAL1 873531.57 -11.645190 0.7175710 -16.22863 3.164692e-59 >> > 5.627771e-55 >> > ... >> > ... >> > >> >> Your lines of code look correct. >> >> > I did not get the normalization factors in the output, whereas it does >> > show up in the results in the vignette, so I am not sure if I missed >> > some step. >> >> I'm confused, which section of the vignette are you referencing? >> >> > >> > Also, would you expect to get meaningful results with just two matched >> > pairs of samples with this analysis? I am wondering if it is better >> > to just report the FC rather than the p-values. >> >> The analysis has 4 samples and 3 parameters to estimate (intercept, >> patient, condition), so you have one residual degree of freedom. So you can >> estimate the dispersion and the p-values should be meaningful. >> >> One drawback with having only two matched pairs is that there are too few >> samples to look for outliers. DESeq2 looks for outliers when there are at >> least 3 samples per condition (actually, per unique combination of factor >> variables in the design). I would take a look at the genes with the highest >> counts and see if these should be filtered manually. >> For example, the first gene in your ordered results table looks to me like >> it just contains an outliers count. However with no replicates, you have to >> perform such filtering manually. >> >> Mike >> >> > >> > Thank you, >> > >> > - Pankaj >> > >> > >> > ________________________________ >> > From: Michael Love [michaelisaiahlove at gmail.com] >> > Sent: Thursday, April 10, 2014 4:08 PM >> > To: Pankaj Agarwal >> > Cc: bioconductor at r-project.org >> > Subject: Re: [BioC] DESeq2: Paired Test >> > >> > hi Pankaj, >> > >> > You should follow the multifactor analysis section of the vignette. >> > You would include the patient effect in the design like so: >> > >> > design(dds) <- ~ patient + condition >> > >> > Where 'patient' and 'condition' are factors which are columns in >> > colData(dds), and condition has levels: normal, tumor (where normal is >> > the base level). >> > >> > Then results(dds) will build a result table for tumor vs normal, >> > controlling for the patient effect. >> > >> > Mike >> > >> > >> > On Thu, Apr 10, 2014 at 3:52 PM, Pankaj Agarwal <p.agarwal at="" duke.edu=""> >> > wrote: >> >> >> >> Hi, >> >> >> >> This query is in reference to the following post >> >> https://stat.ethz.ch/pipermail/bioconductor/2010-April/032869.html >> >> where it was mentioned that "DESeq does not support paired tests >> >> (yet)". >> >> I have a need for performing rna-seq data analysis with matched >> >> tumor/normal pairs of data and was wondering if this capability is >> >> now available in the latest version of DESeq2. I searched the list >> >> site and also the latest vignette, but could not find anything >> >> indicating this capability is now available in DEseq. I would be >> >> grateful if, in case this feature is not yet available, someone could >> >> advice me how to do the analysis. I have samples from two patients, >> >> and from each patient there is a "normal" and a "tumor" >> >> rna-seq sample data. >> >> >> >> Thanks, >> >> >> >> - Pankaj >> >> >> >> [[alternative HTML version deleted]] >> >> >> >> _______________________________________________ >> >> Bioconductor mailing list >> >> Bioconductor at r-project.org >> >> https://stat.ethz.ch/mailman/listinfo/bioconductor >> >> Search the archives: >> >> http://news.gmane.org/gmane.science.biology.informatics.conductor >> > >> >