Variant effects: scoring REF and ALT alleles¶

A PEAgent model reads accessibility straight off the sequence, so the effect of a single-nucleotide variant is simply the difference between the model's output on the REF and ALT windows. peagent_tool builds a 1,344 bp window centred on each SNV, checks the FASTA base against the declared REF allele, and reports the change per cell type on both the probability and the logit scale.

This is the inference step of 08_variant_effect in the PEAgent analysis pipeline. The variants used below are real fine-mapped maize caQTLs with their published posterior inclusion probabilities and effect sizes, so the model output can be compared against measured effects.

In [1]:

import os
from pathlib import Path
from urllib.request import urlretrieve

import numpy as np
import pandas as pd

import peagent_tool

print("peagent_tool", peagent_tool.__version__)

VARIANTS = Path("maize_finemapped_caqtl.tsv")
if not VARIANTS.exists():
    urlretrieve(
        "https://raw.githubusercontent.com/YAOJ-bioin/peagent_tool/main/docs/assets/tutorials/maize_finemapped_caqtl.tsv",
        VARIANTS,
    )
variants = pd.read_csv(VARIANTS, sep="\t")
variants.head()

import os from pathlib import Path from urllib.request import urlretrieve import numpy as np import pandas as pd import peagent_tool print("peagent_tool", peagent_tool.__version__) VARIANTS = Path("maize_finemapped_caqtl.tsv") if not VARIANTS.exists(): urlretrieve( "https://raw.githubusercontent.com/YAOJ-bioin/peagent_tool/main/docs/assets/tutorials/maize_finemapped_caqtl.tsv", VARIANTS, ) variants = pd.read_csv(VARIANTS, sep="\t") variants.head()

peagent_tool 0.1.0

Out[1]:

	chrom	pos	ref	alt	variant_id	acrID	caqtl_celltype	PIP	betahat
0	chr1	119710459	T	C	chr1_119710459_T_C	chr1_119710338_119710488	guard_mother_cell	0.999933	0.410209
1	chr1	166624667	C	A	chr1_166624667_C_A	chr1_166624290_166624599	protodermal_cell	0.999985	1.061403
2	chr1	224726120	C	T	chr1_224726120_C_T	chr1_224726020_224726315	protodermal_cell	0.999979	0.672113
3	chr10	70635320	G	A	chr10_70635320_G_A	chr10_70634962_70635268	bundle_sheath	0.999968	0.472713
4	chr3	92121285	A	G	chr3_92121285_A_G	chr3_92121048_92121238	cortex	0.999997	-0.402201

1. What you need¶

Besides the model assets, FASTA-backed scoring needs the reference genome the model was trained on (maize B73 NAM-5.0) with a pysam-readable index. For a bgzipped FASTA keep the .fai and .gzi files beside it.

pip install "peagent-tool[variant]"
export PEAGENT_GENOME_FASTA=/path/to/Zm-B73-REFERENCE-NAM-5.0.fa.gz

The variant table needs chrom, pos (1-based), ref, alt and optionally variant_id; extra columns are ignored.

In [2]:

GENOME = Path(os.environ.get("PEAGENT_GENOME_FASTA", "Zm-B73-REFERENCE-NAM-5.0.fa.gz"))
print("genome index present:", GENOME.with_suffix(GENOME.suffix + ".fai").exists())
variants[["chrom", "pos", "ref", "alt", "caqtl_celltype", "PIP", "betahat"]]

GENOME = Path(os.environ.get("PEAGENT_GENOME_FASTA", "Zm-B73-REFERENCE-NAM-5.0.fa.gz")) print("genome index present:", GENOME.with_suffix(GENOME.suffix + ".fai").exists()) variants[["chrom", "pos", "ref", "alt", "caqtl_celltype", "PIP", "betahat"]]

genome index present: True

Out[2]:

	chrom	pos	ref	alt	caqtl_celltype	PIP	betahat
