What You Get From RNA-Seq, and How to Analyze It

RNA-seq output comes in three layers. The raw layer is a pair of gzipped FASTQ files per sample, named sample_R1.fastq.gz and sample_R2.fastq.gz. These typically hold 20 to 60 million read pairs at 4 to 12 GB compressed. The alignment layer is a coordinate-sorted BAM with its .bai index, usually 3 to 8 GB, plus a log of mapping rates. The quantification layer is what you will spend your time in. It is a per-transcript abundance table, either quant.sf from Salmon or ReadsPerGene.out.tab from STAR. That table collapses to a gene-by-sample matrix of estimated counts, effective lengths, and TPM. Everything downstream, from normalization to pathway scoring, is derived from that matrix. Analyzing sequencing results means checking QC and choosing a quantifier. You then convert transcript estimates to gene-level counts and normalize correctly. After that, be disciplined about what a single person’s transcriptome can and cannot support.

The file inventory, and what each file is for

The FASTQ files are the only irreplaceable artifact. Each read occupies four lines: identifier, sequence, a + separator, and Phred quality scores in ASCII offset 33. Keep them, checksum them (md5sum *.fastq.gz > manifest.md5), and treat every derived file as reproducible. If your provider delivers a single interleaved FASTQ or unaligned BAM (uBAM), convert with samtools fastq -1 R1.fq.gz -2 R2.fq.gz -0 /dev/null -s /dev/null -n input.bam.

The BAM is an aligned, compressed representation of those same reads. It adds genomic coordinates, CIGAR strings, and mapping quality. You need it for anything positional. That covers splice junction inspection, allele-specific expression, variant calling from RNA and coverage plots in IGV. You do not need it for expression quantification, which is why a lightweight pipeline can skip alignment entirely. If you keep it, keep the genome build and annotation version that produced it, because a BAM against GRCh38 with GENCODE v44 is not interchangeable with one against v29.

The quantification files are small and are where the biology lives. Salmon’s quant.sf has five columns: Name (transcript ID), Length, EffectiveLength, TPM, and NumReads. The effective length correction matters because a 500 bp transcript and a 5,000 bp transcript at the same molar concentration generate very different read counts, and TPM (transcripts per million) is the length-corrected, library-size-normalized unit that fixes this. Alongside these you should have a MultiQC HTML report summarizing FastQC, adapter trimming, duplication, and mapping statistics across every sample.

The pipeline we would run

We would use nf-core/rnaseq with Salmon as the quantifier and STAR alignment retained for inspection. It is versioned, containerized, and produces a consistent directory tree. That matters more for longitudinal data than for a one-off experiment, since you will reprocess your timepoints together. Drift in tool versions or annotation releases is the most common source of fake biology. Community pipeline frameworks exist precisely because reproducibility across multi-omics analyses is otherwise fragile ¹.

A minimal invocation:

nextflow run nf-core/rnaseq -r 3.14.0 \
  --input samplesheet.csv \
  --outdir results \
  --gtf gencode.v44.primary_assembly.annotation.gtf.gz \
  --fasta GRCh38.primary_assembly.genome.fa.gz \
  --pseudo_aligner salmon \
  --extra_salmon_quant_args "--seqBias --gcBias --posBias" \
  -profile docker

If you would rather run the steps yourself, the fast path is fastp -i R1.fq.gz -I R2.fq.gz -o R1.trim.fq.gz -O R2.trim.fq.gz --detect_adapter_for_pe --length_required 25 followed by salmon quant -i salmon_index -l A -1 R1.trim.fq.gz -2 R2.trim.fq.gz --validateMappings --seqBias --gcBias -p 8 -o sample_quant. Build the Salmon index against the transcriptome with the full genome as a decoy (--decoys decoys.txt), which suppresses spurious assignment of intronic and intergenic reads to transcripts. On a laptop with 16 GB of RAM and eight cores, a 40 million read pair sample quantifies in roughly ten minutes.

Then collapse transcripts to genes in R with tximport, which uses effective lengths to produce an offset that DESeq2 consumes directly:

txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM")
dds <- DESeqDataSetFromTximport(txi, colData, ~ timepoint)
dds <- dds[rowSums(counts(dds) >= 10) >= 3, ]
vsd <- vst(dds, blind = FALSE)

The filter is not cosmetic. Of roughly 62,000 GENCODE gene entries, only 12,000 to 15,000 are meaningfully expressed in whole blood, and leaving the rest in inflates multiple-testing burden and destabilizes dispersion estimates. The open-source stack for differential expression and downstream functional analysis is mature and well documented, and we would not write our own ².

Counts, TPM, and the compositional problem

Use raw estimated counts for statistical modeling and TPM only for looking at a gene across samples. DESeq2 and edgeR model counts with a negative binomial distribution because RNA-seq variance exceeds its mean, and feeding them TPM breaks that model. Conversely, comparing raw counts between samples is meaningless without normalization, since library sizes routinely differ twofold.

