Version Compatibility

Reference examples tested with: WASP 0.3.4+, GATK 4.4+, RASQUAL 1.1+, samtools 1.19+, bcftools 1.19+, vcftools 0.1.16+, plink 2.00+, MatrixEQTL 2.3+, QuASAR 0.1+, bowtie2 2.5+, bwa-mem2 2.2.1+.

Verify before use:

• CLI: <tool> --version then <tool> --help to confirm flags
• Python: pip show <package> then help(module.function) to check signatures
• R: packageVersion('<pkg>') then ?function_name to verify parameters

If code throws unexpected errors, introspect the installed package and adapt rather than retrying.

Allele-Specific Accessibility

"Does this heterozygous SNP affect chromatin accessibility on its allele?" -> Count ATAC reads supporting reference vs alternative allele at heterozygous sites in the same individual; significant deviation from 50:50 indicates cis-regulatory effect. Requires careful handling of reference-allele mapping bias (WASP filtering) and within-individual binomial testing.

• CLI: WASP (Geijn 2015) for de-biased reference mapping
• CLI: gatk ASEReadCounter for allele-specific count tables
• CLI: RASQUAL (Kumasaka 2016) for joint cis-mapping with allelic counts
• R: QuASAR (Harvey 2015) for genotype + ASE inference simultaneously

ASE/ASB analysis is fundamentally different from cohort-level differential. Statistical framework is binomial within-individual; sample size is the count of heterozygous SNPs in accessible regions, not the number of individuals.

Algorithmic Taxonomy

Methodology evolves; verify against current Geijn 2015, Kumasaka 2016, Buchkovich 2015 (review) before locking pipelines. Modern caQTL studies (Kumasaka 2019, Anderson 2024) typically combine WASP + RASQUAL at modest N or WASP + GATK + linear caQTL at large N.

Reference Allele Mapping Bias (The Single Most Important Issue)

When aligning reads to the reference genome, reads carrying the reference allele align with 0 mismatches; reads carrying the alternative allele have 1 mismatch and may fail to align (especially with bwa-mem -k or bowtie2 --very-sensitive thresholds). This inflates the apparent reference-allele frequency at every SNP and confounds ASE.

Trigger: Always (mapping bias is universal at heterozygous SNPs).

Mechanism: Aligners are mismatch-penalized. Without correction, ALT reads systematically under-align. Effect size: 1-5% bias toward reference at typical mismatch penalties.

Symptom: Per-SNP REF allele fraction skews above 50% even at sites without true allelic imbalance.

Fix: WASP. Pseudo-alleles every read at heterozygous sites; re-aligns swapped reads; keeps only reads that map identically to both haplotypes. Mandatory for any ASE/ASB analysis.

Goal: Remove reference-allele mapping bias before counting allele-specific ATAC reads.

Approach: Identify reads overlapping heterozygous SNPs, re-align allele-swapped versions, keep only reads consistent across both haplotypes, then count REF/ALT with GATK ASEReadCounter.

# WASP read-correction pipeline (Geijn 2015)
WASP_DIR=/path/to/WASP
PEAKS=peaks.bed                                  # ATAC peaks for filtering
SAMPLE=sample1
OUT=$SAMPLE.wasp_filtered.bam

# 1. Find reads at SNP sites
python $WASP_DIR/mapping/find_intersecting_snps.py \
    --is_paired_end \
    --is_sorted \
    --output_dir wasp_out/ \
    --snp_dir snp_h5/ \
    $SAMPLE.bam

# 2. Re-align flipped-allele reads
bowtie2 -x hg38_idx -1 wasp_out/$SAMPLE.remap.fq1.gz \
                   -2 wasp_out/$SAMPLE.remap.fq2.gz \
                   -S wasp_out/$SAMPLE.remap.sam
samtools view -bS wasp_out/$SAMPLE.remap.sam | samtools sort -o wasp_out/$SAMPLE.remap.bam
samtools index wasp_out/$SAMPLE.remap.bam

# 3. Keep only consistently-mapped reads
python $WASP_DIR/mapping/filter_remapped_reads.py \
    wasp_out/$SAMPLE.to.remap.bam \
    wasp_out/$SAMPLE.remap.bam \
    wasp_out/$SAMPLE.kept.bam

# 4. Merge kept reads with non-overlapping reads
samtools merge $OUT \
    wasp_out/$SAMPLE.kept.bam \
    wasp_out/$SAMPLE.keep.bam

