Visualize DMS scores for different GPC regions

This notebook is the main source of specific DMS score highlights used in the paper.

# Imports
import os
import polyclonal
import polyclonal.alphabets
import itertools
import pandas as pd
import altair as alt
import scipy as sp

# Plotting colors
# re-arranged for plot
tol_muted_adjusted = [
    "#AA4499",
    "#88CCEE",
    "#EE7733",
    "#44AA99",
    "#1f78b4",
    "#CC6677",
    "#117733",
    "#999933",
    "#DDCC77",
    "#CC3311",
    "#882255",
    "#000000",
    "#DDDDDD",
]

# Allow more rows for Altair
_ = alt.data_transformers.disable_max_rows()

# this cell is tagged as `parameters` for papermill parameterization
musculinus_scores = None
callosus_scores = None
dms_scores = None
site_diffs = None

min_times_seen = None
n_selections = None

# Parameters
musculinus_scores = (
    "results/func_effects/averages/293T-TFRC-KO-C-musculinus-TfR1_func_effects.csv"
)
callosus_scores = (
    "results/func_effects/averages/293T-TFRC-KO-C-callosus-TfR1_func_effects.csv"
)
dms_scores = "results/summaries/cell_entry_and_pH_stability.csv"
site_diffs = "results/compare_cell_entry/site_diffs.csv"
min_times_seen = 2
n_selections = 2

# # Uncomment for running interactive
# musculinus_scores = "../results/func_effects/averages/293T-TFRC-KO-C-musculinus-TfR1_func_effects.csv"
# callosus_scores = "../results/func_effects/averages/293T-TFRC-KO-C-callosus-TfR1_func_effects.csv"
# dms_scores = "../results/summaries/cell_entry_and_pH_stability.csv"
# site_diffs = "../results/compare_cell_entry/site_diffs.csv"

# min_times_seen = 2
# n_selections = 2

# Filter functional scores based on min times seen and selection number
musculinus_scores = pd.read_csv(musculinus_scores)
callosus_scores = pd.read_csv(callosus_scores)
all_scores = pd.read_csv(dms_scores)
site_diffs = pd.read_csv(site_diffs)

# Add dummy phenotype column
musculinus_scores["phenotype"] = "functional_effect"
callosus_scores["phenotype"] = "functional_effect"
all_scores["phenotype"] = "functional_effect"

# Plotting settings
alphabet = ['R','K','H','D','E','Q','N','S','T','Y','W','F','A','I','L','M','V','G','P','C']

addtl_tooltip_stats = [
    "n_selections",
    "times_seen",
]

addtl_slider_stats = {
    "times_seen" : min_times_seen,
    "n_selections" : n_selections,
}

Define different GPC regions

# Receptor binding sites
TfR1_sites = [
    # From Abraham et al 2010
    # Motif 1
    111,
    113,
    115,
    117,
    # Motif 2
    114,
    116, # many JUNV isolates have A but Romero has E
    # Motif 3
    219, # MACV 223 -> 219 JUNV
    222, # MACV 226 -> 222 JUNV
    224, # MACV 228 -> 224 JUNV potential difference b/c Abraham et al 2010 has '-'
    # Motif 4
    119,
    122,
    123,
    169,
    171,
    # Motif 5
    97,
    98
]

TfR1_motif_1 = [
    # Motif 1
    111,
    113,
    115,
    117,
]

TfR1_motif_2 = [
    # Motif 2
    114,
    116, # many JUNV isolates have A but Romero has E
]

TfR1_motif_3 = [
    # Motif 3
    219, # MACV 223 -> 219 JUNV
    222, # MACV 226 -> 222 JUNV
    224, # MACV 228 -> 224 JUNV potential difference b/c Abraham et al 2010 has '-'
]

TfR1_motif_4 = [
    # Motif 4
    119,
    122,
    123,
    169,
    171,
]

TfR1_motif_5 = [
    # Motif 5
    97,
    98
]

