How to Build a Bulk RNA-seq Pipeline for Your Own Samples

By the end of this guide you will have a gene-level count matrix and a set of quality-control reports that tell you whether that matrix is trustworthy. You will also have a differential expression table with shrunken log fold changes and adjusted p-values, plus ranked pathway scores. To start you need paired-end FASTQ files, typically 30–60 million read pairs per sample for a standard poly(A) library. You also need a machine with at least 32 GB of RAM and 8 cores, since Salmon indexing of the human gentrome is the memory peak at around 20–25 GB. Budget roughly 200 GB of disk for references and intermediates, and have conda or a container runtime available. Bulk RNA-seq measures the average expression across all cells in the input material, so a whole-blood sample gives you a weighted mixture of neutrophils, lymphocytes, monocytes, and platelets rather than any single cell type. That averaging shapes every interpretation downstream, and we come back to it in step 9.

1. Establish what your library is

Before running anything, write down five facts about the library, because each one changes a flag or a filtering step later. The first is the selection chemistry (poly(A) capture versus ribosomal RNA depletion). The second is whether globin mRNA was depleted (relevant only for whole blood). The others are the strandedness, the read length and whether reads are paired, and the RNA integrity number (RIN) from the prep. Poly(A) libraries concentrate reads on mature mRNA and are the right default for gene-level expression. Ribo-depleted total RNA libraries capture non-polyadenylated transcripts and introns, which raises intronic signal and reduces the effective depth for coding genes at the same read count.

If you do not know the strandedness, infer it rather than guess. Salmon’s -l A will detect the library type and report it in lib_format_counts.json, and a typical dUTP protocol reports ISR for paired-end reverse-stranded data. Getting this wrong silently reduces counts for genes with antisense neighbors.

2. Build the reference

Use GENCODE for human work, take a single release, and pin it. Mixing annotation versions between samples is one of the most common sources of a batch effect that looks biological.

REL=45
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_${REL}/GRCh38.primary_assembly.genome.fa.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_${REL}/gencode.v${REL}.transcripts.fa.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_${REL}/gencode.v${REL}.primary_assembly.annotation.gtf.gz
gunzip *.gz

For Salmon, build a decoy-aware index. Decoys are the genomic sequence itself, added so that reads originating from unannotated or intronic regions can be assigned to the genome instead of being forced onto a transcript. This matters: without decoys, intronic and intergenic reads inflate the expression of transcripts that happen to share short stretches of sequence.

grep "^>" GRCh38.primary_assembly.genome.fa | cut -d " " -f1 | sed 's/>//g' > decoys.txt
cat gencode.v45.transcripts.fa GRCh38.primary_assembly.genome.fa > gentrome.fa
salmon index -t gentrome.fa -d decoys.txt -p 16 -k 31 --gencode -i salmon_v45_decoy

Keep -k 31 for reads of 75 bp or longer. Drop to -k 23 or so only for reads shorter than about 50 bp, where a 31-mer seed becomes too restrictive. The --gencode flag strips GENCODE’s pipe-delimited FASTA headers down to the transcript ID, which saves you a string-parsing step in tximport.

3. Trim and QC the reads

We use fastp because it does adapter detection, quality trimming, poly-G trimming for two-color chemistry, and reporting in one pass.

fastp \
  -i S1_R1.fastq.gz -I S1_R2.fastq.gz \
  -o S1_R1.trim.fastq.gz -O S1_R2.trim.fastq.gz \
  --detect_adapter_for_pe \
  --trim_poly_g \
  --cut_tail --cut_tail_mean_quality 20 \
  --length_required 36 \
  --thread 8 \
  --json S1.fastp.json --html S1.fastp.html

Resist aggressive quality trimming. Selective-alignment and splice-aware aligners tolerate low-quality tails well, and hard trimming to a fixed length mainly reduces mappability. Adapter removal and poly-G removal are the parts that change results. Run FastQC on the trimmed files and aggregate everything with multiqc . at the end of the run so that all per-sample reports sit in one HTML file.

Read three numbers off the MultiQC report. The first is the duplication rate, where high duplication with low library complexity suggests too little input RNA. The second is the GC distribution, where a sharp secondary peak often means adapter dimer or rRNA carryover. The third is the fraction of reads surviving trimming, which should be above roughly 90% for a healthy library.

4. Quantify

