Building an NGS Bioinformatics Pipeline for Your Own Genome and Transcriptome

By the end of this guide you will have a reproducible pipeline that takes your raw sequencing reads and produces a CRAM alignment, a filtered and annotated small-variant VCF, structural variant and repeat-expansion calls, a transcript-level expression matrix, and a validation report showing how accurate your own pipeline is on a benchmark sample. You need the raw FASTQ files from your sequencing provider. For 30x whole-genome sequencing on Illumina, expect roughly 60 to 100 GB gzipped. You also need a machine with at least 32 GB of RAM and about 1 TB of free disk, Docker or Singularity, and either Nextflow or Snakemake. A 16-core workstation will run a 30x genome end to end in roughly a day. On a 64-core cloud instance the alignment and calling steps finish in a few hours. Everything below assumes Illumina short reads, which is what “NGS” usually means in practice, since Illumina’s sequencing-by-synthesis chemistry dominates the installed base. The same pipeline structure applies to Element or Ultima data with different base-quality behavior. Long-read data needs a different aligner and caller entirely.

One framing point before the commands, because it determines what your data can and cannot tell you. A PCR assay amplifies a locus you chose in advance and answers a question you already asked: is this variant present, is this pathogen’s genome in this sample. Sequencing is untargeted within its capture space, so it answers questions you had not thought to ask, at the cost of far more data handling and a real false-positive rate in hard regions. Reviews of genomic detection methods consistently frame this as a sensitivity-and-cost tradeoff against breadth of discovery, with PCR winning on cost per answer and NGS winning on how many answers are available at all.¹ The two combine well: multiplex PCR enrichment followed by sequencing gives you PCR’s sensitivity with the ability to identify what you amplified rather than just detect it.² And depth matters more than people expect. Ultra-deep sequencing has been used to find the handful of postzygotic mutations distinguishing monozygotic twins, differences no standard PCR-based profile can resolve.³

1. Inventory the raw data and check it before you spend compute on it

Start by confirming you received what you paid for, because every downstream problem is cheaper to catch here. Verify checksums against the provider’s manifest. Then confirm read counts, read length, and lane structure directly from the files.

md5sum -c manifest.md5
zcat sample_R1.fastq.gz | head -n 400000 | awk 'NR%4==2 {print length($0)}' | sort -n | uniq -c
seqkit stats -a -j 8 sample_R*.fastq.gz

Then run a quality pass. We use fastp rather than FastQC plus a separate trimmer, because it reports and trims in one pass over the data and writes a JSON summary you can assert on in a test:

fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
  -o trim_R1.fastq.gz -O trim_R2.fastq.gz \
  --detect_adapter_for_pe --cut_tail --cut_tail_mean_quality 20 \
  --length_required 50 --thread 8 \
  --json fastp.json --html fastp.report.html

For whole-genome data, resist heavy trimming. Modern aligners soft-clip adapters well, and aggressive quality trimming biases coverage in low-complexity regions. The numbers we check are duplication rate (above roughly 15% on a 30x library suggests a low-input library amplified too hard), the fraction of bases above Q30 (expect above 90% on a healthy NovaSeq run), and GC distribution (a second mode near 40% GC often means bacterial contamination).

2. Choose a reference and stop changing it

Pick GRCh38 in the no-alt analysis set form: primary chromosomes, unlocalized and unplaced contigs, the EBV decoy, and HLA sequences, with alternate haplotypes excluded. Including ALT contigs without an ALT-aware post-processing step scatters multi-mapping reads and silently drops mapping quality to zero across immune and other polymorphic loci. T2T-CHM13 is the better assembly by every structural metric, but the clinical annotation ecosystem is still centered on GRCh38. That ecosystem includes ClinVar coordinates, gnomAD allele frequencies and most pharmacogenomic star-allele definitions. We therefore align to GRCh38 as the primary analysis and treat CHM13 as a secondary alignment when a specific question involves a region GRCh38 represents badly.

wget https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa
samtools faidx GRCh38_full_analysis_set_plus_decoy_hla.fa
bwa-mem2 index GRCh38_full_analysis_set_plus_decoy_hla.fa   # ~80 GB RAM to build, ~10 GB index

Record the reference MD5 in your pipeline config. Mixing a CRAM written against one reference build with a VCF annotated against another is the single most common way personal genomics analyses go quietly wrong.

3. Align, mark duplicates, and store as CRAM

bwa-mem2 is our default: same output as BWA-MEM within rounding, roughly twice as fast. Set the read group at alignment time, because tools downstream will refuse to run without it.

bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:run1.lane1\tSM:SUBJECT01\tLB:lib1\tPL:ILLUMINA\tPU:run1.lane1' \
  GRCh38_full_analysis_set_plus_decoy_hla.fa trim_R1.fastq.gz trim_R2.fastq.gz \
| samtools fixmate -m -u -@ 4 - - \
| samtools sort -u -@ 8 -m 2G -T /scratch/sort - \
| samtools markdup -@ 8 --reference GRCh38_full_analysis_set_plus_decoy_hla.fa \
  -f markdup.stats.txt - SUBJECT01.cram
samtools index SUBJECT01.cram

The -K 100000000 flag fixes the batch size so output is deterministic regardless of thread count, which you want if you ever need to reproduce a result exactly. -Y uses soft clipping for supplementary alignments, which structural variant callers prefer. Writing CRAM instead of BAM cuts a 30x alignment from roughly 35 GB to 15 to 18 GB with lossless reference-based compression. Skip base quality score recalibration: with DeepVariant it provides no measurable benefit and costs hours.

Now measure coverage, since this is the number that bounds everything you can conclude:

mosdepth --fast-mode --by 1000 --threads 4 \
  --fasta GRCh38_full_analysis_set_plus_decoy_hla.fa SUBJECT01 SUBJECT01.cram

For a 30x genome we expect mean autosomal coverage of 28 to 32x with at least 95% of the callable genome above 10x. Regions below 10x are regions where an absent variant call means nothing.

4. Call small variants

Use DeepVariant. On Illumina 30x data it reaches SNV F1 above 99.5% against the Genome in a Bottle truth sets and beats GATK HaplotypeCaller clearly on indels, without needing a recalibration or joint-genotyping step for a single sample.

docker run -v "$PWD":/data google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \
  --ref=/data/GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --reads=/data/SUBJECT01.cram \
  --output_vcf=/data/SUBJECT01.dv.vcf.gz \
  --output_gvcf=/data/SUBJECT01.dv.g.vcf.gz \
  --num_shards=32

Always keep the gVCF. It carries reference-confidence blocks, so a year from now you can distinguish “reference at this position with 30x support” from “no data here”, which is exactly the distinction that matters when you revisit a locus that has since been reclassified in ClinVar.

Normalize and split multiallelics before annotation, since annotation tools and allele-frequency databases assume left-aligned, biallelic representations:

bcftools norm -m -any -f GRCh38_full_analysis_set_plus_decoy_hla.fa \
  -Ou SUBJECT01.dv.vcf.gz \
| bcftools norm -d exact -Oz -o SUBJECT01.norm.vcf.gz
bcftools index -t SUBJECT01.norm.vcf.gz

5. Call structural variants, copy number, and repeat expansions

Short-read structural variant calling is much less accurate than SNV calling, so run two callers with different evidence models and treat the intersection as your confident set and the union as your hypothesis list.

manta/bin/configManta.py --bam SUBJECT01.cram \
  --referenceFasta GRCh38_full_analysis_set_plus_decoy_hla.fa --runDir manta_run
manta_run/runWorkflow.py -j 16

delly call -g GRCh38_full_analysis_set_plus_decoy_hla.fa \
  -x human.hg38.excl.tsv -o SUBJECT01.delly.bcf SUBJECT01.cram

Manta uses split reads and local assembly, Delly leans on discordant pairs, and they fail differently. For copy number on a single sample without a matched cohort, cnvkit.py batch --method wgs with a flat reference works for events above roughly 50 kb and should not be trusted below that. For tandem repeat loci, run ExpansionHunter against its catalog, keeping in mind that 150 bp reads cannot span long expansions and the estimates saturate:

ExpansionHunter --reads SUBJECT01.cram \
  --reference GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --variant-catalog variant_catalog/grch38/variant_catalog.json \
  --output-prefix SUBJECT01.eh

Any repeat expansion or structural variant call that looks clinically meaningful needs orthogonal confirmation in a certified laboratory and interpretation by a clinical geneticist. Short-read calls in these classes are leads. They are not conclusions.

6. Annotate, then filter in that order

Annotate with Ensembl VEP running offline against a local cache, which keeps the run reproducible and fast.

vep --offline --cache --dir_cache /refs/vep --assembly GRCh38 \
  --fasta GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --input_file SUBJECT01.norm.vcf.gz --vcf --compress_output bgzip \
  --output_file SUBJECT01.vep.vcf.gz \
  --everything --pick_allele_gene --fork 8 \
  --custom /refs/clinvar_20250601.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT \
  --custom /refs/gnomad.v4.1.sites.vcf.gz,gnomADg,vcf,exact,0,AF,AF_popmax \
  --plugin CADD,/refs/CADD_whole_genome_SNVs.tsv.gz \
  --plugin SpliceAI,snv=/refs/spliceai_scores.snv.hg38.vcf.gz,indel=/refs/spliceai_scores.indel.hg38.vcf.gz

