Analyzing Your Own Bulk RNA-Seq Data in Python

Woolf Software
rna-seq python pydeseq2 salmon transcriptomics
A tall glass instrument column where fish-shaped lights swim upstream trailing glowing strands that are combed into a grid of light on a plate.

By the end of this guide you will have a gene-by-sample count matrix with proper length correction, an AnnData object carrying your sample metadata, a differential expression table with shrunken log2 fold changes, and a ranked enrichment result, all produced inside Python with no R interpreter in the loop. AnnData is the annotated-matrix container used by the Python single-cell ecosystem, which keeps the expression values and the per-sample and per-gene annotations together in one object.

To follow along you need a few things in place. You need paired-end FASTQ files from your own samples (typically *_R1_001.fastq.gz and *_R2_001.fastq.gz), where a FASTQ file holds the sequenced reads along with a per-base quality score for each one. You also need a reference transcriptome and genome FASTA plus a GTF annotation from GENCODE, the GTF being the tab-delimited file that says where each gene and transcript sits on the genome. On the hardware side, plan for roughly 32 GB of RAM and 16 cores for the quantification step, and about 50 GB of disk per 10 samples if you keep intermediates. Python 3.11 or 3.12 is the safe range.

If your data came back as a pre-built count matrix rather than FASTQ, skip to step 4. Before you do, ask the provider whether the counts were produced from an aligner plus a counter or from a transcript-level quantifier, because the two require different normalization.

1. Build the environment and pin it

Reproducibility starts with the environment, so it is worth setting this up deliberately before any data moves. Use a conda-forge and bioconda environment for the command-line tools and pip inside it for the Python packages. Mixing pip and conda for the same dependency tree is where most “conda environment wasn’t working” reports come from, so keep the split clean: binaries from conda, Python analysis libraries from pip.

mamba create -n rnaseq -c conda-forge -c bioconda \
  python=3.11 salmon=1.10.3 fastp=0.23.4 multiqc=1.22 samtools=1.20
mamba activate rnaseq
pip install "pytximport==0.11.*" "pydeseq2==0.5.*" "anndata>=0.10" \
  "scanpy>=1.10" "gseapy>=1.1" pyarrow

Pin PyDESeq2 explicitly. Its API changed between the 0.3 and 0.4 lines, where the design_factors argument was replaced by a formula string in design, and tutorials on the web are split across both. Write the pinned versions into a requirements.txt and record the exact reference files you used, including the GENCODE release number. The single most common reason two analyses of the same FASTQ disagree is a different annotation version rather than a different statistical method.

2. Trim, and read the QC before you quantify

Quantification is the expensive step, so it pays to clean the reads first and to look at what the cleaning tells you about the library. Run fastp for adapter and quality trimming. It is fast, emits a JSON report that MultiQC parses into a single cross-sample summary, and does both read-level filtering and per-sample statistics in one pass.

for s in $(cat samples.txt); do
  fastp -i raw/${s}_R1.fastq.gz -I raw/${s}_R2.fastq.gz \
        -o trim/${s}_R1.fq.gz -O trim/${s}_R2.fq.gz \
        --detect_adapter_for_pe --trim_poly_g \
        --length_required 36 --cut_tail --cut_tail_mean_quality 20 \
        --thread 8 --json qc/${s}.fastp.json --html qc/${s}.fastp.html
done
multiqc qc/ -o qc/multiqc

Two numbers deserve your attention before you spend compute on quantification. The first is the duplication rate: for a poly(A) library at 30 million read pairs, a duplication rate above roughly 50 percent usually means low input RNA and PCR over-amplification, which inflates the variance of low-expression genes. The second is the per-base quality drop at the read ends, together with any GC anomaly in the fastp HTML report, which is your early warning of a library prep problem.

Library chemistry also sets expectations. If the library came from whole blood collected in a PAXgene tube without globin depletion, expect hemoglobin transcripts (HBB, HBA1, HBA2) to consume a large fraction of reads, which reduces effective depth for everything else. Check this after quantification and decide whether the remaining depth supports the question you are asking.

3. Quantify with salmon and a decoy-aware index

With trimmed reads in hand, the next step turns sequences into abundances. We recommend transcript-level quantification with salmon over a full spliced alignment for gene-level expression work. It runs in minutes per sample on a laptop-class core count, corrects for fragment GC and sequence-specific bias, and its output plugs directly into the Python import step.

