How to Analyze Raw NGS Data From Your Own Genome

By the end of this guide you will have taken raw sequencer output for one person and produced four things. These are a quality-controlled set of reads, a coordinate-sorted alignment in CRAM, and a small-variant call set in VCF with per-variant quality metrics you can reason about. The last is a transcript-level expression table if you also have RNA. You need the raw files (paired FASTQ or an unaligned BAM, plus checksums) and a Linux machine with at least 16 cores and 64 GB of RAM. You also need about 1 TB of scratch disk for a 30x whole genome, and a package manager that will not fight you. We use micromamba with the bioconda channel, or Docker/Singularity images where the vendor publishes them. Budget six to twelve hours of wall time for a 30x genome on 16 cores, most of it in alignment. Everything below is measurement and interpretation. None of it is diagnosis, and a variant that appears to matter medically needs confirmation in a clinical laboratory and a conversation with a genetics professional before it means anything about your health.

1. Understand what “raw” means, and how sequencing differs from PCR

Illumina instruments write BCL (base call) files, a per-cycle, per-tile binary format that is raw in the strictest sense. Almost nobody works with BCL directly. You convert it once with bcl-convert (the successor to bcl2fastq), which demultiplexes by index sequence and emits FASTQ. Practically, “raw NGS data” means gzipped FASTQ: four lines per read, consisting of an identifier, the base sequence, a + separator, and a Phred quality string where the ASCII character minus 33 gives the probability that the base is wrong. A quality of I is Q40, one error in 10,000.

This is where the common question about PCR belongs. Polymerase chain reaction amplifies a chosen region using primers you design, and it answers a question you already asked: is this sequence present, and how much of it. Sequencing reads out the bases themselves across whatever the library covers, so it answers questions you had not thought to ask. The two are not alternatives. PCR is a step inside most sequencing library preparations, used to amplify the adapter-ligated fragments before loading, which is why duplicate reads exist and why PCR-free library preps are preferred for whole genomes when input DNA allows. Short-read sequencing, as a technology, works by massively parallel imaging of clonal clusters rather than by measuring bulk amplification ¹.

Nanopore and PacBio produce different raw formats (POD5 signal data and unaligned BAM with per-base kinetics, respectively), and the alignment step differs. This guide assumes short reads, 2x150 bp, which is what nearly all consumer and research whole-genome sequencing uses today.

2. Verify the files before you spend compute on them

Every delivery should include checksums. Check them first, because a truncated gzip stream will fail 5 hours into alignment rather than immediately.

md5sum -c delivery.md5
# then confirm the gzip streams are intact and complete
for f in *.fastq.gz; do gzip -t "$f" || echo "CORRUPT: $f"; done

Then get the shape of the data. seqkit is fast and gives you what you need in one pass:

seqkit stats -a -j 8 *.fastq.gz

For a 30x human genome at 2x150 bp you expect roughly 300 Gbp of raw sequence total, which is on the order of 400 million read pairs per lane-set, and compressed FASTQ in the range of 40 to 60 GB for R1 and R2 together. If your total base count divided by 3.1 billion is well under 30, you have less coverage than you were told, and coverage is the single parameter that most determines whether variant calls are trustworthy. Read the identifier line too: Illumina’s format encodes instrument, run, flowcell, lane, and the index sequence, which you will need for read groups in the next step.

3. Trim adapters and filter, but do it conservatively

We use fastp for a single pass that trims adapters, drops unsalvageable reads, and writes an HTML report. Aggressive quality trimming is a mistake for germline variant calling, because modern aligners handle low-quality tails with soft-clipping and base quality score recalibration accounts for the rest. Over-trimming biases coverage at read ends and can remove the very mismatches that signal a real variant.

fastp \
  -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
  -o sample_R1.trim.fastq.gz -O sample_R2.trim.fastq.gz \
  --detect_adapter_for_pe \
  --cut_tail --cut_tail_mean_quality 15 \
  --length_required 50 \
  --thread 16 \
  --html sample.fastp.html --json sample.fastp.json

Read the report rather than filing it. The first number to look at is the duplication rate, where over about 20% on a PCR-free 30x library suggests low input DNA. The next is the proportion of bases at or above Q30, which is typically above 90% on a healthy run. The last is the GC distribution, where a second peak often means bacterial or adapter-dimer contamination. If you also ran FastQC, its per-tile quality module will show you whether a flowcell region failed, which is invisible in aggregate statistics.

4. Align to a reference, and pick the reference deliberately

