How to Analyze Your Own RNA-seq Data, From FASTQ to Gene-Level Interpretation

Oak
rnaseq transcriptomics salmon deseq2 workflow
A glass column of glowing fish-shaped tokens threaded with light feeds a sleek lab machine that sorts them into three lit channels.

By the end of this guide you will have, for each of your RNA-seq timepoints, a per-gene expression matrix in three forms (estimated counts, length-corrected TPM, and variance-stabilized values), a QC report you can defend, an estimate of the blood cell composition driving the signal, and a ranked list of pathways that moved between timepoints. You need paired-end FASTQ files (typically 30 to 60 million read pairs per sample for bulk whole blood), roughly 32 GB of RAM and 8 to 16 cores, about 150 GB of disk, and a working install of fastp, salmon, STAR (optional), samtools, and R 4.3+ with tximport, DESeq2, limma, and fgsea. Everything below assumes human GRCh38 and GENCODE v47 annotation. None of this is a diagnostic procedure, and no output here should be read as a clinical result; expression changes are interpreted by a clinician in the context of your history and standard labs, not by a fold-change threshold.

1. Inventory what you were given and check its provenance

Before any tool runs, write down four facts about each sample: the library chemistry (poly(A) selection or ribosomal RNA depletion), whether globin transcripts were depleted, the strandedness of the protocol, and the sequencing run identifier. These four determine the correct flags downstream, and getting strandedness wrong silently halves your counts for many genes. The standard RNA-seq workflow is preprocessing, alignment or quantification, expression estimation, differential analysis, and functional interpretation, and each of those stages inherits assumptions from the library prep you cannot recover later 1.

Confirm strandedness empirically rather than trusting the metadata. Run salmon quant -l A once (Section 3) and read the lib_format_counts.json it writes: a reverse-stranded Illumina TruSeq library will report the overwhelming majority of fragments as ISR. If the compatible fragment fraction is split roughly evenly between ISF and ISR, the library is unstranded, and you should say so explicitly in your notes.

2. Trim and QC with fastp, then read the numbers that matter

Use fastp for adapter removal and quality trimming. It is fast, writes a single JSON you can parse, and does adapter detection for paired-end data without you supplying sequences.

fastp \
  -i raw/T1_R1.fastq.gz -I raw/T1_R2.fastq.gz \
  -o trim/T1_R1.fq.gz   -O trim/T1_R2.fq.gz \
  --detect_adapter_for_pe \
  --cut_tail --cut_tail_window_size 4 --cut_tail_mean_quality 20 \
  --length_required 36 \
  --correction \
  --thread 8 \
  --json qc/T1.fastp.json --html qc/T1.fastp.html

Do not trim aggressively. Selective-alignment quantifiers tolerate a few low-quality bases far better than they tolerate systematically shortened reads, which bias effective transcript lengths. The four numbers to read from the JSON are the fraction of reads passing filter (expect above 0.93), duplication rate (whole blood poly(A) libraries at 40M pairs usually sit between 0.15 and 0.45), Q30 fraction (expect above 0.90), and insert size peak (a peak below about 150 bp indicates a degraded input RNA and will distort length normalization). Aggregate everything with multiqc qc/ so the timepoints sit side by side, because the comparison you care about is across your own samples over time, not against a published average.

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

We recommend selective alignment with salmon as the primary quantification path and reserve genome alignment for the specific questions that need read coordinates. Salmon on a decoy-aware index gives you transcript-level estimates with bias correction and bootstrap replicates in roughly a tenth of the compute of a STAR run, and the decoy genome suppresses the spurious assignment of intronic and intergenic reads to transcripts they merely resemble.

Build the index once:

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

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

Then quantify each timepoint:

salmon quant -i idx/salmon_gencode_v47 -l A \
  -1 trim/T1_R1.fq.gz -2 trim/T1_R2.fq.gz \
  -p 16 --seqBias --gcBias --posBias \
  --numBootstraps 30 \
  --validateMappings \
  -o quants/T1