Build the index against a “gentrome”, meaning the transcriptome concatenated with the whole genome, with genome sequence names supplied as decoys so that reads from unannotated or intronic regions do not get force-assigned to transcripts.

grep "^>" genome.primary_assembly.fa | cut -d " " -f 1 | sed 's/>//g' > decoys.txt
cat gencode.v46.transcripts.fa genome.primary_assembly.fa > gentrome.fa
salmon index -t gentrome.fa -d decoys.txt -i idx_gencode46_k31 \
  --gencode -k 31 -p 16

Two flags are worth understanding here. The --gencode flag tells salmon to split GENCODE’s pipe-delimited FASTA headers and keep only the transcript ID, which saves you a string-mangling step later. The -k flag sets the k-mer length used for matching: keep -k 31 for read lengths of 75 bp and up, and drop to -k 23 only if your reads are 50 bp or shorter.

for s in $(cat samples.txt); do
  salmon quant -i idx_gencode46_k31 -l A \
    -1 trim/${s}_R1.fq.gz -2 trim/${s}_R2.fq.gz \
    --validateMappings --seqBias --gcBias --posBias \
    --numBootstraps 30 -p 16 -o quants/${s}
done

Once each sample finishes, read quants/${s}/logs/salmon_quant.log for the mapping rate. For human poly(A) data against a decoy-aware GENCODE index, 70 to 90 percent is normal. Below about 60 percent, suspect the wrong species or reference. Other causes are substantial genomic DNA contamination or ribosomal RNA carryover. The -l A flag auto-detects library type, and the log reports what it inferred (ISR for most stranded Illumina TruSeq kits). If the inferred type flips between samples in the same batch, the problem is almost certainly in the metadata rather than the biology.

Salmon is the right default for gene-level expression, but it does not place reads on the genome. If your question is about splicing or novel transcripts or anything else requiring genomic positions, use an aligner instead. Align with STAR to a genome index built with --sjdbOverhang 100, then count with HTSeq. Its 2.0 release added multi-core BAM processing and a Python API you can call directly rather than shelling out to htseq-count.1 HTSeq counts against a GTF are the right tool for exon-level and custom-interval counting, at the cost of an order of magnitude more compute and disk than salmon.

4. Import into Python with pytximport

Salmon reports abundances per transcript, and most downstream analysis happens per gene, so the two have to be reconciled carefully. Transcript abundances need to be aggregated to genes with length correction, because a gene’s effective length changes when isoform usage changes. Summing raw transcript counts without that correction biases differential expression. In R this is tximport; in Python, use pytximport, which reimplements the same aggregation and offset logic and was validated to reproduce R tximport output.2

The aggregation needs a two-column transcript-to-gene map. Generate it from the GTF so that it matches the index exactly:

import re, gzip
import pandas as pd

rows = []
with gzip.open("gencode.v46.annotation.gtf.gz", "rt") as fh:
    for line in fh:
        if line.startswith("#"):
            continue
        f = line.split("\t")
        if f[2] != "transcript":
            continue
        attrs = dict(re.findall(r'(\S+) "([^"]*)";', f[8]))
        rows.append((attrs["transcript_id"], attrs["gene_id"], attrs["gene_name"]))

t2g = pd.DataFrame(rows, columns=["transcript_id", "gene_id", "gene_name"])
t2g.to_csv("tx2gene.v46.tsv", sep="\t", index=False)

Identifier formats are the usual source of trouble at this point. If you used the --gencode flag at index time, salmon’s quant.sf holds bare transcript IDs with version suffixes (ENST00000456328.2), which match the GTF’s transcript_id. If you ever mix a versionless reference with a versioned map, pytximport will report a large fraction of unmapped transcripts, and the right fix is to strip versions on both sides rather than guess.

from pathlib import Path
from pytximport import tximport

samples = [l.strip() for l in open("samples.txt")]
files = [Path(f"quants/{s}/quant.sf") for s in samples]

adata = tximport(
    files,
    data_type="salmon",
    transcript_gene_map="tx2gene.v46.tsv",
    counts_from_abundance="length_scaled_tpm",
    ignore_after_bar=True,
    output_type="anndata",
)
adata.obs_names = samples