# After WASP filtering, GATK ASEReadCounter is safe
gatk ASEReadCounter \
    -I $OUT \
    -V heterozygous_snps.vcf \
    -R hg38.fa \
    -O $SAMPLE.ase_counts.tsv

Per-Tool Failure Modes

GATK ASEReadCounter without WASP -- Reference bias

Trigger: Running ASEReadCounter directly on a standard ATAC BAM without WASP filtering.

Mechanism: Reference allele over-counts due to alignment bias.

Symptom: Per-SNP reference fraction systematically > 0.5; aggregate plots show ~0.51-0.55 instead of 0.5.

Fix: WASP filter first, ALWAYS. There are no exceptions.

Sample size for ASE per SNP

Trigger: Single-individual ATAC; per-SNP heterozygous coverage typically 10-100 reads.

Mechanism: Per-SNP binomial test has limited power; 10 reads at p=0.5 has 95% CI from 0.18 to 0.82 -- effectively no power for moderate effects.

Fix: Aggregate across many SNPs in the same peak (within-peak ASE); aggregate across replicates at same SNP; combine with cis-caQTL across cohort.

RASQUAL -- LD computation requirement

Trigger: Running RASQUAL without pre-computing LD matrix.

Mechanism: RASQUAL needs a per-feature linkage disequilibrium calculation; missing causes the joint model to fail.

Fix: Pre-compute LD with plink --ld-window-kb 5000 per region; supply via --correlation-bias.

WASP -- Phased vs unphased genotypes

Trigger: Using unphased genotypes for ASE.

Mechanism: ASE requires knowing which allele is on which haplotype to assign reads. Unphased het sites ambiguously assign reads.

Fix: Phase genotypes with SHAPEIT5, BEAGLE 5.4, or whatshap (read-based) before running WASP/ASE counter.

Cohort-level caQTL without allelic info

Trigger: MatrixEQTL on peak counts without ASE.

Mechanism: MatrixEQTL maps cohort-level associations; misses cis-mode that ASE captures within individual.

Symptom: Power to detect caQTL is low (typical N=50-100 cohort gives ~hundreds of caQTLs vs ASE-augmented can give thousands).

Fix: Use RASQUAL (joint total + ASE) when N <= 100; or combine MatrixEQTL with separate ASE per individual.

Read-deep peak coverage required

Trigger: Per-peak coverage < 30 reads at SNP site.

Mechanism: Binomial test power at p=0.5, n=30 yields detectable shifts only at |delta_p| >= 0.2.

Fix: Pool replicates if available; or restrict to peaks with sufficient coverage; or aggregate to per-individual peak-level ASE rather than per-SNP.

Decision Tree by Setting

Cohort caQTL Pipeline

Goal: Build cohort-level chromatin QTLs by combining WASP-corrected per-sample counts with cis-genotype association.

Approach: WASP-correct each BAM, build consensus peakset, count reads in peaks per sample, then test cis-genotype association via MatrixEQTL (cohort) or RASQUAL (joint total + allelic).

# 1. WASP-correct each individual's BAM
for sample in $(cat samples.txt); do
    bash wasp_pipeline.sh $sample.bam $sample.vcf
done

# 2. Build consensus peakset (atac-seq/consensus-peakset)
# 3. Count reads in peaks per sample (featureCounts)
featureCounts -F SAF -a consensus.saf -o counts.tsv -p --countReadPairs *.wasp.bam

# 4. Cohort-level caQTL via MatrixEQTL
# (see R script in examples/)

# 5. Per-individual ASE (RASQUAL alternative)
# Parallel: gatk ASEReadCounter per sample, then merge for QuASAR meta-analysis

Within-Peak ASE Aggregation

Goal: Boost per-SNP ASE power by pooling allele counts across heterozygous SNPs in the same peak.

Approach: Map each het SNP to its containing peak, sum REF and ALT counts per peak, run pooled binomial test against 50:50, apply BH FDR, and threshold on effect size.

import pandas as pd, numpy as np
from scipy import stats

# Per-SNP allele counts at heterozygous sites
ase = pd.read_csv('sample.ase_counts.tsv', sep='\t')
ase = ase.rename(columns={'refCount': 'REF', 'altCount': 'ALT'})
ase['totalCount'] = ase['REF'] + ase['ALT']

# Map each SNP to its containing peak
ase['peak'] = map_snps_to_peaks(ase, 'consensus_peaks.bed')

# Aggregate within peak: pooled binomial test
def peak_ase(group):
    ref = group['REF'].sum()
    total = group['totalCount'].sum()
    if total < 30: return pd.Series({'ref_frac': np.nan, 'p_value': np.nan})
    p = stats.binomtest(ref, total, p=0.5).pvalue
    return pd.Series({'ref_frac': ref / total, 'p_value': p, 'snp_count': len(group)})