The --gcBias flag matters more than people assume for longitudinal work, because library prep batches separated by months differ in GC behavior and that difference masquerades as biology. The 30 bootstraps cost little and give you inferential replicates, which is the only way to distinguish technical quantification uncertainty from real change when you have one sample per timepoint. Check the mapping rate in quants/T1/logs/salmon_quant.log: whole blood poly(A) libraries should land between 75% and 90%, and anything below 65% means either the wrong index, heavy rRNA carryover, or genomic DNA contamination.

Run STAR in addition only if you want splice junction counts, fusion detection, or allele-specific expression tied to your genome data:

STAR --genomeDir idx/star_gencode_v47 \
  --readFilesIn trim/T1_R1.fq.gz trim/T1_R2.fq.gz --readFilesCommand zcat \
  --runThreadN 16 --outSAMtype BAM SortedByCoordinate \
  --quantMode GeneCounts \
  --outFilterMultimapNmax 20 --alignSJoverhangMin 8 \
  --waspOutputMode SAMtag --varVCFfile personal.het.vcf --outSAMattributes NH HI AS nM vW \
  --outFileNamePrefix bam/T1.

The WASP flags require a VCF of your heterozygous sites and mark reads that fail the allelic mapping-bias test with vW, which you filter out before counting reference and alternate reads at each site. Without that filter, allele-specific expression estimates are biased toward the reference allele.

4. Collapse transcripts to genes with tximport

Transcript-level estimates are noisier and harder to interpret than gene-level ones, and most downstream gene set resources are gene-keyed. Import with tximport using countsFromAbundance = "lengthScaledTPM", which produces counts that are on a count scale but already corrected for changes in average transcript length between samples.

library(tximport); library(readr)

tx2gene <- read_tsv("gencode.v47.tx2gene.tsv",
                    col_names = c("tx_id", "gene_id", "gene_name"))

files <- file.path("quants", c("T1","T2","T3","T4"), "quant.sf")
names(files) <- c("T1","T2","T3","T4")

txi <- tximport(files, type = "salmon",
                tx2gene = tx2gene[, c("tx_id","gene_id")],
                countsFromAbundance = "lengthScaledTPM",
                ignoreTxVersion = FALSE)

If tximport warns about transcripts missing from tx2gene, your FASTA and GTF are from different GENCODE releases. Fix that rather than setting ignoreTxVersion = TRUE as a workaround, because version mismatches quietly drop genes.

5. Filter, normalize, and stabilize variance across your timepoints

Now build a DESeqDataSet and drop the genes that carry no information. In whole blood, roughly 12,000 to 14,000 of the ~20,000 protein-coding genes are detectable at useful depth, and the rest add multiple-testing burden without adding signal.

library(DESeq2)
coldata <- data.frame(
  row.names = colnames(txi$counts),
  timepoint = factor(c("T1","T2","T3","T4"), levels = c("T1","T2","T3","T4")),
  batch     = factor(c("B1","B1","B2","B2"))
)

dds <- DESeqDataSetFromTximport(txi, coldata, design = ~ batch + timepoint)
keep <- rowSums(counts(dds) >= 10) >= 2
dds  <- dds[keep, ]

dds <- estimateSizeFactors(dds)
vsd <- vst(dds, blind = TRUE)

Inspect sizeFactors(dds). A spread wider than about 0.5 to 2.0 across timepoints means one library is dominated by a small number of transcripts, which in whole blood almost always means globin. Then plot the sample-to-sample distance matrix and a PCA of assay(vsd). If PC1 separates by batch rather than by timepoint, you have a technical structure problem that no downstream test will rescue, and the only clean fixes are to include batch in the design (as above) or to re-prep the libraries together.

6. Ask what changed, with the right model for one person

With a single individual sampled repeatedly, there are no biological replicates in the usual sense, and this is the step where most self-analyses go wrong. Two approaches work and they answer different questions.

