Bio Expression Matrix Counts Ingest

Version Compatibility

Reference examples tested with: pandas 2.2+, numpy 1.26+, scanpy 1.10+, anndata 0.10+, tximport 1.30+, tximeta 1.20+, GenomicFeatures 1.54+, Subread/featureCounts 2.0.6+ (post-v2.0.2 API), STAR 2.7.10+, Salmon 1.10+, kallisto 0.48+, RSEM 1.3.3+, HTSeq 2.0+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• R: packageVersion('<pkg>') then ?function_name to verify parameters
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Count Matrix Ingestion

"Load my featureCounts / Salmon / STAR output into a count matrix" -> Parse the per-sample quantification, strip metadata, choose the correct strandedness or NumReads column, optionally apply length-bias correction via tximport, and return a gene-by-sample matrix appropriate for the downstream DE tool.

The Single Most Important Modern Insight -- lengthScaledTPM is NOT TPM; the naming has fooled many

tximport(..., countsFromAbundance='lengthScaledTPM') returns a count-scale matrix (sum to library size) with the gene-length bias removed, intended as input to DE tools that cannot accept offsets (limma-voom). The word "TPM" in the option name has misled users into reporting these values as normalized abundance -- they are not.

The decision tree for countsFromAbundance (Soneson, Love, Robinson 2015 F1000Res 4:1521):

Two adjacent traps with the same flavor of silent miscount:

1. featureCounts -p since v2.0.2 (March 2021) needs --countReadPairs. Pre-v2.0.2, -p flagged paired-end AND counted pairs as fragments. Post-v2.0.2, -p only flags paired-end -- pairs are NOT counted as one fragment unless --countReadPairs is added. Old scripts run on new installs produce ~2x the expected counts (one count per mate). No warning. Always pass -p --countReadPairs together for paired-end.

2. STAR ReadsPerGene.out.tab column choice depends on library protocol. Column 2 is unstranded; column 3 is forward-stranded; column 4 is reverse-stranded. Illumina TruSeq Stranded (the dominant kit, dUTP-based) is reverse-stranded -- column 4. Reading column 3 for TruSeq throws away ~95% of reads. Reading column 2 conflates antisense expression. Verify with RSeQC infer_experiment.py on a few BAMs before assembling the matrix.

Algorithmic Taxonomy

Decision Tree by Scenario

featureCounts -- The -p and -O Traps

Goal: Run featureCounts correctly on paired-end stranded RNA-seq, then read the output without inheriting the metadata columns.

Approach: CLI invocation with -p --countReadPairs -s <strandedness>, plus optional -O --fraction for overlapping-gene handling; parse with pandas/read.delim stripping the 6 metadata columns.

featureCounts -T 8 -p --countReadPairs -s 2 \
    -t exon -g gene_id \
    -a annotation.gtf \
    -o featurecounts.txt \
    sample1.bam sample2.bam sample3.bam

import pandas as pd

fc = pd.read_csv('featurecounts.txt', sep='\t', comment='#')
counts = fc.set_index('Geneid').iloc[:, 5:]
counts.columns = [c.replace('.bam', '').split('/')[-1] for c in counts.columns]

fc <- read.delim('featurecounts.txt', comment.char = '#', row.names = 1)
counts <- fc[, 6:ncol(fc)]
colnames(counts) <- gsub('.*/|\\.bam$', '', colnames(counts))

STAR --quantMode GeneCounts

Goal: Build a count matrix from STAR's per-sample 4-column output, choosing the correct strandedness column.

Approach: Skip the first 4 summary rows; index column 0 (gene ID); pick column for strandedness; concat across samples.

import pandas as pd
from pathlib import Path

def load_star_genecounts(filepaths, strandedness='reverse'):
    '''Load STAR ReadsPerGene.out.tab files.
    File columns (1-indexed): 1=gene_id, 2=unstranded, 3=forward, 4=reverse.
    After read_csv with index_col=0, the three remaining columns are 0=unstranded, 1=forward, 2=reverse.
    Illumina TruSeq Stranded is 'reverse'.
    '''
    col_map = {'unstranded': 0, 'forward': 1, 'reverse': 2}
    col_idx = col_map[strandedness]
    dfs = {}
    for fp in filepaths:
        sample = Path(fp).name.replace('_ReadsPerGene.out.tab', '')
        df = pd.read_csv(fp, sep='\t', header=None, index_col=0)
        dfs[sample] = df.iloc[4:, col_idx]
    return pd.DataFrame(dfs)

Strandedness verification: for a stranded library, the "wrong" strand column total should be <5% of the "right" strand column total. If comparable, the library was unstranded or there was a kit / config mix-up.

infer_experiment.py -r annotation.bed -i sample.bam

Salmon / kallisto via tximport

Goal: Import transcript-level Salmon or kallisto quantifications to gene-level counts with length-bias correction.

Approach: Build a tx2gene mapping; call tximport() with the right type and countsFromAbundance; hand the result to DESeq2 or limma-voom.

library(tximport)
library(DESeq2)

tx2gene <- read.csv('tx2gene.csv')

files <- file.path('salmon_out', samples$id, 'quant.sf')
names(files) <- samples$id

txi <- tximport(files, type = 'salmon', tx2gene = tx2gene)

dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ condition)

The naive alternative -- summing NumReads to gene level -- is wrong when isoform usage varies across samples. NumReads is normalized against EffectiveLength per transcript; sum-then-DE introduces a length-by-condition bias indistinguishable from differential expression. tximport handles this by carrying the per-sample average transcript length matrix as a DESeq2/edgeR offset.