For gene-level expression we recommend Salmon with selective alignment as the default. It is roughly an order of magnitude faster than genome alignment, handles multi-mapping reads probabilistically through an expectation-maximization step, and gives you transcript-level estimates that can be collapsed to genes. Surveys of RNA-seq practice converge on the same general structure of quality control, quantification, and model-based differential testing regardless of which quantifier you pick.¹

salmon quant \
  -i salmon_v45_decoy \
  -l A \
  -1 S1_R1.trim.fastq.gz -2 S1_R2.trim.fastq.gz \
  -p 16 \
  --gcBias --seqBias --posBias \
  --numBootstraps 30 \
  --validateMappings \
  -o quants/S1

Turn on --gcBias for any dataset where libraries were prepped in more than one batch, since GC-dependent coverage differences are a frequent and correctable driver of apparent expression change. The bootstrap samples give you a technical variance estimate per transcript, which is useful when you later want to know whether a small fold change is even measurable.

Check quants/S1/logs/salmon_quant.log for the mapping rate. For human poly(A) data against a decoy-aware index, we expect 75–90%. Below about 60%, stop and investigate. The usual causes are the wrong species, heavy rRNA contamination, or adapter sequence still in the reads.

Use STAR instead when you need genomic coordinates: variant calling from RNA, allele-specific expression, novel splice junction discovery, or fusion detection. The flow there is STAR --runMode genomeGenerate --sjdbGTFfile gencode.v45.primary_assembly.annotation.gtf --sjdbOverhang 100, then per sample --quantMode TranscriptomeSAM GeneCounts --outSAMtype BAM SortedByCoordinate --twopassMode Basic. Budget about 32 GB of RAM for the human index. If you are already running STAR, featureCounts -p --countReadPairs -s 2 -T 8 -a annotation.gtf -o counts.txt *.bam gives a straightforward gene count matrix.

5. Collapse to gene level

tximport aggregates transcript abundances to genes and produces an offset matrix that carries the average transcript length per gene per sample, so that downstream models correct for changes in isoform usage.

library(tximport); library(GenomicFeatures); library(readr)

txdb  <- makeTxDbFromGFF("gencode.v45.primary_assembly.annotation.gtf")
k     <- keys(txdb, keytype = "TXNAME")
tx2gene <- AnnotationDbi::select(txdb, k, "GENEID", "TXNAME")

files <- file.path("quants", samples$id, "quant.sf")
names(files) <- samples$id

txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "no", ignoreTxVersion = FALSE)

Keep countsFromAbundance = "no" when you are handing the object to DESeq2, because DESeqDataSetFromTximport wants raw-scale counts plus the length offset. Use "lengthScaledTPM" when you want a matrix for limma-voom or for a quick correlation analysis.

6. Interrogate the count matrix before you model it

This is the step that most tutorials skim and where most mistakes are caught. Filter to genes with a reasonable chance of being measurable, apply a variance-stabilizing transform, and look at the sample-sample structure.

library(DESeq2)
dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ condition)
keep <- rowSums(counts(dds) >= 10) >= min(table(samples$condition))
dds <- dds[keep, ]

vsd <- vst(dds, blind = TRUE)
plotPCA(vsd, intgroup = c("condition", "batch"))
sampleDists <- dist(t(assay(vsd)))

Three diagnostics matter. First, library size after filtering: samples with fewer than roughly 10 million assigned read pairs will have noticeably wider confidence intervals and often drive spurious hits. Second, the fraction of counts in the top 50 genes. In whole blood without globin depletion, HBB, HBA1, and HBA2 alone can take 50–70% of reads, which gutted your effective depth. Third, the position of each sample on PC1 and PC2 relative to known technical variables. If batch separates samples more cleanly than the biological variable does, either include batch in the design or remove those samples.

For a first pass with interactive plots, bulkAnalyseR packages the filtering, normalization, dimensionality reduction, and differential expression steps into a Shiny app that generates a shareable app directory from a count matrix.² If you work in Python, OmicVerse covers a similar range and connects bulk results to single-cell references in the same environment.³

7. Fit the differential expression model

DESeq2 is the tool we use for designs with biological replicates. It models counts with a negative binomial distribution, shrinks per-gene dispersion estimates toward a fitted trend, and applies log fold change shrinkage so that low-count genes do not dominate a ranked list.

design(dds) <- ~ batch + condition
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "treated", "control"), alpha = 0.05)
resLFC <- lfcShrink(dds, coef = "condition_treated_vs_control", type = "apeglm")
summary(res)

