How to Analyze Your Own RNA-Seq Data

By the end of this you will have a gene-level expression matrix for your own sample or samples. You will have a set of QC numbers that tell you whether to trust it. You will also have an estimate of the blood cell composition driving most of the signal, and a defensible way to compare one timepoint against another. You need paired-end FASTQ files and their library metadata, about 32 GB of RAM, 100 GB of disk, and a machine with 8 or more cores. The metadata should say whether the library is poly-A selected or ribo-depleted and whether it is stranded or unstranded. It should also give read length and whether globin was depleted. Everything below runs on Linux or macOS with conda, R 4.3+, and about three hours of wall time per sample. If you also have whole-genome data on the same person, keep it handy: the most interesting things you can do with a personal transcriptome come from joining it to your genotypes.

One framing point before the commands. RNA-seq measures relative abundance of transcripts in whatever tissue was sampled, at the moment it was sampled. For a personal profile that tissue is almost always whole blood or PBMCs, and blood transcriptomes are dominated by cell composition. A twofold change in a neutrophil gene usually means you had more neutrophils that morning. It does not usually mean the gene was regulated differently. Plan your analysis around that fact. Nothing here is diagnostic. If a result makes you want to change something about your health, that is a conversation with a physician, and the raw data is the thing to bring.

1. Read the library metadata before touching the data

Three fields determine everything downstream.

Selection chemistry. Poly-A selection captures polyadenylated mRNA and gives you the cleanest per-gene signal per dollar. Ribosomal depletion keeps the non-polyadenylated species, which are lncRNAs, histone mRNAs, pre-mRNA and circRNA precursors. Ribo-depleted libraries typically show 20 to 40 percent intronic reads, which is normal and not a degradation signal by itself. The chemistry also changes what fraction of your reads land in annotated exons, so do not compare a poly-A run to a ribo-depleted run as if they were the same assay. Technology choice has a large effect on which transcript classes you can even see, and short-read versus long-read platforms differ further in isoform resolution ¹.

Strandedness. Most modern kits (dUTP-based, e.g. Illumina Stranded mRNA) are reverse-stranded, which is ISR in Salmon notation and -s 2 in featureCounts. Getting this wrong quietly halves your counts and inverts antisense genes.

Globin. In whole blood without globin depletion, HBB, HBA1, and HBA2 can swallow a majority of your reads, leaving you with an effective depth far below the nominal one. Check it first thing (step 4) and adjust your expectations about low-expression genes accordingly.

Also record input RNA quality if you have it (RIN or DV200). Degraded input shifts coverage toward 3’ ends and inflates the apparent expression of short transcripts.

2. QC and trim the raw reads

conda create -n rnaseq -c bioconda -c conda-forge \
  fastqc multiqc fastp salmon=1.10.3 star samtools subread
conda activate rnaseq

fastqc -t 8 -o qc/ sample_R1.fastq.gz sample_R2.fastq.gz

fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
  -o trim/sample_R1.fq.gz -O trim/sample_R2.fq.gz \
  --detect_adapter_for_pe --trim_poly_g --cut_tail --cut_tail_mean_quality 20 \
  --length_required 36 --thread 8 \
  --html qc/sample.fastp.html --json qc/sample.fastp.json

--trim_poly_g matters on two-color chemistry (NovaSeq, NextSeq), where an absent signal is read as G and creates runs of poly-G at read ends. --detect_adapter_for_pe finds the adapter by overlap rather than requiring you to name it.

Be conservative with quality trimming. Aggressive base trimming shortens reads, reduces mappability, and biases quantification more than the low-quality bases would have. Adapter removal plus a light 3’ trim is the setting we use. The classic pipelines put heavy weight on trimming because they were built for aligners with strict end-to-end matching ². Pseudo-alignment and soft-clipping aligners are less sensitive to it.

Targets to check in the fastp JSON: duplication rate (high duplication with low complexity suggests too little input RNA), Q30 fraction above 85 percent, and adapter content below a few percent after trimming.

3. Quantify with Salmon against a decoy-aware index

