Step 3: SNP Quality Control & VCF Export

After removing low-quality samples (Steps 1–2), this step applies variant-level quality control. Three filters are applied simultaneously: minimum minor allele frequency (MAF), Hardy-Weinberg equilibrium (HWE) test, and per-SNP genotyping rate. These filters remove uninformative, technically artefactual, or poorly genotyped variants before imputation.

Accurate imputation depends critically on having a clean set of typed variants. Low-quality or artefactual typed SNPs propagate errors into imputed genotypes — a phenomenon known as garbage in, garbage out. Pre-imputation QC is therefore more stringent than general analysis QC.

Key Metrics

2a. Minor Allele Frequency (MAF)

MAF is the frequency of the less common allele at a biallelic locus. Variants with very low MAF (<1%) are removed because: (1) they provide little discriminatory power in frequency-based analyses, (2) at very low MAF a single genotyping error dramatically shifts the observed frequency, and (3) imputation accuracy for rare variants in the typed panel is limited.

2b. Hardy-Weinberg Equilibrium (HWE)

Extreme departure from HWE (p < 10⁻⁶) at a given locus almost always indicates a genotyping artefact rather than biology — e.g., miscalling heterozygotes as homozygotes due to poor cluster separation, copy number variation producing excess heterozygosity, or probe-specific technical failures. The lenient 10⁻⁶ threshold avoids removing loci under genuine selection or modest population substructure effects.

2c. Per-SNP Genotyping Rate

Per-SNP missingness measures the fraction of individuals for whom a variant could not be called. The --geno 0.10 threshold requires ≥90% call rate per SNP (≤109 missing out of 1,093). This matches the original winter 2025 pipeline parameters, recovered forensically from the data (max lmiss among retained variants = 0.099270, confirming the 0.10 cutoff).

Threshold Sensitivity (spring 2026 data, 1,093 samples)

Each row shows the combined result of --maf 0.01 --hwe 1e-6 --geno X --chr 1-22 applied to the 620,492 autosomal variants in ConvSK_mind20_dedup.

2d. Autosome Restriction

Distributions of the three QC metrics across all 620,492 autosomal input variants (computed on the 1,093 post-dedup samples). Dashed red lines indicate applied thresholds.

Input

Output — QC-filtered PLINK

Output — Final VCF for imputation

Output — Michigan-ready VCF (after fixref)

Step 3a: Apply combined SNP QC filters

cd /staging/ALSU-analysis/spring2026

plink --bfile ConvSK_mind20_dedup \
  --maf 0.01 \
  --hwe 1e-6 \
  --geno 0.10 \
  --chr 1-22 \
  --make-bed \
  --out ConvSK_mind20_dedup_snpqc

Step 3b: Build VCF-safe SNP list

# 479 variants have I/D alleles (insertion/deletion shorthand)
# These are not valid VCF alleles — must be excluded
# BIM cols: $1=chr  $2=snpID  $3=cM  $4=pos  $5=allele1  $6=allele2
awk '($5!="I" && $5!="D" && $6!="I" && $6!="D") {print $2}' \
  ConvSK_mind20_dedup_snpqc.bim > vcf_ok_snps.txt
wc -l vcf_ok_snps.txt
# 459344

Step 3c: Export per-chromosome VCFs

# --recode vcf bgz: export PLINK binary as bgzip-compressed VCF (.vcf.gz)
# tabix -f -p vcf: create tabix index (-f = overwrite, -p vcf = VCF format preset)
for chr in $(seq 1 22); do
  plink --bfile ConvSK_mind20_dedup_snpqc \
    --chr $chr --extract vcf_ok_snps.txt \
    --recode vcf bgz --out chr${chr}
  tabix -f -p vcf chr${chr}.vcf.gz
done

Step 3d: Concatenate and clean VCF

# bcftools concat: join VCFs from different regions (same samples)
#   -Oz: bgzip-compressed VCF output
bcftools concat -Oz -o impute_in.autosomes.vcf.gz chr{1..22}.vcf.gz
tabix -f -p vcf impute_in.autosomes.vcf.gz
# 459,344 variants

# Keep biallelic SNPs only (no multiallelic found)
# -m2 -M2: require exactly 2 alleles (min 2, max 2 = biallelic)
# -v snps: retain only SNP variant type (exclude indels/MNPs)
bcftools view -m2 -M2 -v snps -Oz \
  -o impute_in.autosomes.clean.vcf.gz impute_in.autosomes.vcf.gz
tabix -f -p vcf impute_in.autosomes.clean.vcf.gz
# 459,344 variants

# Remove duplicate positions (390 removed)
# bcftools norm -d all: remove all types of duplicate records at same position
bcftools norm -d all -Oz \
  -o impute_in.autosomes.clean.nodup.vcf.gz \
  impute_in.autosomes.clean.vcf.gz
tabix -f -p vcf impute_in.autosomes.clean.nodup.vcf.gz
# 458,954 variants

Step 3e: Fix non-ASCII sample IDs

# Auto-detect non-ASCII sample IDs and generate rename mapping
# bcftools query -l: list all sample names from VCF header
# python3 inline: map all Cyrillic homoglyphs → Latin equivalents
bcftools query -l impute_in.autosomes.clean.nodup.vcf.gz | \
  python3 -c "