0	chr1	119710459	T	C	guard_mother_cell	0.999933	0.410209
1	chr1	166624667	C	A	protodermal_cell	0.999985	1.061403
2	chr1	224726120	C	T	protodermal_cell	0.999979	0.672113
3	chr10	70635320	G	A	bundle_sheath	0.999968	0.472713
4	chr3	92121285	A	G	cortex	0.999997	-0.402201
5	chr4	31688757	G	C	xylem_parenchyma	0.999994	0.417374
6	chr5	111057261	G	A	protophloem_SE	0.999967	-0.075214
7	chr5	220111821	C	T	ground_meristem	0.999954	-1.672827
8	chr8	176677455	T	C	protophloem_SE	0.964738	0.514712
9	chr9	89709102	C	T	protodermal_cell	0.999989	0.644270
10	chr9	107204366	C	T	phloem_SE_procambial_precursors	1.000000	0.550715
11	chr9	123753797	T	C	cortex	0.999996	-0.578320

2. One variant, all 94 cell types¶

targets="all" scores every tissue_celltype group. The headline variant below is a fine-mapped caQTL with posterior inclusion probability 0.96.

In [3]:

from peagent_tool import score_variants

headline = variants[variants["variant_id"].str.startswith("chr8_176677455")]
scored = score_variants(
    headline,
    species="maize",
    genome_fasta=GENOME,
    targets="all",
)
print(scored.shape)
scored.sort_values("abs_delta_logit", ascending=False)[
    ["variant_id", "target", "n_cells", "ref_probability", "alt_probability", "delta_probability", "delta_logit"]
].head(10)

from peagent_tool import score_variants headline = variants[variants["variant_id"].str.startswith("chr8_176677455")] scored = score_variants( headline, species="maize", genome_fasta=GENOME, targets="all", ) print(scored.shape) scored.sort_values("abs_delta_logit", ascending=False)[ ["variant_id", "target", "n_cells", "ref_probability", "alt_probability", "delta_probability", "delta_logit"] ].head(10)

(94, 24)

Out[3]:

	variant_id	target	n_cells	ref_probability	alt_probability	delta_probability	delta_logit
61	chr8_176677455_T_C	root_protophloem	212	0.134151	0.435573	0.301422	1.723881
47	chr8_176677455_T_C	leaf_L1_leaf_primordia	91	0.007484	0.036439	0.028955	1.650634
22	chr8_176677455_T_C	tassel_bulliform	339	0.000392	0.001741	0.001349	1.524343
46	chr8_176677455_T_C	leaf_guard_cell	413	0.004562	0.020599	0.016037	1.498205
16	chr8_176677455_T_C	tassel_mesophyll	330	0.000377	0.001637	0.001260	1.492018
56	chr8_176677455_T_C	tassel_bundle_sheath	409	0.000432	0.001808	0.001377	1.483473
87	chr8_176677455_T_C	ear_mesophyll	244	0.000471	0.002006	0.001534	1.469643
55	chr8_176677455_T_C	ear_bulliform	238	0.000482	0.002063	0.001581	1.468121
0	chr8_176677455_T_C	leaf_bulliform	1006	0.003160	0.012864	0.009704	1.464974
23	chr8_176677455_T_C	leaf_mesophyll	1292	0.003528	0.015571	0.012043	1.464575

The window is built from the genome, and the tool verifies that the FASTA base at the variant position matches the declared REF allele before scoring — a mismatch raises rather than silently producing a wrong answer.

In [4]:

scored[["variant_id", "chrom", "pos", "ref", "alt", "fasta_ref_base",
        "seq_start_0based", "seq_end_0based_exclusive", "left_pad", "right_pad"]].head(1)

scored[["variant_id", "chrom", "pos", "ref", "alt", "fasta_ref_base", "seq_start_0based", "seq_end_0based_exclusive", "left_pad", "right_pad"]].head(1)

Out[4]:

	variant_id	chrom	pos	ref	alt	fasta_ref_base	seq_start_0based	seq_end_0based_exclusive	left_pad	right_pad
