Building an RNA-Seq Pipeline for Your Own Blood Transcriptome

Oak
rna-seq transcriptomics bioinformatics pipeline quantification
A small feathered reptile on a black plinth, glowing filaments running through its plumage and fanning into glass cylinders of luminous fluid.

By the end of this guide you will have, for each blood draw you have sequenced, a transcript-level and gene-level expression matrix in TPM and estimated counts, a set of QC metrics you can compare across timepoints, and a splice-junction table you can query for unusual isoform usage. You need the raw reads (paired-end FASTQ files, ideally gzipped, typically 30–60 million read pairs per sample for bulk whole blood), roughly 100 GB of disk per sample during processing, 32 GB of RAM if you plan to align to the genome and 16 GB if you stay with lightweight quantification, and a machine with at least 8 cores. Everything below runs on Linux or macOS with conda or Docker. The commands assume human data and GRCh38. If you have a single timepoint, most of this still applies, but the interpretation section will be thin: a personal transcriptome becomes informative mostly by comparison against your own prior samples, because the between-person variance in gene expression dwarfs the within-person variance for most genes.

A note before the steps. RNA expression in whole blood is dominated by cell composition. Neutrophil fraction alone explains a large share of the variance in a bulk blood transcriptome, so a change you see between January and June may be a change in which cells were in the tube rather than a change in what any cell was doing. The pipeline below produces the numbers. Deciding what those numbers mean about your health is a clinical question, and nothing here is a diagnostic test.

1. Set up a reproducible environment and pin your references

Reproducibility is the single highest-value decision in a longitudinal project, because you will be comparing a sample processed today against one processed two years ago. If the annotation changes underneath you, gene-level counts shift for reasons that have nothing to do with biology. Pipeline frameworks exist precisely to hold this fixed across runs and across machines 1.

Create an environment and pin versions explicitly:

mamba create -n rnaseq -c conda-forge -c bioconda \
  fastp=0.23.4 salmon=1.10.3 star=2.7.11b samtools=1.20 \
  subread=2.0.6 qualimap=2.3 multiqc=1.22 portcullis=1.2.4
mamba activate rnaseq
mamba env export --no-builds > env.rnaseq.yml

Download one reference set and never silently upgrade it. We use GENCODE primary assembly plus the matching comprehensive annotation:

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}.primary_assembly.annotation.gtf.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_${REL}/gencode.v${REL}.transcripts.fa.gz
md5sum *.gz > reference.md5

Record the release number in the output directory name for every run. When you eventually need to move to a newer annotation, reprocess all historical samples with the new one rather than mixing.

2. Trim, and inspect adapter and duplication structure

Run fastp first, both as a trimmer and as your first QC read. It is fast enough to run on raw data without a separate FastQC pass and emits a JSON you can parse into a longitudinal QC table.

fastp \
  -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
  -o sample_R1.trim.fastq.gz -O sample_R2.trim.fastq.gz \
  --detect_adapter_for_pe \
  --trim_poly_g --poly_g_min_len 10 \
  --trim_poly_x --poly_x_min_len 10 \
  --cut_front --cut_tail --cut_mean_quality 20 \
  --length_required 36 \
  --correction \
  --thread 8 \
  --json sample.fastp.json --html sample.fastp.html

--trim_poly_g matters on two-colour Illumina chemistry (NovaSeq, NextSeq), where a dark cycle is read as G and unfilled reads acquire long poly-G tails. Skipping it costs you alignment rate in a way that looks like a library problem. --correction uses overlap between R1 and R2 to fix mismatched bases in the overlapping region, which helps with short inserts.

Read the JSON rather than the HTML for anything you want to track: summary.before_filtering.q30_rate, duplication.rate, adapter_cutting.adapter_trimmed_reads, and insert_size.peak. For a standard poly-A library you want Q30 above about 0.90 and an insert peak comfortably above read length. A duplication rate above 40% in a 30M-pair library usually means low input RNA and over-amplification, not genuinely saturated sequencing.