import sys
# Cyrillic-to-Latin homoglyph map (visually identical characters)
CYR2LAT = {
    '\u0410':'A', '\u0412':'B', '\u0415':'E', '\u041a':'K', '\u041c':'M',
    '\u041d':'H', '\u041e':'O', '\u0420':'P', '\u0421':'C', '\u0422':'T',
    '\u0425':'X', '\u0430':'a', '\u0435':'e', '\u043e':'o', '\u0440':'p',
    '\u0441':'c', '\u0445':'x', '\u043c':'m',
}
count = 0
for line in sys.stdin:
    old = line.rstrip('\n')
    new = ''.join(CYR2LAT.get(c, c) for c in old)
    if new != old:
        remaining = [c for c in new if ord(c) > 127]
        if remaining:
            print(f'WARNING: {old} still has unmapped non-ASCII: {remaining}', file=sys.stderr)
        print(f'{old}\t{new}')
        count += 1
print(f'{count} rename mapping(s) generated', file=sys.stderr)
" > rename_samples.tsv

cat rename_samples.tsv
# Spring 2026 dataset — expected 2 mappings:
#   808_03-25м        808_03-25m
#   1038_08-176Х-00006  1038_08-176X-00006

# Apply rename only if mappings were found
if [ -s rename_samples.tsv ]; then
  # bcftools reheader -s FILE: replace sample names using old→new tsv
  bcftools reheader -s rename_samples.tsv \
    -o impute_in.final.vcf.gz \
    impute_in.autosomes.clean.nodup.vcf.gz
else
  # No non-ASCII IDs found — just rename the file
  mv impute_in.autosomes.clean.nodup.vcf.gz impute_in.final.vcf.gz
fi

# Index and verify
bcftools index impute_in.final.vcf.gz          # create CSI index for random access
bcftools index -n impute_in.final.vcf.gz       # -n: print variant count (quick check)
# 458954

# Confirm no non-ASCII IDs remain
# grep -cP '[^\x00-\x7F]': count lines containing non-ASCII bytes
bcftools query -l impute_in.final.vcf.gz | grep -cP '[^\x00-\x7F]'
# 0

Step 3f: Split back to per-chromosome VCFs

The concatenation + cleanup + reheader steps above produced a single genome-wide VCF. Michigan Imputation Server requires one VCF per chromosome, so we split impute_in.final.vcf.gz back into per-chromosome files:

# bcftools view -r CHR: extract records from region (chromosome) CHR
for chr in $(seq 1 22); do
  bcftools view -r $chr impute_in.final.vcf.gz \
    -Oz -o impute_chr${chr}.vcf.gz
  bcftools index impute_chr${chr}.vcf.gz
done

Step 3g: Chromosome rename, reference strand alignment & palindromic removal

# Create chromosome rename file (1 → chr1, ..., 22 → chr22)
for i in $(seq 1 22); do echo "$i chr$i"; done > chr_rename.tsv

mkdir -p michigan_ready

# For each chromosome: rename → fixref → remove palindromic
for chr in $(seq 1 22); do
  # --rename-chrs FILE: 2-column mapping (old_name  new_name) for chromosome IDs
  bcftools annotate --rename-chrs chr_rename.tsv \
    impute_chr${chr}.vcf.gz -Oz -o tmp_chr${chr}.vcf.gz
  bcftools index tmp_chr${chr}.vcf.gz

  # bcftools +fixref: plugin that aligns alleles to reference genome forward strand
  #   -Ou: uncompressed BCF output (fast for piping, never touches disk)
  #   --: separator between bcftools args and plugin-specific args
  #   -f FILE: reference FASTA for strand determination
  #   -m top: strand matching mode (Illumina TOP/BOT convention)
  # bcftools view -e EXPR: exclude records matching the filter expression
  #   Removes palindromic A/T and C/G SNPs (ambiguous strand, unresolvable)
  # bcftools index -t: create tabix (.tbi) index (vs default CSI)
  bcftools +fixref tmp_chr${chr}.vcf.gz -Ou \
    -- -f /staging/Genomes/Human/chr/GRCh38.fa -m top \
    2>fixref_chr${chr}.log \
  | bcftools view \
    -e 'TYPE="snp" && ((REF="A" && ALT="T") || (REF="T" && ALT="A") || (REF="C" && ALT="G") || (REF="G" && ALT="C"))' \
    -Oz -o michigan_ready/impute_chr${chr}.vcf.gz

  bcftools index -ft michigan_ready/impute_chr${chr}.vcf.gz
done

# Clean up intermediates (tmp_chr = 3g rename intermediates)
rm -f tmp_chr*.vcf.gz tmp_chr*.vcf.gz.csi

Per-chromosome counts (after fixref + palindromic removal)

# Verify PLINK QC variant count
wc -l ConvSK_mind20_dedup_snpqc.bim
# 459823

# --freq: output allele frequencies (.frq)
# --hardy: output Hardy-Weinberg exact test p-values (.hwe)
# --missing: output per-SNP (.lmiss) and per-sample (.imiss) missingness
plink --bfile ConvSK_mind20_dedup_snpqc \
  --freq --hardy --missing \
  --out snpqc_verify

# Check MAF — expect 0
# .frq cols: CHR  SNP  A1  A2  MAF  NCHROBS → $5=MAF
awk 'NR>1 && $5+0 < 0.01' snpqc_verify.frq | wc -l
# 0

# Check HWE — expect 0
# .hwe cols: CHR SNP TEST A1 A2 GENO O(HET) E(HET) P → $9=P
awk 'NR>1 && $9+0 < 1e-6 && $9 != "NA"' snpqc_verify.hwe | wc -l
# 0

# Check per-SNP missingness — expect 0
# .lmiss cols: CHR  SNP  N_MISS  N_GENO  F_MISS → $5=F_MISS
awk 'NR>1 && $5+0 > 0.10' snpqc_verify.lmiss | wc -l
# 0

Per-chromosome variant counts