0	chr8_176677455_T_C	chr8	176677455	T	C	T	176676782	176678126	0	0

3. Where the effect comes from¶

A variant matters when it sits on a base the model is already relying on. Compute the contribution map of the same window and look at the variant position.

In [5]:

import pysam

from peagent_tool import ContributionScorer, plot_attribution

row = headline.iloc[0]
with pysam.FastaFile(str(GENOME)) as fasta:
    start = int(row["pos"]) - 1 - 672
    window = fasta.fetch(row["chrom"], start, start + 1344).upper()

best_target = scored.sort_values("abs_delta_logit", ascending=False).iloc[0]["target"]
scorer = ContributionScorer("maize")
contribution = scorer.compute({row["variant_id"]: window}, target=best_target, normalize="none")

table = contribution.to_dataframe()
print("most affected cell type:", best_target)
table.iloc[668:676]

import pysam from peagent_tool import ContributionScorer, plot_attribution row = headline.iloc[0] with pysam.FastaFile(str(GENOME)) as fasta: start = int(row["pos"]) - 1 - 672 window = fasta.fetch(row["chrom"], start, start + 1344).upper() best_target = scored.sort_values("abs_delta_logit", ascending=False).iloc[0]["target"] scorer = ContributionScorer("maize") contribution = scorer.compute({row["variant_id"]: window}, target=best_target, normalize="none") table = contribution.to_dataframe() print("most affected cell type:", best_target) table.iloc[668:676]

most affected cell type: root_protophloem

Out[5]:

	position	base	contribution	contribution_raw
668	669	G	0.134161	0.134161
669	670	G	0.159659	0.159659
670	671	T	-0.077139	-0.077139
671	672	T	-0.005673	-0.005673
672	673	T	-0.462237	-0.462237
673	674	A	-0.123201	-0.123201
674	675	A	-0.104634	-0.104634
675	676	A	-0.072159	-0.072159

In [6]:

_ = plot_attribution(
    contribution,
    plot_start=643,
    plot_end=703,
    title="%s | %s" % (row["variant_id"], best_target),
)

_ = plot_attribution( contribution, plot_start=643, plot_end=703, title="%s | %s" % (row["variant_id"], best_target), )

4. A batch of fine-mapped caQTLs¶

Scoring the whole table at once and keeping, for each variant, the cell type with the largest absolute effect gives a compact per-variant summary.

In [7]:

all_scored = score_variants(variants, species="maize", genome_fasta=GENOME, targets="all")
per_variant = (
    all_scored.sort_values("abs_delta_logit", ascending=False)
    .groupby("variant_id", as_index=False)
    .head(1)
    .merge(variants[["variant_id", "caqtl_celltype", "PIP", "betahat"]], on="variant_id")
    .sort_values("abs_delta_logit", ascending=False)
    .reset_index(drop=True)
)
per_variant[["variant_id", "caqtl_celltype", "PIP", "betahat", "target", "delta_logit", "abs_delta_logit"]]

all_scored = score_variants(variants, species="maize", genome_fasta=GENOME, targets="all") per_variant = ( all_scored.sort_values("abs_delta_logit", ascending=False) .groupby("variant_id", as_index=False) .head(1) .merge(variants[["variant_id", "caqtl_celltype", "PIP", "betahat"]], on="variant_id") .sort_values("abs_delta_logit", ascending=False) .reset_index(drop=True) ) per_variant[["variant_id", "caqtl_celltype", "PIP", "betahat", "target", "delta_logit", "abs_delta_logit"]]

Out[7]:

	variant_id	caqtl_celltype	PIP	betahat	target	delta_logit	abs_delta_logit
