Cuffdiff or DESeq

[turnersd]

There's a discussion of this going on over at the bioconductor mailing list. I'm quoting some of the recent posts from the last two days:

Dear Bioconductor list,

Sometime ago Simon Anders explained the difference
between DeSeq and Cuffdiff as follows:

"If you have two samples, cuffdiff tests, for each transcript, whether
there is evidence that the concentration of this transcript is not the
same in the two samples.

If you have two different experimental conditions, with replicates for
each condition, DESeq tests, whether, for a given gene, the change in
expression strength between the two conditions is large as compared to
the variation within each replicate group."

Current language on the Cuffdiff site suggests that the current version
of that program tests for whether the change is significant compared to
changes in each condition.

http://cufflinks.cbcb.umd.edu/howitworks.html#hdif

http://cufflinks.cbcb.umd.edu/howitworks.html#reps

Can someone please comment on the relative merits of Cuffdiff and
DeSeq. I ask here because our sequencing core delivers results
based on Cuffdiff and I want to know if I should redo it using
DeSeq,I would greatly appreciate any guidance in this matter.

Thanks and best wishes,
Rich

Reply 1

Not directly relevant to gene-level RNA-seq DE calls, but rather for
exon-level DE,
I found it useful to read this:

http://precedings.nature.com/documents/6837/version/1

In particular, section 4.3 on page 11, and supplementary figures S7 and S8
on page 19.

I was informed by a coworker that since everyone uses
BowTie-TopHat-Cufflinks-Cuffdiff, it is the sensible thing to do.
Conversations with people who know what they are doing (Terry Speed &
BCGSC) suggest the matter is not yet settled.
So I retrieved ~1TB of BAMs, extracted the reads, and started looking into
how that compares to DEXSeq and/or subread.

It would be incredibly informative if the Cufflinks and DEXSeq authors had
time to weigh in on their strengths/weaknesses. DEXSeq & cummeRbund both
offer nice tools for exploring the results; I am curious which pipeline
fits best for my needs.

Thanks for bringing this up.

Reply 2 (from me)

I'm also in the same boat as Rich. I run a new bioinformatics core
here and I'm building a pipeline for RNA-seq. Cufflinks for some time
has supported biological replicates, and I'm also curious about the
relative merits of using
bowtie/tophat-cufflinks-cuffmerge-cuffdiff-cummeRbund versus using
tophat-HTSeq?-customScriptForCreatingMatrix?-DESeq. Cufflinks also
gives me a host of other tests (differential splicing load,
differential TSS usage, differential coding output, etc), which also
seem useful for certain applications.

On a related note, does anyone have a workflow for taking multiple bam
files, running HTSeq-count (or another program), plus some other
program or custom script, to produce a matrix of counts as input to
DESeq?

Stephen

Reply 4:

Dear Stephen,

To your related note, you could have a look at the easyRNASeq package (bioC 2.10) for R (2.15). It reads in your annotation, your bam files and generate a count table for DESeq, all in R. It can actually do the first step of DESeq (estimating size library and dispersion) and give you back a normalized countDataSet object plus some validation plots as described in the DESeq vignettes (the same is true for edgeR). I'm about to push some changes in SVN to correct an issue that prevented the package vignette to be build. The package should be available in a couple of days as binary or you could install it directly from SVN. See http://wiki.fhcrc.org/bioc/SvnHowTo. The package URL is: https://hedgehog.fhcrc.org/bioconduc...cks/easyRNASeq.

Getting the proper set of annotation is definitely the most important step in the whole process and the one that requires most attention, the rest is then pretty straightforward. What I mean by annotation is the description of your feature of interest, be it gene, transcript, exon, enhancers, etc... as genomic loci (chr, start, width, etc...). The main issue in defining the annotation is to avoid counting reads multiple times and the kind of annotation needed does of course depends on your project. If you are interested in looking at isoforms differential expression, you probably want to define synthetic exons (to avoid double counting) and process the obtained count table with DEXseq. If you're looking at gene expression, you would want to create gene models to avoid multiple counting and use these to create your count table. If you are interested in eRNAs, you can define enhancer loci as the count "feature". All this can be done relatively easily in R. Once you have the proper annotation that suits your need, running easyRNASeq is very straightforward. easyRNASeq accepts both RangedData and GRangesList as annotation input, amo ng other formats. easyRNASeq is in addition able to fetch annotations for you from different sources, but most of the time these would need to be post-processed. You can look at my post: "u! sing easyRNASeq examples" from 2 days ago for some examples and comments in addition to the vignette content.

Cheers,

Nico

And finally, the post that Nico mentioned about easyRNAseq

I've fixed a bug in easyRNASeq and made a couple of modifications to deal with novelties introduced in the IRanges and in the DESeq packages. I've just committed that whole bunch, so it will take ~2 days before a binary package is available. You can always install from SVN in the meanwhile.

I've listed below the different examples you sent me in our email exchange last week and I've added some comments which I hope will prove useful to you. I've created a new "thread" to make it easier to read, and I'm confident it will prove useful to others. This leads me to thanking you for your patience, for providing the data and for giving me the opportunity to post such an example

I hope this helps, do not hesitate to contact me if you have further questions and to suggest improvements.

Cheers,

Nico

## load the library
library(easyRNASeq)

## load the chromosome sizes
## Note that you can just provide a named list; using the BSgenome is not necessary (especially as the
Hsapiens ones are huge >800Mb)
library(BSgenome.Hsapiens.UCSC.hg19)
chr.sizes=as.list(seqlengths(Hsapiens))

## set the wd directory to wherever you have your bam files
## that's just for this example, you could give that directory
## directly to the filesDirectory argument of the function
setwd("Desktop/Francesco")