# N-linked glycan sites 
# Mann et al 2025 & Taylor et al 2025 
glycans = [
    95,96,97, # Taylor et al 2025 but only MACV in Mann et al 2025
    105,106,107,
    166,167,168,
    178,179,180,
    357,358,359,
    365,366,367,
    382,383,384,
    387,388,389,
]

# Transmembrane domain 
# 420 - 442 Mann et al 2025
# 24 - 39 & 425 - 439 Taylor et al 2025 & Taylor et al 2026
# Colin suggested 18-41 and 420 - 442
TM = list(range(18,42)) + list(range(420,443))

# Fusion peptides
# N-terminal fusion peptide
# 252 - 258 Mann et al 2025 
NFP = list(range(252, 259))
# Internal fusion loop
# 267 - 291 Mann et al 2025 
IFL = list(range(267,292))

Define specific regions to highlight

# pH stability sites to highlight
pH_sites_1 = [
    10,
    67,
    366,
    368,
]

pH_sites_2 = [
    410,
    411,
    414,
    415,
    416,
    # 419,
    # 420,
    424,
]



pH_sites_3 = [
    427,
    440,
    441,
    443,
    
]

pH_sites_4 = [
    37,
    230,
    282,
    329,
    406,
    431,
    434,
]

prolines = [
    317,
    321,
    322,
    339,
    352,
    371,
    384,
]

Correlation between cell entry and pH effects

subplots = []
for index,region in enumerate(["SSP", "GP1", "GP2", "all"]):

    # Subset data based on region
    subsetted_data = None
    color = None
    if region == "all":
        subsetted_data = all_scores
        color = "black"
    else:
        subsetted_data = all_scores.query("region == @region")
        color = tol_muted_adjusted[index]

    # Remove wildtype residues
    subsetted_data = subsetted_data.loc[subsetted_data["wildtype"] != subsetted_data["mutant"]]

    # Calculate statistics
    r, p = sp.stats.pearsonr(
        subsetted_data.dropna()["entry in C musculinus TfR1 cells"], 
        subsetted_data.dropna()["pH stability"]
    )
    print(f"r correlation for {region}: {r:.2f}")
        
    curr_subplot = alt.Chart(subsetted_data, title=region).mark_point(
        filled=True, 
        color=color, 
        size=50,
        opacity=0.35,
    ).encode(
        alt.X(
            "pH stability",
            axis=alt.Axis(
                title=["pH stability"], 
                values=[-4,-2,0,2,4],
                domainWidth=1,
                domainColor="black",
                tickColor="black",
            ),
            scale=alt.Scale(domain=[-4.5, 4.5])
        ),
        alt.Y(
            "entry in C musculinus TfR1 cells",
            axis=alt.Axis(
                title=["cell entry"], 
                values=[-6,-4,-2,0],
                domainWidth=1,
                domainColor="black",
                tickColor="black",
            ),
            # scale=alt.Scale(domain=[-6.5,1.1])
        ),
        tooltip=[
            "site",
            "wildtype",
            "mutant",
            "entry in C musculinus TfR1 cells",
            "pH stability",
        ], 
    ).properties(
        width=150,
        height=150,
    )
    
    subplots.append(curr_subplot)

# Create final combined plot
pH_vs_func_effects = alt.hconcat(
    subplots[0],
    subplots[1],
    subplots[2],
    subplots[3],
    spacing=5,
    title="cell entry effects vs pH stability",
).configure_axis(
    grid=False,
    labelFontSize=16,
    titleFontSize=16,
    labelFontWeight="normal",
    titleFontWeight="normal",
).configure_title(
    fontSize=24,
).configure_view(
    stroke=None
)

# # Make output dir if doesn't exist
# if not os.path.exists(html_dir):
#     os.mkdir(html_dir)

# print(f"Saving to {html_output}")
# natural_vs_func_effects.save(html_output)

pH_vs_func_effects

r correlation for SSP: 0.10
r correlation for GP1: 0.15
r correlation for GP2: 0.17
r correlation for all: 0.15

Correlation between C musculinus and C callosus entry