Then filter with bcftools +split-vep on the annotated file rather than filtering first. A typical starting query for variants worth a second look is ClinVar pathogenic or likely pathogenic with at least two-star review status and gnomAD popmax allele frequency below 0.001, plus loss-of-function consequences in the same frequency band. Expect on the order of 50 to 100 hits in a healthy genome. Nearly all of them are recessive carrier states, low-penetrance risk alleles, or annotation artifacts. The interpretive step is deciding which of those has any bearing on your health. That is clinical work and belongs with a genetic counselor or medical geneticist rather than with a filtering expression.

Two annotations worth adding separately: pharmacogenomic star alleles via PharmCAT (CYP2D6 needs a dedicated caller because of its structural variation), and HLA typing, which you can derive from either DNA or RNA reads. Pipelines built for cancer immunology have long combined HLA genotyping with expression and neoantigen prediction in one workflow, and the same components work on a healthy individual’s data.⁴

7. Quantify the transcriptome

RNA sequencing tells you what is expressed, which is the layer that changes week to week. For expression levels, use Salmon’s selective alignment against a decoy-aware index built from the transcriptome plus the genome:

salmon index -t gentrome.fa.gz -d decoys.txt -i salmon_idx_v45 -k 31 -p 16

salmon quant -i salmon_idx_v45 -l A \
  -1 rna_R1.fastq.gz -2 rna_R2.fastq.gz \
  --validateMappings --gcBias --seqBias --posBias \
  --numBootstraps 30 -p 16 -o SUBJECT01_rna

Keep the bootstraps. With one biological sample per timepoint, inferential replicates are the only handle you have on technical variance in transcript-level estimates. Import with tximport in R and summarize to gene level. Compare across your own timepoints rather than against population reference values, since batch and library-prep effects between labs swamp most real biological differences. If you also want splicing, fusions, or allele-specific expression, run STAR separately with --quantMode GeneCounts --outSAMtype BAM SortedByCoordinate and feed the alignment to Arriba. Small RNAs need their own library and their own pipeline: the adapters, the read lengths, and the annotation are all different, and the microRNA transcriptome of a purified cell population is only recoverable with a size-selected protocol.⁵ A useful sanity check on your RNA data is the fraction of reads assigned to globin and ribosomal RNA. Above roughly 20% globin in whole blood means depletion did not work and your effective depth is much lower than the read count suggests.

Reads that fail to map to human contain your blood and, in a stool or saliva library, your gut microbial content. Those reads are not garbage. Dedicated microbiome pipelines exist to classify and visualize them rather than making you assemble the steps yourself.⁶

8. Validate the pipeline before you trust it

This is the step almost every personal genomics writeup skips, and it is the one that separates a result from a guess. Run your exact pipeline on a Genome in a Bottle sample with a published truth set, using the same container versions and the same flags. Then measure concordance with hap.py:

docker run -v "$PWD":/data pkrusche/hap.py /opt/hap.py/bin/hap.py \
  /data/HG002_GRCh38_1_22_v4.2.1_benchmark.vcf.gz \
  /data/HG002.mypipeline.vcf.gz \
  -f /data/HG002_GRCh38_1_22_v4.2.1_benchmark_noinconsistent.bed \
  -r /data/GRCh38_full_analysis_set_plus_decoy_hla.fa \
  -o /data/happy_out --engine=vcfeval --threads 16

Read the stratified output and do not stop at the headline F1. Performance in high-confidence unique regions tells you little about segmental duplications, homopolymers longer than 10 bp, or the MHC. Published clinical validation frameworks make exactly this point: a pipeline’s accuracy is a property of the pipeline, the reference, the region class, and the variant type together, and each combination needs its own measured sensitivity and positive predictive value.⁷ Where no truth set covers your situation, simulation lets you inject known variants at known coverage and read length and measure recovery, which is how tools built to compare resequencing pipelines quantify the effect of each parameter choice.⁸ Rerun this validation whenever you bump a tool version. A caller upgrade that improves genome-wide F1 by 0.1% can change your calls in a specific gene you care about.

9. Wrap it in a workflow manager with real provenance

Write the whole thing as a Nextflow or Snakemake workflow with pinned container digests, and for a first pass use the community pipelines rather than hand-rolling: nf-core/sarek for germline WGS and nf-core/rnaseq for expression.

