Import transcript-level quantifications from Salmon/kallisto into R for gene-level analysis with DESeq2/edgeR using tximport or tximeta. Use when importing transcript counts into R for DESeq2/edgeR.
Reference examples tested with: tximport 1.30+, tximeta 1.20+, DESeq2 1.42+, edgeR 4.0+, txdbmaker 1.0+, Salmon 1.10+, kallisto 0.50+
Before using code patterns, verify installed versions match. If versions differ:
packageVersion('<pkg>') then ?function_name to verify parametersIf code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.
"Import Salmon/kallisto results into DESeq2" -> Summarize transcript-level abundance estimates to gene-level counts AND compute a per-gene, per-sample length offset that the DE model consumes.
tximport::tximport(files, type='salmon', tx2gene=tx2gene)tximport does not merely add transcript counts to a gene total. The fragment count a gene produces depends on the average length of the isoforms expressed in that sample, because longer molecules yield more fragments (more start positions). When isoform usage shifts between conditions (differential transcript usage), the gene's average effective length changes, so a naive summed count is length-biased in a condition-correlated way and masquerades as differential expression. tximport corrects this by returning a per-gene, per-sample average-length matrix (txi$length) and passing it as a normalization offset to DESeq2/edgeR. A single per-gene length cannot capture this because the bias is sample-specific (Soneson, Love, Robinson 2015).
Goal: Import transcript-level quantifications into R as gene-level counts plus the length offset for DESeq2 or edgeR.
Approach: Build a transcript-to-gene map, then run tximport over the quant files; the returned txi$counts/txi$length carry both the gene counts and the offset.
library(tximport)
files <- c(sample1 = 'sample1_quant/quant.sf',
sample2 = 'sample2_quant/quant.sf',
sample3 = 'sample3_quant/quant.sf')
tx2gene <- read.csv('tx2gene.csv') # column order: TXNAME, then GENEID
txi <- tximport(files, type = 'salmon', tx2gene = tx2gene)
txi is a list: $abundance (TPM), $counts (estimated counts), $length (the average-length offset source), $countsFromAbundance.
This argument silently determines correctness; nothing errors when it is wrong.
| Mode | What it returns | Use when |
|------|-----------------|----------|
| 'no' (default) | Estimated counts + separate length offset | Full-length library -> DESeq2/edgeR (they consume the offset). The cleanest path. |
| 'lengthScaledTPM' | Counts with the length correction baked in, no separate offset | A tool that cannot take an offset (e.g. limma-voom) |
| 'scaledTPM' | TPM scaled to library size, no length scaling | Transcript-level DTU with txOut=TRUE (DRIMSeq/DEXSeq); the established Love et al. workflow input |
| 'dtuScaledTPM' | Scaled by median isoform length | DTU alternative (tximport >= 1.10), needs tx2gene; helps when isoform lengths within a gene differ widely |
For DTU, scaledTPM is the established default; dtuScaledTPM is the newer purpose-built mode, preferable when a gene's isoforms span very different lengths.
The 3'-tag exception: for 3'-end protocols (10x, QuantSeq, Lexogen) a read count does not scale with transcript length, so there is no length bias to correct, and length-correcting injects one. Do not use the length-scaled modes (lengthScaledTPM/dtuScaledTPM) for tag-seq. Import with the default, but build the DESeqDataSet from the plain counts so the length offset is NOT auto-applied:
# 3'-tag: bypass the length offset that DESeqDataSetFromTximport would otherwise apply
dds <- DESeqDataSetFromMatrix(round(txi$counts), colData = coldata, design = ~ condition)
# Full-length, DESeq2/edgeR (default): keep the offset path
txi <- tximport(files, type = 'salmon', tx2gene = tx2gene)
# Transcript-level for DTU (hand off to alternative-splicing/isoform-switching)
txi_tx <- tximport(files, type = 'salmon', txOut = TRUE,
countsFromAbundance = 'scaledTPM')
The map is a two-column data frame; column ORDER is load-bearing (TXNAME first, GENEID second), names do not matter.
Goal: Map every quantified transcript ID to its gene, with IDs that exactly match the quant files.
Approach: Derive from the annotation that built the index (GTF, ensembldb, biomaRt, or the index t2g); strip version suffixes to match.
# From a GTF: makeTxDbFromGFF moved to txdbmaker in Bioconductor >= 3.19
# (defunct in GenomicFeatures >= 1.61.1; on older Bioconductor use GenomicFeatures::makeTxDbFromGFF)
library(txdbmaker)
txdb <- makeTxDbFromGFF('annotation.gtf')
k <- keys(txdb, keytype = 'TXNAME')
tx2gene <- AnnotationDbi::select(txdb, keys = k, keytype = 'TXNAME',
columns = c('TXNAME', 'GENEID'))
# From biomaRt (useEnsembl; useMart is deprecated)
library(biomaRt)
mart <- useEnsembl(biomart = 'genes', dataset = 'hsapiens_gene_ensembl')
tx2gene <- getBM(attributes = c('ensembl_transcript_id', 'ensembl_gene_id'), mart = mart)
If quant.sf IDs carry version suffixes (ENST00000456328.4) but tx2gene does not (or vice versa), the IDs do not match. Total non-overlap raises an error; partial mismatch silently drops the non-matching transcripts and prints a summary, deflating affected genes toward zero. Fix by stripping versions consistently or with ignoreTxVersion:
txi <- tximport(files, type = 'salmon', tx2gene = tx2gene,
ignoreTxVersion = TRUE, ignoreAfterBar = TRUE)
Goal: Build a DESeqDataSet that uses the tximport length offset without any manual step.
Approach: DESeqDataSetFromTximport stores txi$length as the avgTxLength assay and converts it to per-gene normalization factors inside DESeq().
library(DESeq2)
coldata <- data.frame(condition = factor(c('control', 'control', 'treated', 'treated')),
row.names = names(files))
dds <- DESeqDataSetFromTximport(txi, colData = coldata, design = ~ condition)
dds <- dds[rowSums(counts(dds)) >= 10, ] # light pre-filter (speed); results() does the inferential filter
dds <- DESeq(dds)
res <- results(dds)
Passing a countsFromAbundance='no' txi prints "using counts and average transcript lengths from tximport"; a length-scaled txi prints "using just counts" and applies no offset. Both are handled correctly by the function.
Goal: Carry the length offset into an edgeR DGEList.
Approach: Geometric-mean-center the length matrix, fold in composition-corrected library sizes, log it, attach via scaleOffset.
library(edgeR)
cts <- txi$counts
normMat <- txi$length / exp(rowMeans(log(txi$length))) # center each gene on its geometric mean
normCts <- cts / normMat
eff.lib <- calcNormFactors(normCts) * colSums(normCts)
normMat <- sweep(normMat, 2, eff.lib, '*')
y <- scaleOffset(DGEList(cts), log(normMat))
y <- y[filterByExpr(y, group = coldata$condition), , keep.lib.sizes = FALSE] # group-aware filter
Do not double-apply the offset: if countsFromAbundance='lengthScaledTPM' already baked the correction into the counts, do not also attach a length offset. Use 'no' for the offset path, the scaled modes for the no-offset path, never both.
Gene-level estimates are robust because per-isoform assignment uncertainty cancels on summation. Transcript-level testing must propagate it: edgeR catchSalmon deflates counts by per-transcript overdispersion (differential-expression/edger-basics), swish/fishpond tests across Salmon Gibbs samples (alternative-splicing/isoform-switching), and sleuth uses kallisto bootstraps (expression-matrix/counts-ingest). Generate the replicates at quantification time (rna-quantification/alignment-free-quant).
tximeta hashes the index's reference sequences and looks the digest up against known GENCODE/Ensembl/RefSeq releases, attaching transcript ranges and release metadata automatically, so the exact reference becomes a verified property of the object rather than lab lore.
library(tximeta)
makeLinkedTxome(indexDir = 'salmon_index', source = 'Ensembl', organism = 'Homo sapiens',
release = '110', genome = 'GRCh38', fasta = 'transcripts.fa', gtf = 'annotation.gtf')
coldata <- data.frame(names = names(files), files = files,
condition = c('control', 'control', 'treated', 'treated'))
se <- tximeta(coldata)
gse <- summarizeToGene(se)
dds <- DESeqDataSet(gse, design = ~ condition)
| Symptom | Cause | Fix |
|---------|-------|-----|
| Many genes import as zero or deflated | Transcript-ID version mismatch (partial drop) | ignoreTxVersion = TRUE; or strip \.\d+$ from both sides |
| Error: none of the transcripts present in tx2gene | Total ID mismatch (versions or wrong annotation) | Rebuild tx2gene from the annotation that built the index |
| Summarized at the wrong level, no error | tx2gene columns reversed (GENEID first) | Order as TXNAME, then GENEID |
| Length bias appears in 3'-tag data | DESeqDataSetFromTximport auto-applied the length offset | Build via DESeqDataSetFromMatrix(round(txi$counts), ...) so no offset is applied |
| Fold changes inflated near isoform switches with manual edgeR | Offset double-applied or omitted | One path only: 'no'+offset, or scaled mode without offset |
npx skills add GPTomics/bio-rna-quantification-tximport-workflow下载完整 Skill 目录，包含 SKILL.md 及所有相关文件
Category:science-education