subplots = []
for index,region in enumerate(["SSP", "GP1", "GP2", "all"]):

    # Subset data based on region
    subsetted_data = None
    color = None
    if region == "all":
        subsetted_data = all_scores
        color = "black"
    else:
        subsetted_data = all_scores.query("region == @region")
        color = tol_muted_adjusted[index]

    # Remove wildtype residues
    subsetted_data = subsetted_data.loc[subsetted_data["wildtype"] != subsetted_data["mutant"]]

    # Calculate statistics
    r, p = sp.stats.pearsonr(
        subsetted_data.dropna()["entry in C musculinus TfR1 cells"], 
        subsetted_data.dropna()["entry in C callosus TfR1 cells"]
    )
    print(f"r correlation for {region}: {r:.2f}")
        
    curr_subplot = alt.Chart(subsetted_data, title=region).mark_point(
        filled=True, 
        color=color, 
        size=50,
        opacity=0.35,
    ).encode(
        alt.X(
            "entry in C callosus TfR1 cells",
            axis=alt.Axis(
                title=["C callosus"], 
                values=[-6,-4,-2,0],
                domainWidth=1,
                domainColor="black",
                tickColor="black",
            ),
            scale=alt.Scale(domain=[-6.5, 1.5])
        ),
        alt.Y(
            "entry in C musculinus TfR1 cells",
            axis=alt.Axis(
                title=["C musculinus"], 
                values=[-6,-4,-2,0],
                domainWidth=1,
                domainColor="black",
                tickColor="black",
            ),
            scale=alt.Scale(domain=[-6.5, 1.5])
        ),
        tooltip=[
            "site",
            "wildtype",
            "mutant",
            "entry in C musculinus TfR1 cells",
            "entry in C callosus TfR1 cells",
        ], 
    ).properties(
        width=150,
        height=150,
    )
    
    subplots.append(curr_subplot)

# Create final combined plot
callosus_vs_musculinus = alt.hconcat(
    subplots[0],
    subplots[1],
    subplots[2],
    subplots[3],
    spacing=5,
    title="C callosus vs C musculinus",
).configure_axis(
    grid=False,
    labelFontSize=16,
    titleFontSize=16,
    labelFontWeight="normal",
    titleFontWeight="normal",
).configure_title(
    fontSize=24,
).configure_view(
    stroke=None
)

# # Make output dir if doesn't exist
# if not os.path.exists(html_dir):
#     os.mkdir(html_dir)

# print(f"Saving to {html_output}")
# natural_vs_func_effects.save(html_output)

callosus_vs_musculinus

r correlation for SSP: 0.98
r correlation for GP1: 0.86
r correlation for GP2: 0.97
r correlation for all: 0.92

Divergence between C musculinus and C callosus

color_scale = alt.Scale(
    domain=["SSP", "GP1", "GP2"],
    range=["#AA4499", "#88CCEE", "#EE7733",],
)

base = alt.Chart(site_diffs).encode(
    alt.Y(
        "Jensen-Shannon divergence",
        axis=alt.Axis(
            title=["divergence in effects of mutations", "on entry in C. musculinus vs", "C. callosus TfR1 expressing cells"], 
            values=[0.1,0.2,0.3],
            domainWidth=1,
            domainColor="black",
            tickColor="black",
        ),
        scale=alt.Scale(domain=[-0.005, 0.35])
    ),
    alt.X(
        "site",
        axis=alt.Axis(
            title=["site"], 
            values=[100,200,300,400],
            domainWidth=1,
            domainColor="black",
            tickColor="black",
        ),
        scale=alt.Scale(domain=[0, 485])
    ),
    tooltip=[
        "site",
        "Jensen-Shannon divergence",
    ],
)

line = base.mark_line(strokeWidth=1, opacity=1).encode(
    color=alt.Color("region:N", scale=color_scale, legend=alt.Legend(title="Region"))
)

points = base.mark_point(filled=True).encode(
    color=alt.Color("region:N", scale=color_scale, legend=None)
)

