NGS data analysis pipelines

In short

An NGS data analysis pipeline is the ordered chain of software steps that turns raw sequencer output (millions of short reads in FASTQ files) into biological answers: a list of variants, a gene expression table, an assembled genome or a microbial profile. Every modern genomics project — clinical diagnostics, crop breeding, pathogen surveillance — depends on it, because the sequencer itself produces only base calls with quality scores. In Indian university syllabi it sits in B.Sc/M.Sc Bioinformatics and Genetics papers and in CSIR-UGC NET Life Sciences (Unit XIII), GATE BT/XL and DBT-BET question papers, usually as file formats, quality scores, coverage calculations and the standard workflow order. Employers who hire on this skill include MedGenome, Strand Life Sciences, Eurofins, Neuberg, Mapmygenome, LifeCell, Syngene, Premas Life Sciences, and the services arms of TCS, Wipro and Persistent, plus CSIR-IGIB, NIBMG, BRIC institutes and the INSACOG network.

Explained

Sequencing platforms and what their reads look like

• Second generation (short read): Illumina (NovaSeq, NextSeq, MiSeq) uses sequencing-by-synthesis with reversible dye terminators; reads of 50–300 bp, error rate ~0.1%, errors mostly substitutions. Dominant platform in India.
• Ion Torrent (semiconductor, detects H⁺ release) — fast, small panels; characteristic indel errors in homopolymer runs.
• Third generation (long read): PacBio HiFi (circular consensus, 10–25 kb, high accuracy) and Oxford Nanopore (ionic-current change as DNA passes a protein pore; reads of kb to Mb, higher raw error, real-time output, portable MinION used for field pathogen sequencing).
• Paired-end vs single-end: paired-end reads both ends of a fragment, giving an insert size that helps map repeats, detect structural variants and span gaps. Mate-pair uses long (2–20 kb) inserts for scaffolding.
• Read depth vs read length: depth (coverage) resolves base-level errors; length resolves repeats and phasing. Hybrid assemblies use both.

Stage 1 — Raw data, formats and quality control

• FASTQ: four lines per read — @identifier, sequence, +, ASCII quality string. Modern files use Phred+33 encoding (Sanger/Illumina 1.8+); old Illumina 1.3–1.7 used Phred+64 — misreading the offset shifts every score by 31.
• FastQC / MultiQC: per-base quality boxplots, per-sequence GC, adapter content, duplication level, overrepresented sequences. Quality falls at the 3′ end as phasing errors accumulate.
• Trimming and filtering: Trimmomatic, cutadapt, fastp or Trim Galore remove adapters, low-quality tails (sliding-window below ~Q20), short reads and poly-G artefacts. Over-trimming destroys real signal, so trim by adapter and quality, not by fixed length.
• Demultiplexing: bcl2fastq/BCL Convert splits the flow-cell output into per-sample FASTQ using index (barcode) reads; index hopping on patterned flow cells causes low-level cross-contamination — use unique dual indices.
• Base-call quality (Phred Q) is the central number: Q20 = 1 error in 100, Q30 = 1 in 1,000. A run is usually accepted at >80–90% of bases ≥ Q30.

Stage 2 — Alignment (mapping) and the alignment formats

• Reference-based mapping places each read on a reference genome (GRCh38/hg38 for human; the T2T-CHM13 assembly is now used in research). Burrows–Wheeler transform + FM index aligners: BWA-MEM (DNA, the clinical standard), Bowtie2. Hash/seed-extend: minimap2 for long reads.
• Splice-aware aligners for RNA-seq: STAR, HISAT2 — they allow a read to span an exon–exon junction, which a DNA aligner would clip.
• SAM (text) / BAM (binary, compressed) / CRAM (reference-based, ~30–60% smaller than BAM). Eleven mandatory SAM fields: QNAME, FLAG, RNAME, POS, MAPQ, CIGAR, RNEXT, PNEXT, TLEN, SEQ, QUAL.
• CIGAR string encodes the alignment: 100M = 100 matched/mismatched, S soft-clip, I insertion, D deletion, N skipped region (intron).
• MAPQ = Phred-scaled probability that the read is mis-placed; MAPQ 0 usually means the read maps equally well to several places (repeat/multi-mapper).
• **samtools** does sort, index (.bai), merge, view, flagstat, depth; IGV visualises the BAM.
• Mark duplicates (Picard MarkDuplicates, samtools markdup): PCR and optical duplicates inflate apparent depth and bias allele fractions. In amplicon/UMI data use UMI-aware deduplication (fgbio, UMI-tools) instead — otherwise genuine molecules are discarded.
• BQSR (base quality score recalibration) and indel realignment (legacy) correct systematic base-quality errors before variant calling.

