How to Run a Transcriptomics Workflow on Your Own RNA-Seq Data

By the end of this guide you will have the following for each blood draw.

• a set of quality metrics
• a transcript-level quantification
• a gene-level count matrix assembled across all your timepoints
• a variance-stabilized expression matrix suitable for plotting and clustering
• estimated blood cell-type fractions
• a ranked pathway table describing how one timepoint differs from the rest

You need the raw FASTQ files from your sequencing provider. They are paired-end and gzip-compressed, typically 3 to 8 GB per sample at 30 million read pairs. You also need roughly 100 GB of free disk, 32 GB of RAM, a machine with 8 or more cores, and a working conda or container setup. Everything below runs on a laptop with patience or a single cloud instance in about an hour per sample. The interpretation sections assume you are reading your own data for your own understanding, and anything that looks clinically relevant belongs in a conversation with a physician rather than in a Jupyter notebook.

1. Establish what your library is

Transcriptomics is the measurement of the set of RNA molecules in a sample, and RNA sequencing is the dominant way to do that measurement, so for practical purposes the two terms overlap. The distinction that matters is that transcriptomics is the question and RNA-seq is one instrument, with microarrays, targeted panels, and Nanopore direct RNA sequencing as other instruments answering the same question with different biases ¹. Before you touch a read, write down five properties of your library, because every downstream choice depends on them.

First, the source tissue and collection tube. Whole blood drawn into PAXgene or Tempus tubes stabilizes RNA at the bedside and captures the full leukocyte population. PBMCs isolated by Ficoll drop granulocytes and change the expression profile substantially. Second comes the selection chemistry. Poly-A selection captures polyadenylated mRNA and discards most non-coding RNA. Ribosomal depletion with Ribo-Zero or RiboCop retains long non-coding RNAs and unprocessed transcripts at the cost of a noisier library. Third, read length and pairing, usually 2×100 or 2×150. Fourth, strandedness, which for most modern kits is reverse-stranded (dUTP protocol). Fifth, sequencing depth in read pairs.

For a longitudinal personal profile, we recommend poly-A selected and reverse-stranded 2×150 libraries at 30 to 50 million read pairs per timepoint. Thirty million pairs is enough to quantify roughly the top 12,000 to 14,000 expressed genes in whole blood with usable precision, and going deeper buys mostly low-count transcripts whose variance is dominated by sampling noise. Depth matters less than consistency: if your first three draws were 30M and your fourth is 80M, you have introduced a covariate you will have to model.

2. Quality control the raw reads

Run fastp first, because it does adapter trimming, quality filtering, and reporting in one pass and is fast enough to be unremarkable in the pipeline.

fastp \
  -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
  -o trim_R1.fastq.gz -O trim_R2.fastq.gz \
  --detect_adapter_for_pe \
  --length_required 36 \
  --cut_tail --cut_tail_mean_quality 20 \
  --json fastp.json --html fastp.html \
  --thread 8

We deliberately do not aggressively quality-trim. Pseudoaligners and splice-aware aligners tolerate low-quality tails well, and hard trimming below Q20 removes real signal and shortens reads into multimapping territory. Adapter removal and a minimum length filter are enough.

Then aggregate everything with MultiQC once you have more than one sample:

multiqc -o qc/ fastp_reports/ salmon_out/

The numbers to watch across timepoints: duplication rate, GC distribution, and the fraction of reads surviving filtering. A GC curve with a shoulder at high GC usually means rRNA carryover. A bimodal insert size distribution means your library prep had a size-selection problem. Keep these per-sample metrics in the same table as your expression data, because when a timepoint looks like an outlier six months from now, the explanation is usually in the QC file rather than in biology.

3. Quantify with Salmon against a decoy-aware index

For personal longitudinal profiling, we use Salmon in selective-alignment mode rather than a full genome alignment. It runs in five to ten minutes per sample on eight cores and produces transcript-level estimates with bootstrap uncertainty. It also corrects for GC and sequence-specific bias, which matters when your samples were prepped in different batches months apart. The cost is that you get no BAM file. You cannot inspect splice junctions, call variants, or look at coverage over a specific exon. If any of that is part of your plan, run STAR as well, and keep both.