JS_plot = (line + points).properties(
    width=500,
    height=200,
).configure_axis(
    grid=False,
    labelFontSize=16,
    titleFontSize=16,
    labelFontWeight="normal",
    titleFontWeight="normal",
)

JS_plot

pH site summary

site_pH_mean = (
    all_scores.loc[all_scores["entry in C musculinus TfR1 cells"] > -5]
    .groupby(["site", "wildtype"])
    .aggregate({
        "pH stability" : "mean",
        "region" : "first",
    })
    .reset_index()
)

color_scale = alt.Scale(
    domain=["SSP", "GP1", "GP2"],
    range=["#AA4499", "#88CCEE", "#EE7733",],
)

base = alt.Chart(site_pH_mean).encode(
    alt.Y(
        "pH stability",
        axis=alt.Axis(
            title=["site mean effect", "on pH stability"], 
            values=[-1,-0.5,0,0.5],
            domainWidth=1,
            domainColor="black",
            tickColor="black",
        ),
        scale=alt.Scale(domain=[-1.25, 0.55])
    ),
    alt.X(
        "site",
        axis=alt.Axis(
            title=["site"], 
            values=[100,200,300,400],
            domainWidth=1,
            domainColor="black",
            tickColor="black",
        ),
        scale=alt.Scale(domain=[0, 485])
    ),
    tooltip=[
        "site",
        "pH stability",
    ],
)

line = base.mark_line(strokeWidth=1, opacity=1).encode(
    color=alt.Color("region:N", scale=color_scale, legend=alt.Legend(title="Region"))
)

points = base.mark_point(filled=True).encode(
    color=alt.Color("region:N", scale=color_scale, legend=None)
)

pH_plot = (line + points).properties(
    width=800,
    height=150,
).configure_axis(
    grid=False,
    labelFontSize=16,
    titleFontSize=16,
    labelFontWeight="normal",
    titleFontWeight="normal",
)

pH_plot

# Make copy of df and floor pH stability
all_scores_copy = all_scores.copy()
all_scores_copy["pH stability"] = all_scores_copy["pH stability"].clip(lower=0)

site_positive_pH_mean = (
    all_scores_copy.loc[all_scores_copy["entry in C musculinus TfR1 cells"] > -5]
    .groupby(["site", "wildtype"])
    .aggregate({
        "pH stability" : "mean",
        "region" : "first",
    })
    .reset_index()
)

color_scale = alt.Scale(
    domain=["SSP", "GP1", "GP2"],
    range=["#AA4499", "#88CCEE", "#EE7733",],
)

base = alt.Chart(site_positive_pH_mean).encode(
    alt.Y(
        "pH stability",
        axis=alt.Axis(
            title=["site mean effect on pH stability"], 
            values=[0,0.25,0.5],
            domainWidth=1,
            domainColor="black",
            tickColor="black",
        ),
        scale=alt.Scale(domain=[-0.005, 0.55])
    ),
    alt.X(
        "site",
        axis=alt.Axis(
            title=["site"], 
            values=[100,200,300,400],
            domainWidth=1,
            domainColor="black",
            tickColor="black",
        ),
        scale=alt.Scale(domain=[0, 485])
    ),
    tooltip=[
        "site",
        "pH stability",
    ],
)

line = base.mark_line(strokeWidth=1, opacity=1).encode(
    color=alt.Color("region:N", scale=color_scale, legend=alt.Legend(title="Region"))
)

points = base.mark_point(filled=True).encode(
    color=alt.Color("region:N", scale=color_scale, legend=None)
)

floored_pH_plot = (line + points).properties(
    width=500,
    height=200,
).configure_axis(
    grid=False,
    labelFontSize=16,
    titleFontSize=16,
    labelFontWeight="normal",
    titleFontWeight="normal",
)

floored_pH_plot

pH stability highlights

