varTFBridge is an integrative framework that combines transcription factor (TF) footprinting data with genome-wide association analyses to identify causal noncoding variants and elucidate their regulatory mechanisms in TF-mediated gene regulation.

Common-variant genome-wide association studies have identified thousands of noncoding loci associated with human diseases and complex traits; however, interpreting their functional mechanisms remains a major challenge. varTFBridge addresses this by:

• Leveraging high-resolution FOODIE (single-molecule deaminase footprinting) TF footprints
• Integrating both common and rare variant association analyses
• Employing footprint-gene linking models (ABC-FP-Max)
• Utilizing AlphaGenome for variant effect prediction
• Prioritizing causal noncoding variants that rewire gene regulatory networks

Stage 1: Variant Association

• Common variants (MAF ≥ 0.1%): GWAS followed by genome-wide fine-mapping (GWFM) with SBayesRC
• Rare variants (MAF < 0.1%): Footprint-based burden tests with leave-one-variant-out analysis

Stage 2: Functional Dissection

• VAR2TFBS: Predicts how variants affect TF binding affinity using JASPAR-based position weight matrices
• ABC-FP-Max: Links variants to target genes through Activity-by-Contact scoring adapted for TF footprints
• AlphaGenome: Assesses variant effects across multiple epigenomic layers (histone modifications, TF binding, chromatin accessibility)

Using 490,640 UK Biobank whole-genome sequences across 13 erythroid traits:

• K562 FOODIE footprints show ~70-fold heritability enrichment for erythroid traits (comprising <0.5% of the genome)
• Identified 209 common variants and 18 rare variants linked to TF binding sites and target genes, modulating 207 unique target genes
• Successfully recapitulated the causal variant rs112233623, revealing how disruption of GATA1/TAL1 co-binding alters CCND3 regulation to drive variation in red blood cell count

git clone https://github.com/JasonLinjc/varTFBridge.git
cd varTFBridge

Python packages (install via pip install -r requirements.txt):

• pandas, numpy, pyfaidx, kipoiseq, memelite, tqdm

External tools:

• Python 3.8+
• bedtools (must be on PATH)
• FIMO (from MEME Suite)
• AlphaGenome (v0.5.1)
• REGENIE (v3.3)
• SBayesRC / GCTB
• Bismark (v0.24.2)
• Trim Galore (v0.6.10)

Run the full common or rare variant pipeline end-to-end:

# Common variant pipeline (9 steps)
bash scripts/run_comvar_pipeline.sh

# Rare variant pipeline (4 steps)
bash scripts/run_rarevar_pipeline.sh

Both scripts support resuming from a specific step, running a single step, or skipping AlphaGenome scoring:

bash scripts/run_comvar_pipeline.sh --from 4      # resume from step 4
bash scripts/run_comvar_pipeline.sh --only 5      # run step 5 only
bash scripts/run_rarevar_pipeline.sh --skip-alphag # skip AlphaGenome (requires API key)

AlphaGenome scoring is automatically skipped if ALPHAGENOME_API_KEY is not set.

• FOODIE footprints: K562 and GM12878 cell lines available in this repository
• TF binding motifs: JASPAR 2025 core non-redundant vertebrates
• UK Biobank WGS data: 490,640 participants (requires authorized access)

# See https://github.com/sunneyxie-lab/bulk-foodie-pipeline

Genome-wide fine-mapping (GWFM) uses GCTB 2.5.4 (--impute-summary + --gwfm) to perform Bayesian fine-mapping on GWAS summary statistics.

Input pipeline: REGENIE GWAS (hg38) → rsID mapping → liftover to hg19 → .ma format (MHC excluded)

Output files per trait:

Erythroid traits analysed (13):

The .snpRes files from GCTB use hg19 coordinates. Liftover to hg38 and add Chromosome_hg38, Start_hg38, End_hg38 columns:

python scripts/comvar_liftover_snpRes.py \
    --snpres-dir data/GWFM_erythroids/snpRes \
    --chain data/reference/hg19ToHg38.over.chain.gz \
    --out-dir data/GWFM_erythroids/snpRes_hg38

Since all traits share the same variant set (~13M), the liftover is performed once and applied to all files. Supports multiprocessing with --workers N.

Filter snpRes_hg38 common variants by PIP threshold and annotate with LCS credible set information (PEP_cs, CS_id):

python scripts/comvar_filter_credible_set.py \
    --snpres-dir data/GWFM_erythroids/snpRes_hg38 \
    --lcs-dir data/GWFM_erythroids/lcs \
    --out-dir data/GWFM_erythroids/credible_set_snpRes \
    --pip-threshold 0.1

Output: data/GWFM_erythroids/credible_set_snpRes/{trait}_credible_set_hg38.csv with columns including PEP_cs (credible-set-level PEP from LCS) and CS_id (credible set ID).

Identify GWFM common variants that fall within FOODIE TF footprints. The script accepts both credible set CSV and snpRes (tab-delimited) formats, with auto-detection of delimiter and column name mapping. Outputs per-trait CSVs with PIP, PEP, PEP_cs, and CS_id annotations.