0	chr8_176677455_T_C	protophloem_SE	0.964738	0.514712	root_protophloem	1.723841	1.723841
1	chr3_92121285_A_G	cortex	0.999997	-0.402201	root_protophloem	-0.278256	0.278256
2	chr1_166624667_C_A	protodermal_cell	0.999985	1.061403	root_protophloem	-0.228943	0.228943
3	chr10_70635320_G_A	bundle_sheath	0.999968	0.472713	axillaryBud_LRC	0.124820	0.124820
4	chr5_220111821_C_T	ground_meristem	0.999954	-1.672827	leaf_L1_leaf_primordia	-0.124122	0.124122
5	chr9_123753797_T_C	cortex	0.999996	-0.578320	root_protophloem	-0.118155	0.118155
6	chr1_224726120_C_T	protodermal_cell	0.999979	0.672113	axillaryBud_axillary_meristem_L1_layer	-0.056434	0.056434
7	chr9_107204366_C_T	phloem_SE_procambial_precursors	1.000000	0.550715	tassel_floral_primordia	0.052793	0.052793
8	chr5_111057261_G_A	protophloem_SE	0.999967	-0.075214	root_protophloem	-0.052720	0.052720
9	chr9_89709102_C_T	protodermal_cell	0.999989	0.644270	leaf_subsidiary	0.040077	0.040077
10	chr1_119710459_T_C	guard_mother_cell	0.999933	0.410209	root_protophloem	0.035493	0.035493
11	chr4_31688757_G_C	xylem_parenchyma	0.999994	0.417374	tassel_bulliform	0.026592	0.026592

Do the predicted directions agree with the measured caQTL effect sizes? With a dozen variants this is an illustration rather than a benchmark — the published evaluation uses the full fine-mapped set — but the comparison is one line:

In [8]:

agreement = np.sign(per_variant["delta_logit"]) == np.sign(per_variant["betahat"])
pd.Series({
    "variants": len(per_variant),
    "same direction as betahat": int(agreement.sum()),
    "spearman rho": round(float(per_variant[["delta_logit", "betahat"]].corr(method="spearman").iloc[0, 1]), 3),
})

agreement = np.sign(per_variant["delta_logit"]) == np.sign(per_variant["betahat"]) pd.Series({ "variants": len(per_variant), "same direction as betahat": int(agreement.sum()), "spearman rho": round(float(per_variant[["delta_logit", "betahat"]].corr(method="spearman").iloc[0, 1]), 3), })

Out[8]:

variants                     12.000
same direction as betahat    10.000
spearman rho                  0.315
dtype: float64

5. Scoring sequences you built yourself¶

When the alleles do not come from a genome — designed sequences, haplotypes, edited alleles — pass the REF and ALT windows directly. Both must be the same length and are normalized to 1,344 bp like any other input.

In [9]:

from peagent_tool import score_ref_alt_sequences

ref_window = window
alt_window = window[:672] + ("A" if window[672] != "A" else "G") + window[673:]

designed = score_ref_alt_sequences(
    ref_sequences={"designed_swap": ref_window},
    alt_sequences={"designed_swap": alt_window},
    species="maize",
    targets=[best_target, "global"],
)
designed[["variant_id", "target", "ref_probability", "alt_probability", "delta_probability", "delta_logit"]]

from peagent_tool import score_ref_alt_sequences ref_window = window alt_window = window[:672] + ("A" if window[672] != "A" else "G") + window[673:] designed = score_ref_alt_sequences( ref_sequences={"designed_swap": ref_window}, alt_sequences={"designed_swap": alt_window}, species="maize", targets=[best_target, "global"], ) designed[["variant_id", "target", "ref_probability", "alt_probability", "delta_probability", "delta_logit"]]

Out[9]:

	variant_id	target	ref_probability	alt_probability	delta_probability	delta_logit
0	designed_swap	root_protophloem	0.134151	0.137609	0.003458	0.029522
1	designed_swap	global	0.005906	0.006201	0.000295	0.055159

6. Command line¶

peagent-tool score-variants \
  --species maize \
  --variants maize_finemapped_caqtl.tsv \
  --fasta Zm-B73-REFERENCE-NAM-5.0.fa.gz \
  --target all \
  --out variant_effects.tsv

Next¶

• In-silico synergy — do two motifs act more than additively?
• Contribution maps — the per-base view used in section 3.