How to Run an NGS Pipeline on Your Own Whole-Genome and RNA Data

By the end of this guide you will have five artifacts from one person’s sequencing run: a coordinate-sorted, duplicate-marked CRAM aligned to GRCh38, a normalized and filtered small-variant VCF with per-sample quality metrics, an accuracy estimate for that VCF measured against a public truth set, an annotated shortlist of rare and clinically catalogued variants, and a transcript-level expression table from RNA-seq. You need the raw data (paired FASTQ files, roughly 50–60 GB gzipped for a 30× human genome, plus 5–10 GB per RNA-seq library), a machine with at least 16 cores and 64 GB of RAM, about 1.5 TB of free disk, and Docker or Singularity. Everything below runs on Linux; on a Mac, run the containers under an x86 emulation layer or rent a cloud instance, because several of these tools have no working ARM builds. Budget six to twelve hours of wall time for a 30× genome on 32 cores, most of it in alignment and variant calling.

1. Set up the reference and the toolchain

Every downstream number depends on which reference you align to, so fix that first and record it. We use the GRCh38 “no-alt analysis set” from NCBI, which includes the decoy contigs and excludes the alternate haplotype scaffolds that confuse mapping-quality calculations in most callers. T2T-CHM13 is a better assembly in absolute terms, but nearly all annotation resources (gnomAD, ClinVar, the GIAB benchmark regions) are keyed to GRCh38, so we treat CHM13 as a second alignment you do later for specific questions, not as your primary coordinate system.

mamba create -n ngs -c conda-forge -c bioconda \
  bwa-mem2=2.2.1 samtools=1.21 bcftools=1.21 fastp=0.23.4 \
  mosdepth=0.3.8 multiqc=1.25 salmon=1.10.3 star=2.7.11b ensembl-vep=112
conda activate ngs

REF=GCA_000001405.15_GRCh38_no_alt_analysis_set.fna
curl -O "https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/\
seqs_for_alignment_pipelines.ucsc_ids/${REF}.gz"
gunzip ${REF}.gz
samtools faidx $REF
bwa-mem2 index $REF      # ~1 hour, peaks near 80 GB RAM, writes ~26 GB of index

The index build is the single largest memory spike in this pipeline. If your machine cannot hold it, download a prebuilt index rather than swapping to disk. Note the tool versions and the reference checksum somewhere permanent, because a variant call is only interpretable alongside the reference and caller version that produced it. If you are coming to this without a genomics background, the free curricula assembled for self-directed learners cover the underlying concepts well enough that you do not need a formal course ¹, and published classroom templates for genomic tool training are a reasonable syllabus to follow ².

2. Look at the raw reads before you align them

Run fastp on each library, primarily to get a report rather than to modify reads. Modern aligners soft-clip adapter sequence at read ends, and aggressive trimming tends to cost you more in alignment sensitivity than it gains.

fastp -i r1.fq.gz -I r2.fq.gz -o t1.fq.gz -O t2.fq.gz \
  --detect_adapter_for_pe --disable_quality_filtering \
  --length_required 50 -w 8 -j sample.fastp.json -h sample.fastp.html

Read the report for four things: duplication estimate, adapter content, per-base quality decay in the last 20 cycles, and GC distribution. A single sharp GC peak away from the human mode of about 41 percent usually means contamination or a heavy rRNA fraction in RNA libraries. If adapter content is under a few percent and quality holds above Q30 through most cycles, use the untrimmed FASTQs for alignment and keep the fastp output only as documentation.

3. Align to GRCh38 and write CRAM

Stream alignment, mate fixing, sorting, and duplicate marking through one pipe. Writing intermediate BAMs at each stage triples your I/O for no benefit.

bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:HXXXX.1\tSM:SAMPLE01\tLB:lib1\tPL:ILLUMINA\tPU:HXXXX.1.ATCACG' \
  $REF r1.fq.gz r2.fq.gz \
| samtools fixmate -m -u -@ 4 - - \
| samtools sort -u -@ 8 -m 3G -T /scratch/srt - \
| samtools markdup -@ 8 --reference $REF -O CRAM,level=8 - SAMPLE01.cram
samtools index -@ 8 SAMPLE01.cram

