Introduction

The tximeta package (Love et al. 2019) extends the tximport package (Soneson, Love, and Robinson 2015) for import of transcript-level quantification data into R/Bioconductor. It facilitates addition of annotation data when the RNA-seq data has been quantified with Salmon (Patro et al. 2017) or for scRNA-seq data quantified with Alevin (Srivastava et al. 2019). For more details on these packages – including the motivation for tximeta and description of similar work – consult the References below.

Note: tximeta requires that the entire output directory of Salmon/Alevin is present and unmodified in order to identify the provenance of the reference transcripts. In general, it’s a good idea to not modify or re-arrange the output directory of bioinformatic software as other downstream software rely on and assume a consistent directory structure. For sharing multiple samples, one can use, for example, tar -czf to bundle up a set of Salmon output directories, or to bundle one Alevin output directory.

Analysis starts with sample table

The first step using tximeta is to read in the sample table, which will become the column data, colData, of the final object, a SummarizedExperiment. The sample table should contain all the information we need to identify the Salmon quantification directories. For Alevin quantification, one should point to the quants_mat.gz file that contains the counts for all of the cells.

Here we will use a Salmon quantification file in the tximportData package to demonstrate the usage of tximeta. We do not have a sample table, so we construct one in R. It is recommended to keep a sample table as a CSV or TSV file while working on an RNA-seq project with multiple samples.

## [1] TRUE
##                                                                                        files
## 1 /home/biocbuild/bbs-3.10-bioc/R/library/tximportData/extdata/salmon_dm/SRR1197474/quant.sf
##        names condition
## 1 SRR1197474         A

tximeta expects at least two columns in coldata:

  1. files - a pointer to the quant.sf files
  2. names - the unique names that should be used to identify samples

Running tximeta from a sample table

Normally, we would just run tximeta like so:

However, to avoid downloading remote GTF files during this vignette, we will point to a GTF file saved locally (in the tximportData package). We link the transcriptome of the Salmon index to its locally saved GTF. The standard recommended usage of tximeta would be the code chunk above, or to specify a remote GTF source, not a local one. This following code is therefore not recommended for a typically workflow, but is particular to the vignette code.

## saving linkedTxome in bfc (first time)
## importing quantifications
## reading in files with read_tsv
## 1 
## found matching linked transcriptome:
## [ Ensembl - Drosophila melanogaster - release 98 ]
## building EnsDb with 'ensembldb' package
## Importing GTF file ... OK
## Processing metadata ... OK
## Processing genes ... 
##  Attribute availability:
##   o gene_id ... OK
##   o gene_name ... OK
##   o entrezid ... Nope
##   o gene_biotype ... OK
## OK
## Processing transcripts ... 
##  Attribute availability:
##   o transcript_id ... OK
##   o gene_id ... OK
##   o transcript_biotype ... OK
## OK
## Processing exons ... OK
## Processing chromosomes ... Fetch seqlengths from ensembl ... OK
## Generating index ... OK
##   -------------
## Verifying validity of the information in the database:
## Checking transcripts ... OK
## Checking exons ... OK
## generating transcript ranges
## Warning in checkAssays2Txps(assays, txps): 
## 
## Warning: the annotation is missing some transcripts that were quantified.
## 5 out of 33706 txps were missing from GTF/GFF but were in the indexed FASTA.
## (This occurs sometimes with Ensembl txps on haplotype chromosomes.)
## In order to build a ranged SummarizedExperiment, these txps were removed.
## To keep these txps, and to skip adding ranges, use skipMeta=TRUE
## 
## Example missing txps: [FBtr0307759, FBtr0084079, FBtr0084080, ...]

What happened?

tximeta recognized the hashed checksum of the transcriptome that the files were quantified against, it accessed the GTF file of the transcriptome source, found and attached the transcript ranges, and added the appropriate transcriptome and genome metadata. A remote GTF is only downloaded once, and a local or remote GTF is only parsed to build a TxDb once: if tximeta recognizes that it has seen this Salmon index before, it will use a cached version of the metadata and transcript ranges.

Note the warning above that 5 of the transcripts are missing from the GTF file and so are dropped from the final output. This is a problem coming from the annotation source, and not easily avoided by tximeta.

Pre-computed transcriptome checksums

We plan to support a wide variety of sources and organisms for transcriptomes with pre-computed checksums, though for now the software focuses on predominantly human and mouse transcriptomes (see Next steps below for details). The following checksums are supported in this version of tximeta:

source organism releases
GENCODE Homo sapiens 23-33
GENCODE Mus musculus M6-M24
Ensembl Homo sapiens 76-99
Ensembl Mus musculus 76-99
Ensembl Drosophila melanogaster 79-99
RefSeq Homo sapiens p1-p12
RefSeq Mus musculus p2-p5

For Ensembl transcriptomes, we support the combined protein coding (cDNA) and non-coding (ncRNA) sequences, as well as the protein coding alone (although the former approach combining coding and non-coding transcripts is recommended for more accurate quantification).

tximeta also has functions to support linked transcriptomes, where one or more sources for transcript sequences have been combined or filtered. See the Linked transcriptome section below for a demonstration. (The makeLinkedTxome function was used above to avoid downloading the GTF during the vignette building process.)

SummarizedExperiment output

We, of course, have our coldata from before. Note that we’ve removed files.

## DataFrame with 1 row and 2 columns
##                  names   condition
##            <character> <character>
## SRR1197474  SRR1197474           A

Here we show the three matrices that were imported.

## [1] "counts"    "abundance" "length"

If there were inferential replicates (Gibbs samples or bootstrap samples), these would be imported as additional assays named "infRep1", "infRep2", …

tximeta has imported the correct ranges for the transcripts:

## GRanges object with 33701 ranges and 6 metadata columns:
##               seqnames            ranges strand |       tx_id     tx_biotype
##                  <Rle>         <IRanges>  <Rle> | <character>    <character>
##   FBtr0070129        X     656673-657899      + | FBtr0070129 protein_coding
##   FBtr0070126        X     656356-657899      + | FBtr0070126 protein_coding
##   FBtr0070128        X     656673-657899      + | FBtr0070128 protein_coding
##   FBtr0070124        X     656114-657899      + | FBtr0070124 protein_coding
##   FBtr0070127        X     656356-657899      + | FBtr0070127 protein_coding
##           ...      ...               ...    ... .         ...            ...
##   FBtr0114299       2R 21325218-21325323      + | FBtr0114299         snoRNA
##   FBtr0113582       3R   5598638-5598777      - | FBtr0113582         snoRNA
##   FBtr0091635       3L   1488906-1489045      + | FBtr0091635         snoRNA
##   FBtr0113599       3L     261803-261953      - | FBtr0113599         snoRNA
##   FBtr0113600       3L     831870-832008      - | FBtr0113600         snoRNA
##               tx_cds_seq_start tx_cds_seq_end     gene_id     tx_name
##                      <integer>      <integer> <character> <character>
##   FBtr0070129           657110         657595 FBgn0025637 FBtr0070129
##   FBtr0070126           657110         657595 FBgn0025637 FBtr0070126
##   FBtr0070128           657110         657595 FBgn0025637 FBtr0070128
##   FBtr0070124           657110         657595 FBgn0025637 FBtr0070124
##   FBtr0070127           657110         657595 FBgn0025637 FBtr0070127
##           ...              ...            ...         ...         ...
##   FBtr0114299             <NA>           <NA> FBgn0086023 FBtr0114299
##   FBtr0113582             <NA>           <NA> FBgn0082989 FBtr0113582
##   FBtr0091635             <NA>           <NA> FBgn0086670 FBtr0091635
##   FBtr0113599             <NA>           <NA> FBgn0083014 FBtr0113599
##   FBtr0113600             <NA>           <NA> FBgn0083057 FBtr0113600
##   -------
##   seqinfo: 25 sequences from BDGP6.22 genome

We have appropriate genome information, which prevents us from making bioinformatic mistakes:

## Seqinfo object with 25 sequences from BDGP6.22 genome:
##   seqnames                        seqlengths isCircular   genome
##   211000022278279                      12714       <NA> BDGP6.22
##   211000022278436                       2815       <NA> BDGP6.22
##   211000022278449                       1947       <NA> BDGP6.22
##   211000022278760                       1144       <NA> BDGP6.22
##   211000022279165                       1118       <NA> BDGP6.22
##   ...                                    ...        ...      ...
##   Unmapped_Scaffold_8_D1580_D1567      88768       <NA> BDGP6.22
##   X                                 23542271       <NA> BDGP6.22
##   Y                                  3667352       <NA> BDGP6.22
##   mitochondrion_genome                 19524       <NA> BDGP6.22
##   rDNA                                 76973       <NA> BDGP6.22

Add exons per transcript

Because the SummarizedExperiment maintains all the metadata of its creation, it also keeps a pointer to the necessary database for pulling out additional information, as demonstrated in the following sections.

If necessary, the tximeta package can pull down the remote source to build a TxDb, but given that we’ve already built a TxDb once, it simply loads the cached version. In order to remove the cached TxDb and regenerate, one can remove the relevant entry from the tximeta file cache that resides at the location given by getTximetaBFC().

The se object created by tximeta, has the start, end, and strand information for each transcript. Here, we swap out the transcript GRanges for exons-by-transcript GRangesList (it is a list of GRanges, where each element of the list gives the exons for a particular transcript).

As with the transcript ranges, the exon ranges will be generated once and cached locally. As it takes a non-negligible amount of time to generate the exon-by-transcript GRangesList, this local caching offers substantial time savings for repeated usage of addExons with the same transcriptome.

We have implemented addExons to work only on the transcript-level SummarizedExperiment object. We provide some motivation for this choice in ?addExons. Briefly, if it is desired to know the exons associated with a particular gene, we feel that it makes more sense to pull out the relevant set of exons-by-transcript for the transcripts for this gene, rather than losing the hierarchical structure (exons to transcripts to genes) that would occur with a GRangesList of exons grouped per gene.

Easy summarization to gene-level

Likewise, the tximeta package can make use of the cached TxDb database for the purpose of summarizing transcript-level quantifications and bias corrections to the gene-level. After summarization, the rowRanges reflect the start and end position of the gene, which in Bioconductor are defined by the left-most and right-most genomic coordinates of all the transcripts. As with the transcript and exons, the gene ranges are cached locally for repeated usage.

## loading existing EnsDb created: 2020-03-12 02:59:55
## obtaining transcript-to-gene mapping from TxDb
## generating gene ranges
## summarizing abundance
## summarizing counts
## summarizing length
## GRanges object with 17208 ranges and 6 metadata columns:
##               seqnames            ranges strand |     gene_id      gene_name
##                  <Rle>         <IRanges>  <Rle> | <character>    <character>
##   FBgn0000003       3R   6822498-6822796      + | FBgn0000003 7SLRNA:CR32864
##   FBgn0000008       2R 22136968-22172834      + | FBgn0000008              a
##   FBgn0000014       3R 16807214-16830049      - | FBgn0000014          abd-A
##   FBgn0000015       3R 16927212-16972236      - | FBgn0000015          Abd-B
##   FBgn0000017       3L 16615866-16647882      - | FBgn0000017            Abl
##           ...      ...               ...    ... .         ...            ...
##   FBgn0286199       3R 24279572-24281576      + | FBgn0286199           shps
##   FBgn0286203       2R   5413744-5456095      + | FBgn0286203            stw
##   FBgn0286204       3R   8950246-8963037      - | FBgn0286204            ich
##   FBgn0286213       3L 13023352-13024762      + | FBgn0286213          RpS12
##   FBgn0286222        X   6678424-6681845      + | FBgn0286222           Fum1
##                entrezid   gene_biotype seq_coord_system         symbol
##               <integer>    <character>        <integer>    <character>
##   FBgn0000003      <NA>          ncRNA             <NA> 7SLRNA:CR32864
##   FBgn0000008      <NA> protein_coding             <NA>              a
##   FBgn0000014      <NA> protein_coding             <NA>          abd-A
##   FBgn0000015      <NA> protein_coding             <NA>          Abd-B
##   FBgn0000017      <NA> protein_coding             <NA>            Abl
##           ...       ...            ...              ...            ...
##   FBgn0286199      <NA> protein_coding             <NA>           shps
##   FBgn0286203      <NA> protein_coding             <NA>            stw
##   FBgn0286204      <NA> protein_coding             <NA>            ich
##   FBgn0286213      <NA> protein_coding             <NA>          RpS12
##   FBgn0286222      <NA> protein_coding             <NA>           Fum1
##   -------
##   seqinfo: 25 sequences from BDGP6.22 genome