nextflow run nf-core/sarek -r 3.4.4 -profile docker \
  --input samplesheet.csv --genome GATK.GRCh38 \
  --tools deepvariant,manta,tiddit,vep --outdir results/

The reason to invest here is longitudinal. A one-off analysis can be a shell script. A profile you extend with new timepoints for years needs to record which input files, tool versions, and parameters produced every output file. Institutional NGS platforms built exactly this kind of automated tracking layer because manual bookkeeping fails at scale, and the failure mode is not lost data but data you can no longer interpret.⁹ At minimum, store each run’s nextflow.log, the resolved config, container digests, and the reference MD5 alongside the outputs.

Common problems

Mapping quality zero across whole gene families usually means you aligned to a reference that includes ALT contigs without ALT-aware post-processing. Switch to the no-alt analysis set and re-align.

Indel calls clustering at homopolymers, especially in a single direction, point to PCR stutter from a PCR-plus library. Check the duplication rate in markdup.stats.txt and the strand bias annotations, and treat indels in repeats longer than 10 bp as unresolved without orthogonal confirmation.

A sudden drop in called variants after a tool upgrade is often a normalization mismatch: the new caller emits a different representation of the same indel and your annotation join stops matching. Always run bcftools norm -m -any before any join on position and allele.

Sample identity errors are more common than anyone admits. Run bcftools gtcheck between every new timepoint and your baseline. Discordance above a few percent at common SNPs means a swap somewhere in the chain. The fix is upstream in the lab. It is not in the analysis.

Zero-coverage gaps in genes you care about are not always sequencing failures. Some are genuine reference problems where GRCh38 misrepresents the locus. Check the region against CHM13 before concluding your library was bad.

Finally, RNA expression values that shift dramatically between two timepoints most often reflect collection conditions rather than biology. Time from draw to stabilization, differential white cell counts, and prep batch all move gene-level estimates by amounts that dwarf most regulatory changes. Record collection metadata with the same discipline you record pipeline parameters, and interpret any finding you intend to act on with a clinician.

Footnotes

1. Alfi Sophian, Andhika Harumanto, Abinawanto. Rapid genomic detection of food fraud: A critical review of PCR, CRISPR, next-generation sequencing, and AI-enabled methods for food control applications. Food Control, 2027. https://doi.org/10.1016/j.foodcont.2026.112522 ↩

2. Maksim I. Nadtoka, Anna Y. Bukharina, German V. Roev, et al. Development and Initial Analytical Evaluation of a Multiplex PCR-Based Targeted NGS Assay for Expanded Viral Detection in Clinical Respiratory Specimens. Diagnostics, 2026. https://doi.org/10.3390/diagnostics16182888 ↩

3. Jacqueline Weber-Lehmann, Elmar Schilling, Georg Gradl, et al. Finding the needle in the haystack: Differentiating “identical” twins in paternity testing and forensics by ultra-deep next generation sequencing. Forensic Science International: Genetics, 2014. https://doi.org/10.1016/j.fsigen.2013.10.015 ↩

4. Elias Tappeiner, Francesca Finotello, Pornpimol Charoentong, et al. TIminer: NGS data mining pipeline for cancer immunology and immunotherapy. Bioinformatics, 2017. https://doi.org/10.1093/bioinformatics/btx377 ↩

5. Todd A. Fehniger, Todd Wylie, Elizabeth Germino, et al. Next-generation sequencing identifies the natural killer cell microRNA transcriptome. Genome Research, 2010. https://doi.org/10.1101/gr.107995.110 ↩

6. Teresia M. Buza, Triza Tonui, Francesca Stomeo, et al. iMAP: an integrated bioinformatics and visualization pipeline for microbiome data analysis. BMC Bioinformatics, 2019. https://doi.org/10.1186/s12859-019-2965-4 ↩

7. Somak Roy. Principles and Validation of Bioinformatics Pipeline for Cancer Next-Generation Sequencing. Clinics in Laboratory Medicine, 2022. https://doi.org/10.1016/j.cll.2022.05.006 ↩

8. B. Nevado, M. Perez‐Enciso. Pipeliner: software to evaluate the performance of bioinformatics pipelines for next‐generation resequencing. Molecular Ecology Resources, 2014. https://doi.org/10.1111/1755-0998.12286 ↩

9. Eva Reisinger, Lena Genthner, Jules Kerssemakers, et al. OTP: An automatized system for managing and processing NGS data. Journal of Biotechnology, 2017. https://doi.org/10.1016/j.jbiotec.2017.08.006 ↩