[sudders]

DEseq and edgeR (and other tools that test raw gene-level counts) assume that the features are the same length. However, if you have any differential splicing or isoform switching going on between conditions, this assumption gets you into trouble. Suppose for example, you have an gene with two isoforms, one twice as long as the other. In both conditions being compared, you have 100 reads on the gene. However, in condition A, all the reads are on the shorter one. In B, they're all on the longer one. Comparing the raw gene-level count will show no change, even though the expression of the gene is two-fold higher in condition A! This is a contrived example (the reads wouldn't *all* move like that in a real situation), but it illustrates the point.

I noted this in my original post. Unfortunately, several of the experimental designs I use require more complex designs than just testing group A against group B. Most often this is because of a paired design, but there are several experiments that involve testing interactions between factors etc. which just can't be done in cufflinks. What we can do with the tag based system is test differential expression of exons, and then find some way of selecting an exon to represent a gene, with the understanding that this won't be valid if there is a differential splicing event that changes the usage of that exon.

[turnersd]

Thanks Ethan. Looking forward to trying that one out, or possibly modifying it to work with DEXSeq.

[ETHANol]

This thread was pretty informative to me, so thanks everyone for your comments.

I have been using DESeq. I wrote a little script to wrap up the process, so I thought I would share it. Maybe someone will find it useful.

Anyway, it's here:

http://ethanomics.wordpress.com/2012/04/19/ezdeseq/

[Cole Trapnell]

On the other hand, I'm still to find a solution to the read counting that truly satisfies me. DESeq and edgeR get away without normalising for length because the two entities being compared are the same length. However, if you sum over all exons in a gene then this isn't necessarily the case: different exons might be being used in different samples. .

DEseq and edgeR (and other tools that test raw gene-level counts) assume that the features are the same length. However, if you have any differential splicing or isoform switching going on between conditions, this assumption gets you into trouble. Suppose for example, you have an gene with two isoforms, one twice as long as the other. In both conditions being compared, you have 100 reads on the gene. However, in condition A, all the reads are on the shorter one. In B, they're all on the longer one. Comparing the raw gene-level count will show no change, even though the expression of the gene is two-fold higher in condition A! This is a contrived example (the reads wouldn't *all* move like that in a real situation), but it illustrates the point.

[Cole Trapnell]

Using cuffdiff means you are coming from tophat-cufflinks-cuffmerge workflow. I have found that this workflow is bugged and mixes up gene annotations. It can not be said with 100% certainty if the reported FPMKs indeed belong to the corresponding gene/transcript.

Please view my thread if interested:

http://seqanswers.com/forums/showthread.php?t=18070

I am now using DESeq for my analysis.

The thread you linked to doesn't show a bug, just an incorrect use of Cufflinks, as I just noted in that thread.

Regarding use of Cufflinks or RSEM with count-based tools like DEseq or edgeR - the authors of DEseq have noted on a number of occasions that it is NOT a good idea to use RSEM or Cufflinks with their tool. The reason is that those tools produce (either explicitly or internally during estimation of gene and transcript FPKM) an *estimate* of the number of fragments that came from each transcript. However, that estimate comes with *uncertainty*, and those tools don't account for it when testing for DE. Think of it this way - when you have multiple isoforms, you have to decide to assign your fragments among them. In addition to the biological variability that's present between replicates, there's now also variability in this assignment process. Those tools don't account for it, so using them with the estimated counts from transcript-level tools will lead to false positives.

[epi]

Using cuffdiff means you are coming from tophat-cufflinks-cuffmerge workflow. I have found that this workflow is bugged and mixes up gene annotations. It can not be said with 100% certainty if the reported FPMKs indeed belong to the corresponding gene/transcript.

Please view my thread if interested:

http://seqanswers.com/forums/showthread.php?t=18070

I am now using DESeq for my analysis.

[sudders]

DESeq/edgeR also allow you to use generalized linear modelling to use more complicated statistical models. For example to do paired same analysis (e.g. healthy tissue vs cancer tissue from 10 patients). In this sort of situation I find that cuffdiff finds no genes/isoforms are differentially expressed because the patient-to-patient variance is too high.

On the other hand, I'm still to find a solution to the read counting that truly satisfies me. DESeq and edgeR get away without normalising for length because the two entities being compared are the same length. However, if you sum over all exons in a gene then this isn't necessarily the case: different exons might be being used in different samples. You can test exon-by-exon, but then have to work out how you go from differentially expressed exons to differentially expressed genes. I admit I haven't looked at RSEM, i'll go and do that.

[gaberudy]