The deeper issue is that sequencing measures proportions. It does not measure concentrations. A fixed number of reads is distributed across whatever transcripts are present, so if hemoglobin transcripts rise in one draw, every other gene’s share falls even if its absolute abundance is unchanged. This is the compositional nature of the data, and it is the single most common source of over-interpretation in personal transcriptomics ³. Practical defenses: globin depletion or rRNA depletion at library prep, median-of-ratios or TMM normalization rather than naive TPM comparison, and a centered log-ratio transform when you want to reason about relative change directly. RNA-seq has largely displaced microarrays for exactly the reasons that make it useful here, including dynamic range and the ability to detect isoforms and novel transcripts, but proportionality remains intrinsic to the measurement ⁴.

Failure modes worth checking before you trust anything

Read the MultiQC report first, and check five numbers. Mapping rate below about 80 percent against GRCh38 usually means contamination or an adapter problem, not biology. rRNA content above 10 to 15 percent means depletion failed and your effective depth is much lower than your read count suggests. Duplication above 60 percent in a low-input library signals PCR over-amplification, and duplicates should not be removed in RNA-seq, since highly expressed genes legitimately produce identical fragments. A 3’ bias in gene body coverage points to RNA degradation, quantified as RIN (RNA integrity number) at the bench: below 7, short transcripts and 5’ ends are systematically undercounted. Finally, look at the sex-chromosome genes XIST and RPS4Y1 as a cheap sample-identity check. Confirm genotype concordance between the RNA BAM and your whole-genome VCF when you have both.

For a single person sampled repeatedly, the leading confounder is cell composition. Whole blood expression is dominated by which leukocytes were circulating that morning. Deconvolution with CIBERSORTx or a reference-based method, followed by adjusting for the resulting neutrophil and lymphocyte fractions, often removes most of what looked like a signal. Time of day, fasting state, recent exercise, and acute infection all move interferon and inflammatory modules substantially. Standardize the draw conditions or you will spend your time measuring them.

What you can conclude, and when to involve a clinician

With four or more timepoints you can build a personal reference range per gene and flag deviations, score pathway modules with GSVA or singscore rather than chasing individual genes, and correlate transcript trajectories with your proteomic and glucose series. Integrating layers is where signal that is invisible in any single assay tends to appear, and it is also where statistical care matters most ⁵⁶. What you cannot do is read a diagnosis off an expression matrix. RNA-seq is a measurement. Any finding that looks clinically meaningful, including incidental variants called from the RNA reads, belongs in front of a physician or a certified clinical laboratory. Confirmation has to come first before it means anything about your health.

Questions people also ask

How do I analyze sequencing results? Start with QC using FastQC and MultiQC, then trim adapters with fastp. Quantify with Salmon against a decoy-aware transcriptome index, import to gene level with tximport, and model raw counts in DESeq2. Inspect the BAM in IGV only when a specific gene or junction is in question.

What files will a sequencing provider send? Gzipped paired FASTQ files at minimum, usually with an MD5 manifest. Many providers also return BAM plus index and a gene count matrix. Ask for the reference genome build and annotation version in writing, because the count matrix is unusable without them.

Should I use counts or TPM? Counts for any statistical test, TPM for visualizing one gene across samples or comparing genes within a sample. Never run a negative binomial model on TPM.

Do I need to align reads at all? For expression alone, no. Selective alignment with Salmon or Kallisto is faster and accurate enough. Keep a STAR BAM if you want splice junctions, allele-specific expression, or RNA variant calls.

How deep does sequencing need to be? For bulk gene-level expression in blood, 20 to 30 million read pairs is sufficient. Isoform-level work or low-abundance transcripts benefit from 50 to 100 million, and long reads change the calculus entirely.

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. Moritz Kebschull, Melanie Julia Fittler, Ryan T. Demmer, et al. Differential Expression and Functional Analysis of High-Throughput -Omics Data Using Open Source Tools. Methods in Molecular Biology, 2016. https://doi.org/10.1007/978-1-4939-6685-1_19 ↩

3. Thomas P Quinn, Ionas Erb, Greg Gloor, et al. A field guide for the compositional analysis of any-omics data. GigaScience, 2019. https://doi.org/10.1093/gigascience/giz107 ↩

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

5. Gun Tae Jung, Kwang-Pyo Kim, Kwoneel Kim. How to interpret and integrate multi-omics data at systems level. Animal Cells and Systems, 2020. https://doi.org/10.1080/19768354.2020.1721321 ↩

6. Bonnie Berger, Jian Peng, Mona Singh. Computational solutions for omics data. Nature Reviews Genetics, 2013. https://doi.org/10.1038/nrg3433 ↩