Report lfcShrink output rather than raw log2FoldChange, and use lfcThreshold = 1 inside results() if you want the p-value to test a fold change larger than two, rather than filtering on a threshold after testing. Filtering on effect size after testing on a null of zero is a common and quiet error.

In Python, PyDESeq2 reimplements the DESeq2 model. It reproduces its size factors, dispersion estimates, and Wald statistics closely enough to be used interchangeably for standard two-group designs.⁴

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

dds = DeseqDataSet(counts=counts_df, metadata=meta_df,
                   design_factors=["batch", "condition"], refit_cooks=True)
dds.deseq2()
stat = DeseqStats(dds, contrast=["condition", "treated", "control"])
stat.summary()
stat.lfc_shrink(coef="condition_treated_vs_control")

For longitudinal sampling of one person, which is the case for a personal profile with repeated blood draws, the design changes. With several draws per state you can fit ~ timepoint and use DESeq(dds, test = "LRT", reduced = ~ 1) to find genes whose expression changes across time at all. With paired samples from one subject and a covariate of interest, limma-voom with duplicateCorrelation handles the within-subject correlation properly. What no method can do is give you a valid p-value from a single pair of unreplicated samples. In that situation, report fold changes with the Salmon bootstrap intervals attached. Treat the result as a ranked hypothesis list and do not convert it into a claim.

8. Move to pathway level

Individual gene results are noisy, and interpretation is more stable at the level of coordinated gene sets. We prefer a rank-based approach over a threshold-based overlap test, because it uses the whole ranked list and avoids the arbitrary cutoff. Rank by the Wald statistic. Do not rank by fold change.

library(fgsea); library(msigdbr)

ranks <- res$stat; names(ranks) <- rownames(res)
ranks <- sort(ranks[!is.na(ranks)], decreasing = TRUE)

sets <- msigdbr(species = "Homo sapiens", category = "H")
pathways <- split(sets$gene_symbol, sets$gs_name)

fg <- fgsea(pathways, ranks, minSize = 15, maxSize = 500, eps = 0)
head(fg[order(padj), ], 20)

Start with the 50 Hallmark sets, which are curated to be relatively non-redundant. Only then move to Gene Ontology or Reactome, where thousands of overlapping sets make the multiple-testing correction punishing and the output hard to read. Practical guides to functional analysis are consistent in warning that enrichment output depends heavily on the background gene set and the ranking metric, so record both.⁵

One caution on interpretation. A transcript change in whole blood tells you about the transcriptional state of circulating cells at the moment of the draw. It can shift with a poor night’s sleep, a recent infection, or the time of day. Reviews of transcriptional and multi-omics work in stem cell aging make the same point in a different setting: expression signatures are informative about state and are most convincing when they replicate across timepoints and modalities.⁶ If a finding looks clinically consequential, it belongs in a conversation with a physician. A self-interpreted report is not the place for it.

9. Deconvolve composition if the input is blood or tissue

Because bulk measurement averages over cells, an apparent expression change often reflects a shift in cell proportions rather than a change in per-cell transcription. A rise in neutrophil fraction from 55% to 70% will move hundreds of genes. Estimate composition explicitly with a reference-based deconvolution method. Options are CIBERSORTx with the LM22 signature matrix, EPIC, or granulator in R, which wraps several methods for comparison. Then include the leading composition estimates as covariates in the design, or interpret the results in light of them. Reviews of bulk and single-cell analysis methods treat deconvolution as the standard bridge between the two data types, and a matched single-cell reference gives markedly better estimates than a generic one.⁷

10. Make the run reproducible

Once the pipeline works, replace your shell scripts with nf-core/rnaseq. It implements trimming, Salmon or STAR quantification, and the full QC suite as a versioned Nextflow workflow.

nextflow run nf-core/rnaseq -r 3.14.0 \
  --input samplesheet.csv \
  --outdir results \
  --gtf gencode.v45.primary_assembly.annotation.gtf \
  --fasta GRCh38.primary_assembly.genome.fa \
  --pseudo_aligner salmon \
  --skip_alignment \
  -profile docker

Pin the revision with -r, keep the pipeline version and the GENCODE release in the results directory, and archive software_versions.yml. For longitudinal profiling this matters more than for a one-off study, since a sample collected two years from now has to be comparable to the one you processed today. Integrated multi-omics designs are only as good as the consistency of each layer’s processing. That is why published integration pipelines put so much emphasis on version control and harmonized preprocessing.⁸ When you later want to connect an expression signal to your own genotype, having the RNA and the genome processed against the same reference build makes joint analyses considerably less painful. One such analysis is comparing your expression to published eQTL and GWAS results.⁹