In regards to the question of how to get whole-number counts appropriate for using DESeq, I would recommend using RSEM for quantification: http://deweylab.biostat.wisc.edu/rsem/

RSEM only does quantification of known-transcripts but does so in a very stable manner (comapred to cufflinks) that gives gene-level counts that are the sum of the transcript-level counts for that gene (which cufflinks does not do). It also is more stable from version to version and has less issues in "failing" to produce a count for certain isoforms/genes.

RSEM can take BAM files as input, but expects them to be aligned to the transcriptome (i.e. won't take the output of TopHat that includes both genomic and txome alignments) so you may end up re-doing alignment to use it. But RSEM will happily call bowtie for you if you give it fastq as input.

It will do probablistic assignment of multi-reads giving you counts that can be fractions of whole numbers, but simply rounding these numbers before inputting them to DEseq should give you better results than simply doing a raw counting of reads in genes/exons.

[gringer]

... but I'm sitting at the wrong computer at the moment. I'll attach / embed the useful scripts tomorrow after I get in to work.

I've attached my python script for doing basic contig counting. I use it like this (assuming all experiment IDs start with 'Sample_'):

$ ~/scripts/sam2rawCounts.py $(for x in $(ls -d Sample_*); do \
echo ${x}=$(ls ${x}/*contigs*.bam); done) > contig_reads.txt

In the absence of a consistent directory prefix for sample IDs, I create an 'experiment_dirs.txt' file that identifies the directories that contain the required BAM files:

$ ~/scripts/sam2rawCounts.py $(for x in $(cat experiment_dirs.txt); do \
echo ${x}=$(ls ${x}/*contigs*.bam); done) > contig_reads.txt

I then convert the output to a DESeq-compatible table in R using xtabs:

# get counts from transcriptome mapping
reads.3col.df <- read.table("contig_reads.txt", header = TRUE, stringsAsFactors = FALSE);

## convert to raw count table formatted for DESeq
experiment.counts <- xtabs(Reads ~ Contig + Experiment, data = reads.3col.df);

## assume no biological replicates
experiment.conditions <- sub("^Sample_","",colnames(experiment.counts));

## convert to DESeq CountDataSet object
experiment.cds <- newCountDataSet(unclass(experiment.counts), experiment.conditions);

[lgoff]

Hi all,

For my own purposes, looking for expression differences between isoforms does not make my results anymore meaningful and also makes the types of analysis I am interested in far more difficult.

There are certainly cases when evaluating differences between isoforms is not only the preferred way to conduct the analysis, but is in fact the only way that you would identify a significant change in the expression of a gene. Additionally, biologically significant differential expression of a isoform can be missed when the smallest unit analyzed is a summary model of a 'gene'.

The increased difficulty in interpretation of cuffdiff output is something that I have tried to address with the 'cummeRbund' R/Bioconductor package. The hope is that making the cuffdiff output manageable in a familiar environment, while at the same time maintaining the relationships between genes, isoforms, promoters, etc, will begin to make the interpretation of cuffdiff output more engaging.

[chadn737]

Anyway, enough of my ranting against Bioconductor. I'm still interested in opinions about what is happening under the hood of these methods.

Cuffdiff actually uses the same underlying model as DESeq:

For single isoform genes, Cuffdiff models the variance in fragment counts across replicates using the negative binomial distribution, similar to the method described by Anders and Huber (2010)

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

However since Cufflinks attempts to test differential expression for isoforms, their model gets a lot more complex and they attempt adjust expression values with reads that map to multiple isoforms/genes.

My understanding of statistics is far to rudimentary to actually speak on which performs better. For my own purposes, looking for expression differences between isoforms does not make my results anymore meaningful and also makes the types of analysis I am interested in far more difficult.

Also, in every paper I have seen that compares various methods of testing differential expression, DESeq and EdgeR give very similar results, and DESeq gives fairly similar results to BaySeq, whereas Cufflinks can differ considerably from all of these.

Cufflinks can also take a considerable time to run and even with extra steps, I can finish up using DESeq well before Cufflinks ever completes.

If its not obvious by now, I prefer DESeq and for that matter EdgeR, over Cufflinks. It fits my needs better.

[chadn737]

I have HTseq-count in my pipeline. At the end I just import all the individual count files into R and create the table there.

Or for those who are very unfamiliar with perl, R, or python. Just open the count files that are output by HTseq-count in excel. Then copy and paste to make your table and save it as a csv file. I realize this is labor intensive, but I would expect basically every scientist to at least know excel.

[ETHANol]

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.

Agree. Really, it is one of the few Bioconductor vignettes that make any sense.

Taking note of what else has been mentioned, easyRNASeq (which I hadn't heard about before) sounds like an even better tool for this given that it stays in R. As long as the BAM file reading/processing is quick enough (e.g. considering 5-50GB BAM files for each experiment), there's no reason not to do everything in R.

Actually, I always look for any reason not to use R. R workflows are time consuming in both trying to put together a workflow and using a workflow once you have put it together.

Anyway, enough of my ranting against Bioconductor. I'm still interested in opinions about what is happening under the hood of these methods.

I've heard far too many examples of people using tophat-cufflinks-cuffmerge-cuffdiff for their RNA-seq workflow, and getting extremely different results from one version of cufflinks to the next.

Sometimes people look at data expression data sets that are nearly identical and end up concluding they are very different. For example, you are looking for differentially expressed genes in two data sets. You get a list of differentially expressed genes from software A. Software A gives you a value you use for cutoff say FDR. You make the cutoff and then run the same data set through software B and make the exact same arbitrary cutoff. You look at the list of differentially expressed genes in the two lists and they are quite different because of small variations put some genes above the arbitrary cutoff line and some below. In the conversion from a measurable variable to a nominal variable you lost a lot of information and it makes the data sets look very different when if you look before you made the arbitrary cutoff the data sets looked nearly identical. I'm not saying that's what people have done but it's at the very least possible.

[gringer]

I hacked out a perl script that takes a list of bam files, runs samtools view on them, pipes that output into htseq-count, and drops an output file. It does this for each bam file, dropping a temporary file containing all the features and counts for that particular alignment.
...
I would like to see this discussion continue, and if anyone has any better solution for going from bam alignments to a DESeq-ready counts table, I'm all ears.

There's really no need to use Perl for that if you're also using HTSeq. HTSeq has been designed as a somewhat easy-to-use Python library -- I say 'somewhat' because in some very small points, it didn't quite do what I wanted.

I've got a python script lying around that takes in multiple BAM files (or SAM files, if that's your liking) and spits out counts for each feature in a three-column format (feature, experiment, count) using HTSeq. I'm mapping to a transcriptome for the reads so haven't bothered with GTF files yet, but that's something that can be fairly easily plugged in using the HTSeq library. It would probably also be fairly easy to modify HTSeq-count to support no GFF files and multiple experiments.

Once you have the three column format, it's a somewhat easy process to convert into a count matrix using the xtabs function in R.

... but I'm sitting at the wrong computer at the moment. I'll attach / embed the useful scripts tomorrow after I get in to work.

Taking note of what else has been mentioned, easyRNASeq (which I hadn't heard about before) sounds like an even better tool for this given that it stays in R. As long as the BAM file reading/processing is quick enough (e.g. considering 5-50GB BAM files for each experiment), there's no reason not to do everything in R.

[turnersd]

One of the issues I faced with DESeq was going from alignments (bam files) to the counts table that DESeq requires, e.g.:

      Sample1 Sample2 Sample3 Sample4 Sample5
Gene1    2614    7000    7119    6689    4893
Gene2    1162    2939      72    4988     390
Gene3    5799    8533    9263    8413    5635
Gene4    8129    4023    1524    2270    1816
Gene5    3536    4320    4927    2659    6795

From numerous responses here and elsewhere, I started using HTSeq (specifically, htseq-count) to do this.

htseq-count is dead simple to use. It takes an alignment and a GTF annotation file (e.g. htseq-count [options] <alignment> <gtf>) and produces a list of features and counts of reads that overlap those features (how overlapping reads are counted is specified by the options).

WASH5P 3
WDR18 201
ZNF555 44
ZNF57 82
ZNF77 55
psiTPTE22 0
tAKR 0
no_feature 6321
ambiguous 196
too_low_aQual 0
not_aligned 0
alignment_not_unique 15890

I hacked out a perl script that takes a list of bam files, runs samtools view on them, pipes that output into htseq-count, and drops an output file. It does this for each bam file, dropping a temporary file containing all the features and counts for that particular alignment. I then pipe together a series of join commands to glue together all these temporary count text files that HTSeq produces. This gets me my table that I can use with DESeq.

I've heard far too many examples of people using tophat-cufflinks-cuffmerge-cuffdiff for their RNA-seq workflow, and getting extremely different results from one version of cufflinks to the next. But when it takes some perl hacking around a python script to munge alignments into a format that R can read, it's currently easier IMHO to use the cufflinks solution. In addition, cufflinks gives you other useful information about splicing load, differential TSS usage, differential coding output, etc.

I would like to see this discussion continue, and if anyone has any better solution for going from bam alignments to a DESeq-ready counts table, I'm all ears.