Bio Applied Snp Calling Pipeline

SNP Calling Pipeline

When to Use

• Building an end-to-end single-end DNA SNP pipeline from raw FASTQ to annotated variant report
• Calling variants with samtools mpileup + bcftools call instead of GATK
• Annotating a VCF against dbSNP, RefGene, 1000 Genomes, GWAS Catalog, and ClinVar with ANNOVAR
• Explaining/parsing mpileup/bcftools output (DP4, strand bias, QUAL) or ANNOVAR output files
• Prioritizing (tiering) called variants for clinical/candidate-gene review

Version Compatibility

• Trimmomatic ≥0.39 (pipeline historically used 0.36 — behavior unchanged)
• BWA-MEM2 ≥2.2.1 (preferred for new DNA projects) or HISAT2 ≥2.2.1 run in DNA mode
• SAMtools / BCFtools ≥1.17
• ANNOVAR ≥2020Jun08 release, humandb/ built for the matching genome build (hg19 or hg38 — never mix)
• Python ≥3.10, numpy/pandas for downstream parsing

Prerequisites

conda create -n snp-pipeline python=3.10
conda activate snp-pipeline
conda install -c bioconda trimmomatic bwa-mem2 samtools bcftools
 # ANNOVAR requires manual registration: http://annovar.openbioinformatics.org/

Prior concepts: FASTQ quality encoding (Phred), SAM/BAM basics, VCF format (bio-variant-calling-vcf-basics).

Pipeline Architecture

Raw FASTQ reads
    -> Trimmomatic (TRAILING:20, MINLEN:50)
    -> BWA-MEM2 (or HISAT2 --no-spliced-alignment --no-softclip)
    -> samtools view/sort/index -> idxstats/depth
    -> samtools mpileup -uf ref.fasta | bcftools call -cv -> raw VCF
    -> ANNOVAR convert2annovar.pl -> annotate_variation.pl (dbSNP138, refGene, 1000G, GWAS, ClinVar)
    -> Tiered, annotated variant report

Tool choice note: HISAT2 is a splice-aware RNA-seq aligner; --no-spliced-alignment --no-softclip makes it behave like an end-to-end DNA aligner but BWA-MEM2 is the standard for new WGS/WES projects and is what GATK Best Practices expects. bcftools call -c (consensus caller) is legacy and fast but less accurate than -m (multiallelic) or a GATK/DeepVariant caller — see bio-variant-calling-gatk-variant-calling for production-grade calling.

Goal: turn raw single-end FASTQ into a sorted, indexed BAM and a raw VCF. Approach: trim low-quality 3' bases, align in DNA mode, coordinate-sort, pileup, call.

#!/bin/bash
set -euo pipefail
sample=$1              # e.g. patientA (expects patientA.fastq + Human/patientA.fasta + BWA-MEM2 index)
ref="Human/${sample}.fasta"

 # 1. Quality trim (remove 3' bases below Q20, drop reads <50bp)
java -jar Trimmomatic-0.39.jar SE -phred33 \
    "${sample}.fastq" "${sample}.trimmed.fastq" TRAILING:20 MINLEN:50

 # 2. Align (BWA-MEM2; -R read group required for downstream GATK compatibility)
bwa-mem2 mem -R "@RG\tID:${sample}\tSM:${sample}\tPL:ILLUMINA" \
    "$ref" "${sample}.trimmed.fastq" > "${sample}.sam"

 # 3. SAM -> sorted, indexed BAM
samtools view -b "${sample}.sam" -o "${sample}.bam"
samtools sort "${sample}.bam" -o "${sample}.sorted.bam"
samtools index "${sample}.sorted.bam"
samtools idxstats "${sample}.sorted.bam" > "${sample}.idxstats.txt"
samtools depth "${sample}.sorted.bam" > "${sample}.depth.tsv"

 # 4. Pileup + call (consensus caller; use -m for cohorts/multiallelic sites)
samtools mpileup -uf "$ref" "${sample}.sorted.bam" | \
    bcftools call -cv -o "${sample}.vcf"

 # 5. VCF -> ANNOVAR input, then annotate against 5 databases
perl annovar/convert2annovar.pl -format vcf4 "${sample}.vcf" > "${sample}.avinput"
perl annovar/annotate_variation.pl -filter -out "${sample}.rs" -build hg19 \
    -dbtype snp138 "${sample}.avinput" annovar/humandb/
perl annovar/annotate_variation.pl -out "${sample}.refgene" -build hg19 \
    "${sample}.avinput" annovar/humandb/
perl annovar/annotate_variation.pl -filter -out "${sample}.1000g" -buildver hg19 \
    -dbtype 1000g2014oct_all "${sample}.avinput" annovar/humandb/
