Bio Chipseq Allele Specific Binding

Version Compatibility

Reference examples tested with: WASP 0.3.4+, RASQUAL 1.1+, BaalChIP 1.30+ (Bioconductor), AlleleSeq 2.0+, samtools 1.19+, bcftools 1.19+, GATK 4.5+, pysam 0.22+.

Allele-Specific Binding (ASB)

"Identify variants that affect transcription factor or histone modification binding in cis" -> Compare ChIP-seq read counts at the reference and alternate alleles of heterozygous variants in a single sample. Differential read counts (ALT vs REF at hetSNPs in peaks) reveal allele-specific binding.

• CLI (mandatory bias filter): WASP mapping pipeline to remove reference-allele mapping bias
• CLI (joint association): RASQUAL for cis-QTL + ASB (genotype VCF piped in via tabix)
• R (Bayesian beta-binomial): BaalChIP with copy-number-aware overdispersion
• CLI (personalized genome): AlleleSeq with phased diploid genome
• Statistical test: beta-binomial likelihood ratio or chi-squared on count tables

ASB analysis has three universal pitfalls: reference-allele mapping bias (universal across short-read aligners), imprinted loci (constitutively allele-skewed by biology), and copy-number variation (changes effective allele dose). All three must be addressed or results are unreliable.

Method Taxonomy

Universal First Step: WASP Reference-Bias Filter

Goal: Remove reads that show reference-allele mapping bias before any ASB testing.

Approach: Align reads, identify those overlapping heterozygous SNPs, swap alleles and re-align; reads that don't map consistently to the same position with both alleles are discarded. The output is a bias-corrected BAM at the cost of 22-31% read loss.

Reference-allele mapping bias is systematic: reads with the reference allele align more readily because the reference is the alignment target. This inflates REF allele frequency by 1-5% genome-wide. WASP fixes this:

# WASP mapping pipeline
# 1. Initial alignment
bowtie2 -x hg38 -1 R1.fq -2 R2.fq -S step1.sam
samtools view -bS step1.sam | samtools sort -o step1.bam
samtools index step1.bam

# 2. Identify reads overlapping hetSNPs; swap alleles; re-map
python /path/to/WASP/mapping/find_intersecting_snps.py \
    --is_paired_end \
    --is_sorted \
    --output_dir wasp_out/ \
    --snp_tab snps_tab.h5 \
    --snp_index snps_index.h5 \
    --haplotype haplotypes.h5 \
    --samples sample_list.txt \
    step1.bam

# 3. Re-map swapped reads
bowtie2 -x hg38 -1 wasp_out/step1.remap.fq.gz -S step2.sam
# (process step2.sam similarly)

# 4. Filter reads that don't map back consistently
python /path/to/WASP/mapping/filter_remapped_reads.py \
    step1.to.remap.bam step2.bam step1.keep.bam

# 5. Final WASP-filtered BAM (use this for all downstream ASB analysis)
samtools sort -o step1.wasp.bam step1.keep.bam
samtools index step1.wasp.bam

WASP always drops 22-31% of reads. This is the cost of bias correction; downstream power is reduced but ASB calls are trustworthy.

Alternative to WASP filter: RASQUAL's phi parameter models bias within the test rather than filtering reads. More sophisticated but assumes binomial bias structure.

Workflow: BaalChIP (Recommended for Cancer / Copy-Number-Imbalanced Samples)

library(BaalChIP)
library(BSgenome.Hsapiens.UCSC.hg38)

# Sample metadata
samples <- data.frame(
    SampleID = c('HCC1395_FOXA1_rep1', 'HCC1395_FOXA1_rep2'),
    Tissue = 'TNBC',
    Target = 'FOXA1',
    BAM = c('rep1.wasp.bam', 'rep2.wasp.bam'),
    Peaks = c('rep1_peaks.bed', 'rep2_peaks.bed'),
    Group = 'HCC1395'
)

# hetSNP file: VCF or BED with chrom, pos, ref, alt, allele frequencies
hetSNPs <- 'het_snps.bed'

# CNV file for copy-number-aware overdispersion (critical for cancer)
cnvs <- 'cnvs.bed'

# Initialize BaalChIP object
res <- BaalChIP(samplesheet = samples, hets = hetSNPs)

# Run filters and Bayesian test
res <- alleleCounts(res, min_base_quality = 10, min_mapq = 15)
res <- QCfilter(res, RegionsToFilter = list(blacklist = rtracklayer::import('blacklist_v2.bed')))
res <- mergePerGroup(res)
res <- filter1allele(res)
res <- getASB(res, Iter = 5000, conf_level = 0.95)
# Verify parameter names against the installed BaalChIP version (`?getASB`); some releases
# use `nIter` instead of `Iter`.

# Results
asb_table <- BaalChIP.report(res)
head(asb_table)

BaalChIP outputs per-hetSNP: allelic ratio (AR), bias-corrected ratio (Corrected.AR), a Bayesian credible interval (Bayes_lower/Bayes_upper), and the ASB call (isASB).

