Getting raw counts needed for Deseq/EdgeR from TCGA RSEM files

[dnet]

Hi,

I wish to run DE analysis using DESeq or EdgeR on RNA-seq data downloaded from TCGA. I would like to use the Level 3 RNA-Seq data, which is already processed using RSEM.

I wonder if I can use the column named "raw counts" in the RSEM un-normalized output as the raw read counts needed for the input for DESeq and EdgeR.

For example, the column marked in bold in the file :

Filename:
unc.edu__IlluminaHiSeq_RNASeqV2__TCGA-A1-A0SB-01A-11R-A144-07__expression_rsem_gene.txt


barcode gene_id raw_count scaled_estimate transcript_id
TCGA-A1-A0SB-01A-11R-A144-07 ?|100130426 0 0 uc011lsn.1
TCGA-A1-A0SB-01A-11R-A144-07 ?|100133144 34.05 1.23812E-06 uc010unu.1,uc010uoa.1


Thanks !
D. N.

[mmuurr]

i don't have an answer, but essentially am curious about the same point.
i believe the TCGA Level 3 RNASeqv2 "unnormalized" data represents the 'raw' RSEM counts, and thus piping this input into edgeR would be fine... but perhaps i'm mistaken.
one spot where i've seen conflicting opinions is on how edgeR handles non-integer based counts, which will be the case with the RSEM output.

in a few test cases i've run, i haven't encountered any glaring errors, though i found a HUGE number of differentially expressed genes when comparing prostate cancer samples (both unmatched cases using exact tests and using a subset of matched cases using a GLM approach to handle the paired samples).
at an FDR threshold of 0.05 (using B-H correction), nearly half the genome qualified as differentially expressed, which -- at first glance -- seemed high to me.

[Simon Anders]

Are you sure that you used counts per gene, and not counts per transcript, as input? RSEM outputs the latter by default, but these are unsuitable (even in principle) for downstream analysis for differential expression testing. (If you don't know why, see my earlier posts on the subject.)

[mmuurr]

yes, the counts are represented at the gene level for the publicly available TCGA RNASeqv2 (unnormalized) data.
while i can't find the verbose output for the execution of their RNASeqv2 pipeline, my guess is that the RSEM mappings to transcripts are collapsed to the individual gene level by summing counts.
so, there are ~20,000 genes represented in their "Level 3" files (lower levels, representing increasingly raw data -- e.g. the reads themselves -- are not all publicly available).

as for the large number of differentially expressed genes, additional reading lends me to believe the non-integer counts do need to be rounded prior to edgeR analyses.
(though even this rounding step has been the focus of some debate on the R/Bioconductor-help mailing list.)

[GenePool]

Hi All,

For what it's worth, we're committed to making this sort of data more freely available and usable by the community. In that spirt, we've included a freely available reference library of genomics data in our product, GenePool. This library happens to include the RNASeqV2 gene-level counts computed by the UNC pipeline that leveraged RSEM. We've also taken the time to extract and curate the sample-level metadata and make it easily available to researchers to subset the samples, and analyze the data accordingly. For more advanced users, you can easily just export out the counts and sample level metadata and get into more high-powered statistical analyses, that hopefully some day we just roll right back into the GenePool platform :-)

Incidentally, we've also included the isoform-, splice-junction-, and exon-level counts as part of GenePool's premium content.

If you're interested in learning more, please check out GenePool's growing genomics library, check out the following threads:

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

We'd love to have your feedback on this effort.

------------------------------
GenePool is making genomics data management, analysis, and sharing easier!
Products @ www.stationxinc.com