Do not aggressively quality-trim beyond the settings above. Quantifiers handle mismatches fine, and hard trimming shortens reads in a way that increases multimapping ambiguity. The distinction between reads that are lost to library artifacts and reads that are lost to over-processing is exactly where most usable data goes missing in practice 2.

3. Quantify with selective alignment against a decoy-aware index

For a personal profile whose primary output is expression levels, we recommend salmon in selective-alignment mode as the main quantifier, with genome alignment as a second, optional track. The reasoning is straightforward: selective alignment gives transcript-level estimates with explicit handling of GC and positional bias, it runs in roughly ten minutes per sample on 8 cores, and its error modes are well characterized. Pseudoalignment-class methods were the change that made routine transcript-level quantification tractable at all 3.

Build a decoy-aware index once. The decoys are the genome sequence, which prevents reads from unannotated genomic regions being forced onto a transcript:

grep "^>" <(gunzip -c GRCh38.primary_assembly.genome.fa.gz) \
  | cut -d " " -f1 | sed 's/>//g' > decoys.txt
cat gencode.v45.transcripts.fa.gz GRCh38.primary_assembly.genome.fa.gz > gentrome.fa.gz

salmon index \
  -t gentrome.fa.gz -d decoys.txt \
  -p 16 -i salmon_idx_gencode45 --gencode -k 31

--gencode strips the pipe-delimited GENCODE FASTA headers down to the transcript ID, which saves you a string-parsing step later. Keep -k 31 unless your reads are shorter than 50 bp, in which case drop to 21 or 23.

Then quantify:

salmon quant \
  -i salmon_idx_gencode45 -l A \
  -1 sample_R1.trim.fastq.gz -2 sample_R2.trim.fastq.gz \
  --validateMappings \
  --seqBias --gcBias --posBias \
  --numGibbsSamples 30 \
  --threads 16 \
  -o quant/sample

-l A infers library type, which you should then verify in quant/sample/lib_format_counts.json: a stranded dUTP protocol should report ISR with over 95% of fragments consistent. --gcBias corrects fragment-level GC bias and is worth the extra runtime when comparing samples from different library preps, which is exactly your situation over years of draws. --numGibbsSamples 30 gives you posterior samples so you can quantify inferential uncertainty per transcript, which matters enormously for genes with many similar isoforms.

Collapse to gene level in R with tximport, which correctly handles the average-transcript-length offset rather than naively summing TPM:

library(tximport); library(readr)
tx2gene <- read_tsv("tx2gene.gencode45.tsv")  # transcript_id, gene_id
files <- setNames(file.path("quant", samples, "quant.sf"), samples)
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM")

Use lengthScaledTPM counts when you plan to feed a downstream tool that expects counts but you also want length correction baked in.

4. Align to the genome for QC and junction-level work

Quantification alone will not tell you whether your library is contaminated with genomic DNA, how much of it is intronic, or whether a specific splice junction is behaving oddly. For that you need a genome alignment. STAR is the practical choice for human data with 32 GB of RAM.

STAR --runMode genomeGenerate \
  --genomeDir star_idx_gencode45 \
  --genomeFastaFiles GRCh38.primary_assembly.genome.fa \
  --sjdbGTFfile gencode.v45.primary_assembly.annotation.gtf \
  --sjdbOverhang 99 --runThreadN 16

STAR --genomeDir star_idx_gencode45 \
  --readFilesIn sample_R1.trim.fastq.gz sample_R2.trim.fastq.gz \
  --readFilesCommand zcat \
  --outSAMtype BAM SortedByCoordinate \
  --outSAMattributes NH HI AS nM NM MD \
  --outFilterMultimapNmax 20 \
  --outFilterMismatchNoverLmax 0.04 \
  --alignSJoverhangMin 8 --alignSJDBoverhangMin 1 \
  --alignIntronMin 20 --alignIntronMax 1000000 \
  --quantMode GeneCounts \
  --runThreadN 16 \
  --outFileNamePrefix star/sample.
samtools index star/sample.Aligned.sortedByCoord.out.bam

Set --sjdbOverhang to read length minus one. If your reads are 150 bp, use 149 and rebuild the index. --quantMode GeneCounts gives you a free ReadsPerGene.out.tab with unstranded, forward, and reverse columns, and comparing those three columns is the fastest way to confirm your library’s strandedness independently of salmon’s inference.

