Loading [MathJax]/extensions/Safe.js

Validation neutralization assays versus polyclonal fits

Compare actual measured neutralization values for specific mutants to the polyclonal fits.

# Imports
import os
import warnings
import neutcurve
import numpy as np
import pandas as pd
import scipy as sp
import seaborn as sns
from matplotlib import pyplot as plt
from matplotlib import ticker as mticker

# Plotting colors
tol_muted_adjusted = [
    "#000000",
    "#CC6677", 
    "#1f78b4", 
    "#DDCC77", 
    "#117733", 
    "#882255", 
    "#88CCEE",
    "#44AA99", 
    "#999933", 
    "#AA4499", 
    "#DDDDDD",
]

# Seaborn style settings
sns.set(rc={
    "figure.dpi":300, 
    "savefig.dpi":300,
    "svg.fonttype":"none",
})
sns.set_style("ticks")
sns.set_palette(tol_muted_adjusted)

# Suppress warnings
warnings.simplefilter("ignore")

File paths for data:

# this cell is tagged as `parameters` for papermill parameterization
fraction_infected_89F_path = None
validation_neuts_89F_path = None
model_predictions_path = None
out_dir = None
neuts_image_path = None
corr_image_path = None

# Parameters
fraction_infected_89F_path = "data/validation_frac_infected_89F.csv"
validation_neuts_89F_path = "data/validation_neuts_89F.csv"
model_predictions_path = "results/antibody_escape/averages/89F_mut_effect.csv"
out_dir = "results/validation_plots/"
neuts_image_path = "results/validation_plots/validation_neut_curves_89F.svg"
corr_image_path = "results/validation_plots/89F_validation_correlation.svg"

# # Uncomment for running interactive
# fraction_infected_89F_path = "../data/validation_frac_infected_89F.csv"
# validation_neuts_89F_path = "../data/validation_neuts_89F.csv"
# model_predictions_path = "../results/antibody_escape/averages/89F_mut_effect.csv"
# out_dir = "../results/validation_plots/"
# neuts_image_path = "../results/validation_plots/validation_neut_curves_89F.svg"
# corr_image_path = "../results/validation_plots/89F_validation_correlation.svg"

Read validation assay measurements and calculate log2 fold change IC50 values

# Read nuetralization data
frac_infected_89F = pd.read_csv(fraction_infected_89F_path, na_filter=None)
validation_neuts_89F = pd.read_csv(validation_neuts_89F_path, na_filter=None)

Create nuetralization curves for 8.9F validation mutations

# Rename viruses from plasmid number to mutation
rename_dict = {
    "3017" : "Unmutated",
    "3893" : "N89D",
    "3894" : "N119S",
    "3895" : "K125L",
    "3896" : "K126L",
    "3897" : "Y129L",
    "3898" : "S135L",
    "3899" : "S138N",
    "3900" : "N148R",
    "3901" : "Q149H",
}
frac_infected_89F["virus"] = frac_infected_89F["virus"].replace(rename_dict)

# Drop VSVG
frac_infected_89F = frac_infected_89F.loc[frac_infected_89F["virus"] != "VSVG"]

# Fit hill curves using neutcurve
fits = neutcurve.curvefits.CurveFits(
    data=frac_infected_89F,
    fixbottom=0,
    fixtop=1,
)

# IC values to calculate
fitParams = fits.fitParams(ics=[50, 80, 90, 95, 97, 98, 99])

# Markers
markers = [
    "o",
    "o",
    "o",
    "o",
    "o",
    "o",
    "o",
    "o",
    "o",
    "o",
    "o",
]