We use counts_from_abundance="length_scaled_tpm", which produces gene counts scaled so that library size and average transcript length are already accounted for. Those counts go straight into a count-based model without a separate offset matrix, which keeps the downstream code simple. The alternative is passing raw counts plus an average transcript length offset. It is slightly more statistically faithful but requires the model to accept an offset, and PyDESeq2 does not expose one as cleanly.

Attach metadata now, while you still remember what the samples were:

meta = pd.read_csv("metadata.tsv", sep="\t", index_col="sample").loc[samples]
adata.obs = meta          # e.g. columns: timepoint, collection_date, batch, rin
adata.var["gene_name"] = t2g.drop_duplicates("gene_id").set_index("gene_id") \
                            .reindex(adata.var_names)["gene_name"].values

5. Sanity-check the matrix before modeling

Fitting a model to a matrix you have not inspected is how subtle problems become published problems, so look at the data first. Three checks catch most of what goes wrong.

The first is library composition: compute the fraction of counts taken by the top handful of genes. In whole blood without globin depletion, HBB alone can exceed 20 percent of counts, and that ceiling shapes what you can detect.

import numpy as np, scanpy as sc

X = adata.to_df()
cpm = X.div(X.sum(axis=1), axis=0) * 1e6
top = cpm.mean(axis=0).sort_values(ascending=False).head(15)
print(pd.DataFrame({"gene": adata.var["gene_name"].reindex(top.index).values,
                    "mean_cpm": top.values}))

The second is the number of detected genes. With 30 million read pairs from blood, expect roughly 13,000 to 16,000 genes at 1 CPM or more. CPM is counts per million, a simple depth normalization. A sample well below that range is usually degraded or shallow. The third check concerns sample relationships: run a principal component analysis on log-transformed, variance-stabilized-ish values and confirm that the dominant axis tracks the biology you care about rather than batch or collection date.

adata.layers["counts"] = adata.X.copy()
keep = (cpm > 1).sum(axis=0) >= max(2, int(0.2 * adata.n_obs))
adata = adata[:, keep.values].copy()
sc.pp.log1p(adata)
sc.pp.pca(adata, n_comps=min(10, adata.n_obs - 1))
print(np.round(adata.uns["pca"]["variance_ratio"], 3))

The result tells you what to do next. If PC1 separates sequencing batches, include batch in the design. If PC1 separates by RIN, your fold changes will be contaminated by degradation artifacts, and no amount of modeling fully fixes that. RIN is the RNA integrity number, a 1-to-10 measure of degradation.

6. Differential expression with PyDESeq2

With a clean matrix and sensible metadata, you can fit the model. Yes, you can run DESeq2 in Python. PyDESeq2 is a reimplementation of the DESeq2 workflow rather than a wrapper around the R package, so you do not need R or rpy2 installed. Feed it the raw (length-scaled) counts layer, never the log-transformed values.

from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

counts = pd.DataFrame(adata.layers["counts"], index=adata.obs_names,
                      columns=adata.var_names).round().astype(int)

dds = DeseqDataSet(
    counts=counts,
    metadata=adata.obs,
    design="~batch + timepoint",
    refit_cooks=True,
    n_cpus=8,
)
dds.deseq2()

res = DeseqStats(dds, contrast=["timepoint", "week12", "week0"], n_cpus=8)
res.summary()
res.lfc_shrink(coeff="timepoint[T.week12]")
tab = res.results_df.join(adata.var["gene_name"])
tab = tab.sort_values("padj").query("padj < 0.05 and abs(log2FoldChange) > 0.5")

Always call lfc_shrink before you rank genes by effect size. Unshrunken log2 fold changes for low-count genes are dominated by noise and will put uninformative genes at the top of your list.

The hard constraint in personal transcriptomics is sample size. A single individual’s before-and-after comparison with one sample per timepoint has no replication, and DESeq2 cannot estimate dispersion from it. There are two reasonable responses:

  • Collect a longitudinal series. Take five or more timepoints per condition and treat timepoints within a condition as replicates, while including a term for anything that changes systematically. Collection time of day matters for blood, where circadian genes swing severalfold.
  • Drop per-gene hypothesis testing. Work instead with pathway-level scores or ranked correlations across time, which are far more stable at small n.