Then get the QC picture:

qualimap rnaseq \
  -bam star/sample.Aligned.sortedByCoord.out.bam \
  -gtf gencode.v45.primary_assembly.annotation.gtf \
  -pe -s --java-mem-size=16G \
  -outdir qc/qualimap/sample
multiqc . -o qc/multiqc

The numbers we watch per sample: exonic fraction (expect 65–85% for a poly-A blood library), intronic fraction (a jump above roughly 30% suggests degraded RNA or nuclear contamination), intergenic fraction (above 10% points to genomic DNA carryover), and 5’–3’ coverage bias. If the median 5’ bias drops well below 0.7 relative to 3’, the RNA was degraded before library prep and gene-length-dependent artifacts will contaminate every comparison you make with that sample.

For junction-level analysis, STAR’s SJ.out.tab is a starting point but contains a substantial number of spurious junctions from misalignment near repeats. Filtering junctions on alignment-derived features rather than raw read count is considerably more accurate, which is what Portcullis does 4:

portcullis full -t 16 -o portcullis_out \
  GRCh38.primary_assembly.genome.fa \
  star/sample.Aligned.sortedByCoord.out.bam

The filtered BED it produces is what you should query if you ever want to ask whether a particular gene is using an unusual isoform.

5. Normalize across your own timepoints

This is the step where a personal longitudinal dataset differs most from a standard two-group experiment. You are not testing a treatment against a control. You are asking whether gene X is different in you today than in you six months ago, against a background of technical variation and cell-composition drift.

Start with a variance-stabilizing transform on the gene-level counts and look at the sample-sample correlation structure:

library(DESeq2)
dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ 1)
keep <- rowSums(counts(dds) >= 10) >= ceiling(0.5 * ncol(dds))
dds <- dds[keep, ]
vsd <- vst(dds, blind = TRUE)
plotPCA(vsd, intgroup = c("batch", "draw_date"))

If PC1 separates your samples by sequencing batch rather than by date, you have a batch effect large enough to swamp biology, and you should either include batch in the design or process future samples in the same prep. Filtering low-count genes before any transformation is not cosmetic: unfiltered low-count genes inflate multiple-testing burden and destabilize dispersion estimates, and standardized frameworks built for regulatory-grade reporting make this filter an explicit, documented step 5.

Then estimate cell composition and regress it out or model it explicitly. Deconvolution tools such as CIBERSORTx or the granulator Bioconductor package with the LM22 signature matrix will give you fractions for the major leukocyte types from bulk TPM. Add neutrophil and lymphocyte fractions as covariates:

design(dds) <- ~ neutrophil_frac + lymphocyte_frac + condition

Without this, a seasonal infection two days before a draw will show up as hundreds of “differentially expressed” genes that are a shift in cell counts. Any interpretation of what a persistent expression change means for you individually belongs with a clinician who can see your full picture, including a CBC drawn at the same time.

6. Move from gene lists to something interpretable

A ranked list of genes is not yet knowledge. Run gene set enrichment on the full ranked statistic rather than a thresholded list, because thresholding throws away the coordinated small changes that pathway analysis is best at detecting:

library(fgsea)
ranks <- setNames(res$stat, res$symbol)
ranks <- ranks[!is.na(ranks)]
fg <- fgseaMultilevel(pathways = hallmark, stats = ranks,
                      minSize = 15, maxSize = 500)

Use MSigDB Hallmark first (50 well-curated sets, low redundancy), then Reactome if you need finer resolution. The main interpretive trap is redundancy: fifteen overlapping interferon sets all reaching significance is one signal, not fifteen. Collapse with collapsePathways() before you read anything into the result. The broader pattern of choices here, from quantifier to enrichment method to how functional results are reported, is laid out well in recent practical reviews 6.

If you also have proteomics from the same draw, correlate transcript and protein for the genes where you have both. Expect a Spearman correlation in the 0.4–0.6 range genome-wide, higher for abundant secreted proteins and poor for anything post-translationally regulated. Where the two disagree persistently, the disagreement is itself information, and integrating layers formally rather than eyeballing them is what multi-omics pipelines are built for 7.