Stage 3 — Variant calling, filtering and annotation (DNA-seq / WES / WGS)

1. Germline calling: GATK HaplotypeCaller (local de novo haplotype assembly over active regions), bcftools mpileup/call, FreeBayes, DeepVariant (CNN on pileup images). GVCF mode allows joint genotyping of cohorts (GenomicsDBImport → GenotypeGVCFs).
2. Somatic (tumour) calling: Mutect2, VarScan2, Strelka2 — tumour-vs-normal pairs, low VAF calls, needs higher depth (500–1,000× for panels).
3. VCF format: header ## lines, then #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT sample; FORMAT keys GT (0/0, 0/1, 1/1), DP (depth), AD (allele depths), GQ (genotype quality). bgzip + tabix for indexed access; bcftools/vcftools for filtering and merging.
4. Filtering: VQSR (machine-learning model on annotation profiles, needs large cohorts) or hard filters (QD, FS, MQ, SOR, depth, strand bias) for small sets.
5. Structural and copy-number variants: Manta, DELLY, Lumpy (SV); CNVkit, GATK gCNV, ExomeDepth (CNV from read depth); long reads call SVs far better.
6. Annotation: VEP, ANNOVAR, SnpEff add gene, transcript, consequence (missense, nonsense, frameshift, splice), population frequency (gnomAD; IndiGenomes for Indian allele frequencies), in-silico pathogenicity (SIFT, PolyPhen-2, CADD, REVEL) and clinical databases (ClinVar, OMIM, HGMD).
7. Clinical interpretation follows the ACMG/AMP 2015 five-tier system — pathogenic, likely pathogenic, VUS, likely benign, benign — with evidence codes (PVS1, PS1–4, PM1–6, PP1–5, BA1, BS1–4, BP1–7).

Stage 4 — Pipelines for the other assay types

• RNA-seq (expression): QC → splice-aware alignment (STAR/HISAT2) → counting (featureCounts, HTSeq-count against a GTF/GFF3) → normalisation → differential expression with DESeq2 or edgeR (negative binomial model, shrinkage of dispersion), then multiple-testing correction. Alignment-free quantifiers salmon and kallisto estimate transcript abundance directly and are much faster.
• Normalisation units: raw counts for DE tools; CPM, RPKM/FPKM and TPM for comparing genes within a sample; TPM is preferred because its total is fixed at 1,000,000 per sample.
• Single-cell RNA-seq: Cell Ranger/STARsolo → barcode and UMI counting → QC (mitochondrial %, gene count) → normalisation, PCA, graph clustering, UMAP, marker genes (Seurat in R, Scanpy in Python).
• ChIP-seq / ATAC-seq: align → filter duplicates and blacklist regions → peak calling with MACS2/MACS3 against an input control → BED files → motif analysis (HOMER, MEME) and annotation (ChIPseeker).
• De novo assembly: de Bruijn graph assemblers for short reads (SPAdes for bacteria, MEGAHIT, SOAPdenovo), OLC/string-graph for long reads (Canu, Flye, hifiasm); scaffolding, gap filling, polishing (Pilon, Racon, Medaka); quality by QUAST (N50, L50, contig count) and completeness by BUSCO.
• Metagenomics: 16S rRNA amplicon (QIIME2, DADA2 → ASVs → diversity) versus shotgun metagenomics (Kraken2/Bracken for taxonomy, MetaPhlAn, HUMAnN for function, MAG binning with MetaBAT2).
• Genome surveillance: amplicon tiling (ARTIC protocol) → consensus with iVar/Medaka → lineage assignment (Pangolin, Nextclade) — the workflow INSACOG laboratories used for SARS-CoV-2 and now for other pathogens.

Reproducibility, infrastructure and Indian practice

