Import and summarize transcript-level abundance estimates for transcript- and gene-level analysis with Bioconductor packages, such as edgeR, DESeq2, and limma-voom. The motivation and methods for the functions provided by the tximport package are described in the following article (Soneson, Love, and Robinson 2015):

Charlotte Soneson, Michael I. Love, Mark D. Robinson (2015): Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research http://dx.doi.org/10.12688/f1000research.7563.1

In particular, the tximport pipeline offers the following benefits: (i) this approach corrects for potential changes in gene length across samples (e.g. from differential isoform usage) (Trapnell et al. 2013), (ii) some of the upstream quantification methods (Salmon, Sailfish, kallisto) are substantially faster and require less memory and disk usage compared to alignment-based methods that require creation and storage of BAM files, and (iii) it is possible to avoid discarding those fragments that can align to multiple genes with homologous sequence, thus increasing sensitivity (Robert and Watson 2015).

Import transcript-level estimates

We begin by locating some prepared files that contain transcript abundance estimates for six samples, from the tximportData package. The tximport pipeline will be nearly identical for various quantification tools, usually only requiring one change the type argument. We begin with quantification files generated by the Salmon software, and later show the use of tximport with any of:

First, we locate the directory containing the files. (Here we use system.file to locate the package directory, but for a typical use, we would just provide a path, e.g. "/path/to/dir".)

library(tximportData)
dir <- system.file("extdata", package = "tximportData")
list.files(dir)
##  [1] "cufflinks"            "kallisto"             "kallisto_boot"       
##  [4] "rsem"                 "sailfish"             "salmon"              
##  [7] "salmon_gibbs"         "samples.txt"          "samples_extended.txt"
## [10] "tx2gene.csv"

Next, we create a named vector pointing to the quantification files. We will create a vector of filenames first by reading in a table that contains the sample IDs, and then combining this with dir and "quant.sf".

samples <- read.table(file.path(dir, "samples.txt"), header = TRUE)
samples
##   pop center                assay    sample experiment       run
## 1 TSI  UNIGE NA20503.1.M_111124_5 ERS185497  ERX163094 ERR188297
## 2 TSI  UNIGE NA20504.1.M_111124_7 ERS185242  ERX162972 ERR188088
## 3 TSI  UNIGE NA20505.1.M_111124_6 ERS185048  ERX163009 ERR188329
## 4 TSI  UNIGE NA20507.1.M_111124_7 ERS185412  ERX163158 ERR188288
## 5 TSI  UNIGE NA20508.1.M_111124_2 ERS185362  ERX163159 ERR188021
## 6 TSI  UNIGE NA20514.1.M_111124_4 ERS185217  ERX163062 ERR188356
files <- file.path(dir, "salmon", samples$run, "quant.sf")
names(files) <- paste0("sample", 1:6)
all(file.exists(files))
## [1] TRUE

Transcripts need to be associated with gene IDs for gene-level summarization. If that information is present in the files, we can skip this step. For Salmon, Sailfish, and kallisto the files only provide the transcript ID. We first make a data.frame called tx2gene with two columns: 1) transcript ID and 2) gene ID. The column names do not matter but this column order must be used. The transcript ID must be the same one used in the abundance files.

Creating this tx2gene data.frame can be accomplished from a TxDb object and the select function from the AnnotationDbi package. The following code could be used to construct such a table:

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
k <- keys(txdb, keytype = "GENEID")
df <- select(txdb, keys = k, keytype = "GENEID", columns = "TXNAME")
tx2gene <- df[, 2:1]  # tx ID, then gene ID

Note: if you are using an Ensembl transcriptome, the easiest way to create the tx2gene data.frame is to use the ensembldb packages. The annotation packages can be found by version number, and use the pattern EnsDb.Hsapiens.vXX. The transcripts function can be used with return.type="DataFrame", in order to obtain something like the df object constructed in the code chunk above. See the ensembldb package vignette for more details.

Here we read in a pre-constructed tx2gene table:

tx2gene <- read.csv(file.path(dir, "tx2gene.csv"))
head(tx2gene)
##         TXNAME   GENEID
## 1    NM_130786     A1BG
## 2    NR_015380 A1BG-AS1
## 3 NM_001198818     A1CF
## 4 NM_001198819     A1CF
## 5 NM_001198820     A1CF
## 6    NM_014576     A1CF

