How to Build an NGS Analysis Pipeline for Your Own Genome

By the end of this guide you will have a reproducible pipeline that takes paired-end FASTQ files and produces a set of outputs. These are a coordinate-sorted CRAM aligned to GRCh38, a small-variant VCF with per-sample genotype quality, and a structural variant and copy-number call set. You also get an annotated variant table joined to gnomAD frequencies and ClinVar assertions, transcript-level RNA quantifications, and a set of QC metrics that tell you whether any of it can be trusted. You need a Linux machine with at least 16 cores and 64 GB of RAM (DeepVariant and STAR are the memory-hungry steps). You also need roughly 1 TB of free disk for a single 30x genome plus RNA sample. That is FASTQ at about 100 GB gzipped, an intermediate BAM around 90 GB, and a final CRAM closer to 18 GB. The reference bundles and annotation caches take another 80 GB or so. Storage and compute, rather than algorithms, are usually the binding constraint when people first move sequencing analysis in-house ¹. Budget a day of wall-clock time per genome on 16 cores, less if you have 32 or more.

1. Set up a reproducible environment and verify your raw data

Start by pinning every tool version, because a pipeline whose outputs you cannot regenerate two years from now is not a pipeline. We use a conda environment for the command-line tools and Docker for DeepVariant, which ships models tied to specific container tags.

mamba create -n ngs -c conda-forge -c bioconda \
  fastqc=0.12.1 fastp=0.23.4 multiqc=1.21 \
  bwa-mem2=2.2.1 samtools=1.19 bcftools=1.19 gatk4=4.5.0.0 \
  mosdepth=0.3.6 somalier=0.2.19 verifybamid2=2.0.1 \
  manta=1.6.0 smoove=0.2.8 expansionhunter=5.0.0 \
  salmon=1.10.2 star=2.7.11b ensembl-vep=111.0 nextflow=24.04.2
mamba activate ngs

Before anything else, check the integrity of what you were given. Sequencing providers deliver FASTQ (or unaligned BAM) with checksums, and a truncated transfer is the single most common cause of a run that fails four hours in.

md5sum -c checksums.md5
seqkit stats -a -j 8 *.fastq.gz      # read counts, length distribution, N50

For a 30x human genome on 2x150 chemistry you should see roughly 600 million read pairs. If the count is far below that, the depth you paid for is not in the files.

2. Trim, QC, and decide what to throw away

Read-level QC is where you catch adapter read-through, quality collapse at the end of a run, and library complexity problems. We run fastp because it handles adapter detection, polyG trimming and QC reporting in a single pass over the data. PolyG trimming is needed for two-colour chemistry, where a dark base is called as G.

fastp -i R1.fq.gz -I R2.fq.gz -o R1.trim.fq.gz -O R2.trim.fq.gz \
  --detect_adapter_for_pe --trim_poly_g --cut_tail --cut_tail_mean_quality 20 \
  --length_required 50 --thread 8 \
  --json fastp.json --html fastp.html

Be conservative with quality trimming for DNA. Modern aligners soft-clip poor tails on their own, and aggressive trimming shortens reads into repetitive regions and biases indel calling. Our rule is to remove adapters and polyG, drop reads below 50 bp, and leave the rest alone. Then aggregate everything with multiqc . and look at three numbers: duplication rate (above 20% on a PCR-free WGS library means the library was over-amplified or under-loaded), GC distribution (a second peak suggests contamination), and per-base quality at cycle 150.

3. Align to GRCh38 with the right reference build

The reference you choose propagates into every downstream coordinate and every annotation lookup, so make it an explicit decision. We use the GRCh38 full analysis set with decoys and HLA contigs (GRCh38_full_analysis_set_plus_decoy_hla.fa), because the decoy sequences absorb reads from sequences absent in the primary assembly and reduce false positives in segmental duplications. If you plan to compare against a dataset built on hs37d5 or on T2T-CHM13, you will need liftover, and liftover of indels near assembly discrepancies is lossy.

bwa-mem2 index GRCh38_full_analysis_set_plus_decoy_hla.fa

bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:run1\tSM:sample1\tLB:lib1\tPL:ILLUMINA' \
  GRCh38_full_analysis_set_plus_decoy_hla.fa R1.trim.fq.gz R2.trim.fq.gz \
| samtools sort -@ 8 -m 4G -o sample1.bam -
samtools index -@ 8 sample1.bam