We recommend Salmon with a decoy-aware transcriptome index over a full genome alignment for gene-level expression. It runs in about 10 minutes per sample on 16 threads, corrects for GC and sequence bias, and gives per-transcript estimates with uncertainty. The tradeoff: you get no BAM. That means no coverage plots, no splice junction inspection and no variant calling from RNA. If you want those, run STAR too (step 7).

Build the index once:

wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_47/gencode.v47.transcripts.fa.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_47/GRCh38.primary_assembly.genome.fa.gz

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 -i idx_gencode47 \
  -k 31 -p 16 --gencode

The decoys are the genome sequence. Without them, reads from intronic or intergenic regions that happen to resemble a transcript get assigned to that transcript, and you see phantom expression in genes you do not express. The --gencode flag strips the pipe-delimited GENCODE FASTA headers down to transcript IDs. Use -k 31 for reads of 75 bp or longer, -k 23 if you are stuck with 50 bp reads.

Quantify:

salmon quant -i idx_gencode47 -l A \
  -1 trim/sample_R1.fq.gz -2 trim/sample_R2.fq.gz \
  --seqBias --gcBias --posBias \
  --numBootstraps 30 --validateMappings \
  -p 16 -o quant/sample

-l A infers library type. Check quant/sample/lib_format_counts.json afterward and confirm it inferred ISR (or whatever your kit is). If it reports IU for a stranded kit, something is wrong with your read pairing or your metadata. --gcBias is worth the runtime on any library where fragment GC could differ across batches, which is every longitudinal profile. The bootstraps give you per-transcript inferential replicates, which matter for isoform-level work where reads are ambiguous between transcripts of the same gene ³.

4. Collapse to genes and run the first sanity checks

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

# tx2gene from the GENCODE GTF, or parse the transcript FASTA headers
tx2gene <- read_tsv("tx2gene_gencode_v47.tsv",
                    col_names = c("tx_id", "gene_id", "gene_name"))

files <- c(t1 = "quant/t1/quant.sf", t2 = "quant/t2/quant.sf")
txi <- tximport(files, type = "salmon", tx2gene = tx2gene[,1:2],
                countsFromAbundance = "lengthScaledTPM")

cts <- txi$counts
tpm <- txi$abundance

Use lengthScaledTPM if the counts are headed for a count-based model like DESeq2 without a separate offset. Use raw counts plus txi$length as an offset if you prefer DESeq2’s DESeqDataSetFromTximport, which handles the length correction internally.

Now the checks that catch problems:

# 1. Globin burden (whole blood)
globin <- c("HBB","HBA1","HBA2","HBD")
sum(tpm[tx2gene$gene_name[match(rownames(tpm), tx2gene$gene_id)] %in% globin, ]) / 1e6

# 2. Mitochondrial fraction: high values suggest degraded or stressed input
# 3. Genes detected at TPM > 1
colSums(tpm > 1)

Rules of thumb for whole blood, poly-A, 30M read pairs start with 12,000 to 15,000 genes at TPM greater than 1. Globin is under 10 percent if depleted and possibly over 50 percent if not. Mitochondrial genes are under 15 percent. Mapping rate from Salmon should be above 80 percent for poly-A libraries. Below 60 percent, suspect contamination, rRNA carryover, or the wrong index.

On depth: 20 to 30 million read pairs is enough for gene-level expression of moderately and highly expressed genes. Isoform-level work, fusion detection, or low-abundance transcripts want 50 to 100 million. Depth and library prep both change which genes clear the detection threshold. In small-sample designs they change the conclusions you draw, and the effect does not stop at the error bars ⁴.

5. Decide what question a single sample can answer

With one sample you cannot do differential expression. There is no dispersion to estimate and no null to compare to. What you can do:

Rank yourself against a reference. Convert your TPMs to ranks and compare against GTEx whole blood medians for the same tissue. Percentile ranks are far more stable across protocols than absolute TPMs. Anything you find is a hypothesis for the next timepoint. It is not a finding.

Estimate cell composition. This explains most of the variance in blood and is worth doing before anything else. CIBERSORTx with the LM22 signature matrix on your TPM table, or xCell in R, will give fractions for 22 leukocyte subsets. Run it, then ask whether the genes that looked interesting are simply markers of a subset whose fraction moved.