The tximport package has a single function for importing transcript-level estimates. The type argument is used to specify what software was used for estimation (“kallisto”, “salmon”, “sailfish”, and “rsem” are implemented). A simple list with matrices, “abundance”, “counts”, and “length”, is returned, where the transcript level information is summarized to the gene-level. The “length” matrix can be used to generate an offset matrix for downstream gene-level differential analysis of count matrices, as shown below.

Note: While tximport works without any dependencies, it is significantly faster to read in files using the readr package. If tximport detects that readr is installed, then it will use the readr::read_tsv function by default. A change from version 1.2 to 1.4 is that the reader is not specified by the user anymore, but chosen automatically based on the availability of the readr package. Advanced users can still customize the import of files using the importer argument.

library(tximport)
library(readr)
txi <- tximport(files, type = "salmon", tx2gene = tx2gene)
names(txi)
## [1] "abundance"           "counts"              "length"             
## [4] "countsFromAbundance"
head(txi$counts)
##             sample1   sample2    sample3   sample4   sample5   sample6
## A1BG     109.232000 316.22400 110.638000 116.00000  86.38430  76.91630
## A1BG-AS1  83.969700 138.44900 119.274000 151.08300 123.98500 103.25100
## A1CF       9.030691  10.01847   5.019242  13.01820  25.21914  25.07356
## A2M       24.000000   2.00000  21.000000   6.00000  38.00000   8.00000
## A2M-AS1    1.000000   1.00000   1.000000   1.00000   0.00000   0.00000
## A2ML1      3.047950   1.02987   4.076160   1.04945   3.07761   5.12409

We could alternatively generate counts from abundances, using the argument countsFromAbundance, scaled to library size, "scaledTPM", or additionally scaled using the average transcript length, averaged over samples and to library size, "lengthScaledTPM". Using either of these approaches, the counts are not correlated with length, and so the length matrix should not be provided as an offset for downstream analysis packages. For more details on these approaches, see the article listed under citation("tximport").

We can avoid gene-level summarization by setting txOut=TRUE, giving the original transcript level estimates as a list of matrices.

txi.tx <- tximport(files, type = "salmon", txOut = TRUE, tx2gene = tx2gene)

These matrices can then be summarized afterwards using the function summarizeToGene. This then gives the identical list of matrices as using txOut=FALSE (default) in the first tximport call.

txi.sum <- summarizeToGene(txi.tx, tx2gene)
all.equal(txi$counts, txi.sum$counts)
## [1] TRUE

Salmon / Sailfish

Salmon or Sailfish quant.sf files can be imported by setting type to "salmon" or "sailfish".

files <- file.path(dir, "salmon", samples$run, "quant.sf")
names(files) <- paste0("sample", 1:6)
txi.salmon <- tximport(files, type = "salmon", tx2gene = tx2gene)
head(txi.salmon$counts)
##             sample1   sample2    sample3   sample4   sample5   sample6
## A1BG     109.232000 316.22400 110.638000 116.00000  86.38430  76.91630
## A1BG-AS1  83.969700 138.44900 119.274000 151.08300 123.98500 103.25100
## A1CF       9.030691  10.01847   5.019242  13.01820  25.21914  25.07356
## A2M       24.000000   2.00000  21.000000   6.00000  38.00000   8.00000
## A2M-AS1    1.000000   1.00000   1.000000   1.00000   0.00000   0.00000
## A2ML1      3.047950   1.02987   4.076160   1.04945   3.07761   5.12409
files <- file.path(dir, "sailfish", samples$run, "quant.sf")
names(files) <- paste0("sample", 1:6)
txi.sailfish <- tximport(files, type = "sailfish", tx2gene = tx2gene)
head(txi.sailfish$counts)
##             sample1   sample2    sample3   sample4   sample5   sample6
## A1BG     109.165000 316.13800 110.525000 116.00000  86.26030  76.75400
## A1BG-AS1  85.582100 141.15800 120.811000 153.49000 127.03400 105.32600
## A1CF       9.034322  10.02056   5.020884  13.02060  25.23152  25.07919
## A2M       24.000000   2.00000  21.000000   6.00000  38.00000   8.00000
## A2M-AS1    1.000000   1.00000   1.000000   1.00000   0.00000   0.00000
## A2ML1      3.046160   1.02936   4.074320   1.04871   3.07782   5.11940