Using credible set files (filtered by PIP/PEP thresholds):

python scripts/comvar_overlap_foodie_footprints.py \
    --snp-dir data/GWFM_erythroids/credible_set_snpRes \
    --footprint-dir data/FOODIE_footprints \
    --lcs-dir data/GWFM_erythroids/lcs \
    --out-dir results/comvar_footprint_overlap

Using full snpRes_hg38 files (all ~13M common variants per trait):

python scripts/comvar_overlap_foodie_footprints.py \
    --snp-dir data/GWFM_erythroids/snpRes_hg38 \
    --snp-suffix .snpRes \
    --footprint-dir data/FOODIE_footprints \
    --lcs-dir data/GWFM_erythroids/lcs \
    --out-dir results/comvar_footprint_overlap_snpRes

Output: per-trait CSVs ({out_dir}/{footprint}/{trait}_{footprint}.csv) with columns SNP, Chromosome, Start, End, A1, A2, freq, PIP, PEP, PEP_cs, CS_id, footprint_region, plus combined BED files per footprint.

Trait-agnostic: takes the merged BED from Step 1 (all unique variants across traits) and predicts variant effects on TF binding using FIMO motif scanning against JASPAR PWMs. Classifies TF binding changes as Create, Disrupt, Increase, Decrease, or Unchange.

python scripts/comvar_var2tfbs.py \
    --input-bed results/comvar_footprint_overlap_credible/GWFM_variants_in_K562.merged.hg38.bed \
    --allele-src results/comvar_footprint_overlap_credible/K562.merged.hg38 \
    --ref-genome data/reference/hg38.fa \
    --jaspar-meme data/JASPAR_MEME/JASPAR2024_CORE_vertebrates_non-redundant_pfms_meme.txt \
    --out-dir results/comvar_var2tfbs_results

Output: {out_dir}/{cell}_var2tfbs.csv with ref/alt FIMO hits, TF change classification (Create/Disrupt/Increase/Decrease/Unchange), and FASTA files in {out_dir}/fasta/.

Identifies driver rare variants from footprint-based burden test leave-one-out (LOO) analysis and predicts their effects on TF binding. For each significant footprint (Bonferroni-corrected p < 0.05/N_footprints), the driver variant is the one whose removal causes the largest increase in burden test p-value. Footprints where no variant has more than 30 carriers (MAC > 30) are excluded.

python scripts/rarevar_var2tfbs.py \
    --burden-dir data/burdentest_erythroids \
    --loo-file data/leaveoneout_results/K562.leave_one_out.all_traits.20251120.csv \
    --ref-genome data/reference/hg38.fa \
    --jaspar-meme data/JASPAR_MEME/JASPAR2024_CORE_vertebrates_non-redundant_pfms_meme.txt \
    --out-dir results/rarevar_var2tfbs_results

Output: driver_variants_summary.csv (driver variants per trait-footprint) and K562_rarevar_var2tfbs.csv (TF binding effect predictions).

Adapted from the ABC model to link TF footprints to target genes using chromatin accessibility (ATAC-seq), Hi-C contact frequency, and footprint activity scores. See ABC-FP/README.md for license details.

1. Configure biosamples — edit ABC-FP/config/config_FOODIE_ATAC.tsv with paths to your input files:

2. Configure parameters — edit ABC-FP/config/config.yaml to set results_dir and reference file paths.

3. Run the pipeline:

cd ABC-FP
snakemake -j <cores> --use-conda

Output: {results_dir}/{biosample}/Predictions/EnhancerPredictionsAllPutative.tsv.gz containing ABC scores for all putative enhancer-gene links.

Links variants to target genes by bridging: variant → footprint → enhancer → gene using ABC-FP-Max scores. For each variant, the target gene is the one with the highest ABC-FP score. Also annotates with TF binding changes and cell-type-specific TF RNA expression.

python scripts/link_var2gene.py \
    --var2tfbs results/comvar_var2tfbs_results/K562_var2tfbs.csv \
    --footprint-bed data/FOODIE_footprints/K562.merged.hg38.bed \
    --enhancer-bed data/ABC_FP_results/K562_FOODIE_ATAC/Neighborhoods/EnhancerList.bed \
    --abc-predictions data/ABC_FP_results/K562_FOODIE_ATAC/Predictions/EnhancerPredictionsAllPutative.tsv.gz \
    --tf-expr data/gene_expr/K562_ENCFF485RIA_gene.tsv \
    --cell K562 \
    --out-dir results/var2gene_results

Output: {cell}_{prefix}_ABC-FP-Full.csv (variant-TF-gene table with rsID, TF, TF_change, TF expression, TargetGene, ABC.Score.FP, distance) and {cell}_{prefix}_ABC-FP-Max.csv (one row per variant with top ABC-FP-Max gene).

AlphaGenome (Google DeepMind) predicts variant effects across gene expression, splicing, chromatin features, and contact maps at single-bp resolution from up to 1Mbp DNA sequences. Used as orthogonal validation of VAR2TFBS motif-based predictions.

Setup (requires an API key, free for non-commercial use):

