How to analyse NGS data, start to finish

Most guides to this start at the FASTQ files and walk you through commands. That is the wrong place to start, because two of the decisions that most affect whether your analysis works at all are made before the sequencer runs, and no amount of careful processing recovers from getting them wrong.

So this goes in the order the work actually happens. RNA-seq is the worked example, since it is the most common thing a lab outsources, and there is a note at the end on where other assay types diverge.

Before you sequence: the two decisions that matter

How many replicates

This is the single highest-leverage number in the whole project, and it is fixed before any data exists.

Biological replicates, not technical ones. The benchmark most people cite recommends at least six per condition, rising to twelve if you need to detect genes across all fold changes rather than only the large ones. Three per group is extremely common in published work and is genuinely underpowered for anything subtle.

If your budget forces a choice between more reads per sample and more samples, choose more samples almost every time. Depth beyond roughly 20 to 30 million reads buys very little for standard differential expression, while replicates buy statistical power directly.

What you will actually receive

Ask the provider for raw FASTQ files, not just a processed count matrix and a PDF of figures. A count matrix is somebody else's analysis decisions baked in, and you cannot revisit them two years later when a reviewer asks. Also ask which reference genome and annotation version they used, because you will need to match it.

Step 1: Do not build the pipeline yourself

The most useful advice in this entire post: use nf-core/rnaseq, currently at 3.26.0. It is a community-maintained, versioned, containerised pipeline that runs QC, alignment, quantification and a combined report, and it has been tested by more people than your hand-rolled bash script ever will be.

nextflow run nf-core/rnaseq -r 3.26.0 \
  --input samplesheet.csv \
  --outdir results \
  --genome GRCh38 \
  -profile docker

The samplesheet is a CSV listing sample names and FASTQ paths. That command is most of the job.

The rest of this post explains what that pipeline is doing, because running something you do not understand is how people end up defending results they cannot explain. The goal is comprehension, not reimplementation.

Step 2: Quality control on the raw reads

FastQC reports per-base quality, adapter content, duplicate levels and GC distribution for each file. MultiQC collapses all of those into a single report, which is the only practical way to look at more than about six samples.

What you are looking for is not perfection, it is outliers:

One sample whose quality collapses while the others are fine, which usually means a library or lane problem.
Unexpected adapter content, which tells you the insert size was shorter than the read length.
A GC distribution with a second peak, which often means contamination.
Wildly different read counts between samples, which will bite you at normalisation.

Step 3: Trimming, which you probably should not do

This surprises people, so it is worth being direct. For gene-level quantification you usually do not need to trim.

Modern aligners soft-clip: they locally ignore the parts of a read that do not match rather than requiring a full-length match. A systematic evaluation in NAR Genomics and Bioinformatics found that read trimming is not required for mapping and quantification at the gene level, and that aligners rescue many low-quality bases a trimmer would have discarded. Aggressive trimming can actively degrade differential expression results.

Trim when you have a specific reason: heavy adapter contamination visible in QC, or an application like assembly or variant calling where a mismatched base is a different kind of problem. Trimming by reflex because a tutorial said to is the common error.

Step 4: Alignment, or pseudoalignment

Two routes, and the choice is less consequential than it looks.

Alignment with STAR or HISAT2 places each read on the genome. Slower and memory-hungry, but you get a BAM you can look at in a browser, which matters if you care about novel isoforms, splice junctions, or anything positional.
Pseudoalignment with Salmon or kallisto assigns reads to transcripts without a full genomic placement. Dramatically faster, less memory, and for gene-level differential expression the results are very close.

The nf-core default runs STAR then Salmon, which gives you both the genomic QC metrics and the quantification. If you are unsure, that is the right answer.

Step 5: Quality control on the alignment

This is the step people skip, and it is where the real problems surface. Raw read quality being fine tells you very little.

Mapping rate. Below about 70 percent for a well-annotated organism means something is wrong: wrong reference, contamination, or degraded RNA.
Where reads land. A high fraction in intergenic or intronic regions can mean DNA contamination or degradation.
Duplication. High duplication with low complexity suggests too little input material and too much amplification.
rRNA fraction. If depletion failed, most of your expensive reads are ribosomal.

Step 6: Counts to results

Now you have a matrix of genes by samples. Load it into DESeq2 or edgeR. Both are well validated and, below about twelve replicates, both outperform the alternatives.

Two things to do before you look at any gene:

Plot a PCA of your samples. If they cluster by batch, processing date, or the person who prepped them rather than by condition, you have a confound, and that must go in your model rather than be ignored.
Check the dispersion plot. It should show the characteristic downward trend. If it looks strange, your design or your counts are.

Feed raw counts, never TPM or FPKM. These tools model count noise directly and normalise internally. Handing them pre-normalised values breaks their statistics silently rather than loudly.

Step 7: Reading the output honestly

Use adjusted p-values. You tested twenty thousand genes, so an unadjusted 0.05 gives you a thousand false positives by construction.
Filter on effect size as well as significance. With enough replicates a 1.02-fold change becomes significant and means nothing biologically.
Treat enrichment analysis as hypothesis generation. Pathway databases are biased toward well-studied biology, so an immune signature partly reflects immunology having been studied more.
Validate something orthogonally. A handful of qPCR reactions is cheap relative to building a project on a bioinformatics artefact.

Where other assay types diverge

The shape holds, the middle changes:

Variant calling aligns with BWA-MEM, then follows GATK best practices. Trimming matters more here, since a mismatched base is a candidate variant.
ATAC-seq and ChIP-seq replace quantification with peak calling, usually MACS. Fragment size distribution becomes a key QC metric.
Single-cell is a different pipeline entirely, with cell calling, doublet detection and ambient RNA correction before you get anywhere near a comparison.

How this usually goes wrong

In rough order of how often it happens: too few replicates, discovered after sequencing when it cannot be fixed. A batch effect perfectly confounded with the condition, so the two can never be separated. Mismatched genome and annotation versions producing quietly wrong counts. And running a pipeline nobody understood, then being unable to answer a reviewer asking why a parameter was set the way it was.

None of those are compute problems. Three of the four are decided before the data exists, which is why this post started where it did.