# Bio Applied Snp Calling Pipeline

> SNP calling pipeline: Trimmomatic trim, BWA-MEM2/HISAT2 align, samtools mpileup + bcftools call, ANNOVAR annotate (dbSNP, RefGene, 1000G, ClinVar). Use for FASTQ-to-VCF pipelines or ANNOVAR variant annotation.

- Skill: `pavel-kravchenko/bio-applied-snp-calling-pipeline` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-snp-calling-pipeline`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-snp-calling-pipeline/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: DevOps & Infra
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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

```bash
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

```text
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.

```bash
#!/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.

```python
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