pH_sites_1_chart = polyclonal.plot.lineplot_and_heatmap(
    data_df=all_scores,
    stat_col="pH stability",
    category_col="phenotype",
    alphabet=alphabet,
    # addtl_tooltip_stats=addtl_tooltip_stats,
    addtl_slider_stats={"entry in C musculinus TfR1 cells" : -5},
    init_floor_at_zero=False,
    init_site_statistic="mean",
    show_zoombar=False,
    show_lineplot=False,
    sites=pH_sites_1,
    plot_title="Sites 10, 67, 366, & 368",
    heatmap_max_fixed=3,
    heatmap_min_fixed=-3,
    heatmap_color_scheme="pinkyellowgreen",
)

pH_sites_1_chart

pH_sites_2_chart = polyclonal.plot.lineplot_and_heatmap(
    data_df=all_scores,
    stat_col="pH stability",
    category_col="phenotype",
    alphabet=alphabet,
    # addtl_tooltip_stats=addtl_tooltip_stats,
    addtl_slider_stats={"entry in C musculinus TfR1 cells" : -5},
    init_floor_at_zero=False,
    init_site_statistic="mean",
    show_zoombar=False,
    show_lineplot=False,
    sites=pH_sites_2,
    plot_title="pH core",
    heatmap_max_fixed=3,
    heatmap_min_fixed=-3,
    heatmap_color_scheme="pinkyellowgreen",
)

pH_sites_2_chart

pH_sites_3_chart = polyclonal.plot.lineplot_and_heatmap(
    data_df=all_scores,
    stat_col="pH stability",
    category_col="phenotype",
    alphabet=alphabet,
    # addtl_tooltip_stats=addtl_tooltip_stats,
    addtl_slider_stats={"entry in C musculinus TfR1 cells" : -5},
    init_floor_at_zero=False,
    init_site_statistic="mean",
    show_zoombar=False,
    show_lineplot=False,
    sites=pH_sites_3,
    plot_title="427, 440, 441, & 443",
    heatmap_max_fixed=3,
    heatmap_min_fixed=-3,
    heatmap_color_scheme="pinkyellowgreen",
)

pH_sites_3_chart

pH_sites_4_chart = polyclonal.plot.lineplot_and_heatmap(
    data_df=all_scores,
    stat_col="pH stability",
    category_col="phenotype",
    alphabet=alphabet,
    # addtl_tooltip_stats=addtl_tooltip_stats,
    addtl_slider_stats={"entry in C musculinus TfR1 cells" : -5},
    init_floor_at_zero=True,
    init_site_statistic="mean",
    show_zoombar=False,
    show_lineplot=False,
    sites=pH_sites_4,
    plot_title="pH stabilizing candidates",
    heatmap_max_fixed=3,
    heatmap_min_fixed=-3,
    heatmap_color_scheme="pinkyellowgreen",
)

pH_sites_4_chart

TfR1 binding sites

musculinus_TfR1_chart = polyclonal.plot.lineplot_and_heatmap(
    data_df=musculinus_scores,
    stat_col="effect",
    category_col="phenotype",
    alphabet=alphabet,
    addtl_tooltip_stats=addtl_tooltip_stats,
    addtl_slider_stats=addtl_slider_stats,
    init_floor_at_zero=False,
    init_site_statistic="mean",
    show_zoombar=False,
    show_lineplot=False,
    sites=TfR1_sites,
    plot_title="C musculinus TfR1 binding sites",
    heatmap_max_at_least=2,
    heatmap_min_at_least=-1,
)

musculinus_TfR1_chart

individual_motif_charts = []
for motif in [TfR1_motif_1,TfR1_motif_2,TfR1_motif_3,TfR1_motif_4,TfR1_motif_5]:
    chart = polyclonal.plot.lineplot_and_heatmap(
        data_df=musculinus_scores,
        stat_col="effect",
        category_col="phenotype",
        alphabet=alphabet,
        addtl_tooltip_stats=addtl_tooltip_stats,
        addtl_slider_stats=addtl_slider_stats,
        init_floor_at_zero=False,
        init_site_statistic="mean",
        show_zoombar=False,
        show_lineplot=False,
        sites=motif,
        plot_title="C musculinus TfR1 binding sites",
        heatmap_max_at_least=2,
        heatmap_min_at_least=-1,
    )

    individual_motif_charts.append(chart)