Three flags matter here. -K 100000000 fixes the batch size so that output is identical regardless of thread count, which you want when you are comparing two runs. -Y uses soft clipping for supplementary alignments, which structural-variant callers need in order to recover breakpoint sequence. The @RG line is not optional: DeepVariant and GATK both refuse to run without a sample name, and the platform unit encodes flowcell and lane so that you can trace a batch effect later. CRAM at level 8 with the reference available costs roughly 40 percent less disk than BAM, landing near 16–20 GB for a 30× genome, and samtools reads it transparently everywhere a BAM is accepted.

4. Measure coverage, duplication, and contamination

Before calling a single variant, decide whether this library is worth calling on. Three tools give you the answer in under twenty minutes.

samtools stats -@ 8 --reference $REF SAMPLE01.cram > SAMPLE01.stats
mosdepth -t 4 -n --fast-mode --by 1000 SAMPLE01 SAMPLE01.cram
multiqc .

The thresholds we hold to for a 30× PCR-free human genome: mean autosomal depth at or above 30×, at least 95 percent of the genome covered at 20× or more, duplicate rate under 10 percent (PCR-free libraries usually come in nearer 2 percent), insert-size mode between 300 and 450 bp, and an error rate from samtools stats below 0.5 percent. A library that hits 28× mean depth but only 88 percent of the genome at 20× has a coverage uniformity problem, typically from PCR bias, and its indel calls in GC-extreme regions will be unreliable.

Contamination is the failure mode that quietly corrupts everything downstream, because a 3 percent foreign DNA fraction shows up as a diffuse excess of low-allele-fraction heterozygous calls that look like real variants.

VerifyBamID --SVDPrefix resource/1000g.phase3.100k.b38.vcf.gz.dat \
  --Reference $REF --BamFile SAMPLE01.cram

Treat a FREEMIX estimate above 0.02 as a reason to stop and resequence rather than a number to correct for.

5. Call small variants

For a single germline human genome from Illumina short reads, we use DeepVariant and do not reach for GATK’s HaplotypeCaller unless we need joint calling across a cohort. DeepVariant needs no base-quality recalibration step, has fewer tunable knobs to get wrong, and performs better on indels in homopolymers out of the box.

docker run -v "$PWD":/data google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \
  --ref=/data/$REF \
  --reads=/data/SAMPLE01.cram \
  --output_vcf=/data/SAMPLE01.vcf.gz \
  --output_gvcf=/data/SAMPLE01.g.vcf.gz \
  --num_shards=32 \
  --intermediate_results_dir=/data/dv_tmp

Keep the gVCF. It records confidence at reference positions, which is the only way to distinguish “this position matches the reference” from “this position had no usable coverage,” and it is what you need if you later add family members or a second timepoint.

Normalize before you do anything else with the VCF. Left-alignment and allele splitting make variants comparable between call sets and between annotation runs.

bcftools norm -f $REF -m -any -Oz -o SAMPLE01.norm.vcf.gz SAMPLE01.vcf.gz
bcftools index -t SAMPLE01.norm.vcf.gz
bcftools stats -F $REF SAMPLE01.norm.vcf.gz > SAMPLE01.bcfstats

Sanity-check the aggregate numbers against population expectations: roughly 4.5 to 5 million variants against GRCh38, a genome-wide transition/transversion ratio near 2.0 to 2.1, and a heterozygous-to-homozygous ratio between about 1.5 and 2.0. A Ti/Tv below 1.8 genome-wide means false positives are leaking in. Ancestry shifts total variant count by several hundred thousand, so compare against a reference sample of similar background rather than a single global number.

6. Benchmark your call set against a truth set

This is the step most tutorials skip, and it is what turns a pipeline from something that runs into something you can trust. Genome in a Bottle publishes high-confidence variant calls and confident-region BED files for several samples. If you sequenced a GIAB sample as a control, compare directly. If you did not, run the same pipeline on public HG002 FASTQs once and treat the result as your pipeline’s accuracy profile.

docker run -v "$PWD":/data jmcdani20/hap.py:v0.3.12 /opt/hap.py/bin/hap.py \
  /data/HG002_GRCh38_1_22_v4.2.1_benchmark.vcf.gz \
  /data/HG002.norm.vcf.gz \
  -f /data/HG002_GRCh38_1_22_v4.2.1_benchmark_noinconsistent.bed \
  -r /data/$REF -o /data/happy --engine=vcfeval --threads 32