Check identity. If you also have WGS, call genotypes at common SNPs from an RNA BAM and compare to your DNA VCF with bcftools gtcheck. Discordance above a few percent means a sample swap. This is the single highest-value five-minute check in the whole pipeline.

6. Do the comparison you can do: two or more timepoints on yourself

With three or more timepoints, you have a real design. Paired within-person comparisons remove the largest confounder in transcriptomics, which is between-person variation.

library(DESeq2)
dds <- DESeqDataSetFromTximport(txi, colData = meta,
                                design = ~ neutrophil_frac + condition)
keep <- rowSums(counts(dds) >= 10) >= 3
dds <- dds[keep, ]
dds <- DESeq(dds)
res <- lfcShrink(dds, coef = "condition_post_vs_pre", type = "apeglm")

Put estimated cell fractions in the design as covariates. Without them, a shift in neutrophil percentage reads out as hundreds of “differentially expressed” genes. Filter low counts before testing: genes with fewer than 10 counts in fewer than the smallest group size contribute noise and cost you multiple-testing power. Use apeglm shrinkage so that low-count genes do not top your fold-change ranking.

Exploratory view first, always:

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

If PC1 separates by sequencing batch or by collection date rather than by your variable of interest, stop and deal with that before interpreting anything. Batch and biology confounded in a personal time series is the most common way these analyses go wrong ⁵.

7. When to run STAR instead of, or alongside, Salmon

Run an alignment when you need coordinates:

STAR --runThreadN 16 --genomeDir star_idx_gencode47 \
  --readFilesIn trim/sample_R1.fq.gz trim/sample_R2.fq.gz \
  --readFilesCommand zcat \
  --outSAMtype BAM SortedByCoordinate \
  --outSAMattributes NH HI AS nM NM MD \
  --quantMode GeneCounts \
  --outFilterMultimapNmax 20 --alignSJoverhangMin 8 \
  --outFileNamePrefix star/sample.

Build the index with --sjdbOverhang set to read length minus one (99 for 100 bp reads). STAR needs roughly 32 GB of RAM for GRCh38.

The BAM gets you: allele-specific expression at heterozygous sites from your WGS (GATK ASEReadCounter, with WASP filtering to avoid reference-mapping bias), splice junction counts for checking whether a splice-region variant in your genome has a transcript-level consequence, fusion calling with Arriba or STAR-Fusion, and coverage plots to eyeball whether a gene’s signal is real or a single pile-up artifact. A reference-guided alignment is also the starting point if you want to assemble novel isoforms rather than only quantify the annotation ⁶.

8. Functional interpretation without fooling yourself

Use rank-based GSEA over a fixed gene set collection rather than a hypergeometric test on a thresholded gene list. Thresholded overrepresentation throws away the effect sizes and is sensitive to where you set the cutoff.

library(fgsea); library(msigdbr)
h <- msigdbr(species = "Homo sapiens", collection = "H")
pathways <- split(h$gene_symbol, h$gs_name)
stats <- res$log2FoldChange; names(stats) <- res$gene_name
stats <- sort(stats[!is.na(stats)], decreasing = TRUE)
fg <- fgsea(pathways, stats, minSize = 15, maxSize = 500, nPermSimple = 10000)

Start with the 50 Hallmark sets. They are curated to be non-redundant, which keeps you from reporting fifteen versions of “inflammation.” Reactome and GO are useful next, with the caveat that GO terms overlap heavily and produce correlated hits that look like independent confirmation and are not. Enrichment results are downstream of every choice you made in quantification and filtering, so record the versions of the annotation and the gene set collection with the result ⁷. Note also that pathway databases are biased toward well-studied genes, which means your enrichment output partly reflects the history of biomedical research rather than your biology.

If you have other omics layers on the same person, the integration is where the value is. Joint models across genotype, expression, and protein data recover relationships that no single layer shows ⁸.

Common problems