individual_motif_charts[0]

individual_motif_charts[1]

individual_motif_charts[2]

individual_motif_charts[3]

individual_motif_charts[4]

Plot functional effects for N-linked glycans

glycan_musculinus = polyclonal.plot.lineplot_and_heatmap(
    data_df=musculinus_scores,
    stat_col="effect",
    category_col="phenotype",
    alphabet=alphabet,
    addtl_tooltip_stats=addtl_tooltip_stats,
    addtl_slider_stats=addtl_slider_stats,
    init_floor_at_zero=False,
    init_site_statistic="mean",
    show_zoombar=False,
    show_lineplot=False,
    sites=glycans,
    plot_title="C musculinus N-linked glycans",
    heatmap_max_at_least=2,
    heatmap_min_at_least=-1,
)

glycan_musculinus

Plot functional effects for fusion peptides

NFP_musculinus = polyclonal.plot.lineplot_and_heatmap(
    data_df=musculinus_scores,
    stat_col="effect",
    category_col="phenotype",
    alphabet=alphabet,
    addtl_tooltip_stats=addtl_tooltip_stats,
    addtl_slider_stats=addtl_slider_stats,
    init_floor_at_zero=False,
    init_site_statistic="mean",
    show_zoombar=False,
    show_lineplot=False,
    sites=NFP,
    plot_title="C musculinus N-terminal fusion peptide",
    heatmap_max_at_least=2,
    heatmap_min_at_least=-1,
)

NFP_musculinus

IFL_musculinus = polyclonal.plot.lineplot_and_heatmap(
    data_df=musculinus_scores,
    stat_col="effect",
    category_col="phenotype",
    alphabet=alphabet,
    addtl_tooltip_stats=addtl_tooltip_stats,
    addtl_slider_stats=addtl_slider_stats,
    init_floor_at_zero=False,
    init_site_statistic="mean",
    show_zoombar=False,
    show_lineplot=False,
    sites=IFL,
    plot_title="C musculinus internal fusion loop",
    heatmap_max_at_least=2,
    heatmap_min_at_least=-1,
)

IFL_musculinus

Candidate prefusion stabilizing mutations

prefusion = polyclonal.plot.lineplot_and_heatmap(
    data_df=musculinus_scores,
    stat_col="effect",
    category_col="phenotype",
    alphabet=alphabet,
    addtl_tooltip_stats=addtl_tooltip_stats,
    addtl_slider_stats=addtl_slider_stats,
    init_floor_at_zero=False,
    init_site_statistic="mean",
    show_zoombar=False,
    show_lineplot=False,
    sites=prolines,
    plot_title="prefusion stabilizing candidates",
    heatmap_max_at_least=2,
    heatmap_min_at_least=-1,
)

prefusion 

Distributions of functional scores for different regions

Plot the distribution of averaged functional scores per site (except stop codons) stratified by region of GPC to show more mutationally tolerant and intolerant regions.

# Load data and remove wildtype
functional_scores = pd.read_csv(dms_scores)
functional_scores = (
    functional_scores.loc[functional_scores["wildtype"] != functional_scores["mutant"]]
)

# Label regions as SSP, GP1, or GP2
functional_scores["region"] = (
    functional_scores.apply(
        lambda x: "SSP" if x["site"] <= 58 else ("GP1" if x["site"] <= 251 else "GP2"), axis=1
    )
)


# Add GP2 transmembrane domain 
# 420 - 442 Mann et al 2025
# 425 - 439 Taylor et al 2025
TM_sites = (
    functional_scores.loc[functional_scores["site"].isin(TM)].copy()
)
TM_sites["region"] = "TM"
functional_scores = (
    pd.concat([
        functional_scores, 
        TM_sites,
    ], ignore_index = True)
)