Note: for previous version of Salmon or Sailfish, in which the quant.sf files start with comment lines, it is recommended to specify the importer argument as a function which reads in the lines beginning with the header. For example:

txi <- tximport("quant.sf", type = "none", txOut = TRUE, txIdCol = "Name", abundanceCol = "TPM", 
    countsCol = "NumReads", lengthCol = "Length", importer = function(x) read_tsv(x, 
        skip = 8))

Salmon / Sailfish with inferential replicates

If inferential replicates (Gibbs or bootstrap samples) are present in expected locations relative to the quant.sf file, tximport will import these as well, if txOut=TRUE. tximport will not summarize inferential replicate information to the gene-level. Here we demonstrate using Salmon, run with only 5 Gibbs replicates (usually more Gibbs samples would be useful for estimating variability).

files <- file.path(dir, "salmon_gibbs", samples$run, "quant.sf")
names(files) <- paste0("sample", 1:6)
txi.inf.rep <- tximport(files, type = "salmon", txOut = TRUE)
names(txi.inf.rep)
## [1] "abundance"           "counts"              "infReps"            
## [4] "length"              "countsFromAbundance"
names(txi.inf.rep$infReps)
## [1] "sample1" "sample2" "sample3" "sample4" "sample5" "sample6"
dim(txi.inf.rep$infReps$sample1)
## [1] 178136      5

The tximport arguments varReduce and dropInfReps can be used to summarize the inferential replicates into a single variance per transcript and per sample, or to not import inferential replicates, respectively.

kallisto

kallisto abundance.h5 files can be imported by setting type to "kallisto". Note that this requires that you have the Bioconductor package rhdf5 installed. (Here we only demonstrate reading in transcript-level information.)

files <- file.path(dir, "kallisto_boot", samples$run, "abundance.h5")
names(files) <- paste0("sample", 1:6)
txi.kallisto <- tximport(files, type = "kallisto", txOut = TRUE)
head(txi.kallisto$counts)
##                   sample1 sample2 sample3 sample4 sample5 sample6
## ENST00000448914.1       0       0       0       0       0       0
## ENST00000631435.1       0       0       0       0       0       0
## ENST00000632684.1       0       0       0       0       0       0
## ENST00000434970.2       0       0       0       0       0       0
## ENST00000415118.1       0       0       0       0       0       0
## ENST00000633010.1       0       0       0       0       0       0

kallisto with inferential replicates

Because the kallisto_boot directory also has inferential replicate information, it was imported as well (and because txOut=TRUE). As with Salmon, inferential replicate information will not be summarized to the gene level.

names(txi.kallisto)
## [1] "abundance"           "counts"              "infReps"            
## [4] "length"              "countsFromAbundance"
names(txi.kallisto$infReps)
## [1] "sample1" "sample2" "sample3" "sample4" "sample5" "sample6"
dim(txi.kallisto$infReps$sample1)
## [1] 178136      5

kallisto with TSV files

kallisto abundance.tsv files can be imported as well, but this is typically slower than the approach above.

files <- file.path(dir, "kallisto", samples$run, "abundance.tsv")
names(files) <- paste0("sample", 1:6)
txi.kallisto.tsv <- tximport(files, type = "kallisto", tx2gene = tx2gene)
head(txi.kallisto.tsv$counts)
##             sample1   sample2    sample3   sample4   sample5  sample6
## A1BG     108.581000 314.42400 110.450000 116.00000  85.80300 75.91360
## A1BG-AS1  86.163600 140.10700 129.994000 146.40800 136.92800 97.29540
## A1CF       9.003863  12.01096   3.005232  15.01082  24.01285 22.01611
## A2M       24.000000   2.00000  21.000000   6.00000  38.00000  8.00000
## A2M-AS1    1.000000   1.00000   1.000000   1.00000   0.00000  0.00000
## A2ML1      3.012760   1.01650   3.049480   2.04965   2.02477  3.04483

RSEM

RSEM sample.genes.results files can be imported by setting type to "rsem".

files <- file.path(dir, "rsem", samples$run, paste0(samples$run, ".genes.results"))
names(files) <- paste0("sample", 1:6)
txi.rsem <- tximport(files, type = "rsem")
head(txi.rsem$counts)
##          sample1 sample2 sample3 sample4 sample5 sample6
## A1BG       94.64  278.03   94.07   96.00   55.00   64.03
## A1BG-AS1   64.28  114.08   98.88  109.05   95.32   73.11
## A1CF        0.00    2.00    1.00    1.00    0.00    1.00
## A2M        24.00    2.00   18.00    4.00   35.00    8.00
## A2M-AS1     1.00    1.00    1.00    0.00    0.00    0.00
## A2ML1       0.84    2.89    0.00    1.00    2.00    3.11