For limma-voom (no offset mechanism):

txi <- tximport(files, type = 'salmon', tx2gene = tx2gene,
                countsFromAbundance = 'lengthScaledTPM')
y <- DGEList(counts = txi$counts)
v <- voom(y, design, plot = TRUE)

For DTU (DRIMSeq, DEXSeq), use dtuScaledTPM with txOut=TRUE.

For 3'-tagged libraries (10x Chromium 3', QuantSeq), length bias is negligible; use countsFromAbundance='no' WITHOUT the length offset (tximeta's argument or manual disable).

tximeta -- Automatic Provenance

Goal: Import Salmon/kallisto with automatic linkage to the exact annotation release used in the index.

Approach: tximeta (Love, Soneson, Hickey et al. 2020 PLoS Comp Biol 16:e1007664) inspects the Salmon index hash and pulls matching annotation and metadata; the resulting SummarizedExperiment carries the provenance.

library(tximeta)

coldata <- data.frame(
    names = samples$id,
    files = file.path('salmon_out', samples$id, 'quant.sf'),
    condition = samples$condition
)

se <- tximeta(coldata)
gse <- summarizeToGene(se)

library(DESeq2)
dds <- DESeqDataSet(gse, design = ~ condition)

For new projects, tximeta is the modern preference over hand-managed tx2gene -- it eliminates a class of "wrong annotation" bugs.

HTSeq-count

import pandas as pd
from pathlib import Path

def load_htseq_counts(filepaths):
    '''Load HTSeq count files, dropping the __no_feature etc. summary rows.'''
    dfs = {}
    for fp in filepaths:
        sample = Path(fp).stem.replace('_counts', '')
        df = pd.read_csv(fp, sep='\t', header=None, index_col=0,
                         names=['gene', 'count'])
        df = df[~df.index.str.startswith('__')]
        dfs[sample] = df['count']
    return pd.DataFrame(dfs)

HTSeq overlap modes (Anders, Pyl, Huber 2015 Bioinformatics 31:166):

-s yes|no|reverse: TruSeq Stranded is reverse. The --mode and -s flags must match the library; mismatches are silent miscounts.

The __no_feature and __ambiguous totals are QC signals. If __no_feature > 30%: annotation incomplete, wrong reference, or strandedness misspecified. If __ambiguous > 10%: many overlapping gene annotations (common with comprehensive GTFs); consider intersection-nonempty.

RSEM expected_count

library(tximport)
library(DESeq2)

files <- file.path('rsem_out', paste0(samples$id, '.genes.results'))
names(files) <- samples$id
txi <- tximport(files, type = 'rsem', txIn = FALSE, txOut = FALSE)

txi$length[txi$length == 0] <- 1

dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ condition)

RSEM expected_count is non-integer (EM-derived). DESeq2 requires integers, but DESeqDataSetFromTximport rounds appropriately. The txi$length == 0 substitution handles the "Error: all(lengths > 0) is not TRUE" panic from rRNA-filtered or zero-length transcripts.

Avoid round() + DESeqDataSetFromMatrix -- that path loses the length correction.

kallisto + sleuth (Uncertainty-Aware DE)

kallisto quant -i index -o sample1_out -b 100 -t 8 read1_1.fq.gz read1_2.fq.gz

library(sleuth)

s2c <- data.frame(sample = samples$id,
                  condition = samples$condition,
                  path = file.path('kallisto_out', samples$id))

so <- sleuth_prep(s2c, ~ condition,
                  target_mapping = tx2gene,
                  aggregation_column = 'gene_id',
                  gene_mode = TRUE)
so <- sleuth_fit(so, ~ condition, 'full')
so <- sleuth_fit(so, ~ 1, 'reduced')
so <- sleuth_lrt(so, 'reduced', 'full')
results <- sleuth_results(so, 'reduced:full')

Sleuth (Pimentel et al. 2017 Nat Methods 14:687) uses the 100 bootstrap replicates from kallisto to decouple BIOLOGICAL variance (between-replicate) from INFERENTIAL variance (within-replicate, from EM uncertainty). For transcript-level analyses where quantification uncertainty matters (similar isoforms, low coverage), sleuth's response error linear model is more conservative than DESeq2 on plain counts. For well-quantified gene-level DGE, sleuth and DESeq2-via-tximport converge.

Salmon Selective Alignment

Selective alignment (Srivastava 2020 Genome Biol 21:239) -- which combines fast pseudo-mapping with traditional alignment of seed extensions to improve quantification accuracy for transcripts with sequence similarity -- has been the default since Salmon 1.0.0. The historical --validateMappings flag that explicitly requested it is now a deprecated no-op; passing it does nothing.

salmon quant -i index -l A \
    -1 read_1.fq.gz -2 read_2.fq.gz \
    --gcBias --seqBias \
    -p 8 -o sample_out

--gcBias corrects fragment GC bias; --seqBias corrects random hexamer priming bias. Both are recommended for any Salmon run from RNA-seq with biological condition variation. They affect EffectiveLength estimates and propagate through tximport.

10X Genomics

import scanpy as sc

adata = sc.read_10x_mtx('filtered_feature_bc_matrix/')
adata = sc.read_10x_h5('filtered_feature_bc_matrix.h5')

library(Seurat)
mat <- Read10X(data.dir = 'filtered_feature_bc_matrix/')
mat <- Read10X_h5('filtered_feature_bc_matrix.h5')

10X convention: cells in rows (AnnData) or cells in columns (Seurat after Read10X). The R/Python convention differs by transpose.