Common problems

Low mapping rate with high-quality reads usually means rRNA. Quantify it directly by checking counts on rRNA gene biotypes in the GTF, or run a quick bbduk.sh ref=ribokmers.fa screen. Above about 20% rRNA, the library is depth-limited and re-sequencing is cheaper than fighting the analysis.

A batch effect on PC1 that correlates with sequencing run is correctable in the model (~ batch + condition) but not if batch and condition are confounded. If every treated sample ran in one batch, the design is not identifiable and no software will save it. Randomize sample-to-batch assignment at collection.

Duplicate gene symbols after collapsing from Ensembl IDs break downstream joins. Keep Ensembl gene IDs as the primary key throughout and attach symbols only for display. Also decide once whether to strip version suffixes (ENSG00000141510.17 to ENSG00000141510) and apply it everywhere.

Thousands of significant genes at an adjusted p-value below 0.05 with tiny fold changes usually indicates an unmodeled covariate or a composition shift. A dramatic biological effect is not the usual cause. Check the PCA, check the deconvolution estimates, and consider adding a lfcThreshold.

Zero significant genes with a clear biological expectation often reflects depth or replicate count. Three replicates per group is the practical floor for detecting moderate effects, and additional biological replicates buy more power than additional reads once you are past roughly 20 million assigned pairs per sample.

On cost, library prep plus 30–50 million paired reads per sample typically lands in the low hundreds of dollars at an academic core or commercial provider, with the prep often exceeding the sequencing cost at this depth. Compute is minor by comparison: a full Salmon run on a 40 million pair sample takes a few minutes on 16 cores.

Footnotes

1. Ana Conesa, Pedro Madrigal, Sonia Tarazona, et al. A survey of best practices for RNA-seq data analysis. Genome Biology, 2016. https://doi.org/10.1186/s13059-016-0881-8 ↩

2. Ilias Moutsopoulos, Eleanor C Williams, Irina I Mohorianu. bulkAnalyseR: an accessible, interactive pipeline for analysing and sharing bulk multi-modal sequencing data. Briefings in Bioinformatics, 2022. https://doi.org/10.1093/bib/bbac591 ↩

3. Zehua Zeng, Yuqing Ma, Lei Hu, et al. OmicVerse: a framework for bridging and deepening insights across bulk and single-cell sequencing. Nature Communications, 2024. https://doi.org/10.1038/s41467-024-50194-3 ↩

4. Boris Muzellec, Maria Teleńczuk, Vincent Cabeli, et al. PyDESeq2: a python package for bulk RNA-seq differential expression analysis. Bioinformatics, 2023. https://doi.org/10.1093/bioinformatics/btad547 ↩

5. Jiung-Wen Chen, Lisa Shrestha, George Green, et al. The hitchhikers’ guide to RNA sequencing and functional analysis. Briefings in Bioinformatics, 2023. https://doi.org/10.1093/bib/bbac529 ↩

6. Bongsoo Park, Hagai Yanai, Jun Ding, et al. Hematopoietic stem cell aging: a review of transcriptional and multi-omics insights and potential paths for AI integration. Experimental & Molecular Medicine, 2026. https://doi.org/10.1038/s12276-026-01805-0 ↩

7. Jorge A. Tzec‐Interián, Daianna González‐Padilla, Elsa B. Góngora‐Castillo. Bioinformatics perspectives on transcriptomics: A comprehensive review of bulk and single‐cell RNA sequencing analyses. Quantitative Biology, 2025. https://doi.org/10.1002/qub2.78 ↩

8. Esraa Hamdi Abdelaziz, Rasha Ismail, Mai S. Mabrouk, et al. Multi-omics data integration and analysis pipeline for precision medicine: Systematic review. Computational Biology and Chemistry, 2024. https://doi.org/10.1016/j.compbiolchem.2024.108254 ↩

9. Mingkai Gong, Xian Zhao, Qingze Li, et al. Integrated Analysis of Bulk RNA Sequencing, eQTL, GWAS, and Single‐Cell RNA Sequencing Reveals Key Genes in Hepatocellular Carcinoma. Journal of Cellular and Molecular Medicine, 2025. https://doi.org/10.1111/jcmm.70359 ↩