# Add NFP and IFL
NFP_sites = (
    functional_scores.loc[functional_scores["site"].isin(NFP)].copy()
)
NFP_sites["region"] = "NFP"
functional_scores = (
    pd.concat([
        functional_scores, 
        NFP_sites,
    ], ignore_index = True)
)
IFL_sites = (
    functional_scores.loc[functional_scores["site"].isin(IFL)].copy()
)
IFL_sites["region"] = "IFL"
functional_scores = (
    pd.concat([
        functional_scores, 
        IFL_sites,
    ], ignore_index = True)
)

# Add TfR1 binding sites
TfR1_binding_sites = (
    functional_scores.loc[functional_scores["site"].isin(TfR1_sites)].copy()
)
TfR1_binding_sites["region"] = "TfR1 binding sites"
functional_scores = (
    pd.concat([
        functional_scores, 
        TfR1_binding_sites,
    ], ignore_index = True)
)

# Add glycans
glycan_sites = (
    functional_scores.loc[functional_scores["site"].isin(glycans)].copy()
)
glycan_sites["region"] = "N-glycosylation sites"
functional_scores = (
    pd.concat([
        functional_scores, 
        glycan_sites,
    ], ignore_index = True)
)

def plot_GPC_region_effects(stat_col, df):

    # Plot score distrbutions for each region
    distribution_plot = alt.Chart(
            df, title=f"{stat_col}"
        ).mark_circle(opacity=0.15, size=75).encode(
        y=alt.Y(
            "region:N",
            title="GPC region",
            sort=None,
            axis=alt.Axis(
                domainWidth=1,
                domainColor="black",
                tickColor="black",
            ),
        ),
        x=alt.X(
            stat_col+":Q",
            title=stat_col,
            axis=alt.Axis(
                values=[-5,-4,-3,-2,-1,0,1],
                domainWidth=1,
                domainColor="black",
                tickColor="black",
            ),
        ),
        yOffset="jitter:Q",
        color=alt.Color(
            "region:N",
            scale=alt.Scale(
                domain=functional_scores["region"].unique().tolist(), 
                range=tol_muted_adjusted
            ),
        ).legend(None),
        tooltip=[
            "site",
            "wildtype",
            "mutant",
            alt.Tooltip(
                stat_col, format=".2f", title="effect on cell entry"
            ),
        ],
    ).transform_calculate(
        # Generate Gaussian jitter with a Box-Muller transform
        jitter="sqrt(-2*log(random()))*cos(2*PI*random())"
    )
    
    # Plot median line
    median_plot = alt.Chart(df).mark_tick(size=45, thickness=3, color="#000000").encode(
        y=alt.Y(
            "region:N",
            title="GPC region",
            sort=None,
            axis=alt.Axis(
                domainWidth=1,
            ),
        ),
        x=alt.X(
            "median("+stat_col+"):Q",
            axis=alt.Axis(
                values=[-5,-4,-3,-2,-1,0,1],
                domainWidth=1,
            ),
        ),
        tooltip=[
            "region:N",
            alt.Tooltip(
                "median("+stat_col+"):Q", format=".2f", title="median effect on cell entry",
            ),
        ],
    )
    
    # Combine striplot and median
    combined_plot = (
        (distribution_plot + median_plot)
        .configure_axis(
            grid=False,
            labelFontSize=16,
            titleFontSize=16,
            labelFontWeight="normal",
            titleFontWeight="normal",
        )
        .properties(
            width=400, 
            height=400,
        )
        .configure_title(
            fontSize=24,
        ).configure_view(
            stroke=None
        )
    )
    
    # # Make output dir if doesn't exist
    # if not os.path.exists(html_dir):
    #     os.mkdir(html_dir)
    
    # print(f"Saving to {html_output}")
    # combined_plot.save(html_output)
    
    return combined_plot

effects = [
    "entry in C musculinus TfR1 cells",
    "entry in C callosus TfR1 cells",
]
GPC_region_plots = []

for effect in effects:
    GPC_region_plots.append(plot_GPC_region_effects(effect, functional_scores))

GPC_region_plots[0]

GPC_region_plots[1]