On 30× Illumina data with this pipeline, expect SNV F1 above 99.5 percent inside the benchmark regions, with indel F1 a little lower and concentrated losses in homopolymers longer than about eight bases. Those numbers apply only inside the high-confidence regions, which exclude segmental duplications and much of the repetitive genome. Roughly 5 to 10 percent of the genome remains inaccessible to short reads, and any gene of interest sitting in a segmental duplication (for example the paralog-rich regions where pseudogenes mirror a real gene) needs orthogonal confirmation. Clinical laboratories handle this with targeted follow-up, and the analytical validation frameworks they use are described in the clinical-genetics bioinformatics literature ³.

7. Annotate and build a shortlist

Annotation turns five million coordinates into a few hundred you can read. We use Ensembl VEP offline with a local cache, plus custom annotation from gnomAD and ClinVar.

vep -i SAMPLE01.norm.vcf.gz -o SAMPLE01.vep.vcf.gz \
  --vcf --compress_output bgzip --fork 16 \
  --cache --offline --dir_cache "$HOME/.vep" --assembly GRCh38 --fasta $REF \
  --everything --pick_allele_gene \
  --custom file=gnomad.genomes.v4.1.sites.vcf.bgz,short_name=gnomADg,\
format=vcf,type=exact,fields=AF_grpmax \
  --custom file=clinvar_20250601.vcf.gz,short_name=ClinVar,\
format=vcf,type=exact,fields=CLNSIG%CLNREVSTAT%CLNDN

Then filter down in two passes. First, keep variants with a maximum population frequency under 0.001 that are predicted loss-of-function or missense with high deleteriousness scores. Second, separately, keep everything ClinVar labels pathogenic or likely pathogenic with at least two-star review status, regardless of predicted consequence.

bcftools view -i 'INFO/ClinVar_CLNSIG~"athogenic"' SAMPLE01.vep.vcf.gz \
  | bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%INFO/ClinVar_CLNDN\t[%GT]\n'

Most of what survives this filter in any healthy genome is recessive carrier status, meaning a single copy of a variant that causes disease only when both copies are affected. A low single-digit percentage of people carry a variant in the small set of genes that professional bodies consider medically actionable. Research-grade sequencing is not a diagnostic test: a pathogenic call here needs confirmation in a CLIA-certified laboratory and interpretation with a clinician or genetic counselor before it means anything about your health, and that confirmation is still routinely done by Sanger sequencing, which remains the reference method for reading a single locus at very high per-base accuracy ⁴. The same tiering-then-confirmation structure is what clinical oncology pipelines use for somatic variant review ⁵.

8. Quantify the RNA

RNA-seq answers a different question from DNA: which transcripts were present, in what quantity, in the tissue you sampled at the moment you sampled it. For expression levels we use Salmon with a decoy-aware index, which is fast and handles multi-mapping reads by probabilistic assignment rather than discarding them.

# gentrome.fa.gz = transcripts + full genome; decoys.txt = genome contig names
salmon index -t gentrome.fa.gz -d decoys.txt -i salmon_idx --gencode -k 31 -p 16

salmon quant -i salmon_idx -l A -1 r1.fq.gz -2 r2.fq.gz \
  --validateMappings --gcBias --seqBias --numBootstraps 30 \
  -p 16 -o quant_SAMPLE01

Check the mapping rate in quant_SAMPLE01/logs/salmon_quant.log. Under about 75 percent for a polyA-selected human library usually means rRNA carryover, globin dominance in whole blood, or degraded input. The bootstrap replicates give you a technical variance estimate per transcript, which matters enormously when your design is one person over time rather than two groups of twenty. For a longitudinal single-subject profile, work in log2 ratios against that person’s own earlier timepoints and ignore any transcript whose bootstrap spread is comparable to the change you are looking at.

Salmon gives you abundance but no alignments, so add a STAR run when you need splice junctions, fusion detection, or isoform structure. Purpose-built pipelines exist for the harder RNA questions, including isoform diversity and RNA-editing analysis from standard RNA-seq input ⁶ and circular RNA detection, which requires backsplice-junction-aware tooling that ordinary quantifiers do not provide ⁷.