Build the index once, using the genome as a decoy so that reads from unannotated genomic regions are not forced into transcripts:

# GENCODE v46 human transcripts + primary assembly genome
grep "^>" GRCh38.primary_assembly.genome.fa | cut -d " " -f 1 | sed 's/>//g' > decoys.txt
cat gencode.v46.transcripts.fa.gz GRCh38.primary_assembly.genome.fa.gz > gentrome.fa.gz

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

The -k 31 default is right for 100 bp and longer reads. Drop to -k 23 only if your reads are 50 bp or shorter. Then quantify:

salmon quant \
  -i salmon_index_gencode46 \
  -l A \
  -1 trim_R1.fastq.gz -2 trim_R2.fastq.gz \
  --validateMappings --gcBias --seqBias \
  --numBootstraps 30 \
  -p 16 -o salmon_out/T1

-l A infers the library type. Check salmon_out/T1/lib_format_counts.json afterward and confirm it reports ISR for a standard reverse-stranded paired library. If it reports IU (unstranded) when you expected stranded, your provider used a different kit than you assumed, and you should fix the assumption rather than the flag.

Mapping rate is the first real signal about sample quality. Whole blood poly-A libraries typically land between 75% and 90% against a decoy-aware index. Below 60%, suspect rRNA contamination, adapter dimer, or degraded RNA.

4. Assemble the gene-level matrix

Salmon gives you transcript-level estimates. Most interpretation happens at the gene level, and tximport handles the aggregation correctly by carrying the effective-length information into an offset that DESeq2 uses, which is more accurate than summing TPMs.

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

txdb <- makeTxDbFromGFF("gencode.v46.annotation.gtf.gz")
k <- keys(txdb, keytype = "TXNAME")
tx2gene <- select(txdb, k, "GENEID", "TXNAME")

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

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

The single most common failure here is a version-suffix mismatch: Salmon reports ENST00000456328.2 and your tx2gene table has ENST00000456328. Set ignoreTxVersion = TRUE if you see a large fraction of transcripts dropped, and verify the overlap explicitly with mean(rownames(txi$counts) %in% tx2gene$GENEID) before proceeding. Pin the annotation version in your project notes. Gene models change between GENCODE releases, and a count matrix built from v44 is not directly comparable to one built from v46.

5. Normalize with a longitudinal single-subject design in mind

This is where a personal profile departs from a standard experiment. A conventional RNA-seq study compares groups of individuals, and the statistical machinery assumes biological replicates drawn from a population ². You have one individual measured repeatedly, so your replicates are timepoints, and the variance you are modeling is within-person temporal variance plus technical variance. That is a legitimate design, but it changes what a p-value means. You are asking whether a gene moved more than it usually moves in you. You are not asking whether it differs between people.

library(DESeq2)

dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ batch + condition)
keep <- rowSums(counts(dds) >= 10) >= 3
dds <- dds[keep, ]

vsd <- vst(dds, blind = FALSE)
plotPCA(vsd, intgroup = c("condition", "batch"))

Filter before normalizing. Requiring at least 10 counts in at least 3 samples typically leaves 13,000 to 16,000 genes from whole blood and removes most of the noise that inflates multiple-testing burden. Use vst rather than rlog once you have more than about a dozen samples, on speed grounds. Then look at the PCA. If PC1 separates your samples by sequencing batch rather than by anything biological, you have a batch effect large enough to dominate every downstream result. You must either model it explicitly in the design or remove it for visualization with limma::removeBatchEffect(assay(vsd), batch = vsd$batch). Never remove it for visualization and then ignore it in the model.

6. Estimate cell composition before interpreting any gene

To a first approximation, whole blood expression is a weighted average of the expression profiles of neutrophils, lymphocytes, monocytes, and platelets. The weights are the abundances of those cell types. If your neutrophil fraction rises from 55% to 72% between two draws, hundreds of genes will appear to change, and none of that change tells you anything about regulation inside a cell. Estimate composition first and treat it as the primary explanation until you have ruled it out.

# CIBERSORTx (web or docker) or immunedeconv wrapping several methods
library(immunedeconv)
tpm <- txi$abundance[keep, ]
frac <- deconvolute(tpm, method = "quantiseq")

