Dear BiocCommunity,

i have imported from the R package curatedCRCdata, an RNA-Seq dataset regarding a specific type of cancer. In detail, after contacting the maintainers of the package, they kindely provided the information, that the data are IlluminaHiSeq_RNASeqV2 Level 3 data, which were quantified using RSEM (https://wiki.nci.nih.gov/display/tcga/rnaseq+version+2). When I load the dataset from the package:

library(curatedCRCData)

data(TCGA.RNASeqV2_eset)

TCGA.RNASeqV2_eset

ExpressionSet (storageMode: lockedEnvironment)

assayData: 20502 features, 195 samples

  element names: exprs

protocolData: none

phenoData

  sampleNames: TCGA.AA.3662 TCGA.A6.4105 ... TCGA.A6.6652 (195 total)

  varLabels: unique_patient_ID alt_sample_name ... uncurated_author_metadata (59 total)

  varMetadata: labelDescription

featureData

  featureNames: ? A1BG ... ZZZ3 (20502 total)

  fvarLabels: probeset gene

  fvarMetadata: labelDescription

experimentData: use 'experimentData(object)'

  pubMedIds: 22810696

Annotation: NA

class(TCGA.RNASeqV2_eset)

[1] "ExpressionSet"

attr(,"package")

[1] "Biobase"

head(exprs(TCGA.RNASeqV2_eset)) # a small output

         TCGA.AA.3662 TCGA.A6.4105 TCGA.F4.6463 TCGA.F4.6806 TCGA.A6.6650 TCGA.AZ.6600

?         9.282712     9.933779   10.0443941     9.910121    10.088809     9.875006     9.893715

A1BG      6.027692     4.707000    3.6559242     5.592373     3.253914     5.622860     3.563683

A1CF      8.273718     7.445153    7.0927571     7.118704     7.855290     8.577953     7.774680

range(exprs(TCGA.RNASeqV2_eset))

[1]  0.00000 20.34961

My main questions are the following:

Because from a relative search in other posts/papers, the RSEM does not provide "essentially raw counts", but estimated counts (which are also not rounded). Thus:

1. Is it possible to use for downstream analysis for RNA-Seq data like edgeR ? Or because the counts have e to be integer, other methodologies/packages are eligible for a simple differential expression analysis (for instance a two-group comparison) ?
2. Secondly, also from the values above, it seems that the  counts are also somehow normalized or transformed (perhaps it is the output of "rsem.genes.normalized_results" in the above link of wiki NCI). Unfortunately, i could not find further information about the above level 3 transformation (for how the rsem counts are normalized).  Again, a proper normalization methodology (like TMM) should be necessary for any downstream analysis ?

Please excuse me for any naive questions, but im relatively new to RNA-Seq and at this point any suggestions or opinions would be essential !!

I support the use of rounded estimated counts from software like RSEM (also Sailfish, Salmon, kallisto), and we have a pipeline built for this using the Bioconductor package tximport. See the software vignette which has examples for all of these software including RSEM's *.genes.results file, and downstream example use of edgeR, DESeq2, and limma.

What you should never do is use *normalized* counts with the count-based tools.

Yes, another question is what you have in that matrix. Those do not look like estimated counts, because they max out at 20.3.

Have you looked at the relevant man page:

?TCGA.RNASeqV2_eset

edgeR and limma can both accept RSEM expected counts, and the counts do not need to be rounded to an integer. Both edgeR and limma work with fractional counts.

It is also possible to analyze the (approximately) log2(count+1) values directly in limma, without converting back to counts. You would normalize these with normalizeQuantiles and then use the eBayes(trend=TRUE) pipeline (without voom). Here is an example run to illustrate this simple approach. Here I compare tumour to adjacent normal samples. I have kept genes that have about 10 counts or more in at least 14 samples:

> library(curatedCRCData)
> data(TCGA.RNASeqV2_eset)
> targets <- pData(TCGA.RNASeqV2_eset)
> y <- normalizeQuantiles(exprs(TCGA.RNASeqV2_eset))
> type <- factor(targets$sample_type)
> table(type)
type
adjacentnormal          tumor 
            14            181 
> keep <- rowSums(y > log2(11)) >= 14
> table(keep)
keep
FALSE  TRUE 
 4209 16293 
> y2 <- y[keep,]
> design <- model.matrix(~type)
> fit <- lmFit(y2,design)
> fit <- eBayes(fit,robust=TRUE,trend=TRUE)
> topTable(fit,coef=2)
        logFC AveExpr     t  P.Value adj.P.Val     B
CDH3     6.56  10.294  25.2 1.24e-63  2.02e-59 133.7
KRT80    7.48   8.973  24.0 2.08e-60  1.70e-56 126.5
OTOP2   -9.11   1.504 -22.7 6.43e-57  3.49e-53 118.6
ETV4     6.06  10.469  21.6 1.06e-53  4.31e-50 111.3
ESM1     5.97   6.688  19.2 7.10e-47  2.31e-43  95.9
MTHFD1L  2.65   9.659  17.8 6.66e-43  1.81e-39  86.9
LGI1    -4.51   0.957 -17.8 9.26e-43  2.15e-39  86.5
KRT24   -5.90   0.817 -17.5 5.43e-42  1.11e-38  84.8
GRIN2D   6.07   8.724  17.2 4.46e-41  8.07e-38  82.7
JUB      3.24   8.707  16.9 3.15e-40  5.13e-37  80.8

I note that

plotMDS(y2, label=type)

shows a perfect separation of tumour from normal samples. Also

plotSA(fit)

shows the mean-variance relationship typical of RNA-seq data.

I can offer a couple things -

1) the survival data are in pData(TCGA.RNASeqV2_eset)[, c("days_to_death", "vital_status")] not in the exprs() - apologies if the man page makes that unclear. 

2) the level-III TCGA RNA-seq data are provided as-is, BUT are log2-transformed if the maximum value is greater than 50 (from attach_PROCESSED.R in the pipeline):

if(max(exprs(eset),na.rm=TRUE) > 50){
  ##Do log-transformation
  if(min(exprs(eset),na.rm=TRUE) <= 0){
    exprs(eset) <- exprs(eset) - min(exprs(eset),na.rm=TRUE) + 1
  }
  exprs(eset) <- log2(exprs(eset))
}

FYI the script used to download TCGA RNA-seq data is here (see download_TCGA-RNASeqV2-rectum.sh and download_TCGA-RNASeqV2.sh). These download URLs are broken now that TCGA data have been moved to the GDC, but it they do show the original filenames if you'd like to check them out from GDC:

• unc.edu_COAD.IlluminaHiSeq_RNASeqV2.Level_3.1.4.0.tar.gz
• unc.edu_READ.IlluminaHiSeq_RNASeqV2.Level_3.1.5.0.tar.gz