We prefer the first when a plan allows repeated sampling, and the second when it does not.

7. Ranked enrichment instead of threshold lists

That second option deserves its own treatment, because it is the one most personal datasets end up needing. With a small design, a list of genes passing an adjusted p-value cutoff is often empty or tiny. Rank-based enrichment uses the whole gene ranking and does not depend on a cutoff. Use gseapy with a signed ranking statistic.

import gseapy as gp

full = res.results_df.join(adata.var["gene_name"]).dropna(subset=["stat"])
rnk = (full.assign(gene=full["gene_name"])
           .groupby("gene")["stat"].mean()
           .sort_values(ascending=False))

pre = gp.prerank(rnk=rnk, gene_sets=["GO_Biological_Process_2023",
                                     "Reactome_2022"],
                 permutation_num=1000, min_size=15, max_size=500,
                 threads=8, seed=42, outdir=None)
print(pre.res2d.sort_values("FDR q-val").head(20)
        [["Term", "NES", "FDR q-val", "Tag %"]])

Interpret the normalized enrichment scores (NES) with the gene set size in view. Sets smaller than about 15 genes produce unstable scores, and sets above 500 are usually too diffuse to mean anything specific. Read the leading-edge genes, which are the ones driving each score, rather than stopping at the term name.

8. Put your sample in context with public data

A single person’s expression values are hard to interpret without a reference distribution. Before you conclude that a gene is high or low, find out what its normal range looks like in the same tissue and assay. Recount3, ARCHS4, and GTEx give you bulk expression compendia you can download and compare against. Browser-based tools now assemble harmonized bulk and single-cell references with per-cell-type expression views. That is a quick way to check whether a signal you see in whole blood is plausibly from a specific leukocyte population.3

Comparisons across studies need care. Match the tissue and library chemistry as closely as you can, and be aware that cross-study comparison of absolute TPM (transcripts per million, the length-normalized abundance unit) is unreliable. Ranks and within-sample comparisons travel better than absolute values.

One final point on interpretation. Nothing in an expression table is a diagnosis. Transcript abundances in blood move with time of day and recent exercise. They also move with infection, sleep, and the leukocyte mix of that particular draw. If something in your results looks clinically relevant, take the raw data and the full methods to a physician or a clinical genetics service rather than acting on it yourself.

Common problems

A handful of failure modes account for most of the time lost on this pipeline, and each has a characteristic signature.

A low mapping rate with a correct reference usually means rRNA or genomic DNA in the library. Check the top-expressed features for RNA45S, RN7SL1, or mitochondrial genes. Then check the fastp insert-size distribution for a long tail that suggests DNA carryover.

A pytximport warning that many transcripts are absent from the map almost always means an ID version mismatch or a reference version skew. Verify it by intersecting the first column of a quant.sf with your map directly, and regenerate the map from the same GENCODE release as the index rather than patching IDs by hand.

PyDESeq2 raising an error during dispersion fitting typically means too few samples, or a design matrix that is not full rank. The latter happens when a covariate is perfectly collinear with the condition, for example when every week-0 sample was sequenced in batch 1. Drop the redundant term or redesign the sampling.

All-zero or implausibly small adjusted p-values across thousands of genes point to counts that were transformed before modeling, or to a metadata frame whose row order no longer matches the count matrix. Align by index explicitly, never by position.

Finally, gene symbols are unstable across annotation releases, and multiple Ensembl gene IDs map to the same symbol. Keep Ensembl IDs as the primary key through the whole pipeline and attach symbols only at the reporting layer.

Woolf Software builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. Givanna H Putri, Simon Anders, Paul Theodor Pyl, et al. Analysing high-throughput sequencing data in Python with HTSeq 2.0. Bioinformatics, 2022. https://doi.org/10.1093/bioinformatics/btac166 ↩

  2. Malte Kuehl, Milagros N Wong, Nicola Wanner, et al. Gene count estimation with pytximport enables reproducible analysis of bulk RNA sequencing data in Python. Bioinformatics, 2024. https://doi.org/10.1093/bioinformatics/btae700 ↩

  3. Linh Truong, Thao Truong, Huy Nguyen. OmnibusX: A unified platform for accessible multi-omics analysis. PLOS Computational Biology, 2025. https://doi.org/10.1371/journal.pcbi.1013480 ↩