Add different identifiers

We would like to add support to easily map transcript or gene identifiers from one annotation to another. This is just a prototype function, but we show how we can easily add alternate IDs given that we know the organism and the source of the transcriptome. (This function currently only works for GENCODE and Ensembl gene or transcript IDs but could be extended to work for arbitrary sources.)

## Loading required package: AnnotationDbi
## 
## mapping to new IDs using 'org.Dm.eg.db' data package
## if all matching IDs are desired, and '1:many mappings' are reported,
## set multiVals='list' to obtain all the matching IDs
## 'select()' returned 1:many mapping between keys and columns
## DataFrame with 17208 rows and 7 columns
##                 gene_id      gene_name  entrezid   gene_biotype
##             <character>    <character> <integer>    <character>
## FBgn0000003 FBgn0000003 7SLRNA:CR32864        NA          ncRNA
## FBgn0000008 FBgn0000008              a        NA protein_coding
## FBgn0000014 FBgn0000014          abd-A        NA protein_coding
## FBgn0000015 FBgn0000015          Abd-B        NA protein_coding
## FBgn0000017 FBgn0000017            Abl        NA protein_coding
## ...                 ...            ...       ...            ...
## FBgn0286199 FBgn0286199           shps        NA protein_coding
## FBgn0286203 FBgn0286203            stw        NA protein_coding
## FBgn0286204 FBgn0286204            ich        NA protein_coding
## FBgn0286213 FBgn0286213          RpS12        NA protein_coding
## FBgn0286222 FBgn0286222           Fum1        NA protein_coding
##             seq_coord_system         symbol       REFSEQ
##                    <integer>    <character>  <character>
## FBgn0000003               NA 7SLRNA:CR32864           NA
## FBgn0000008               NA              a NM_001014543
## FBgn0000014               NA          abd-A NM_001170161
## FBgn0000015               NA          Abd-B NM_001275719
## FBgn0000017               NA            Abl NM_001104153
## ...                      ...            ...          ...
## FBgn0286199               NA           shps    NM_142982
## FBgn0286203               NA            stw NM_001144134
## FBgn0286204               NA            ich NM_001275464
## FBgn0286213               NA          RpS12    NM_168534
## FBgn0286222               NA           Fum1    NM_132111

Run a differential expression analysis

The following code chunk demonstrates how to build a DESeqDataSet and begin a differential expression analysis.

## using counts and average transcript lengths from tximeta

The following code chunk demonstrates how to build a DGEList object, for use with edgeR.

The following code chunk demonstrates how one could use the Swish method in the fishpond Bioconductor package. Here we use the transcript-level object se. This dataset only has a single sample and no inferential replicates, but the analysis would begin with such code. See the Swish vignette in the fishpond package for a complete example:

For limma with voom transformation we recommend, as in the tximport vignette to generate counts-from-abundance instead of providing an offset for average transcript length.

## loading existing EnsDb created: 2020-03-12 02:59:55
## obtaining transcript-to-gene mapping from TxDb
## loading existing gene ranges created: 2020-03-12 03:00:17
## summarizing abundance
## summarizing counts
## summarizing length

Above we generated counts-from-abundance when calling summarizeToGene. The counts-from-abundance status is then stored in the metadata:

## [1] "lengthScaledTPM"

Additional metadata

The following information is attached to the SummarizedExperiment by tximeta:

## [1] "tximetaInfo"         "quantInfo"           "countsFromAbundance"
## [4] "level"               "txomeInfo"           "txdbInfo"
## List of 31
##  $ salmon_version                                        : chr "0.14.1"
##  $ samp_type                                             : chr "none"
##  $ opt_type                                              : chr "vb"
##  $ quant_errors                                          :List of 1
##   ..$ : list()
##  $ num_libraries                                         : int 1
##  $ library_types                                         : chr "ISR"
##  $ frag_dist_length                                      : int 1001
##  $ seq_bias_correct                                      : logi FALSE
##  $ gc_bias_correct                                       : logi TRUE
##  $ num_bias_bins                                         : int 4096
##  $ mapping_type                                          : chr "mapping"
##  $ num_valid_targets                                     : int 33706
##  $ num_decoy_targets                                     : int 0
##  $ num_eq_classes                                        : int 70718
##  $ serialized_eq_classes                                 : logi FALSE
##  $ length_classes                                        : int [1:5, 1] 867 1533 2379 3854 71382
##  $ index_seq_hash                                        : chr "7ba5e9597796ea86cf11ccf6635ca88fbc37c2848d38083c23986aa2c6a21eae"
##  $ index_name_hash                                       : chr "b6426061057bba9b7afb4dc76fa68238414cf35b4190c95ca6fc44280d4ca87c"
##  $ index_seq_hash512                                     : chr "05f111abcda1efd2e489ace6324128cdaaa311712a28ed716d957fdfd8706ec41ca9177ebf12f54e99c2a89582d06f31c5e09dc1dce2d13"| __truncated__
##  $ index_name_hash512                                    : chr "ccdf58f23e48c8c53cd122b5f5990b5adce9fec87ddf8bd88153afbe93296d87b818fba89d12dbc20c882f7d98353840394c5040fea7432"| __truncated__
##  $ num_bootstraps                                        : int 0
##  $ num_processed                                         : int 42422337
##  $ num_mapped                                            : int 34098209
##  $ num_decoy_fragments                                   : int 0
##  $ num_dovetail_fragments                                : int 2048810
##  $ num_fragments_filtered_vm                             : int 989383
##  $ num_alignments_below_threshold_for_mapped_fragments_vm: int 267540
##  $ percent_mapped                                        : num 80.4
##  $ call                                                  : chr "quant"
##  $ start_time                                            : chr "Sat Oct 12 13:55:01 2019"
##  $ end_time                                              : chr "Sat Oct 12 14:08:11 2019"
## List of 8
##  $ index   : chr "Dm.BDGP6.22.98_salmon-0.14.1"
##  $ source  : chr "Ensembl"
##  $ organism: chr "Drosophila melanogaster"
##  $ release : chr "98"
##  $ genome  : chr "BDGP6.22"
##  $ fasta   :List of 1
##   ..$ : chr [1:2] "ftp://ftp.ensembl.org/pub/release-98/fasta/drosophila_melanogaster/cdna/Drosophila_melanogaster.BDGP6.22.cdna.all.fa.gz" "ftp://ftp.ensembl.org/pub/release-98/fasta/drosophila_melanogaster/ncrna/Drosophila_melanogaster.BDGP6.22.ncrna.fa.gz"
##  $ gtf     : chr "/home/biocbuild/bbs-3.10-bioc/R/library/tximportData/extdata/salmon_dm/Drosophila_melanogaster.BDGP6.22.98.gtf.gz"
##  $ sha256  : chr "7ba5e9597796ea86cf11ccf6635ca88fbc37c2848d38083c23986aa2c6a21eae"
## $version
## [1] '1.4.5'
## 
## $importTime
## [1] "2020-03-11 22:59:54 EDT"
##                                      Db type 
##                                      "EnsDb" 
##                              Type of Gene ID 
##                            "Ensembl Gene ID" 
##                           Supporting package 
##                                  "ensembldb" 
##                                Db created by 
##        "ensembldb package from Bioconductor" 
##                               script_version 
##                                      "0.0.1" 
##                                Creation time 
##                   "Wed Mar 11 23:00:03 2020" 
##                              ensembl_version 
##                                         "98" 
##                                 ensembl_host 
##                                    "unknown" 
##                                     Organism 
##                    "Drosophila_melanogaster" 
##                                 genome_build 
##                                   "BDGP6.22" 
##                              DBSCHEMAVERSION 
##                                        "1.0" 
##                                  source_file 
## "Drosophila_melanogaster.BDGP6.22.98.gtf.gz"

Transcriptomes not known to tximeta

tximeta automatically imports relevant metadata when the transcriptome matches a known source – known in the sense that it is in the set of pre-computed hashed checksums in tximeta (GENCODE, Ensembl, and RefSeq for human and mouse). tximeta also facilitates the linking of transcriptomes used in building the Salmon index with relevant public sources, in the case that these are not part of this pre-computed set known to tximeta. The linking of the transcriptome source with the quantification files is important in the case that the transcript sequence no longer matches a known source (uniquely combined or filtered FASTA files), or if the source is not known to tximeta. Combinations of coding and non-coding human, mouse, and fruit fly Ensembl transcripts should be automatically recognized by tximeta and does not require making a linkedTxome. As the package is further developed, we plan to roll out support for all common transcriptomes, from all sources. Below we demonstrate how to make a linkedTxome and how to share and load a linkedTxome.

We point to a Salmon quantification file which was quantified against a transcriptome that included the coding and non-coding Drosophila melanogaster transcripts, as well as an artificial transcript of 960 bp (for demonstration purposes only).

## [1] TRUE

Trying to import the files gives a message that tximeta couldn’t find a matching transcriptome, so it returns an non-ranged SummarizedExperiment.

## importing quantifications
## reading in files with read_tsv
## 1 
## couldn't find matching transcriptome, returning non-ranged SummarizedExperiment

Linked transcriptomes

If the transcriptome used to generate the Salmon index does not match any transcriptomes from known sources (e.g. from combining or filtering known transcriptome files), there is not much that can be done to automatically populate the metadata during quantification import. However, we can facilitate the following two cases:

  1. the transcriptome was created locally and has been linked to its public source(s)
  2. the transcriptome was produced by another group, and they have produced and shared a file that links the transcriptome to public source(s)

tximeta offers functionality to assist reproducible analysis in both of these cases.

To make this quantification reproducible, we make a linkedTxome which records key information about the sources of the transcript FASTA files, and the location of the relevant GTF file. It also records the checksum of the transcriptome that was computed by Salmon during the index step.

Multiple GTF/GFF files: linkedTxome and tximeta do not currently support multiple GTF/GFF files, which is a more complicated case than multiple FASTA, which is supported. Currently, we recommend that users should add or combine GTF/GFF files themselves to create a single GTF/GFF file that contains all features used in quantification, and then upload such a file to Zenodo, which can then be linked as shown below. Feel free to contact the developers on the Bioconductor support site or GitHub Issue page for further details or feature requests.

By default, linkedTxome will write out a JSON file which can be shared with others, linking the checksum of the index with the other metadata, including FASTA and GTF sources. By default, it will write out to a file with the same name as the indexDir, but with a .json extension added. This can be prevented with write=FALSE, and the file location can be changed with jsonFile.

First we specify the path where the Salmon index is located.

Typically you would not use system.file and file.path to locate this directory, but simply define indexDir to be the path of the Salmon directory on your machine. Here we use system.file and file.path because we have included parts of a Salmon index directory in the tximeta package itself for demonstration of functionality in this vignette.

Now we provide the location of the FASTA files and the GTF file for this transcriptome.

Note: the basename for the GTF file is used as a unique identifier for the cached versions of the TxDb and the transcript ranges, which are stored on the user’s behalf via BiocFileCache. This is not an issue, as GENCODE, Ensembl, and RefSeq all provide GTF files which are uniquely identified by their filename, e.g. Drosophila_melanogaster.BDGP6.22.98.gtf.gz.

The recommended usage of tximeta would be to specify a remote GTF source, as seen in the commented-out line below:

Instead of the above commented-out FTP location for the GTF file, we specify a location within an R package. This step is just to avoid downloading from a remote FTP during vignette building. This use of file.path to point to a file in an R package is specific to this vignette and should not be used in a typical workflow. The following GTF file is a modified version of the release 98 from Ensembl, which includes description of a one transcript, one exon artificial gene which was inserted into the transcriptome (for demonstration purposes only).

Finally, we create a linkedTxome. In this vignette, we point to a temporary directory for the JSON file, but a more typical workflow would write the JSON file to the same location as the Salmon index by not specifying jsonFile.

makeLinkedTxome performs two operation: (1) it creates a new entry in an internal table that links the transcriptome used in the Salmon index to its sources, and (2) it creates a JSON file such that this linkedTxome can be shared.

## writing linkedTxome to /tmp/Rtmp11ig3j/Dm.BDGP6.22.98.plus_salmon-0.14.1.json
## saving linkedTxome in bfc

After running makeLinkedTxome, the connection between this Salmon index (and its checksum) with the sources is saved for persistent usage. Note that because we added a single transcript of 960bp to the FASTA file used for quantification, tximeta could tell that this was not quantified against release 98 of the Ensembl transcripts for Drosophila melanogaster. Only when the correct set of transcripts were specified does tximeta recognize and import the correct metadata.

With use of tximeta and a linkedTxome, the software figures out if the remote GTF has been accessed and compiled into a TxDb before, and on future calls, it will simply load the pre-computed metadata and transcript ranges.

Note the warning that 5 of the transcripts are missing from the GTF file and so are dropped from the final output. This is a problem coming from the annotation source, and not easily avoided by tximeta.

## importing quantifications
## reading in files with read_tsv
## 1 
## found matching linked transcriptome:
## [ Ensembl - Drosophila melanogaster - release 98 ]
## building EnsDb with 'ensembldb' package
## Importing GTF file ... OK
## Processing metadata ... OK
## Processing genes ... 
##  Attribute availability:
##   o gene_id ... OK
##   o gene_name ... OK
##   o entrezid ... Nope
##   o gene_biotype ... OK
## OK
## Processing transcripts ... 
##  Attribute availability:
##   o transcript_id ... OK
##   o gene_id ... OK
##   o transcript_biotype ... OK
## OK
## Processing exons ... OK
## Processing chromosomes ... Fetch seqlengths from ensembl ... OK
## Generating index ... OK
##   -------------
## Verifying validity of the information in the database:
## Checking transcripts ... OK
## Checking exons ... OK
## generating transcript ranges
## Warning in checkAssays2Txps(assays, txps): 
## 
## Warning: the annotation is missing some transcripts that were quantified.
## 5 out of 33707 txps were missing from GTF/GFF but were in the indexed FASTA.
## (This occurs sometimes with Ensembl txps on haplotype chromosomes.)
## In order to build a ranged SummarizedExperiment, these txps were removed.
## To keep these txps, and to skip adding ranges, use skipMeta=TRUE
## 
## Example missing txps: [FBtr0307759, FBtr0084079, FBtr0084080, ...]

We can see that the appropriate metadata and transcript ranges are attached.

## GRanges object with 33702 ranges and 6 metadata columns:
##               seqnames            ranges strand |       tx_id     tx_biotype
##                  <Rle>         <IRanges>  <Rle> | <character>    <character>
##       Newgene       3R             1-960      + |     Newgene protein_coding
##   FBtr0070129        X     656673-657899      + | FBtr0070129 protein_coding
##   FBtr0070126        X     656356-657899      + | FBtr0070126 protein_coding
##   FBtr0070128        X     656673-657899      + | FBtr0070128 protein_coding
##   FBtr0070124        X     656114-657899      + | FBtr0070124 protein_coding
##           ...      ...               ...    ... .         ...            ...
##   FBtr0114299       2R 21325218-21325323      + | FBtr0114299         snoRNA
##   FBtr0113582       3R   5598638-5598777      - | FBtr0113582         snoRNA
##   FBtr0091635       3L   1488906-1489045      + | FBtr0091635         snoRNA
##   FBtr0113599       3L     261803-261953      - | FBtr0113599         snoRNA
##   FBtr0113600       3L     831870-832008      - | FBtr0113600         snoRNA
##               tx_cds_seq_start tx_cds_seq_end     gene_id     tx_name
##                      <integer>      <integer> <character> <character>
##       Newgene             <NA>           <NA>     Newgene     Newgene
##   FBtr0070129           657110         657595 FBgn0025637 FBtr0070129
##   FBtr0070126           657110         657595 FBgn0025637 FBtr0070126
##   FBtr0070128           657110         657595 FBgn0025637 FBtr0070128
##   FBtr0070124           657110         657595 FBgn0025637 FBtr0070124
##           ...              ...            ...         ...         ...
##   FBtr0114299             <NA>           <NA> FBgn0086023 FBtr0114299
##   FBtr0113582             <NA>           <NA> FBgn0082989 FBtr0113582
##   FBtr0091635             <NA>           <NA> FBgn0086670 FBtr0091635
##   FBtr0113599             <NA>           <NA> FBgn0083014 FBtr0113599
##   FBtr0113600             <NA>           <NA> FBgn0083057 FBtr0113600
##   -------
##   seqinfo: 25 sequences from BDGP6.22 genome
## Seqinfo object with 25 sequences from BDGP6.22 genome:
##   seqnames                        seqlengths isCircular   genome
##   211000022278279                      12714       <NA> BDGP6.22
##   211000022278436                       2815       <NA> BDGP6.22
##   211000022278449                       1947       <NA> BDGP6.22
##   211000022278760                       1144       <NA> BDGP6.22
##   211000022279165                       1118       <NA> BDGP6.22
##   ...                                    ...        ...      ...
##   Unmapped_Scaffold_8_D1580_D1567      88768       <NA> BDGP6.22
##   X                                 23542271       <NA> BDGP6.22
##   Y                                  3667352       <NA> BDGP6.22
##   mitochondrion_genome                 19524       <NA> BDGP6.22
##   rDNA                                 76973       <NA> BDGP6.22

Clear linkedTxomes

The following code removes the entire table with information about the linkedTxomes. This is just for demonstration, so that we can show how to load a JSON file below.