• Workflow managers: Nextflow (with nf-core community pipelines: nf-core/sarek, nf-core/rnaseq), Snakemake, WDL/Cromwell, CWL. They give resumable runs, parallel scheduling on HPC/cloud and a record of every command.
• Environment control: conda/mamba, Docker and Singularity/Apptainer containers pin tool versions; a pipeline without pinned versions is not reproducible.
• Compute: WGS at 30× is ~90–120 GB of FASTQ per sample; alignment needs 16–64 GB RAM, STAR needs ~32 GB for the human index, a bacterial SPAdes assembly runs on a laptop. Indian academic users commonly use NSM supercomputers (PARAM series), institute clusters or cloud.
• Indian data submission and governance:
  • IBDC (Indian Biological Data Centre, BRIC-RCB, Faridabad) is India's national repository — INSDC-linked, with INSDC-style submission of raw reads (BioProject/BioSample/SRA-equivalent) and the IndiGen/IndiGenomes variant resource from CSIR-IGIB.
  • GenomeIndia (DBT) completed whole-genome sequencing of ~10,000 individuals across Indian populations, giving a reference set of Indian variation; IndiGen (CSIR) sequenced 1,029 genomes.
  • Research on human participants follows the ICMR National Ethical Guidelines for Biomedical and Health Research Involving Human Participants with Institutional Ethics Committee approval; access to Indian biological resources is governed by the Biological Diversity Act, 2002 as amended in 2023 (NBA/SBB approvals).
  • Diagnostic NGS labs work to NABL (ISO 15189) accreditation; NGS-based IVD kits are regulated by CDSCO under the Medical Devices Rules, 2017.
  • Personal data: genomic data is personal data under the DPDP Act, 2023. Note the sequencing of duties: under the DPDP Rules, 2025 the Data Fiduciary obligations (notice, consent, breach intimation, erasure) and the omission of IT Act s.43A take effect only from 13 May 2027; until then IT Act s.43A and the SPDI Rules, 2011 apply.

Key points

• Pipeline order - QC → trim → align → sort/dedup → recalibrate → call variants → filter → annotate; the order is itself an exam answer.
• FASTQ - read name, sequence, +, quality; Phred+33 encoding in all current data.
• Phred score - Q = −10 log₁₀(P of wrong base); Q30 = 99.9% accuracy.
• Coverage (depth) - average times a base is read; ~30× for germline WGS, 100× for WES, 500–1,000× for somatic panels, 50–100× for de novo bacterial assembly.
• SAM/BAM/CRAM - same alignment information as text, binary and reference-compressed; BAM must be coordinate-sorted and indexed.
• CIGAR and MAPQ - how the read aligns, and how confident the placement is; MAPQ 0 = multi-mapping.
• Duplicate marking - removes PCR/optical duplicates; never deduplicate UMI or amplicon data by position alone.
• GATK Best Practices - the reference germline workflow (BWA-MEM → MarkDuplicates → BQSR → HaplotypeCaller GVCF → joint genotyping → VQSR).
• VCF - one row per variant site; GT/DP/AD/GQ carry the genotype evidence.
• Annotation ≠ interpretation - VEP/ANNOVAR add facts; ACMG/AMP criteria assign clinical class.
• RNA-seq - counts for DESeq2/edgeR, TPM for cross-gene comparison; salmon/kallisto skip alignment.
• Peak calling - MACS2 for ChIP-seq/ATAC-seq, always against an input control.
• Assembly metrics - N50, L50, contig count, BUSCO completeness; de Bruijn for short reads, OLC for long reads.
• Long reads - resolve repeats, phase haplotypes and call structural variants that short reads miss.
• Batch effect and reference version - never mix hg19 and hg38 coordinates or sequencing batches without correction.
• Multiple testing - thousands of genes tested, so report adjusted p (FDR) by Benjamini–Hochberg, not raw p.
• Reproducibility - Nextflow/Snakemake plus containers and pinned versions; log every parameter.
• Indian infrastructure - IBDC for deposition, GenomeIndia and IndiGenomes for Indian allele frequencies, NABL/CDSCO for clinical use.

Formulas and examples

Example 1 — Coverage of a human WGS run. 400 million 150 bp paired-end reads (i.e. 400 million read *pairs*? here take 400 million single reads of 150 bp), G = 3.1 × 10⁹ bp.