If you have several timepoints within a stable condition and several within a changed condition (before and after a sustained intervention, for instance), treat the timepoints within each condition as replicates and run a standard DESeq2 contrast. The dispersion estimated this way includes your own week-to-week biological variability, which is the correct null for asking whether something is different.

dds$condition <- factor(c("pre","pre","post","post"))
design(dds) <- ~ condition
dds <- DESeq(dds)
res <- lfcShrink(dds, coef = "condition_post_vs_pre", type = "apeglm")
sum(res$padj < 0.05 & abs(res$log2FoldChange) > log2(1.5), na.rm = TRUE)

If instead you have one sample per condition, do not run a two-sample differential test. Use the salmon bootstraps through fishpond/swish to get quantification uncertainty, and express each gene as a deviation from your own running baseline: a z-score computed on assay(vsd) across all prior timepoints. A gene that sits three standard deviations from your personal baseline is a hypothesis worth following into the next draw, not a result. This personal-baseline framing is the whole point of longitudinal molecular data, and it makes the interpretation step depend on your own history rather than on a population reference range.

7. Correct for cell composition before you believe any pathway result

Whole blood RNA is a weighted average of neutrophil, monocyte, lymphocyte, and platelet transcriptomes, and the weights move with sleep, exercise, acute infection, and time of day. A large fraction of apparent “immune activation” signatures in longitudinal blood RNA-seq is a shift in the neutrophil-to-lymphocyte ratio, not a change in per-cell transcription.

Estimate composition from the expression matrix itself using a reference signature matrix (LM22 with CIBERSORTx, or the granulator Bioconductor package with dtangle or nnls deconvolution) run on TPMs, not on variance-stabilized values. If you have a complete blood count with differential drawn at the same time, use the measured proportions instead, since they are more accurate than any deconvolution. Then add the two or three dominant proportions as covariates:

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

Fit that model and compare the gene list to the uncorrected one. The genes that survive are the ones worth spending attention on.

8. Interpret with ranked gene set enrichment, not a threshold list

Cut-off-based enrichment throws away the ranking information you paid for. Use fgsea on the full shrunken log fold change vector against Reactome or Hallmark sets.

library(fgsea); library(msigdbr)
sets <- split(msigdbr(species = "Homo sapiens", collection = "H")$gene_symbol,
              msigdbr(species = "Homo sapiens", collection = "H")$gs_name)
ranks <- setNames(res$log2FoldChange, symbols)
ranks <- sort(ranks[!is.na(ranks)], decreasing = TRUE)
fg <- fgseaMultilevel(sets, ranks, minSize = 15, maxSize = 500)

Report normalized enrichment score, adjusted p-value, and leading-edge genes together. A pathway with an NES of 2.1 driven by four leading-edge genes that are all in the same operon-like cluster is a weaker result than one with NES 1.5 spread across forty genes. Web-based suites such as ExpressAnalyst within the Analyst software family will run the same enrichment and network views without local installation if you prefer a graphical path, and they integrate transcriptome results with metabolomic and proteomic layers in one session 2. Interactive workflows like BIOMEX are useful at this stage for exploring bulk and single-cell expression matrices visually before committing to a statistical claim 3. Older but still sound reviews of omics interpretation workflows are worth reading for how to keep the biological question in front of the enrichment machinery rather than behind it 4.

9. Make the run reproducible and store it so future-you can query it

Pin every version: GENCODE release, salmon version and index hash (quants/T1/cmd_info.json records both), R package versions via renv::snapshot(). Six months from now you will add a timepoint, and if you quantify it against a different index the new sample will separate from the old ones on PC1 for purely technical reasons. Re-quantify all timepoints together whenever the index changes.