fig, axes = fits.plotSera(
    colors=tol_muted_adjusted,
    markers=markers,
    max_viruses_per_subplot=11,
    xlabel="",
    ylabel="",
    attempt_shared_legend=False,
    despine=True,
)
axes[0,0].set_title(
    "antibody 8.9F", 
    weight="bold",
    fontsize=8,
    # color="#117733",
)
axes[0,0].set_xlabel(
    "concentration (\u03BCg/mL)", 
    fontsize=8,
    # weight="bold",
)
axes[0,0].set_ylabel(
    "fraction infectivity", 
    fontsize=8,
    # weight="bold",
)
axes[0,0].set_ylim(-0.1, 1.3)
axes[0,0].set_yticks([0, 0.5, 1.0])
axes[0,0].set_yticklabels(labels=[0, 0.5, 1.0], fontsize=8)
axes[0,0].set_xlim(0.0005, 12.5)
axes[0,0].set_xticks([0.001, 0.01, 0.1, 1, 10])
axes[0,0].set_xticklabels(labels=["$10^{-3}$", "$10^{-2}$", "$10^{-1}$", "$10^0$", "$10^1$"], fontsize=8)
plt.setp(axes[0,0].collections, alpha=0.8, linewidths=0.5, colors="black") # for vertical error bar segment
plt.setp(axes[0,0].lines, alpha=0.8, markeredgewidth=0.5, markeredgecolor="black", linewidth=1) # for the lines and markers
sns.move_legend(
    axes[0,0], 
    bbox_to_anchor=(1.05, 1), 
    loc="upper left",
    borderaxespad=0,
    frameon=False,
    fontsize=8,
    title="amino acid\nsubstitutions",
    title_fontproperties={"weight" : "bold", "size" : 8},
    alignment="left"
)

# Add edges to legend markers to match scatter plot
for ha in axes[0,0].legend_.legendHandles:
    ha.set_markeredgecolor("black")
    ha.set_markeredgewidth(0.5)
    ha.set_linewidth(0)
    
# Change all spines
for axis in ["top", "bottom", "left", "right"]:
    axes[0,0].spines[axis].set_linewidth(1)
axes[0,0].tick_params(axis="both", length=4, width=1)

width = 2.5
height = 2.5
fig.set_size_inches(width, height)

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

# Save fig
plt.savefig(neuts_image_path)

Now get the predictions by the averaged polyclonal model fits:

# Read model predictions
model_prediction_df = pd.read_csv(model_predictions_path)

# Create new column to match single mutant viruses
model_prediction_df["aa_substitutions"] = (
    model_prediction_df["wildtype"] + 
    model_prediction_df["site"].astype(str) + 
    model_prediction_df["mutant"]
)
model_prediction_df = (
    model_prediction_df.drop(columns=[
        "wildtype",
        "site",
        "mutant",
    ])
)

# Merge model predictions with measured ICs
validation_vs_prediction = (
    validation_neuts_89F.merge(
        model_prediction_df,
        how="left",
        on=["aa_substitutions"],
        validate="one_to_one",
    )
    .rename(columns={"aa_substitutions" : "amino acid substitutions", "ic50_bound" : "lower bound"})
    .fillna(0)
)

# Log10 of IC50 value for accurate calculation of r value
validation_vs_prediction["log10 ic50"] = np.log10(validation_vs_prediction["ic50"])

Calculate the Pearson correlation coefficient r and R^2 between the predicted IC50s from our models and the IC50s measured in validation assays for the antibody 8.9F.

# Calculate correlation between predicted and measured
r, p = sp.stats.pearsonr(
    x=validation_vs_prediction["log10 ic50"], 
    y=validation_vs_prediction["escape_median"]
)
print(f"R={r}")
print(f"R^2={r**2}")

# Plot predicted vs measured
fig, ax = plt.subplots(figsize=(2.25, 2.25))
corr_chart = sns.scatterplot(
    data=validation_vs_prediction,
    x="ic50", 
    y="escape_median",
    hue="amino acid substitutions",
    # style="lower bound",
    # markers=["o", "s"],
    alpha=0.8,
    edgecolor="black",
    ax=ax,
    s=40,
)
corr_chart.set_title(
    "antibody 8.9F",
    fontsize=8,
    weight="bold",
    horizontalalignment="center",
    loc="center",
)
# Change all spines
for axis in ["top", "bottom", "left", "right"]:
    corr_chart.spines[axis].set_linewidth(1)