• Total bases = 400 × 10⁶ × 150 = 6.0 × 10¹⁰ bp = 60 Gb.
• C = 6.0 × 10¹⁰ / 3.1 × 10⁹ ≈ 19.4× — below the 30× clinical germline norm, so about 9.3 × 10¹⁰ bases (93 Gb) would be needed for 30×.
• Uncovered bases at 19.4× = 3.1 × 10⁹ × e^−19.4 ≈ 3.1 × 10⁹ × 3.75 × 10⁻⁹ ≈ 12 bases (random model only; real gaps are far more, because coverage is not Poisson in repeats and GC extremes).

Example 2 — Data required for a clinical exome. Target 50 Mb, desired mean depth 100×, on-target 70%, duplicates 10%, mapped 98%.

• Useful bases needed = 50 × 10⁶ × 100 = 5.0 × 10⁹ bp.
• Raw bases = 5.0 × 10⁹ / (0.70 × 0.90 × 0.98) = 5.0 × 10⁹ / 0.6174 ≈ 8.1 × 10⁹ bp (≈8.1 Gb).
• As 150 bp PE reads: 8.1 × 10⁹ / 150 ≈ 54 million reads ≈ 27 million read pairs.

Example 3 — Phred scores in practice.

• Q30 → P = 10⁻³ → 1 wrong base in 1,000. A 150 bp read entirely at Q30 carries 150 × 0.001 = 0.15 expected errors.
• Q20 → P = 0.01 → the same read carries 1.5 expected errors; this is why 3′ tails below Q20 are trimmed.
• A site with depth 30 and a true heterozygous variant: chance that all 30 reads are the reference allele = 0.5³⁰ ≈ 9.3 × 10⁻¹⁰ — negligible, which is why 30× suffices for germline SNV calling but not for a 5% somatic VAF.

Example 4 — TPM from counts.

Scaling factor = ΣRPK / 10⁶ = 1,400 / 10⁶ = 0.0014. TPM_A = 500 / 0.0014 = 357,142.9; TPM_B = 600 / 0.0014 = 428,571.4; TPM_C = 300 / 0.0014 = 214,285.7. Note gene B has the fewest reads but the highest expression per transcript — the point of length normalisation.

Example 5 — N50 and L50 of an assembly. Contigs: 500, 400, 300, 200, 100 kb. Total = 1,500 kb; half = 750 kb.

• Cumulative (largest first): 500 (<750), 900 (≥750) → N50 = 400 kb, L50 = 2.
• If the same 1,500 kb came as 15 contigs of 100 kb each, N50 = 100 kb and L50 = 8 — identical total, far worse contiguity.

Example 6 — VAF interpretation at depth 100. AD = 72,28 → VAF = 28/100 = 0.28. Too low for a clean germline heterozygote (0.5): consistent with a somatic/subclonal variant, tumour admixture, a CNV, or allele dropout — flag for orthogonal confirmation.

Where it is used

• Clinical genomics and diagnostics: carrier screening, hereditary cancer and rare-disease exomes, NIPT, oncology hotspot panels and liquid biopsy at MedGenome, Strand Life Sciences, Neuberg, Eurofins, LifeCell, Mapmygenome and hospital genetics labs (AIIMS, Tata Memorial, CMC Vellore) — roles as NGS bioinformatics analyst and variant curation scientist.
• Public health genomics: INSACOG/NCDC pathogen genome surveillance, AMR gene detection from bacterial isolates, outbreak phylogenomics in ICMR institutes and state labs.
• Agri-genomics and animal breeding: ICAR institutes, NBPGR and seed companies use GBS/WGS pipelines for marker-assisted selection, QTL mapping and genomic prediction in rice, wheat, cotton, cattle and buffalo.
• Pharma and biotech R&D: Syngene, Biocon, Dr Reddy's and CROs use RNA-seq and single-cell pipelines for target discovery, biomarker work, cell-line genome verification and biosimilar characterisation.
• IT services and product teams: TCS, Wipro, Persistent, Infosys and Illumina/Agilent/Thermo Fisher India build and validate cloud pipelines, LIMS integration and reporting software — roles for pipeline/DevOps engineers with Nextflow, Docker, AWS/Azure and Python/R.
• Academia and national resources: CSIR-IGIB, NIBMG Kalyani, BRIC institutes, IISc, IITs and DBT centres; data deposition to IBDC and the use of IndiGenomes/GenomeIndia Indian allele frequencies to filter common local variants out of candidate lists.