9. Wrap it in a workflow manager

Once the pipeline above works end to end, stop running it by hand. We use Nextflow with the nf-core pipelines: nf-core/sarek covers alignment through variant calling and annotation, and nf-core/rnaseq covers the RNA side.

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

The value is caching and provenance. -resume restarts from the last completed process rather than the beginning, containers pin every tool version, and the run report records exactly which parameters produced which file. Writing down the procedure so that a second person, or you in a year, can reproduce a result is the part of analysis that decays fastest without discipline, and the standard-operating-procedure conventions developed for shared workflow platforms are a good model for how much to document ⁸.

Common problems

Reference mismatch is the most frequent silent failure. A VCF called against GRCh37 coordinates and annotated against a GRCh38 cache produces plausible-looking gene names attached to the wrong positions. Check the contig lines in your VCF header against your FASTA index before every annotation run, and prefer chr1-style UCSC naming throughout rather than converting midway.

Sample swaps are the second. Any time you have both DNA and RNA from the same person, genotype the RNA at common SNP sites with bcftools mpileup restricted to a 1000 Genomes site list and check concordance against the WGS calls. Genuine same-person concordance runs above 95 percent at well-covered sites, and anything near 50 percent means two different people.

Memory exhaustion during samtools sort shows up as a pipeline that dies hours in. The -m flag is per thread, so -@ 8 -m 3G reserves 24 GB before the aligner’s own footprint. Give the sort a -T path on a fast disk with at least 200 GB free, since it spills thousands of temporary files.

Multi-allelic sites and unnormalized indels break joins between call sets. If two VCFs disagree about a variant you can see in IGV, run both through bcftools norm -m -any -f $REF and compare again before concluding anything.

Finally, expect what short reads cannot do. Structural variants called by Manta or smoove from 150 bp reads recover only part of what long reads find, with recall worst for insertions and for anything inside tandem repeats. Repeat expansions, methylation, and phasing across long distances all need a different assay. That limitation is intrinsic to the read length, not to the software, and no amount of pipeline tuning fixes it ^{9 10}.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

1. David B. Searls. An Online Bioinformatics Curriculum. PLoS Computational Biology, 2012. https://doi.org/10.1371/journal.pcbi.1002632 ↩

2. Scott Hotaling, Brittany L. Slabach, David W. Weisrock. Next-generation teaching: a template for bringing genomic and bioinformatic tools into the classroom. Journal of Biological Education, 2017. https://doi.org/10.1080/00219266.2017.1357650 ↩

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. Jerzy K. Kulski. Next-Generation Sequencing — An Overview of the History, Tools, and “Omic” Applications. Next Generation Sequencing - Advances, Applications and Challenges, 2016. https://doi.org/10.5772/61964 ↩

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

6. Noel-Marie Plonski, Emily Johnson, Madeline Frederick, et al. Automated Isoform Diversity Detector (AIDD): a pipeline for investigating transcriptome diversity of RNA-seq data. BMC Bioinformatics, 2020. https://doi.org/10.1186/s12859-020-03888-6 ↩

7. Edward A. Salinas, Yvonne J. K. Edwards. Circular RNA Identification and Characterization with CircRNAFlow: A Bioinformatics Approach. Advances in Experimental Medicine and Biology, 2025. https://doi.org/10.1007/978-981-96-9428-0_6 ↩

8. Sona Charles, Kiran K. Telukunta, Darshana Joshi, et al. Standard Operating Procedures for Effective Galaxy Workflows. Next-Generation Sequencing, 2025. https://doi.org/10.1201/9781003354062-15 ↩

9. Chelliah Ramachandran, Eric Banan-Mwine Daliri, Elahi Fazle, et al. Impact of Sequencing and Bioinformatics Tools in Food Microbiology. Sequencing Technologies in Microbial Food Safety and Quality, 2021. https://doi.org/10.1201/9780429329869-8 ↩

10. Sheng-Yong Niu, Jinyu Yang, Adam McDermaid, et al. Bioinformatics tools for quantitative and functional metagenome and metatranscriptome data analysis in microbes. Briefings in Bioinformatics, 2017. https://doi.org/10.1093/bib/bbx051 ↩