For storage, a plain directory of Parquet files plus a samples.tsv manifest will carry you a long way, and for larger longitudinal collections a graph-backed store lets you query relationships among genes, timepoints, and enriched terms rather than re-deriving them from flat files each time 5. If you would rather not maintain infrastructure, the nf-core/rnaseq pipeline pinned to a specific revision (nextflow run nf-core/rnaseq -r 3.14.0 -profile docker) gives you the same salmon and STAR steps with provenance recorded automatically, and cloud analysis platforms provide comparable end-to-end handling for large multi-omic collections 6.

Common problems

Globin transcripts dominating the library is the most frequent failure in whole blood. Without globin depletion, HBB, HBA1, and HBA2 can absorb over half of all reads, which collapses your effective depth for everything else. Check quant.sf for those three genes as a fraction of total TPM before anything else. If they are high and you cannot re-prep, remove them and their pseudogenes prior to estimateSizeFactors so normalization is not driven by them, and accept that lowly expressed genes are underpowered in that sample.

A mapping rate below 65% usually has one of three causes: rRNA carryover in a ribo-depletion library that underperformed (check RNA45SN1, RNA18SN5), genomic DNA contamination (intronic reads, visible as a low exonic fraction in qualimap rnaseq or Picard CollectRnaSeqMetrics), or the wrong organism or annotation build. Run Picard’s RnaSeqMetrics and read PCT_INTRONIC and PCT_RIBOSOMAL_BASES before re-running anything.

Degraded input RNA produces 3’ bias that inflates expression estimates for short transcripts and deflates them for long ones. CollectRnaSeqMetrics reports MEDIAN_5PRIME_TO_3PRIME_BIAS; values far from 1.0 mean gene-length-correlated artifacts, and comparing such a sample to a high-quality one will generate a long list of false positives that all share a length distribution. Plot log fold change against transcript length as a diagnostic whenever a result surprises you.

Batch confounded with time is the structural trap of longitudinal single-person data, because later timepoints were necessarily prepped later. Mitigate it by banking RNA and prepping in balanced batches, or by including a shared reference aliquot in every batch that you can use to estimate and remove batch effects with limma::removeBatchEffect for visualization (never on the counts fed to DESeq2).

Finally, resist over-interpreting a single flagged gene. Transcript abundance in blood is noisy, responsive to the previous night’s sleep, and only loosely coupled to protein levels for many genes. The value of the workflow above is the trajectory across many timepoints, and any finding that looks clinically meaningful belongs in a conversation with a physician, supported by the appropriate validated assay, rather than in a spreadsheet by itself 1.

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. In Seok Yang, Sangwoo Kim. Analysis of Whole Transcriptome Sequencing Data: Workflow and Software. Genomics & Informatics, 2015. https://doi.org/10.5808/gi.2015.13.4.119 ↩ ↩2

  2. Jessica D. Ewald, Guangyan Zhou, Yao Lu, et al. Web-based multi-omics integration using the Analyst software suite. Nature Protocols, 2024. https://doi.org/10.1038/s41596-023-00950-4 ↩

  3. Federico Taverna, Jermaine Goveia, Tobias K Karakach, et al. BIOMEX: an interactive workflow for (single cell) omics data interpretation and visualization. Nucleic Acids Research, 2020. https://doi.org/10.1093/nar/gkaa332 ↩

  4. Irmgard Mühlberger, Julia Wilflingseder, Andreas Bernthaler, et al. Computational Analysis Workflows for Omics Data Interpretation. Methods in Molecular Biology, 2011. https://doi.org/10.1007/978-1-61779-027-0_17 ↩

  5. Raquel L. Costa, Luiz Gadelha, Marcelo Ribeiro-Alves, et al. GeNNet: an integrated platform for unifying scientific workflows and graph databases for transcriptome data analysis. PeerJ, 2017. https://doi.org/10.7717/peerj.3509 ↩

  6. Wen Li, Zhining Zhang, Bo Xie, et al. HiOmics: A cloud-based one-stop platform for the comprehensive analysis of large-scale omics data. Computational and Structural Biotechnology Journal, 2024. https://doi.org/10.1016/j.csbj.2024.01.002 ↩