sns.despine()
corr_chart.tick_params(axis="both", length=4, width=1)
corr_chart.set_xlabel(
    "IC50 (\u03BCg/mL) measured\nby pseuodvirus neutralization", 
    # weight="bold",
    fontsize=8,
)
corr_chart.set_ylabel(
    "escape score measured\nby DMS (arbitrary units)", 
    # weight="bold",
    fontsize=8,
)
corr_chart.set_xscale("log")
corr_chart.set_xlim(0.02, 50)
corr_chart.set_ylim(-1.25, 5.25)
corr_chart.xaxis.set_minor_locator(mticker.NullLocator())  # no minor ticks
corr_chart.set_xticks([0.1, 1, 10])
corr_chart.set_xticklabels(["$10^{-1}$", "$10^{0}$", "$10^{1}$"], size=8)
corr_chart.set_yticks([-1, 0, 1, 2, 3, 4, 5])
corr_chart.set_yticklabels(corr_chart.get_yticks(), size=8)
sns.move_legend(
    corr_chart, 
    "upper left", 
    bbox_to_anchor=(1.2, 1),
    fontsize=8,
    markerscale=1,
    handletextpad=0.1,
    frameon=False,
    borderaxespad=0.75,
    title="amino acid\nsubstitutions",
    title_fontproperties = {
        "size" : 8, 
        "weight" : "bold",
    },
)
corr_chart.text(
    0.03, 
    4, 
    f"r={r:.2f}", 
    horizontalalignment="left",  
    weight="bold",
    fontsize=8,
)

# Label points
for i in range(0, validation_vs_prediction.shape[0]):
    x_pos = validation_vs_prediction.at[i, "ic50"]
    y_pos = validation_vs_prediction.at[i, "escape_median"]
    name = validation_vs_prediction.at[i, "amino acid substitutions"]
    if name == "Q149H":
        corr_chart.text(
            x_pos, 
            y_pos-0.425, 
            f"{name}", 
            horizontalalignment="left",  
            # weight="bold",
            fontsize=6,
        )
    elif name == "N89D":
        corr_chart.text(
            x_pos+0.025, 
            y_pos-0.15, 
            f"{name}", 
            horizontalalignment="left",  
            # weight="bold",
            fontsize=6,
        )
    elif name == "N148R":
        corr_chart.text(
            x_pos+0.0095, 
            y_pos-0.15, 
            f"{name}", 
            horizontalalignment="left",  
            # weight="bold",
            fontsize=6,
        )
    elif name == "N119S":
        corr_chart.text(
            x_pos+8, 
            y_pos-0.2, 
            f"{name}", 
            horizontalalignment="left",  
            # weight="bold",
            fontsize=6,
        )

    elif name == "Y129L":
        corr_chart.text(
            x_pos+8, 
            y_pos+0.1, 
            f"{name}", 
            horizontalalignment="left",  
            # weight="bold",
            fontsize=6,
        )
    elif name == "S138N":
        corr_chart.text(
            x_pos+8, 
            y_pos, 
            f"{name}", 
            horizontalalignment="left",  
            # weight="bold",
            fontsize=6,
        )
    elif name == "unmutated":
        corr_chart.text(
            x_pos, 
            y_pos+0.25, 
            f"{name}       ", 
            horizontalalignment="center",  
            # weight="bold",
            fontsize=6,
        )
    else:
        corr_chart.text(
            x_pos, 
            y_pos+0.25, 
            f"{name}", 
            horizontalalignment="left",  
            # weight="bold",
            fontsize=6,
        )

# Add edges to legend markers to match scatter plot
for ha in corr_chart.legend_.legendHandles:
    ha.set_edgecolor("black")
    ha.set_linewidths(0.5)
        
# Set square ratio
corr_chart.set_box_aspect(1)

# Add vertical line for maximum concentration
corr_chart.axvline(
    x=20,
    color="#000000",
    ls="--",
    lw=0.5,
    alpha=0.5,
)
corr_chart.text(
    22, 
    -0.75, 
    "max\nantibody\nconcentration\ntested", 
    horizontalalignment="left",  
    # weight="bold",
    fontsize=6,
)

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

# Save fig
plt.savefig(corr_image_path)

R=0.904837749168338
R^2=0.8187313523200243