Workflow: RASQUAL (Joint cis-QTL + ASB)

# Prepare input
# - BAM filtered by WASP
# - Genotype VCF (phased)
# - Peak BED

# Run RASQUAL. The genotype VCF is piped in from tabix (there is NO --vcf flag);
# options are single-dash. Inputs -y/-k/-x are BINARY files built by RASQUAL's
# txt2bin utilities (not .txt). One run per feature: -j selects the feature row,
# -l = number of test (cis) SNPs, -m = number of feature SNPs in the window.
tabix genotypes.vcf.gz chr:start-end | \
  rasqual -y Y.bin -k K.bin -x X.bin \
    -n N_samples -j FEATURE_INDEX -l N_TEST_SNPS -m N_FEATURE_SNPS \
    -s EXON_STARTS -e EXON_ENDS -f PEAK_ID \
    > rasqual_results.txt

# Output columns: chrom, peak_id, n_RSNPs, n_FSNPs, n_imputed, summarized_phi,
#                 summarized_overdispersion, summarized_pi, beta, log10_BF, ...

RASQUAL's phi parameter is the per-feature bias estimate; pi is the allelic ratio.

Workflow: AlleleSeq (Personalized Diploid Genome)

# Build personalized diploid genome from phased VCF
java -jar vcf2diploid.jar -id SAMPLE -chr hg38.fa -vcf SAMPLE.phased.vcf   # per-haplotype FASTAs written to CWD (no -outDir option)

# Align reads to both maternal and paternal copies
bowtie2-build personalized/maternal.fa maternal_index
bowtie2-build personalized/paternal.fa paternal_index
bowtie2 -x maternal_index -1 R1.fq -2 R2.fq -S maternal.sam
bowtie2 -x paternal_index -1 R1.fq -2 R2.fq -S paternal.sam

# AlleleSeq2 pipeline (Makefile-based; there is no AlleleSeq2.pl entry point)
make -f PIPELINE.mk PGENOME_DIR=personalized/ REFGENOME_VERSION=GRCh38 ALIGNMENT_MODE=ASB NTHR=8

# Output: per-hetSNP allelic counts and binomial test

Personalized genome avoids reference bias by construction. Cost: per-sample diploid genome generation and indexing.

Three Universal Pitfalls

Pitfall 1: Imprinted Loci Are Constitutively Skewed

Imprinted loci (H19, IGF2, MEG3, MEG8, KCNQ1OT1, etc.) show extreme allele bias by biology, not from differential binding.

# Filter imprinted loci before ASB analysis
# Imprinted-gene coordinates are derived from a catalog (geneimprint.com or the
# Otago Imprinted Gene Catalogue, igc.otago.ac.nz) mapped to hg38; there is no
# canonical hosted hg38 BED. Given imprinted_loci_hg38.bed:
bedtools intersect -v -a hetSNPs.bed -b imprinted_loci_hg38.bed > hetSNPs.non_imprinted.bed

Pitfall 2: X-Inactivation in Females

In female samples, X-linked genes show extreme allele skew because each cell silences one X chromosome. This appears as ASB at every X-linked hetSNP.

# Filter chrX in female samples
awk '$1 != "chrX"' hetSNPs.bed > hetSNPs.autosomal.bed
# Or analyze chrX separately with imprinting-aware methods

Pitfall 3: Copy-Number Imbalance (Cancer Genomes)

In cancer cells, copy-number gain of one allele alters effective allele dose; raw allelic ratios mix dose and binding effects. BaalChIP's copy-number-aware overdispersion handles this; other methods require pre-filtering CN-altered regions.

# Use ASCAT / Sequenza / FACETS to call allele-specific CNVs
# Exclude CN-altered regions from ASB analysis OR use BaalChIP

Per-Tool Failure Modes

WASP -- Reference panel mismatch

Trigger: Using a WASP SNP file from a different population than the sample.

Mechanism: WASP swaps alleles at known hetSNPs; if the variant isn't in the SNP file, no swap happens; reads retain reference bias.

Symptom: Sample-specific hetSNPs (not in 1KG) still show reference bias after WASP.

Fix: Build WASP SNP file from the sample's own genotype VCF, not a population panel; OR use RASQUAL which handles novel hetSNPs.

WASP -- Excessive read loss

Trigger: WASP filter removes >40% of reads.

Mechanism: Many reads span multiple hetSNPs; each must re-map consistently after every allele swap; combinatorial loss.

Fix: Accept the loss (genuine bias correction) OR switch to AlleleSeq (personalized genome avoids the swap-and-remap step) OR RASQUAL (no read filtering).

RASQUAL -- Convergence failure

Trigger: Sparse data (few hetSNPs per peak); strong copy-number imbalance.

Mechanism: EM convergence requires enough hetSNPs per feature; sparse data underspecifies the model.

Fix: Require well-imputed SNPs (--imputation-quality-fsnp); combine replicates; or switch to BaalChIP for sparse-data robustness.