Haplotype Phasing

## ================================================================
## SHAPEIT4 statistical phasing
## Prerequisites: SHAPEIT4 v4.2+, bcftools, tabix
## Dataset: Rice 3K chr04 multi-sample VCF
## ================================================================

## Step 0: Pre-flight checks
## SHAPEIT4 requires:
##   1. Multi-sample VCF (biallelic SNPs, genotyped calls -- not missing)
##   2. Genetic map in PLINK format (chr pos_bp cM pos_bp)
##   3. Optionally: a pre-phased reference panel for scaffold phasing

VCF="rice3k_chr04_filtered.vcf.gz"
MAP="chr04_genetic_map_plink.txt"   ## 4-column: chr snpID cM bp
OUTVCF="rice3k_chr04_phased.vcf.gz"

## Step 1: Filter to high-quality biallelic SNPs with <10% missing
bcftools view ${VCF} \
  --min-alleles 2 --max-alleles 2 \
  --type snps \
  --min-ac 1:minor | \
bcftools filter \
  -e 'F_MISSING > 0.1' \
  -o rice3k_chr04_forPhasing.vcf.gz \
  -O z
tabix -p vcf rice3k_chr04_forPhasing.vcf.gz

## Step 2: Run SHAPEIT4
shapeit4.2 \
  --input rice3k_chr04_forPhasing.vcf.gz \
  --map ${MAP} \
  --region chr04 \
  --output ${OUTVCF} \
  --thread 16 \
  --mcmc-iterations 5b,1p,1b,1p,1b,1p,5m \  ## burn-in+pruning+main iterations
  --pbwt-depth 4 \                             ## number of PBWT neighbours to condition on (4-8 typical)
  --pbwt-mac 2 \                               ## min MAC for PBWT conditioning
  --log shapeit4_chr04.log

tabix -p vcf ${OUTVCF}

## Step 3: Verify phasing -- check PS (phase set) tags and switch to | notation
bcftools view -h ${OUTVCF} | grep "FORMAT"
## Phased genotypes are coded as 0|1 or 1|0 (pipe = phased)
## Unphased remaining genotypes are 0/1 (slash = unphased)

## Count phased vs unphased sites
bcftools stats ${OUTVCF} | grep "^PSC"  ## per-sample phasing stats
bcftools view ${OUTVCF} | grep -v "^#" | \
  awk '{for(i=10;i<=NF;i++) if($i~/\|/) phased++; else unphased++}
  END {print "Phased:", phased, "Unphased:", unphased,
       "Rate:", phased/(phased+unphased)}'

## Step 4: Scaffold with reference panel (optional -- improves rare variant phasing)
## Use 1000G phased reference or in-house phased panel as scaffold
shapeit4.2 \
  --input rice3k_chr04_forPhasing.vcf.gz \
  --map ${MAP} \
  --region chr04 \
  --reference 1kg_phased_reference_chr04.vcf.gz \  ## pre-phased reference
  --output rice3k_chr04_phased_scaffolded.vcf.gz \
  --thread 16

## ================================================================
## WhatsHap read-backed phasing
## Works with short reads, PacBio HiFi, Nanopore, and trio data
## ================================================================

BAM="sample01_sorted_markdup.bam"
VCF="sample01_genotyped.vcf.gz"
REF="IRGSP1.0_genome.fasta"
OUT="sample01_phased.vcf.gz"

## Step 1: Basic read-backed phasing (single individual)
whatshap phase \
  --output ${OUT} \
  --reference ${REF} \
  --indels \            ## also phase insertions/deletions
  --distrust-genotypes \ ## re-genotype using read evidence (improves accuracy)
  ${VCF} ${BAM}
bgzip ${OUT}; tabix -p vcf ${OUT}.gz

## Step 2: Phasing statistics
whatshap stats \
  --gtf phase_blocks.gtf \   ## optional: output phase blocks as GTF for IGV
  ${OUT}.gz

## Step 3: Trio phasing (parents remove switch errors completely)
## Requires: proband VCF, father BAM, mother BAM
whatshap phase \
  --ped pedigree.ped \       ## pedigree file: child father mother sex
  --output trio_phased.vcf.gz \
  --reference ${REF} \
  child.vcf.gz child.bam father.bam mother.bam

## Step 4: Haplotag BAM (tag each read with its inferred haplotype HP:i:1 or HP:i:2)
## Required for haplotype-specific methylation, allele-specific expression
whatshap haplotag \
  --output sample01_haplotagged.bam \
  --reference ${REF} \
  ${OUT}.gz ${BAM}
samtools index sample01_haplotagged.bam
## Now reads carry HP:1 or HP:2 tags -- usable in methylation and ASE analyses

## Step 5: Combine statistical + read-backed phasing
## Best practice: run SHAPEIT4 first (population-level), then refine with WhatsHap
## This gives chromosome-scale blocks with read-level accuracy within blocks
shapeit4.2 --input sample01.vcf.gz --map chr04.map \
  --region chr04 --output sample01_stat_phased.vcf.gz --thread 8

whatshap phase \
  --reference ${REF} \
  --chromosome chr04 \
  --output sample01_combined_phased.vcf.gz \
  sample01_stat_phased.vcf.gz ${BAM}  ## WhatsHap refines the SHAPEIT4 phase

## Step 6: Long-read phasing with HiFi (PacBio)
## No genetic map needed; reads are long enough to bridge most het sites
whatshap phase \
  --output sample01_hifi_phased.vcf.gz \
  --reference ${REF} \
  --indels \
  sample01.vcf.gz sample01_hifi.bam   ## HiFi BAM replaces short-read BAM

## ================================================================
## Haplotype block detection using PLINK and LDBlockShow
## ================================================================
library(dplyr); library(ggplot2)

## Step 1: Compute pairwise LD with PLINK2
## plink2 --vcf phased.vcf.gz --r2-phased --ld-window-kb 500 \
##   --ld-window-r2 0.05 --out chr04_ld_phased --allow-extra-chr --double-id

## Step 2: Define blocks using Gabriel 2002 method in PLINK1.9
## plink --vcf phased.vcf.gz --blocks no-pheno-req --blocks-max-kb 500 \
##   --out chr04_blocks --allow-extra-chr --double-id
## Output: chr04_blocks.blocks.det -- block start, end, n_SNPs, kb

## Step 3: Visualise block structure in R
blocks <- read.table("chr04_blocks.blocks.det", header=TRUE)
blocks$size_kb <- (blocks$BP2 - blocks$BP1) / 1000

## Block size distribution
ggplot(blocks, aes(size_kb)) +
  geom_histogram(bins=60, fill="#2171b5", colour="white", linewidth=0.2) +
  scale_x_log10() +
  labs(x="Block size (kb, log scale)", y="Count",
       title=paste0("Haplotype block size distribution (n=", nrow(blocks), " blocks)")) +
  theme_bw(base_size=12)
ggsave("block_size_distribution.png", width=7, height=5, dpi=200)

cat("Block summary:\n")
cat("  N blocks:", nrow(blocks), "\n")
cat("  Median size (kb):", round(median(blocks$size_kb),1), "\n")
cat("  Max size (kb):", round(max(blocks$size_kb),1), "\n")
cat("  Total covered (Mb):", round(sum(blocks$size_kb)/1000, 1), "\n")

## Step 4: LD heatmap for a candidate region (e.g. Fst peak)
## Use LDheatmap package for focused visualisation
## library(LDheatmap); library(genetics)
## genos <- read.table("candidate_region_genotypes.txt", header=TRUE)
## ldh <- LDheatmap(genos, distances=positions, color=heat.colors(20))