The -K 100000000 flag fixes the number of bases processed per batch, which makes output deterministic regardless of thread count. Without it, two runs of the same data on different core counts produce byte-different BAMs and slightly different mapping quality in ambiguous regions. -Y uses soft clipping for supplementary alignments, which structural variant callers need.

Mark duplicates next, and set the optical duplicate distance to match your flow cell. Patterned flow cells (NovaSeq, NextSeq) produce ExAmp duplicates that sit far apart in pixel space.

gatk MarkDuplicates -I sample1.bam -O sample1.md.bam -M sample1.dup.txt \
  --OPTICAL_DUPLICATE_PIXEL_DISTANCE 2500 --CREATE_INDEX true

samtools view -@ 8 -C -T GRCh38_full_analysis_set_plus_decoy_hla.fa \
  -o sample1.cram sample1.md.bam && samtools index sample1.cram

We skip GATK base quality score recalibration when calling with DeepVariant, which learns its own error model from pileup images and gains nothing from BQSR. If you are calling with HaplotypeCaller, run BQSR.

4. Measure coverage and confirm the sample is what you think it is

Do this before variant calling, because the failure modes here invalidate everything downstream. Three checks matter.

mosdepth --by 1000 --fast-mode -t 8 sample1 sample1.cram
verifybamid2 --SVDPrefix resource/1000g.phase3.100k.b38.vcf.gz.dat \
  --Reference GRCh38_full_analysis_set_plus_decoy_hla.fa --BamFile sample1.cram
somalier extract -d extracted/ --sites sites.hg38.vcf.gz \
  -f GRCh38_full_analysis_set_plus_decoy_hla.fa sample1.cram