perl annovar/annotate_variation.pl -regionanno -out "${sample}.gwas" -build hg19 \
    -dbtype gwasCatalog "${sample}.avinput" annovar/humandb/
perl annovar/annotate_variation.pl -filter -out "${sample}.clinvar" -buildver hg19 \
    -dbtype clinvar_20221231 "${sample}.avinput" annovar/humandb/

Goal: parse a real VCF produced by bcftools call and compute per-variant QC metrics. Approach: read DP4 (ref-fwd, ref-rev, alt-fwd, alt-rev) from the INFO field to flag strand bias, without any external VCF library.

import re


def parse_vcf_records(vcf_path):
    """Parse a bcftools-style VCF into a list of dicts with CHROM/POS/REF/ALT/QUAL/INFO fields."""
    records = []
    with open(vcf_path) as fh:
        for line in fh:
            if line.startswith("#"):
                continue
            fields = line.rstrip("\n").split("\t")
            chrom, pos, _id, ref, alt, qual, filt, info = fields[:8]
            info_dict = dict(
                kv.split("=", 1) if "=" in kv else (kv, True)
                for kv in info.split(";")
            )
            records.append({
                "chrom": chrom, "pos": int(pos), "ref": ref, "alt": alt,
                "qual": float(qual), "filter": filt, "info": info_dict,
            })
    return records


def strand_bias_ratio(dp4_str):
    """Compute the alt-allele forward/total ratio from a DP4 string 'rf,rr,af,ar'.

    Returns None if there are no alt-supporting reads; values far from 0.5
    indicate the variant is only seen on one strand (likely an artifact).
    """
    rf, rr, af, ar = (int(x) for x in dp4_str.split(","))
    alt_total = af + ar
    if alt_total == 0:
        return None
    return af / alt_total


def flag_variants(records, min_qual=30, bias_low=0.2, bias_high=0.8):
    """Annotate each record with pass/fail QC flags for QUAL and strand bias."""
    flagged = []
    for r in records:
        dp4 = r["info"].get("DP4")
        bias = strand_bias_ratio(dp4) if dp4 else None
        flagged.append({
            **r,
            "qual_pass": r["qual"] >= min_qual,
            "strand_bias": bias,
            "bias_pass": bias is None or (bias_low <= bias <= bias_high),
        })
    return flagged


if __name__ == "__main__":
    # Self-check with a two-line inline VCF (no file I/O needed)
    demo = [
        "##fileformat=VCFv4.1",
        "chr1\t925952\t.\tG\tA\t222\t.\tDP=30;DP4=0,0,15,15;MQ=60",
        "chr1\t930000\t.\tT\tC\t40\t.\tDP=20;DP4=9,0,10,1;MQ=55",
    ]
    open("/tmp/demo.vcf", "w").write("\n".join(demo) + "\n")
    recs = parse_vcf_records("/tmp/demo.vcf")
    flagged = flag_variants(recs)
    assert flagged[0]["bias_pass"] is True   # balanced DP4 -> no bias
    assert flagged[1]["bias_pass"] is False  # 10/11 alt reads forward -> biased
    print("OK:", [(f["pos"], f["strand_bias"], f["bias_pass"]) for f in flagged])

Pitfalls

• Coordinate systems: BED is 0-based half-open; VCF/GFF are 1-based inclusive — off-by-one errors are the most common bug when cross-referencing.
• HISAT2 for DNA: must disable spliced alignment and soft-clipping (--no-spliced-alignment --no-softclip), otherwise false split-read alignments appear as spurious indels.
• BWA/HISAT2 read groups: always set -R/--rg — GATK and most downstream tools require a read group to identify the sample.
• Genome build mismatch: ANNOVAR humandb/ files, the reference FASTA, and the alignment index must all be the same build (hg19 or hg38) — silent wrong annotations otherwise.
• bcftools call -c is legacy: fine for a quick single-sample survey, but has no multiallelic support and is less accurate than GATK HaplotypeCaller/DeepVariant for cohorts or clinical use.
• No duplicate marking or BQSR in this pipeline — for WGS/WES, add samtools markdup/Picard MarkDuplicates before calling.

See Also

• bio-variant-calling-gatk-variant-calling — production-grade germline calling with HaplotypeCaller
• bio-variant-calling-vcf-basics — VCF format, INFO/FORMAT fields, genotype decoding
• bio-variant-calling-variant-annotation — VEP/SnpEff alternatives to ANNOVAR
• bio-read-alignment-bwa-alignment — BWA-MEM2 alignment details and read-group flags