Use GRCh38 with alt contigs and the decoy sequences (the “full analysis set”), or the T2T-CHM13v2.0 assembly if you care about centromeres, segmental duplications, and the short arms of acrocentric chromosomes. GRCh38 remains the pragmatic default because every annotation resource is keyed to it. Do not mix references between steps, and do not mix chromosome naming conventions: chr1 from UCSC and 1 from Ensembl will silently produce empty intersections downstream.

bwa-mem2 index -p GRCh38 GCA_000001405.15_GRCh38_full_analysis_set.fna

bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:HXXXX.1\tSM:subject01\tLB:lib1\tPL:ILLUMINA\tPU:HXXXX.1.ATCACG' \
  GRCh38 sample_R1.trim.fastq.gz sample_R2.trim.fastq.gz \
| samtools sort -@ 8 -m 4G -o sample.sorted.bam -

The -R read group is not optional. Variant callers use SM to name the sample and duplicate markers use LB to decide what counts as a library duplicate. -K 100000000 fixes the chunk size so the output is deterministic across thread counts, which matters when you want a rerun to reproduce a previous result exactly. -Y uses soft clipping for supplementary alignments, which keeps sequence available for structural variant callers later.

Mark duplicates, then convert to CRAM. CRAM stores bases by reference difference and typically cuts file size by roughly half against BAM, at the cost of needing the exact reference FASTA present to read it.

samtools markdup -@ 8 --write-index sample.sorted.bam sample.markdup.bam
samtools view -@ 8 -T GRCh38.fa -C -o sample.cram --write-index sample.markdup.bam

Now measure the alignment rather than trusting it. mosdepth gives coverage distribution in seconds:

mosdepth --by 1000 --fast-mode -t 4 sample sample.cram
samtools flagstat -@ 8 sample.cram

You want above 95% of reads mapped, above 90% properly paired, and a mean autosomal depth close to what you paid for, with at least 90% of the callable genome at 20x or better. The alignment and its quality metrics, not the FASTQ, are the substrate for everything that follows ¹.

5. Call small variants

For a single germline genome, we use DeepVariant. It is a convolutional network over pileup images, it requires no per-sample tuning, and it consistently beats hand-built hard filters on the GIAB benchmark samples for SNVs. GATK HaplotypeCaller followed by VQSR remains the right answer for large joint-called cohorts, where the variant quality model has thousands of samples to learn from. For n=1 it has nothing to learn from, and hard filtering is fragile. The field has produced a long list of callers with meaningfully different error profiles, and the choice matters more than most people assume ².

singularity run -B /data:/data docker://google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \
  --ref=/data/GRCh38.fa \
  --reads=/data/sample.cram \
  --output_vcf=/data/sample.vcf.gz \
  --output_gvcf=/data/sample.g.vcf.gz \
  --num_shards=32

Keep the gVCF. It records confidence at every reference position, which is what lets you later say “this position was covered at 35x and no variant was seen” instead of “this variant was absent from a file,” a distinction that matters whenever you check whether a specific reported variant is present.

Structural variants need separate tools. Run manta for breakpoint-based calls and smoove/lumpy for read-pair and split-read evidence, then take the intersection as your higher-confidence set. Short reads miss a large fraction of structural variation, especially insertions and anything inside repeats, and long-read sequencing is the fix for that class of question.

Sanity-check the call set before interpretation:

bcftools stats -F GRCh38.fa sample.vcf.gz > sample.bcfstats

A whole genome from one person should yield roughly 4 to 5 million small variants, a transition/transversion ratio near 2.0 to 2.1 genome-wide, and a heterozygous-to-homozygous-alternate ratio around 1.5 to 2.0. A Ti/Tv of 1.5 means your filtering let noise through. A het/hom ratio far below 1 suggests either a reference or contamination problem. Run somalier or verifyBamID2 to estimate cross-sample contamination and to confirm reported sex from X and Y coverage.

6. Annotate, and restrict what you look at

Annotation is where a sound-looking file becomes a misleading one if you are careless. Use Ensembl VEP with the offline cache, pick a single consistent transcript set, and record the version of every database.

vep -i sample.vcf.gz -o sample.vep.vcf.gz \
  --cache --offline --dir_cache /data/vep --assembly GRCh38 \
  --vcf --compress_output bgzip --fork 16 \
  --everything --pick_allele_gene \
  --plugin CADD,/data/cadd/whole_genome_SNVs.tsv.gz \
  --custom /data/clinvar.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT

Then filter hard. Of the millions of variants in your file, the interpretable fraction is tiny. An interpretable variant is rare, meaning a gnomAD allele frequency below roughly 0.1% for dominant conditions. It sits in a gene with established disease association, and it is classified by multiple submitters in ClinVar with the review status to match. We start from ClinVar entries with two or more stars and ignore everything else on a first pass.