conda create -n alphagenome_env python=3.12 -y
conda activate alphagenome_env
pip install ./alphagenome

Example — score a variant:

from alphagenome.data import genome
from alphagenome.models import dna_client

model = dna_client.create('YOUR_API_KEY')

variant = genome.Variant(
    chromosome='chr6', position=41957259,
    reference_bases='C', alternate_bases='T',
)
interval = genome.Interval(chromosome='chr6', start=41432972, end=42481548)

outputs = model.predict_variant(
    interval=interval, variant=variant,
    ontology_terms=['CL:0000038'],  # erythroid progenitor cell
    requested_outputs=[dna_client.OutputType.CHIP_SEQ],
)

See alphagenome/colabs/ for tutorials: quick_start.ipynb, batch_variant_scoring.ipynb, example_analysis_workflow.ipynb, and visualization_modality_tour.ipynb.

Score all VAR2TFBS variants with AlphaGenome's recommended CenterMaskScorer configurations (11 output types, 19 scorers). Supports both common and rare variant ID formats, with checkpoint/resume for large batches.

# Rare variants (no allele-src needed, ~48s for 19 variants)
python scripts/alphag_score_variants.py \
    --var2tfbs results/rarevar_var2tfbs_results/K562_rarevar_var2tfbs.csv \
    --api-key $ALPHAGENOME_API_KEY \
    --cell K562 \
    --prefix rarevar \
    --out-dir results/alphag_scores

# Common variants (requires allele-src and ref-genome for coordinate resolution)
python scripts/alphag_score_variants.py \
    --var2tfbs results/comvar_var2tfbs_results/K562_var2tfbs.csv \
    --allele-src results/comvar_footprint_overlap_credible/K562.merged.hg38 \
    --ref-genome data/reference/hg38.fa \
    --api-key $ALPHAGENOME_API_KEY \
    --cell K562 \
    --prefix comvar \
    --out-dir results/alphag_scores

Output: {cell}_{prefix}_alphag_scores.tsv — tidy scores with columns: rsID, variant_id, output_type, variant_scorer, track_name, raw_score, quantile_score, transcription_factor, biosample_name, etc.

Merge all pipeline outputs into a single linkage table per variant type, showing the full chain: variant → TF binding change → target gene → trait association. Each row represents one rsID × TF × trait combination.

Common variants — merges snpRes overlap, VAR2TFBS, ABC-FP-Max gene linkage, TF expression, and AlphaGenome scores:

python scripts/merge_comvar_linkage_table.py --project-root .

Output: results/K562_comvar2grn.csv — 2,977 rows, 224 variants, 342 TFs, 208 genes, 13 traits. Key columns: rsID, Chromosome, Position, trait, PIP, PEP, TF, TF_change, TF_K562_rna_tpm, TargetGene, ABC.Score.FP, alphag_H3K27ac_score, alphag_ATAC_score.

Rare variants — merges burden test driver variants, VAR2TFBS, ABC-FP-Max gene linkage, TF expression, and AlphaGenome scores:

python scripts/merge_rarevar_linkage_table.py --project-root .

Output: results/K562_rarevar2grn.csv — 609 rows, 19 variants, 99 TFs, 15 genes, 13 traits. Key columns: rsID, Chromosome, Position, trait, burden_p, driver_loo_p, driver_MAC, TF, TF_change, TF_K562_rna_tpm, TargetGene, ABC.Score.FP, alphag_H3K27ac_score, alphag_ATAC_score.

Annotate common variants with PIP > 0.7 with AlphaGenome TF ChIP DIFF_LOG2_SUM scores. For co-binding TFs (e.g. GATA1::TAL1), each component TF is expanded into a separate row with its own AlphaGenome score and K562 RNA expression.

python scripts/annotate_comvar_alphag_tf_chip.py --project-root .

Output: results/K562_comvar_pip70_alphag_tf_chip.csv — one row per variant × component TF. Key columns: rsID, max_PIP, TF_motif, TF_change, TF_alphag, TF_alphag_K562_rna_tpm, alphag_TF_chip_score, alphag_TF_chip_quantile.

All pipeline outputs are stored under results/:

Variant-footprint overlap results from Step 1. Contains per-cell-type subdirectories with per-trait CSVs and merged BED files.

Common variant TF binding effect predictions from Step 2.

Rare variant driver identification and TF binding effect predictions.

Variant-to-gene linking via ABC-FP-Max scores for both common and rare variants.

AlphaGenome variant effect prediction scores across 11 output types (CHIP_TF, CHIP_HISTONE, DNASE, ATAC, RNA_SEQ, CAGE, PROCAP, CONTACT_MAPS, SPLICE_JUNCTIONS, SPLICE_SITES, SPLICE_SITE_USAGE).

Integrated variant→TF→gene→trait linkage tables merging all pipeline outputs.

This project includes the following third-party components, each retaining its original license (see THIRD_PARTY_NOTICES.md):

This project is licensed under the MIT License - see the LICENSE file for details.

For questions and feedback, please open an issue on GitHub or contact Jiecong Lin ([email protected]) and Yajie Zhao ([email protected]).