Use with downstream Bioconductor differential expression packages

Note: there are two suggested ways of importing estimates for use with gene-level differential expression methods. The first method, which we show below for edgeR and for DESeq2, is to use the estimated counts from the quantification tools, and additionally to use the transcript-level abundance estimates to calculate an offset that corrects for changes to the average transcript length across samples. The code examples below accomplish these steps for you, keeping track of appropriate matrices and calculating these offsets.

The second method is to use the tximport argument countsFromAbundance="lengthScaledTPM" or "scaledTPM", and then to use the count matrix txi$counts directly as you would a regular count matrix with these software.

edgeR

An example of creating a DGEList for use with edgeR (M. D. Robinson, McCarthy, and Smyth 2010):

library(edgeR)
cts <- txi$counts
normMat <- txi$length
normMat <- normMat/exp(rowMeans(log(normMat)))
library(edgeR)
o <- log(calcNormFactors(cts/normMat)) + log(colSums(cts/normMat))
y <- DGEList(cts)
y$offset <- t(t(log(normMat)) + o)
# y is now ready for estimate dispersion functions see edgeR User's Guide

DESeq2

An example of creating a DESeqDataSet for use with DESeq2 (Love, Huber, and Anders 2014):

library(DESeq2)

The user should make sure the rownames of sampleTable align with the colnames of txi$counts, if there are colnames. The best practice is to read sampleTable from a CSV file, and to construct files from a column of sampleTable, as was shown in the tximport examples above.

sampleTable <- data.frame(condition = factor(rep(c("A", "B"), each = 3)))
rownames(sampleTable) <- colnames(txi$counts)
dds <- DESeqDataSetFromTximport(txi, sampleTable, ~condition)
# dds is now ready for DESeq() see DESeq2 vignette

limma-voom

An example of creating a data object for use with limma-voom (Law et al. 2014). Because limma-voom does not use the offset matrix stored in y$offset, we recommend using the scaled counts generated from abundances, either "scaledTPM" or "lengthScaledTPM":

files <- file.path(dir, "salmon", samples$run, "quant.sf")
names(files) <- paste0("sample", 1:6)
txi <- tximport(files, type = "salmon", tx2gene = tx2gene, countsFromAbundance = "lengthScaledTPM")
library(limma)
y <- DGEList(txi$counts)
y <- calcNormFactors(y)
design <- model.matrix(~condition, data = sampleTable)
v <- voom(y, design)
# v is now ready for lmFit() see limma User's Guide

Acknowledgments

The development of tximport has benefited from contributions and suggestions from:

Rob Patro (inferential replicates import), Andrew Parker Morgan (RHDF5 support), Ryan C. Thompson (RHDF5 support), Matt Shirley (ignoreTxVersion), Stephen Turner, Richard Smith-Unna, Rory Kirchner, Martin Morgan,

Session info