Compare the estimates against your CBC differential from the same draw if you have one. Agreement within a few percentage points for the major populations is a good sign that both measurements are sound. Where they disagree badly, trust the CBC and treat the deconvolution as a soft prior. Once you have fractions, you can add them as covariates: design = ~ batch + neutrophil_frac + condition. This is the single change that most improves the interpretability of longitudinal blood transcriptomics.

7. Test for change across timepoints

With covariates in place, run the contrast you care about.

dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "post", "baseline"),
               alpha = 0.05, lfcThreshold = 0.5)
resLFC <- lfcShrink(dds, coef = "condition_post_vs_baseline", type = "apeglm")
summary(res)

We set lfcThreshold = 0.5 deliberately. Testing against zero in a well-powered design returns long lists of genes with a 10% change that are statistically real and biologically uninterpretable. Testing against a fold-change floor asks a question you can act on. Apply shrinkage before ranking genes by effect size, because low-count genes otherwise produce enormous, meaningless log fold changes.

If you have many timepoints and want to model continuous trends rather than discrete contrasts, switch to limma-voom with a spline basis on time, or use DESeq2’s likelihood ratio test with a reduced model that drops the time term. Either approach tells you which genes have any time-dependent structure, which is usually a better first question than which genes differ between two arbitrary draws.

8. Move to pathways, because single genes are noisy

Individual transcript measurements in blood are volatile. Pathway-level summaries aggregate across dozens of genes and are considerably more stable across draws, which is why published transcriptomic workflows consistently route from differential expression into enrichment and network analysis rather than stopping at a gene list ³. Use a ranked-list method rather than a threshold-based overlap test, because you keep the full ordering and avoid the arbitrariness of a cutoff.

library(fgsea); library(msigdbr)

stats <- resLFC$log2FoldChange * -log10(resLFC$pvalue)
names(stats) <- rownames(resLFC)
stats <- sort(stats[is.finite(stats)])

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

fg <- fgsea(pathways, stats, minSize = 15, maxSize = 500, nPermSimple = 10000)
head(fg[order(padj)], 20)

Start with the Hallmark collection (50 curated sets) before Reactome or GO, because Hallmark sets are non-redundant and each one maps to a process you can name. Interferon response, inflammatory response, oxidative phosphorylation, and cholesterol homeostasis are the sets that move most often in longitudinal blood. Transcriptome-wide profiling is used the same way in drug and biomarker discovery, where the readout of interest is a coordinated program rather than any one transcript ⁴.

9. Connect the transcriptome to your other layers

A transcriptome read in isolation answers fewer questions than the same transcriptome read against a genome, a proteome, and a set of blood biomarkers from the same draw. Integration across molecular layers is the direction the field has moved, precisely because each layer constrains the interpretation of the others ⁵. Three concrete joins are worth doing.

The first is allele-specific expression, which requires the BAM file from a STAR alignment plus your phased genotypes, and tells you whether a heterozygous variant produces imbalanced expression of the two alleles. The second is the protein-transcript comparison: correlate your measured plasma protein levels against the corresponding transcript counts and note that for most genes the correlation is weak, which is information about post-transcriptional regulation rather than an error. The third is simply plotting a pathway score alongside a biomarker trajectory over time. Multi-omic resources and integration methods for exactly this kind of cross-layer analysis are now well catalogued ⁶, and machine learning methods for combining the layers are an active area ⁷.

10. Make the whole thing reproducible

Run the pipeline as a workflow rather than as a shell history. We use nf-core/rnaseq, which wraps trimming, Salmon or STAR quantification, and QC aggregation with container-pinned tool versions:

nextflow run nf-core/rnaseq -r 3.14.0 \
  --input samplesheet.csv \
  --outdir results/ \
  --genome GRCh38 \
  --pseudo_aligner salmon \
  --skip_alignment \
  -profile docker

Then record the annotation release, the Salmon index checksum, and the container digests. Also keep an renv.lock for the R side. The reason is that you will add a timepoint next year and need the new sample processed identically to the old ones, otherwise the pipeline version becomes a batch effect. Reproducibility and extensibility are the design goals of current multi-omics workflow frameworks for exactly this reason ⁸.