Salmon mapping rate of 40 percent with a clean FASTQ. Almost always the wrong index (mouse instead of human, or a cDNA-only index without decoys) or heavy rRNA content in a ribo-depleted library that underperformed. Check meta_info.json for percent_mapped and run one million reads through a contamination screen before rebuilding anything.

Library type inferred as IU when you ordered a stranded kit. Look at lib_format_counts.json. If the ISF and ISR counts are roughly equal, the strandedness was lost, and you should quantify with -l IU and accept that antisense genes and overlapping loci will be unreliable.

Thousands of significant genes between two of your own timepoints. Check cell fractions first. A three-point swing in neutrophil percentage will do this. Add composition covariates and rerun before believing anything.

All fold changes concentrated in short or long genes. Sign of degraded input or missing GC bias correction. Plot log fold change against transcript length and against GC content. A trend means you are looking at a technical artifact.

Different gene counts every time you rerun. You changed annotation versions. GENCODE adds and retires genes at every release and changes transcript sets more often than that. Pin the release in your pipeline and record it alongside the matrix.

Version drift generally. Record tool versions, index checksums, and parameters with every quantification run. Public omics repositories store processed data next to full pipeline metadata for exactly this reason, and a personal time series is useless if timepoint 1 and timepoint 4 went through different code ⁹. Keep the trimming and quantification steps in a single script and rerun the whole series when you upgrade anything ¹⁰.

Reads that map nowhere in high volume. Run a small sample through a taxonomic classifier. Bacterial or phage sequence in blood RNA usually means a collection or prep issue. It is not usually a finding.

Footnotes

1. Maria Koutra, Paschalis Evangelidis, Apostolia Papalexandri, et al. RNA Sequencing Technologies in Acute Lymphoblastic Leukemia: A Comparative Technical Review. Current Issues in Molecular Biology, 2026. https://doi.org/10.3390/cimb48080768 ↩

2. M. Lohse, A. M. Bolger, A. Nagel, et al. RobiNA: a user-friendly, integrated software solution for RNA-Seq-based transcriptomics. Nucleic Acids Research, 2012. https://doi.org/10.1093/nar/gks540 ↩

3. Jochen B. W. Wolf. Principles of transcriptome analysis and gene expression quantification: an RNA ‐seq tutorial. Molecular Ecology Resources, 2013. https://doi.org/10.1111/1755-0998.12109 ↩

4. Dongying Li, Binsheng Gong, Joshua Xu, et al. Impact of Sequencing Depth and Library Preparation on Toxicological Interpretation of RNA-Seq Data in a “Three-Sample” Scenario. Chemical Research in Toxicology, 2020. https://doi.org/10.1021/acs.chemrestox.0c00368 ↩

5. Michele Araújo Pereira, Eddie Luidy Imada, Rafael Lucas Muniz Guedes. RNA‐seq: Applications and Best Practices. Applications of RNA-Seq and Omics Strategies - From Microorganisms to Human Health, 2017. https://doi.org/10.5772/intechopen.69250 ↩

6. Pierre De Wit, Melissa H. Pespeni, Jason T. Ladner, et al. The simple fool’s guide to population genomics via RNA‐Seq: an introduction to high‐throughput sequencing data analysis. Molecular Ecology Resources, 2012. https://doi.org/10.1111/1755-0998.12003 ↩

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

8. Chen Wang, Zhijie You, Siqi Chen, et al. AI-enabled multi-omics pharmacogenomic modeling guides resistance-aware multitarget optimization of venetoclax therapy in acute myeloid leukemia. npj Precision Oncology, 2026. https://doi.org/10.1038/s41698-026-01608-6 ↩

9. Shayoni Ray, Samrawit Gebre, Homer Fogle, et al. GeneLab: Omics database for spaceflight experiments. Bioinformatics, 2018. https://doi.org/10.1093/bioinformatics/bty884 ↩

10. Rowayna Shouib, Gary Eitzen, Rineke Steenbergen. A Guide to Basic RNA Sequencing Data Processing and Transcriptomic Analysis. BIO-PROTOCOL, 2025. https://doi.org/10.21769/bioprotoc.5295 ↩