bio-applied-snp-calling-pipeline

name	bio-applied-snp-calling-pipeline
description	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.
tool_type	bash
primary_tool	bcftools

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 

 
perl annovar/convert2annovar.pl -format vcf4  > 
perl annovar/annotate_variation.pl -filter -out  -build hg19 \
    -dbtype snp138  annovar/humandb/
perl annovar/annotate_variation.pl -out  -build hg19 \
     annovar/humandb/
perl annovar/annotate_variation.pl -filter -out  -buildver hg19 \
    -dbtype 1000g2014oct_all  annovar/humandb/
perl annovar/annotate_variation.pl -regionanno -out  -build hg19 \
    -dbtype gwasCatalog  annovar/humandb/
perl annovar/annotate_variation.pl -filter -out  -buildver hg19 \
    -dbtype clinvar_20221231  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  dp4_str.split())
    alt_total = af + ar
     alt_total == :
         
     af / alt_total


 ():
    
    flagged = []
     r  records:
        dp4 = r[].get()
        bias = strand_bias_ratio(dp4)  dp4  
        flagged.append({
            **r,
            : r[] >= min_qual,
            : bias,
            : bias    (bias_low <= bias <= bias_high),
        })
     flagged


 __name__ == :
    
    demo = [
        ,
        ,
        ,
    ]
    (, ).write(.join(demo) + )
    recs = parse_vcf_records()
    flagged = flag_variants(recs)
     flagged[][]     
     flagged[][]    
    (, [(f[], f[], f[])  f  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