Note: Running this code will clear any information about linkedTxomes. Don’t run this unless you really want to clear this table!

## Loading required package: dbplyr
## # A tibble: 6 x 10
##   rid   rname create_time access_time rpath rtype fpath last_modified_t… etag 
##   <chr> <chr> <chr>       <chr>       <chr> <chr> <chr>            <dbl> <chr>
## 1 BFC1  link… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## 2 BFC2  Dros… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## 3 BFC3  txpR… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## 4 BFC4  gene… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## 5 BFC5  Dros… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## 6 BFC6  txpR… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## # … with 1 more variable: expires <dbl>
## # A tibble: 5 x 10
##   rid   rname create_time access_time rpath rtype fpath last_modified_t… etag 
##   <chr> <chr> <chr>       <chr>       <chr> <chr> <chr>            <dbl> <chr>
## 1 BFC2  Dros… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## 2 BFC3  txpR… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## 3 BFC4  gene… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## 4 BFC5  Dros… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## 5 BFC6  txpR… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## # … with 1 more variable: expires <dbl>

Loading linkedTxome JSON files

If a collaborator or the Suppmentary Files for a publication shares a linkedTxome JSON file, we can likewise use tximeta to automatically assemble the relevant metadata and transcript ranges. This implies that the other person has used tximeta with the function makeLinkedTxome demonstrated above, pointing to their Salmon index and to the FASTA and GTF source(s).

We point to the JSON file and use loadLinkedTxome and then the relevant metadata is saved for persistent usage. In this case, we saved the JSON file in a temporary directory.

## saving linkedTxome in bfc (first time)

Again, using tximeta figures out whether it needs to access the remote GTF or not, and assembles the appropriate object on the user’s behalf.

## importing quantifications
## reading in files with read_tsv
## 1 
## found matching linked transcriptome:
## [ Ensembl - Drosophila melanogaster - release 98 ]
## loading existing EnsDb created: 2020-03-12 03:00:23
## loading existing transcript ranges created: 2020-03-12 03:00:38
## Warning in checkAssays2Txps(assays, txps): 
## 
## Warning: the annotation is missing some transcripts that were quantified.
## 5 out of 33707 txps were missing from GTF/GFF but were in the indexed FASTA.
## (This occurs sometimes with Ensembl txps on haplotype chromosomes.)
## In order to build a ranged SummarizedExperiment, these txps were removed.
## To keep these txps, and to skip adding ranges, use skipMeta=TRUE
## 
## Example missing txps: [FBtr0307759, FBtr0084079, FBtr0084080, ...]

Clear linkedTxomes

Finally, we clear the linkedTxomes table again so that the above examples will work. This is just for the vignette code and not part of a typical workflow.

Note: Running this code will clear any information about linkedTxomes. Don’t run this unless you really want to clear this table!

## # A tibble: 6 x 10
##   rid   rname create_time access_time rpath rtype fpath last_modified_t… etag 
##   <chr> <chr> <chr>       <chr>       <chr> <chr> <chr>            <dbl> <chr>
## 1 BFC2  Dros… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## 2 BFC3  txpR… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## 3 BFC4  gene… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## 4 BFC5  Dros… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## 5 BFC6  txpR… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## 6 BFC7  link… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## # … with 1 more variable: expires <dbl>
## # A tibble: 5 x 10
##   rid   rname create_time access_time rpath rtype fpath last_modified_t… etag 
##   <chr> <chr> <chr>       <chr>       <chr> <chr> <chr>            <dbl> <chr>
## 1 BFC2  Dros… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## 2 BFC3  txpR… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## 3 BFC4  gene… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## 4 BFC5  Dros… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## 5 BFC6  txpR… 2020-03-12… 2020-03-12… /tmp… rela… 7301…               NA <NA> 
## # … with 1 more variable: expires <dbl>

Other quantifiers

tximeta can import the output from any quantifiers that are supported by tximport, and if these are not Salmon, Alevin, or Sailfish output, it will simply return a non-ranged SummarizedExperiment. We are working to allow manually passing of the hash value of the transcriptome, the cDNA sequences of which can be hashed with FastaDigest (can be installed with pip install fasta_digest).

Acknowledgments

The development of tximeta has benefited from suggestions from these and other individuals in the community:

  • Vincent Carey
  • Lori Shepherd
  • Martin Morgan
  • Koen Van den Berge

Next steps

Integration with GA4GH / refget API

  • We are collaborating and in communication with GA4GH working groups to build out functionality to perform lookup on collections of transcripts, i.e. transcriptomes, provided by Ensembl. This will greatly expand the applicability of tximeta. The current version of tximeta relies on hashing of common transcriptomes, which are stored within the package’s extdata directory in a CSV file, or on the use of linkedTxome for any additional transcriptomes. However, with the GA4GH / refget API in place, tximeta will be able to identify and access transcriptome metadata for vastly more reference transcriptomes (collections of transcripts).

Facilitate plots and summaries

  • Basic plots across samples: abundances, mapping rates, rich bias model parameters
  • Time summaries: when quantified? when imported?

Extended functionality

  • Facilitate functional annotation, either with vignettes/workflow or with additional functionality. E.g.: housekeeping genes, arbitrary gene sets, genes expressed in GTEx tissues
  • liftOver is clunky and doesn’t integrate with GenomeInfoDb. It requires user input and there’s a chance to mis-annotate. Ideally this should all be automated.

Session info