peak_ase_df = ase.groupby('peak').apply(peak_ase)
peak_ase_df['adj_p'] = stats.false_discovery_control(peak_ase_df['p_value'].dropna())
sig_ase = peak_ase_df[(peak_ase_df['adj_p'] < 0.05) & (abs(peak_ase_df['ref_frac'] - 0.5) >= 0.2)]

|ref_frac - 0.5| >= 0.2 is a 30:70 effect; smaller imbalances are detectable but biologically minor. Per-peak SNP count >= 2 strengthens the call.

RASQUAL Joint Modeling

RASQUAL uses a non-standard CLI: it reads the VCF from stdin via tabix and uses single-letter flags. The canonical invocation pattern is:

Goal: Combine total accessibility counts and allele-specific counts into one joint cis-caQTL test per peak.

Approach: Pre-build binary count and offset files via rasqualTools, then iterate per-feature, streaming the cis-window VCF through tabix into RASQUAL with feature coordinates and SNP counts.

# 1. Pre-compute genotype offsets and binary count files (rasqualTools R package)
# (see rasqualTools::saveRasqualMatricesAsBinary; produces .bin files for -y, -k)

# 2. Per-feature (peak) RASQUAL call: pipe tabix VCF in via stdin
# Per-line meaning: feature name, chromosome, start, end, n_total_SNPs, n_test_SNPs
while IFS=$'\t' read -r name chr start end n_all n_test; do
    tabix cohort.vcf.gz $chr:$((start-500000))-$((end+500000)) | \
        rasqual -y counts.bin \
                -k offsets.bin \
                -n $N_SAMPLES \
                -j $FEATURE_INDEX \
                -l $n_all -m $n_test \
                -s $start -e $end \
                -f $name \
                > $name.rasqual.txt
done < features.tsv

-y is the binary count file; -k is the binary size-factor / offset file (both produced by rasqualTools::saveRasqualMatricesAsBinary from R); -j is the row index of the feature; -l/-m are total / test SNP counts in the feature window. RASQUAL does NOT use --features, --counts, --vcf flags; those are common in newer caQTL wrappers but not in stock RASQUAL.

For modern usage, the rasqualTools R wrapper (Kumasaka GitHub) handles this orchestration. Verify against rasqual --help because the flag set is unusual.

RASQUAL output includes joint p-values, total-only and ASE-only sub-tests; the joint test typically gains 1.5-3x power vs MatrixEQTL alone.

Reconciliation

Operational rule for high-confidence reporting: A cis-regulatory variant must show (a) WASP-filtered allelic imbalance with adjusted p < 0.05 and effect size >= 0.2, AND (b) cohort-level caQTL p < 1e-5 (or RASQUAL joint p < 1e-5), AND (c) accessibility peak overlap. Validation against MPRA or CRISPRi-FlowFISH increases confidence.

Common Errors

References

• van de Geijn B et al 2015 Nat Methods 12:1061 (WASP; reference allele bias correction)
• Castel SE et al 2015 Genome Biol 16:195 (GATK ASEReadCounter; ASE framework)
• Kumasaka N et al 2016 Nat Genet 48:206 (RASQUAL; joint total + ASE caQTL)
• Harvey CT et al 2015 Bioinformatics 31:1235 (QuASAR)
• Buchkovich ML et al 2015 Genome Biol 16:185 (mapping bias review)
• (Modern caQTL workflow citations: consult current literature; earlier placeholder removed pending verified citation)
• Browning SR & Browning BL 2007 Am J Hum Genet 81:1084 (BEAGLE phasing)
• Patterson M et al 2015 J Comput Biol 22:498 (whatshap; read-based phasing)

Related Skills

• atac-seq/atac-peak-calling - Generate peaks for ASE within-peak aggregation
• atac-seq/consensus-peakset - Cohort consensus peakset
• atac-seq/differential-accessibility - Cohort-level (cis + trans) differential
• atac-seq/deep-learning-atac - Predicted variant effects vs observed allelic imbalance
• atac-seq/enhancer-gene-linking - Map ASE-supported variants to target genes
• variant-calling/vcf-basics - VCF input
• variant-calling/joint-calling - Cohort genotype inputs
• phasing-imputation/haplotype-phasing - Phasing before WASP
• causal-genomics/fine-mapping - Use caQTL for fine-mapping
• population-genetics/association-testing - GWAS context