bcftools view -i 'INFO/ClinVar_CLNSIG ~ "Pathogenic" && INFO/ClinVar_CLNREVSTAT ~ "multiple_submitters"' \
  sample.vep.vcf.gz | bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%INFO/CSQ\n'

Two limits deserve stating plainly. First, most rare variants you find will be variants of uncertain significance, and interpretation depends on population frequency data, functional evidence, and family context that a single genome cannot supply ³. Second, a pathogenic-looking call from a research pipeline is a hypothesis. It is not a result. Short-read calls in homologous regions such as PMS2, SMN1, and the HLA locus are unreliable by construction, and clinical confirmation by an orthogonal method in an accredited laboratory is the only basis for medical action. Sequencing your genome also surfaces findings you did not ask for, including carrier status and risk alleles that carry implications for relatives, which is worth deciding how to handle before you look rather than after ⁴.

7. Quantify RNA, if you have it

RNA-seq answers a different question: which transcripts are present and at what abundance in the tissue sampled, usually whole blood or PBMCs. For expression quantification we use salmon in selective-alignment mode, which is fast, accounts for GC and positional bias, and does not require you to keep a genome BAM.

salmon index -t gencode.v45.transcripts.fa.gz -i salmon_idx --gencode -k 31 -p 16

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

Check the mapping rate in rna_quant/logs/salmon_quant.log. Below about 70% for a poly-A blood library usually means rRNA or globin carryover rather than a pipeline error. Aggregate transcripts to genes with tximport in R, and remember that a single RNA sample gives you relative abundance at one moment under one set of conditions. Expression without a longitudinal baseline or a comparison group is difficult to interpret. If you want to map expression back to your own variants, the use is checking whether a candidate splice variant shows the predicted junction, which requires a genome-aligned BAM from STAR and inspection in IGV rather than a transcript abundance table.

Common problems

Reference mismatches are the most frequent silent failure. A CRAM cannot be read without the exact reference it was written against, and an md5-mismatched FASTA will produce garbage rather than an error. Keep the reference next to the data and record its checksum.

Coverage that looks adequate on average can be inadequate where it matters. Mean 30x hides the fact that GC-extreme first exons, and any region with a segmental duplication, may sit at 5x or zero. Always check per-gene coverage with mosdepth --by targets.bed before concluding a gene is variant-free.

Contamination and sample swaps are more common than people expect in any multi-sample workflow. Run somalier relate across every sample you handle and confirm that genotype concordance and inferred sex match what you expect. A 2% contamination fraction is enough to generate spurious low-allele-fraction calls that look like mosaicism.

Duplicate rates above 25% on a PCR-based library mean your effective coverage is far below nominal, because duplicates are removed before calling. The fix is more input DNA or a PCR-free prep. More sequencing of the same library does not help.

Allele-specific dropout produces false homozygous calls when a variant under a primer or probe site suppresses one allele. This is a targeted-panel problem more than a whole-genome one, but it appears in whole genomes at sites with dense variation.

Finally, treat the data as identifying. A genome is the most persistent identifier a person has, sufficient on its own to link a sample to named relatives through genealogy databases, and it does not expire or rotate ⁵. Encrypt the volume, keep the raw FASTQ off shared machines, and think carefully before uploading to any analysis service whose data retention terms you have not read.

Footnotes

1. Rachel L Goldfeder, Dennis P Wall, Muin J Khoury, et al. Human Genome Sequencing at the Population Scale: A Primer on High-Throughput DNA Sequencing and Analysis. American Journal of Epidemiology, 2017. https://doi.org/10.1093/aje/kww224 ↩ ↩²

2. S. Pabinger, A. Dander, M. Fischer, et al. A survey of tools for variant analysis of next-generation genome sequencing data. Briefings in Bioinformatics, 2013. https://doi.org/10.1093/bib/bbs086 ↩

3. Tuuli Lappalainen, Alexandra J. Scott, Margot Brandt, et al. Genomic Analysis in the Age of Human Genome Sequencing. Cell, 2019. https://doi.org/10.1016/j.cell.2019.02.032 ↩

4. Wim Pinxten, Heidi Carmen Howard. Ethical issues raised by whole genome sequencing. Best Practice & Research Clinical Gastroenterology, 2014. https://doi.org/10.1016/j.bpg.2014.02.004 ↩

5. Andreas Tillmar, Peter Sjölund, Bo Lundqvist, et al. Whole-genome sequencing of human remains to enable genealogy DNA database searches – A case report. Forensic Science International: Genetics, 2020. https://doi.org/10.1016/j.fsigen.2020.102233 ↩