Interview questions

• Q: What are the four lines of a FASTQ record? A: A header beginning with @, the base sequence, a + separator, and an ASCII quality string of the same length as the sequence.
• Q: What does Q30 mean, and why is a run judged by %Q30? A: Q30 means a 1-in-1,000 chance that the base call is wrong (99.9% accuracy); the fraction of bases at Q30 or above is the single best summary of run quality, with >80–90% expected on Illumina.
• Q: Difference between a BAM and a CRAM file? A: Both store alignments; BAM is BGZF-compressed binary, while CRAM compresses against the reference sequence and so is substantially smaller but needs that exact reference to decode.
• Q: Why must duplicates be marked before variant calling, and when must they not be removed? A: PCR/optical duplicates are copies of one molecule, so they falsely inflate depth and allele fractions; in UMI-tagged or amplicon data reads share start positions by design, so deduplicate by UMI, not by coordinates.
• Q: What is MAPQ 0 telling you? A: The read maps equally well to more than one location (typically a repeat or segmental duplication), so its placement is unreliable and most callers ignore it.
• Q: Why is a splice-aware aligner needed for RNA-seq? A: Reads crossing exon–exon junctions have no contiguous match on the genome; STAR/HISAT2 split them with an N (skipped) CIGAR operation, whereas BWA would soft-clip or mis-map them.
• Q: Why does DESeq2 take raw counts rather than TPM? A: Its negative binomial model estimates dispersion from count variance, which length/library-normalised values destroy; DESeq2 does its own median-of-ratios size-factor normalisation.
• Q: TPM versus FPKM — which would you use to compare two genes in one sample and why? A: TPM, because TPM values sum to one million in every sample, so the proportion of transcripts is directly comparable; FPKM totals vary between samples.
• Q: What is N50 and what is its weakness as a metric? A: The contig size above which half the assembly lies; it rewards a few long contigs and can be inflated by mis-joins, so it must be read with contig count, total length, and BUSCO completeness.
• Q: 30× is standard for germline WGS but a tumour panel is run at 1,000× — why? A: A germline heterozygote appears at ~50% VAF so modest depth suffices; somatic variants in impure tumours may be at 1–5% VAF, needing hundreds of reads to distinguish from sequencing noise.
• Q: How would you make an NGS pipeline reproducible for an audit? A: Encode it in Nextflow/Snakemake, pin tool versions in a container, record the reference build and all parameters in a versioned config, and keep per-sample run logs and checksums.
• Q: What legal and ethical requirements apply to human genomic data in India? A: IEC approval and informed consent under the ICMR National Ethical Guidelines, NBA/SBB clearance for Indian biological resources under the Biological Diversity Act as amended in 2023, NABL (ISO 15189) accreditation and CDSCO regulation for diagnostic use, and the DPDP Act, 2023 — with the Data Fiduciary duties under the DPDP Rules, 2025 effective from 13 May 2027, the IT Act s.43A/SPDI Rules applying until then.

Common mistakes

• Confusing per-base quality (Q) with MAPQ or with coverage depth — three different confidences.
• Mixing reference builds: calling on hg38 but annotating with hg19 coordinates, without liftover, which silently mislabels every variant.
• Quoting "coverage" without saying whether it is mean target depth or the percentage of target above a threshold (e.g. %bases ≥20×), which is what clinical reports actually need.
• Treating annotation output as a diagnosis — a VEP consequence or a high CADD score is not a pathogenicity classification; ACMG/AMP evidence weighing is.
• Reporting raw p-values from thousands of genes instead of FDR-adjusted values, and calling fold-change significant without any statistical test or replicates.
• Using TPM/FPKM as input to DESeq2/edgeR, or comparing counts across samples without library-size normalisation.
• Aggressive uniform trimming and filtering (fixed length cuts, deduplicating amplicons, discarding all multi-mappers) that removes genuine biological signal; and running pipelines with unpinned tool versions, so results cannot be reproduced.

Related topics

Sequence alignment: BLAST and multiple alignment · Genomics and proteomics · Biological databases: NCBI, UniProt and PDB · Bioinformatics basics · Phylogenetic analysis · Molecular biology techniques: PCR, cloning and sequencing