Common problems

Globin transcripts dominate whole blood libraries. HBB, HBA1, and HBA2 can consume half or more of your reads in a non-depleted poly-A library from whole blood. That effectively halves your usable depth and compresses the dynamic range for everything else. Either request globin depletion at library prep (GLOBINclear or an equivalent) or budget extra depth. Computationally, you can remove globin genes from the matrix before normalization, but that does not recover the reads you lost.

Degraded RNA produces 3’ bias. If the RNA integrity number (RIN) is below 7, coverage concentrates at transcript 3’ ends and effective lengths are wrong. Long transcripts are then systematically underestimated. Ask your provider for the RIN or DV200 value on every sample and treat it as a covariate if it varies.

Time of day is a real covariate. Blood transcript abundance for a substantial fraction of genes varies with the circadian cycle, and a 7 a.m. draw is not comparable to a 4 p.m. draw. Standardize your draw time and fasting state across timepoints. This single procedural choice removes more variance than any analytical correction.

Acute illness overwhelms everything. An unremarkable viral infection produces an interferon signature that dominates PC1 and makes the sample nearly useless for tracking slower changes. Record symptoms at each draw and be willing to exclude a timepoint.

Spatial transcriptomics is a different measurement and in no sense an upgrade. It quantifies RNA while preserving the physical position of each measurement within a tissue section, using either spatially barcoded capture arrays or imaging-based in situ hybridization. Its analysis pipeline adds segmentation, registration, and neighborhood statistics on top of the steps above ⁹. It requires a tissue biopsy, so it is not part of a blood-based personal profile, and searching for a spatial workflow when you have blood FASTQs will send you down the wrong path.

Finally, resist the pull of a single dramatic gene. A twofold change in one transcript in one draw is within the range of ordinary variation for many blood genes. Believe a result when it appears as a coherent pathway, replicates across draws, and agrees with something measured on a different platform. And if a pattern in your data seems to bear on your health, take the data to a clinician rather than resolving it yourself.

Footnotes

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

2. Laetitia B. B. Martin, Zhangjun Fei, James J. Giovannoni, et al. Catalyzing plant science research with RNA-seq. Frontiers in Plant Science, 2013. https://doi.org/10.3389/fpls.2013.00066 ↩

3. Amadeo Muñoz García, Lars M. Eijssen, Martina Kutmon, et al. A bioinformatics workflow to decipher transcriptomic data from vitamin D studies. The Journal of Steroid Biochemistry and Molecular Biology, 2019. https://doi.org/10.1016/j.jsbmb.2019.01.003 ↩

4. Xiaonan Yang, Ling Kui, Min Tang, et al. High-Throughput Transcriptome Profiling in Drug and Biomarker Discovery. Frontiers in Genetics, 2020. https://doi.org/10.3389/fgene.2020.00019 ↩

5. Biswapriya B Misra, Carl Langefeld, Michael Olivier, et al. Integrated omics: tools, advances and future approaches. Journal of Molecular Endocrinology, 2019. https://doi.org/10.1530/jme-18-0055 ↩

6. Xianjun Dong, Chunyu Liu, Mikhail Dozmorov. Review of multi-omics data resources and integrative analysis for human brain disorders. Briefings in Functional Genomics, 2021. https://doi.org/10.1093/bfgp/elab024 ↩

7. Somayah Albaradei. AI-Driven Approaches to Utilization of Multi-Omics Data for Personalized Diagnosis and Treatment of Cancer: A Comprehensive Review. Computer Modeling in Engineering & Sciences, 2025. https://doi.org/10.32604/cmes.2025.072584 ↩

8. Zeping Gui, Huidong Su, Yuhao Chen, et al. Ragomics‐Meta: An open‐source, scalable, reproducible, and extensible multi‐omics analysis workflow. iMetaOmics, 2026. https://doi.org/10.1002/imo2.70121 ↩

9. Xiao Li, Ximo Pechuan-Jorge, Tyler Risom, et al. SPEX: A modular end-to-end platform for high-plex tissue spatial omics analysis. GigaScience, 2025. https://doi.org/10.1093/gigascience/giaf090 ↩