Common problems

Low alignment rate with normal Q30. Check for poly-G first (rerun fastp with --trim_poly_g and compare), then check for rRNA. Whole blood libraries prepared without globin depletion will also have 50–70% of reads consumed by HBA1, HBA2, and HBB. Salmon will happily quantify them, but your effective depth for everything else collapses. Count reads assigned to the globin genes in quant.sf and compute your usable-read fraction before concluding anything about library quality.

Strandedness disagreement. If STAR’s ReadsPerGene.out.tab column 3 and column 4 are roughly equal, the library is unstranded regardless of what the kit claimed. Forcing -l ISR on unstranded data silently halves your counts. Always verify empirically.

TPM that will not compare across timepoints. TPM is compositional: it sums to one million, so a single gene rising forces every other gene down. If one sample has a large induction of a few transcripts, TPM comparisons elsewhere become misleading. Use DESeq2’s median-of-ratios normalization on counts for comparisons, and reserve TPM for within-sample ranking. For absolute comparisons, a spike-in gives you an external anchor that composition-based normalization cannot provide 8.

Annotation drift. Between GENCODE releases, gene IDs are retired and transcript models change. If you notice a gene “disappearing” between two timepoints, check whether the gene ID still exists in the newer annotation before looking for a biological explanation. This is the argument for reprocessing everything when you upgrade.

Unstable transcript-level estimates. Genes with many highly similar isoforms (immunoglobulins, HLA, long non-coding RNAs) will have high inferential variance. Check the Gibbs posterior spread from --numGibbsSamples before believing an isoform-switch story, and do gene-level analysis for these families.

Too few replicates. With one sample per timepoint, you have no within-timepoint variance estimate, and any single-sample “differential expression” is a point estimate with unquantified error. The fix is time: five or six draws give you an empirical within-person distribution against which a new sample can be scored. Single-cell approaches sidestep the composition problem entirely by measuring each cell separately, at higher cost and with different technical artifacts 9.

Oak 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. Kathleen M. Fisch, Tobias Meißner, Louis Gioia, et al. Omics Pipe: a community-based framework for reproducible multi-omics data analysis. Bioinformatics, 2015. https://doi.org/10.1093/bioinformatics/btv061 ↩

  2. Felix Pförtner, Eva Briem, Wolfgang Enard, et al. Increasing usable reads in RNA-seq protocols. iScience, 2026. https://doi.org/10.1016/j.isci.2026.116984 ↩

  3. Paul A McGettigan. Transcriptomics in the RNA-seq era. Current Opinion in Chemical Biology, 2013. https://doi.org/10.1016/j.cbpa.2012.12.008 ↩

  4. Daniel Mapleson, Luca Venturini, Gemy Kaithakottil, et al. Efficient and accurate detection of splice junctions from RNA-seq with Portcullis. GigaScience, 2018. https://doi.org/10.1093/gigascience/giy131 ↩

  5. M.C. Verheijen, T.W. Gant, W. Tong, et al. R-ODAF: Omics data analysis framework for regulatory application. Toxicology Letters, 2021. https://doi.org/10.1016/s0378-4274(21)00539-7 ↩

  6. 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 ↩

  7. Bianka Alexandra Pasat, Eleftherios Pilalis, Katarzyna Mnich, et al. MultiOmicsIntegrator: a nextflow pipeline for integrated omics analyses. Bioinformatics Advances, 2024. https://doi.org/10.1093/bioadv/vbae175 ↩

  8. Sabrina Tan, Jia Min Lee, Yan Lin Li, et al. A biological spike-in enables cost-effective multi-species RNA-seq and reveals global transcriptional collapse under dark stress. 2026. https://doi.org/10.64898/2026.06.16.732598 ↩

  9. Yanxiang Deng, Amanda Finck, Rong Fan. Single-Cell Omics Analyses Enabled by Microchip Technologies. Annual Review of Biomedical Engineering, 2019. https://doi.org/10.1146/annurev-bioeng-060418-052538 ↩