## get the bam filenames
bamfiles=dir(getwd(),pattern="*\\.bam$")

## run easyRNASeq to get an RNAseq object.
## we use biomaRt to fetch the annotation directly from ensembl
## 
rnaSeq <- easyRNASeq(filesDirectory=getwd(),
                     organism="Hsapiens",
                     chr.sizes=chr.sizes,
                     readLength=100L,
                     annotationMethod="biomaRt",
                     format="bam",
                     count="exons",
                     filenames=bamfiles[1],
                     outputFormat="RNAseq"
                     )

## Fetching the annotation that way is time consuming
## Using the RNAseq object you can extract the genic information
## and save them as an rda for a faster processing time
## later on
gAnnot <- genomicAnnotation(rnaSeq)

## There are 244 "chromosomes" in that annotation, let's keep only what we need
gAnnot <- gAnnot[space(gAnnot) %in% paste("chr",c(1:22,"X","Y","M"),sep=""),]
save(gAnnot,file="gAnnot.rda")

## you could use it this way
countTable <- easyRNASeq(filesDirectory=getwd(),
                         organism="Hsapiens",
                         chr.sizes=chr.sizes,
                         readLength=100L,
                         annotationMethod="rda",
                         annotationFile="gAnnot.rda",
                         format="bam",
                         count="exons",
                         filenames=bamfiles[1]
                         )

## using different annotation (makeTranscriptDB)
library(GenomicFeatures)
hg19.tx <- makeTranscriptDbFromUCSC(
                                    genome="hg19",
                                    tablename="refGene")

## easyRNASeq can deal with GRangesList object, so no need to modify it much, i.e. no need to convert it to a RangedData
gAnnot <- exons(hg19.tx)
## change the metadata column name to suit easyRNASeq
colnames(elementMetadata(gAnnot)) <- "exon"
## finally turn it into a GrangesList
gAnnot <- split(gAnnot,seqnames(gAnnot))

## you could save it as an rda object, you can as well use it directly
## by using the annotationMethod "env" and the annotationObject arguments.
## In addition we will be selecting for a single chromosome: chr1
## as the corresponding name in your bam file is 1, that's what we'll use
countTable <- easyRNASeq(filesDirectory=getwd(),
                         organism="Hsapiens",
                         chr.sizes=chr.sizes,
                         readLength=100L,
                         annotationMethod="env",
                         annotationObject=gAnnot,
                         format="bam",
                         count="exons",
                         filenames=bamfiles[1],
                         chr.sel="1"
                         )

## applying DESeq
## this will not yield very sensitive results as we have no replicates (biological)
## These are important for DESeq to accurately model the technical and biological variance
## With no replicates for every condition, the dispersion will be based on a "pooled" estimate
## making the differential expression call lose sensitivity. In addition, DESeq is in such case
## using a conservative approach (which is good) so you'd get even less significant results.

## Defining the conditions
conditions <- c("A","B")
names(conditions) <- bamfiles

## running DESeq
## NOTE that the two last arguments are for the DESeq estimateDispersions method. They are just
transferred through.
## There are necessary as we have only 1 replicate per condition; see the DESeq vignette for more details.
countDataSet <- easyRNASeq(filesDirectory=getwd(),
                           organism="Hsapiens",
                           chr.sizes=chr.sizes,
                           readLength=100L,
                           annotationMethod="env",
                           annotationObject=gAnnot,
                           format="bam",
                           count="exons",
                           filenames=bamfiles,
                           chr.sel="1",
                           outputFormat="DESeq",
                           normalize=TRUE,
                           conditions=conditions,
                           fitType="local",
                           method="blind"
                         )

## finally my session info
> sessionInfo()
R Under development (unstable) (2012-02-07 r58290)
Platform: x86_64-apple-darwin10.8.0 (64-bit)

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] GenomicFeatures_1.7.16 AnnotationDbi_1.17.15  easyRNASeq_1.1.6      
 [4] ShortRead_1.13.12      latticeExtra_0.6-19    RColorBrewer_1.0-5    
 [7] Rsamtools_1.7.29       DESeq_1.7.6            locfit_1.5-6          
[10] lattice_0.20-0         akima_0.5-7            Biobase_2.15.3        
[13] BSgenome_1.23.2        GenomicRanges_1.7.24   Biostrings_2.23.6     
[16] IRanges_1.13.24        genomeIntervals_1.11.0 intervals_0.13.3      
[19] edgeR_2.5.9            limma_3.11.11          biomaRt_2.11.1        
[22] BiocGenerics_0.1.4    

loaded via a namespace (and not attached):
 [1] annotate_1.33.2    bitops_1.0-4.1     DBI_0.2-5          genefilter_1.37.0 
 [5] geneplotter_1.33.1 grid_2.15.0        hwriter_1.3        RCurl_1.91-1      
 [9] RSQLite_0.11.1     rtracklayer_1.15.7 splines_2.15.0     survival_2.36-12  
[13] tools_2.15.0       XML_3.9-4          xtable_1.6-0       zlibbioc_1.1.1

[chadn737]

The vignette for DESeq is exceptionally well written and makes it easy for a biologist with no R experience to follow it. I say this having at one time being that biologist with NO understanding of R.

[ETHANol]

Without starting some kind of pointless back and forth, I am wondering what are the advantages of using Cuffdiff or DESeq to identify differentially expressed genes from RNA-seq data. The statistical side is a little over my head at the moment but is also what spurs my interest in this question. I find Simon Anders posts very informative so I am inclined to use DESeq but everything R and Bioconductor is a little biologist unfriendly. The Cuffdiff workflow is easier to follow but am I sacrificing anything by using Cuffdiff?