## Loading required package: usethis
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 3.6.3 (2020-02-29)
##  os       Ubuntu 18.04.4 LTS          
##  system   x86_64, linux-gnu           
##  ui       X11                         
##  language (EN)                        
##  collate  C                           
##  ctype    en_US.UTF-8                 
##  tz       America/New_York            
##  date     2020-03-11                  
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package              * version  date       lib source        
##  acepack                1.4.1    2016-10-29 [2] CRAN (R 3.6.3)
##  annotate               1.64.0   2020-03-11 [2] Bioconductor  
##  AnnotationDbi        * 1.48.0   2020-03-11 [2] Bioconductor  
##  AnnotationFilter       1.10.0   2020-03-11 [2] Bioconductor  
##  askpass                1.1      2019-01-13 [2] CRAN (R 3.6.3)
##  assertthat             0.2.1    2019-03-21 [2] CRAN (R 3.6.3)
##  backports              1.1.5    2019-10-02 [2] CRAN (R 3.6.3)
##  base64enc              0.1-3    2015-07-28 [2] CRAN (R 3.6.3)
##  Biobase              * 2.46.0   2020-03-11 [2] Bioconductor  
##  BiocFileCache        * 1.10.2   2020-03-11 [2] Bioconductor  
##  BiocGenerics         * 0.32.0   2020-03-11 [2] Bioconductor  
##  BiocParallel         * 1.20.1   2020-03-11 [2] Bioconductor  
##  biomaRt                2.42.0   2020-03-11 [2] Bioconductor  
##  Biostrings             2.54.0   2020-03-11 [2] Bioconductor  
##  bit                    1.1-15.2 2020-02-10 [2] CRAN (R 3.6.3)
##  bit64                  0.9-7    2017-05-08 [2] CRAN (R 3.6.3)
##  bitops                 1.0-6    2013-08-17 [2] CRAN (R 3.6.3)
##  blob                   1.2.1    2020-01-20 [2] CRAN (R 3.6.3)
##  callr                  3.4.2    2020-02-12 [2] CRAN (R 3.6.3)
##  checkmate              2.0.0    2020-02-06 [2] CRAN (R 3.6.3)
##  cli                    2.0.2    2020-02-28 [2] CRAN (R 3.6.3)
##  cluster                2.1.0    2019-06-19 [2] CRAN (R 3.6.3)
##  colorspace             1.4-1    2019-03-18 [2] CRAN (R 3.6.3)
##  crayon                 1.3.4    2017-09-16 [2] CRAN (R 3.6.3)
##  curl                   4.3      2019-12-02 [2] CRAN (R 3.6.3)
##  data.table             1.12.8   2019-12-09 [2] CRAN (R 3.6.3)
##  DBI                    1.1.0    2019-12-15 [2] CRAN (R 3.6.3)
##  dbplyr               * 1.4.2    2019-06-17 [2] CRAN (R 3.6.3)
##  DelayedArray         * 0.12.2   2020-03-11 [2] Bioconductor  
##  desc                   1.2.0    2018-05-01 [2] CRAN (R 3.6.3)
##  DESeq2               * 1.26.0   2020-03-11 [2] Bioconductor  
##  devtools             * 2.2.2    2020-02-17 [2] CRAN (R 3.6.3)
##  digest                 0.6.25   2020-02-23 [2] CRAN (R 3.6.3)
##  dplyr                  0.8.5    2020-03-07 [2] CRAN (R 3.6.3)
##  edgeR                * 3.28.1   2020-03-11 [2] Bioconductor  
##  ellipsis               0.3.0    2019-09-20 [2] CRAN (R 3.6.3)
##  ensembldb              2.10.2   2020-03-11 [2] Bioconductor  
##  evaluate               0.14     2019-05-28 [2] CRAN (R 3.6.3)
##  fansi                  0.4.1    2020-01-08 [2] CRAN (R 3.6.3)
##  foreign                0.8-76   2020-03-03 [2] CRAN (R 3.6.3)
##  Formula                1.2-3    2018-05-03 [2] CRAN (R 3.6.3)
##  fs                     1.3.2    2020-03-05 [2] CRAN (R 3.6.3)
##  genefilter             1.68.0   2020-03-11 [2] Bioconductor  
##  geneplotter            1.64.0   2020-03-11 [2] Bioconductor  
##  GenomeInfoDb         * 1.22.0   2020-03-11 [2] Bioconductor  
##  GenomeInfoDbData       1.2.2    2020-03-05 [2] Bioconductor  
##  GenomicAlignments      1.22.1   2020-03-11 [2] Bioconductor  
##  GenomicFeatures        1.38.2   2020-03-11 [2] Bioconductor  
##  GenomicRanges        * 1.38.0   2020-03-11 [2] Bioconductor  
##  ggplot2                3.3.0    2020-03-05 [2] CRAN (R 3.6.3)
##  glue                   1.3.1    2019-03-12 [2] CRAN (R 3.6.3)
##  gridExtra              2.3      2017-09-09 [2] CRAN (R 3.6.3)
##  gtable                 0.3.0    2019-03-25 [2] CRAN (R 3.6.3)
##  highr                  0.8      2019-03-20 [2] CRAN (R 3.6.3)
##  Hmisc                  4.3-1    2020-02-07 [2] CRAN (R 3.6.3)
##  hms                    0.5.3    2020-01-08 [2] CRAN (R 3.6.3)
##  htmlTable              1.13.3   2019-12-04 [2] CRAN (R 3.6.3)
##  htmltools              0.4.0    2019-10-04 [2] CRAN (R 3.6.3)
##  htmlwidgets            1.5.1    2019-10-08 [2] CRAN (R 3.6.3)
##  httr                   1.4.1    2019-08-05 [2] CRAN (R 3.6.3)
##  IRanges              * 2.20.2   2020-03-11 [2] Bioconductor  
##  jpeg                   0.1-8.1  2019-10-24 [2] CRAN (R 3.6.3)
##  jsonlite               1.6.1    2020-02-02 [2] CRAN (R 3.6.3)
##  knitr                  1.28     2020-02-06 [2] CRAN (R 3.6.3)
##  lattice                0.20-40  2020-02-19 [2] CRAN (R 3.6.3)
##  latticeExtra           0.6-29   2019-12-19 [2] CRAN (R 3.6.3)
##  lazyeval               0.2.2    2019-03-15 [2] CRAN (R 3.6.3)
##  lifecycle              0.2.0    2020-03-06 [2] CRAN (R 3.6.3)
##  limma                * 3.42.2   2020-03-11 [2] Bioconductor  
##  locfit                 1.5-9.1  2013-04-20 [2] CRAN (R 3.6.3)
##  magrittr               1.5      2014-11-22 [2] CRAN (R 3.6.3)
##  Matrix                 1.2-18   2019-11-27 [2] CRAN (R 3.6.3)
##  matrixStats          * 0.55.0   2019-09-07 [2] CRAN (R 3.6.3)
##  memoise                1.1.0    2017-04-21 [2] CRAN (R 3.6.3)
##  munsell                0.5.0    2018-06-12 [2] CRAN (R 3.6.3)
##  nnet                   7.3-13   2020-02-25 [2] CRAN (R 3.6.3)
##  openssl                1.4.1    2019-07-18 [2] CRAN (R 3.6.3)
##  org.Dm.eg.db         * 3.10.0   2020-03-05 [2] Bioconductor  
##  pillar                 1.4.3    2019-12-20 [2] CRAN (R 3.6.3)
##  pkgbuild               1.0.6    2019-10-09 [2] CRAN (R 3.6.3)
##  pkgconfig              2.0.3    2019-09-22 [2] CRAN (R 3.6.3)
##  pkgload                1.0.2    2018-10-29 [2] CRAN (R 3.6.3)
##  png                    0.1-7    2013-12-03 [2] CRAN (R 3.6.3)
##  prettyunits            1.1.1    2020-01-24 [2] CRAN (R 3.6.3)
##  processx               3.4.2    2020-02-09 [2] CRAN (R 3.6.3)
##  progress               1.2.2    2019-05-16 [2] CRAN (R 3.6.3)
##  ProtGenerics           1.18.0   2020-03-11 [2] Bioconductor  
##  ps                     1.3.2    2020-02-13 [2] CRAN (R 3.6.3)
##  purrr                  0.3.3    2019-10-18 [2] CRAN (R 3.6.3)
##  R6                     2.4.1    2019-11-12 [2] CRAN (R 3.6.3)
##  rappdirs               0.3.1    2016-03-28 [2] CRAN (R 3.6.3)
##  RColorBrewer           1.1-2    2014-12-07 [2] CRAN (R 3.6.3)
##  Rcpp                   1.0.3    2019-11-08 [2] CRAN (R 3.6.3)
##  RCurl                  1.98-1.1 2020-01-19 [2] CRAN (R 3.6.3)
##  readr                  1.3.1    2018-12-21 [2] CRAN (R 3.6.3)
##  remotes                2.1.1    2020-02-15 [2] CRAN (R 3.6.3)
##  rlang                  0.4.5    2020-03-01 [2] CRAN (R 3.6.3)
##  rmarkdown              2.1      2020-01-20 [2] CRAN (R 3.6.3)
##  rpart                  4.1-15   2019-04-12 [2] CRAN (R 3.6.3)
##  rprojroot              1.3-2    2018-01-03 [2] CRAN (R 3.6.3)
##  Rsamtools              2.2.3    2020-03-11 [2] Bioconductor  
##  RSQLite                2.2.0    2020-01-07 [2] CRAN (R 3.6.3)
##  rstudioapi             0.11     2020-02-07 [2] CRAN (R 3.6.3)
##  rtracklayer            1.46.0   2020-03-11 [2] Bioconductor  
##  S4Vectors            * 0.24.3   2020-03-11 [2] Bioconductor  
##  scales                 1.1.0    2019-11-18 [2] CRAN (R 3.6.3)
##  sessioninfo            1.1.1    2018-11-05 [2] CRAN (R 3.6.3)
##  stringi                1.4.6    2020-02-17 [2] CRAN (R 3.6.3)
##  stringr                1.4.0    2019-02-10 [2] CRAN (R 3.6.3)
##  SummarizedExperiment * 1.16.1   2020-03-11 [2] Bioconductor  
##  survival               3.1-11   2020-03-07 [2] CRAN (R 3.6.3)
##  testthat               2.3.2    2020-03-02 [2] CRAN (R 3.6.3)
##  tibble                 2.1.3    2019-06-06 [2] CRAN (R 3.6.3)
##  tidyselect             1.0.0    2020-01-27 [2] CRAN (R 3.6.3)
##  tximeta              * 1.4.5    2020-03-11 [1] Bioconductor  
##  tximport               1.14.0   2020-03-11 [2] Bioconductor  
##  usethis              * 1.5.1    2019-07-04 [2] CRAN (R 3.6.3)
##  utf8                   1.1.4    2018-05-24 [2] CRAN (R 3.6.3)
##  vctrs                  0.2.4    2020-03-10 [2] CRAN (R 3.6.3)
##  withr                  2.1.2    2018-03-15 [2] CRAN (R 3.6.3)
##  xfun                   0.12     2020-01-13 [2] CRAN (R 3.6.3)
##  XML                    3.99-0.3 2020-01-20 [2] CRAN (R 3.6.3)
##  xtable                 1.8-4    2019-04-21 [2] CRAN (R 3.6.3)
##  XVector                0.26.0   2020-03-11 [2] Bioconductor  
##  yaml                   2.2.1    2020-02-01 [2] CRAN (R 3.6.3)
##  zlibbioc               1.32.0   2020-03-11 [2] Bioconductor  
## 
## [1] /tmp/RtmpDuFmiq/Rinst685d458b0396
## [2] /home/biocbuild/bbs-3.10-bioc/R/library

References

Love, Michael I., Charlotte Soneson, Peter F. Hickey, Lisa K. Johnson, N. Tessa Pierce, Lori Shepherd, Martin Morgan, and Rob Patro. 2019. “Tximeta: reference sequence checksums for provenance identification in RNA-seq.” bioRxiv. https://doi.org/10.1101/777888.

Patro, Rob, Geet Duggal, Michael I. Love, Rafael A. Irizarry, and Carl Kingsford. 2017. “Salmon Provides Fast and Bias-Aware Quantification of Transcript Expression.” Nature Methods. https://doi.org/10.1038/nmeth.4197.

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

Srivastava, Avi, Laraib Malik, Tom Sean Smith, Ian Sudbery, and Rob Patro. 2019. “Alevin efficiently estimates accurate gene abundances from dscRNA-seq data.” Genome Biology 20 (65). https://doi.org/10.1186/s13059-019-1670-y.