sessionInfo()
## R version 3.4.0 (2017-04-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.2 LTS
## 
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.5-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.5-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] DESeq2_1.16.0                          
##  [2] SummarizedExperiment_1.6.0             
##  [3] DelayedArray_0.2.0                     
##  [4] matrixStats_0.52.2                     
##  [5] edgeR_3.18.0                           
##  [6] limma_3.32.0                           
##  [7] readr_1.1.0                            
##  [8] tximport_1.4.0                         
##  [9] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
## [10] GenomicFeatures_1.28.0                 
## [11] AnnotationDbi_1.38.0                   
## [12] Biobase_2.36.0                         
## [13] GenomicRanges_1.28.0                   
## [14] GenomeInfoDb_1.12.0                    
## [15] IRanges_2.10.0                         
## [16] S4Vectors_0.14.0                       
## [17] BiocGenerics_0.22.0                    
## [18] tximportData_1.3.2                     
## [19] knitr_1.15.1                           
## 
## loaded via a namespace (and not attached):
##  [1] splines_3.4.0            Formula_1.2-1           
##  [3] latticeExtra_0.6-28      GenomeInfoDbData_0.99.0 
##  [5] Rsamtools_1.28.0         yaml_2.1.14             
##  [7] RSQLite_1.1-2            backports_1.0.5         
##  [9] lattice_0.20-35          digest_0.6.12           
## [11] RColorBrewer_1.1-2       XVector_0.16.0          
## [13] checkmate_1.8.2          colorspace_1.3-2        
## [15] htmltools_0.3.5          Matrix_1.2-9            
## [17] plyr_1.8.4               XML_3.98-1.6            
## [19] biomaRt_2.32.0           genefilter_1.58.0       
## [21] zlibbioc_1.22.0          xtable_1.8-2            
## [23] scales_0.4.1             BiocParallel_1.10.0     
## [25] annotate_1.54.0          tibble_1.3.0            
## [27] htmlTable_1.9            ggplot2_2.2.1           
## [29] nnet_7.3-12              lazyeval_0.2.0          
## [31] survival_2.41-3          magrittr_1.5            
## [33] memoise_1.1.0            evaluate_0.10           
## [35] foreign_0.8-67           data.table_1.10.4       
## [37] tools_3.4.0              hms_0.3                 
## [39] formatR_1.4              stringr_1.2.0           
## [41] munsell_0.4.3            locfit_1.5-9.1          
## [43] cluster_2.0.6            Biostrings_2.44.0       
## [45] compiler_3.4.0           rhdf5_2.20.0            
## [47] grid_3.4.0               RCurl_1.95-4.8          
## [49] htmlwidgets_0.8          rjson_0.2.15            
## [51] bitops_1.0-6             base64enc_0.1-3         
## [53] rmarkdown_1.4            gtable_0.2.0            
## [55] DBI_0.6-1                R6_2.2.0                
## [57] GenomicAlignments_1.12.0 gridExtra_2.2.1         
## [59] rtracklayer_1.36.0       Hmisc_4.0-2             
## [61] rprojroot_1.2            stringi_1.1.5           
## [63] Rcpp_0.12.10             geneplotter_1.54.0      
## [65] rpart_4.1-11             acepack_1.4.1

References

Bray, Nicolas, Harold Pimentel, Pall Melsted, and Lior Pachter. 2016. “Near-Optimal Probabilistic Rna-Seq Quantification.” Nature Biotechnology 34: 525–27. http://dx.doi.org/10.1038/nbt.3519.

Law, Charity W., Yunshun Chen, Wei Shi, and Gordon K. Smyth. 2014. “voom: precision weights unlock linear model analysis tools for RNA-seq read counts.” Genome Biology 15 (2): 29. http://dx.doi.org/10.1186/gb-2014-15-2-r29.

Li, Bo, and Colin N. Dewey. 2011. “RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome.” BMC Bioinformatics 12: 323+. doi:10.1186/1471-2105-12-3231.

Love, Michael I., Wolfgang Huber, and Simon Anders. 2014. “Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.” Genome Biology 15 (12): 550. http://dx.doi.org/10.1186/s13059-014-0550-8.

Patro, Rob, Geet Duggal, Michael I. Love, Rafael A. Irizarry, and Carl Kingsford. 2016. “Salmon Provides Accurate, Fast, and Bias-Aware Transcript Expression Estimates Using Dual-Phase Inference.” BioRxiv. http://biorxiv.org/content/early/2016/08/30/021592.

Patro, Rob, Stephen M. Mount, and Carl Kingsford. 2014. “Sailfish enables alignment-free isoform quantification from RNA-seq reads using lightweight algorithms.” Nature Biotechnology 32: 462–64. http://dx.doi.org/10.1038/nbt.2862.

Robert, Christelle, and Mick Watson. 2015. “Errors in RNA-Seq quantification affect genes of relevance to human disease.” Genome Biology. doi:10.1186/s13059-015-0734-x.

Robinson, Mark D., Davis J. McCarthy, and Gordon K. Smyth. 2010. “edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.” Bioinformatics 26 (1): 139. http://dx.doi.org/10.1093/bioinformatics/btp616.

Soneson, Charlotte, Michael I. Love, and Mark Robinson. 2015. “Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences.” F1000Research 4 (1521). http://dx.doi.org/10.12688/f1000research.7563.1.

Trapnell, Cole, David G Hendrickson, Martin Sauvageau, Loyal Goff, John L Rinn, and Lior Pachter. 2013. “Differential analysis of gene regulation at transcript resolution with RNA-seq.” Nature Biotechnology. doi:10.1038/nbt.2450.