somalier relate extracted/*.somalier

Mean autosomal coverage should be within 10% of what you ordered, and the fraction of the callable genome at 20x or better should exceed 90% for a PCR-free 30x library. VerifyBamID2 reports a FREEMIX contamination estimate. Above 0.02 (2% foreign DNA) you will see an excess of low-allele-fraction heterozygous calls, and above 0.05 you should re-sequence rather than analyze. Somalier gives you genotype fingerprints, which is how you confirm later that your RNA sample and your DNA sample came from the same person. Sample identity and contamination checks are the parts of clinical pipeline validation that individual users most often skip and most often need ².

5. Call small variants

For germline SNVs and short indels on Illumina WGS we use DeepVariant. It consistently outperforms the traditional haplotype-based callers on indels in homopolymers, and it requires no tuning beyond picking the right model type.

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/sample1.cram \
  --output_vcf=/data/sample1.vcf.gz \
  --output_gvcf=/data/sample1.g.vcf.gz \
  --num_shards=32 --intermediate_results_dir=/data/tmp

Emit the gVCF even if you only want the VCF. The gVCF records reference confidence at every position, which is the only way to distinguish “homozygous reference” from “no data here” when you revisit a locus later. A 30x genome typically yields 4.5 to 5 million variants, of which about 4 million are SNVs, and the transition/transversion ratio should land near 2.0 genome-wide. A Ti/Tv far below 2.0 signals a false positive problem.

bcftools stats -F GRCh38_full_analysis_set_plus_decoy_hla.fa sample1.vcf.gz > sample1.stats

6. Call structural and repeat variation separately

Small-variant callers miss everything larger than about 50 bp, which includes many of the variants with the largest functional effect. Run at least two orthogonal SV callers and intersect, because single-caller precision on deletions and duplications is poor.

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

smoove call --name sample1 --fasta GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --outdir smoove/ --exclude exclude.cnvnator_100bp.GRCh38.20170403.bed -p 8 sample1.cram

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

Manta uses split reads and discordant pairs and is fast. The smoove tool wraps lumpy plus duphold, which adds depth-based confirmation to each call. ExpansionHunter genotypes known tandem repeat loci, which short-read SV callers otherwise ignore entirely. Expect 6,000 to 10,000 SV calls per genome before filtering, and expect most singleton calls to be artifacts near assembly gaps.

7. Validate the pipeline against a truth set before you trust it on yourself

This is the step that separates a pipeline from a pile of commands. Download the Genome in a Bottle HG002 FASTQ files, run them through the exact pipeline you just built, and compare against the GIAB v4.2.1 benchmark VCF and high-confidence BED.

hap.py HG002_GRCh38_1_22_v4.2.1_benchmark.vcf.gz hg002.mypipeline.vcf.gz \
  -f HG002_GRCh38_1_22_v4.2.1_benchmark_noinconsistent.bed \
  -r GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --engine=vcfeval --threads 16 -o happy_hg002 \
  --stratification stratification-regions/GRCh38/GRCh38-all-stratifications.tsv

On a 30x PCR-free library with bwa-mem2 and DeepVariant you should see SNV F1 above 0.995 and indel F1 above 0.99 inside the high-confidence regions. If your numbers are materially worse, something in your configuration is wrong, and you now know it before interpreting your own data. The stratification files break performance down by region class, and the numbers in low-mappability regions and in homopolymers longer than 10 bp are dramatically worse than the genome-wide average. Benchmarking against reference material with known genotypes is the accepted method for establishing analytical sensitivity and specificity of a sequencing pipeline, and it is what clinical laboratories are required to do before reporting results ^{2 3}. Keep the hap.py output alongside your pipeline code as its performance specification.

8. Annotate, then filter with population frequency

Annotation is where a VCF becomes something you can reason about. We use Ensembl VEP offline with a local cache, plus custom overlays for ClinVar and gnomAD, so that the whole process is reproducible and does not depend on a web service that changes underneath you.

vep -i sample1.vcf.gz --vcf --compress_output bgzip -o sample1.vep.vcf.gz \
  --cache --offline --dir_cache $VEP_CACHE --assembly GRCh38 \
  --fasta GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --everything --pick_allele_gene --fork 8 \
  --plugin CADD,whole_genome_SNVs.tsv.gz \
  --custom clinvar_20240609.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT \
  --custom gnomad.genomes.v4.0.sites.chr.vcf.gz,gnomADg,vcf,exact,0,AF,AF_popmax

Then flatten to a table and filter. The single most useful filter is population allele frequency. A variant present at 1% or more in gnomAD is not by itself a candidate for a rare high-penetrance effect.

bcftools +split-vep sample1.vep.vcf.gz -f '%CHROM\t%POS\t%REF\t%ALT\t%SYMBOL\t%Consequence\t%gnomADg_AF\t%ClinVar_CLNSIG\t%CADD_PHRED[\t%GT\t%GQ]\n' \
  -d -A tab | awk -F'\t' '($7=="." || $7<0.01) && $9>20' > candidates.tsv

Two cautions. First, ClinVar assertions vary enormously in quality, and the review status field (CLNREVSTAT) matters as much as the significance field: a single-submitter “likely pathogenic” with no assertion criteria is weak evidence. Second, annotation databases and transcript sets update continuously. A variant’s interpretation can change without the underlying data changing at all, so record the cache and database versions with every analysis ^{3 4}. If anything in this table looks clinically relevant to you, the correct next step is confirmation in an accredited clinical laboratory and a conversation with a physician or genetic counselor. Research-grade calls from a self-built pipeline are not a diagnosis and cannot support any treatment decision.

9. Add the RNA layer

RNA sequencing tells you what is being transcribed, which the genome cannot. For quantification we use Salmon with selective alignment against a decoy-aware index. It is fast, accounts for multi-mapping across isoforms, and corrects sequence and GC bias.

# build gentrome: transcriptome + genome as decoy
cat gencode.v45.transcripts.fa GRCh38.primary_assembly.genome.fa > gentrome.fa
grep "^>" GRCh38.primary_assembly.genome.fa | cut -d " " -f1 | sed 's/>//g' > decoys.txt
salmon index -t gentrome.fa -d decoys.txt -i salmon_idx -k 31 -p 16 --gencode

salmon quant -i salmon_idx -l A -1 rna_R1.fq.gz -2 rna_R2.fq.gz \
  --gcBias --seqBias --numGibbsSamples 20 -p 16 -o rna_quant/

The Gibbs samples give you posterior uncertainty per transcript, which matters because isoform-level estimates for genes with many similar transcripts are far less certain than the point estimate suggests. Check the mapping rate in rna_quant/logs/salmon_quant.log: below 70% usually means rRNA or globin carryover, or the wrong strandedness assumption (-l A infers it, and the inferred library type is printed in lib_format_counts.json).

With a single sample there is no differential expression analysis to run, since there is no contrast. What you can do is compare your TPMs against a reference distribution such as GTEx for the matched tissue, and treat the comparison as approximate because library prep and batch differences move expression estimates substantially. If you have repeated sampling over time, the same person measured at several timepoints is a much stronger comparison than any external cohort. Also run somalier on the RNA alignment and confirm the fingerprint matches your DNA sample.

10. Wrap the whole thing in a workflow manager

Once each step works, stop running them by hand. We use Nextflow with containers. It gives you resumable execution, per-process resource declarations, and a provenance record of exactly which container digest produced which file.

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

nextflow run nf-core/rnaseq -r 3.14.0 -profile docker \
  --input rna_samplesheet.csv --outdir rna_results/ \
  --genome GRCh38 --pseudo_aligner salmon -resume

Starting from nf-core is our recommendation over writing your own DSL2 workflow, because the community pipelines already encode the parameter choices and QC steps discussed above and are continuously tested. Workflow systems exist precisely because manual orchestration of multi-step sequencing analysis does not scale and does not reproduce ^{5 6}. Keep the nextflow.config, the run report, and the pipeline_info/ directory with your results.

Common problems

Chromosome naming mismatches break joins silently. A VCF with chr1 will not intersect a BED with 1, and bedtools will report zero overlaps rather than an error. Normalize everything to the UCSC-style chr prefix at the start and check with bcftools view -h | grep contig.

Unnormalized indels cause false negatives in annotation. The same deletion can be represented several ways, and ClinVar lookups will miss unless both sides are left-aligned and parsimonious. Always run bcftools norm -f ref.fa -m -any -c w before comparison or annotation.

Coverage looks fine on average but fails locally. Mean coverage hides GC-extreme dropout, and first exons of GC-rich genes are routinely under-covered even in a 30x genome. Use the mosdepth per-region output against your genes of interest rather than trusting the genome-wide mean.

Different pipelines give different answers on the same reads. This is expected. Choice of aligner, caller, and filtering thresholds changes call sets measurably. The effect is largest in repetitive and low-complexity regions, so a variant that appears with one tool chain and not another should be treated as uncertain until confirmed ^{7 8}. This is also why step 7 is not optional.

RNA quantification drifts between runs. Library prep batch, RNA integrity, and the fraction of reads consumed by rRNA all shift TPM estimates. When comparing timepoints, process them with identical tool versions and index builds, and include the RNA integrity number from the lab in your metadata.

Everything runs out of disk at 3 a.m. Intermediate BAMs and Nextflow work directories are the culprits. Pipe alignment directly into samtools sort instead of writing an unsorted BAM, convert to CRAM once the BAM is validated, and clean work/ with nextflow clean -f -before <run_name> after you have archived the results.

Footnotes

1. Rama R. Gullapalli, Ketaki V. Desai, Lucas Santana-Santos, et al. Next generation sequencing in clinical medicine: Challenges and lessons for pathology and biomedical informatics. Journal of Pathology Informatics, 2012. https://doi.org/10.4103/2153-3539.103013 ↩

2. Jeffrey A SoRelle, Megan Wachsmann, Brandi L. Cantarel. Assembling and Validating Bioinformatic Pipelines for Next-Generation Sequencing Clinical Assays. Archives of Pathology & Laboratory Medicine, 2020. https://doi.org/10.5858/arpa.2019-0476-ra ↩ ↩²

3. Rute Pereira, Jorge Oliveira, Mário Sousa. Bioinformatics and Computational Tools for Next-Generation Sequencing Analysis in Clinical Genetics. Journal of Clinical Medicine, 2020. https://doi.org/10.3390/jcm9010132 ↩ ↩²

4. Simon Cabello-Aguilar, Julie A. Vendrell, Jérôme Solassol. A Bioinformatics Toolkit for Next-Generation Sequencing in Clinical Oncology. Current Issues in Molecular Biology, 2023. https://doi.org/10.3390/cimb45120608 ↩

5. GunHwan Ko, Pan-Gyu Kim, Jongcheol Yoon, et al. Closha: bioinformatics workflow system for the analysis of massive sequencing data. BMC Bioinformatics, 2018. https://doi.org/10.1186/s12859-018-2019-3 ↩

6. Darlene D. Wagner, Rachel L. Marine, Edward Ramos, et al. VPipe: an Automated Bioinformatics Platform for Assembly and Management of Viral Next-Generation Sequencing Data. Microbiology Spectrum, 2022. https://doi.org/10.1128/spectrum.02564-21 ↩

7. Imane Allali, Jason W. Arnold, Jeffrey Roach, et al. A comparison of sequencing platforms and bioinformatics pipelines for compositional analysis of the gut microbiome. BMC Microbiology, 2017. https://doi.org/10.1186/s12866-017-1101-8 ↩

8. D. S. Horner, G. Pavesi, T. Castrignano, et al. Bioinformatics approaches for genomics and post genomics applications of next-generation sequencing. Briefings in Bioinformatics, 2009. https://doi.org/10.1093/bib/bbp046 ↩