Compare DMS data to natural sequences¶

In [1]:
# Imports
import os
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
from Bio import SeqIO, AlignIO 


# Plotting colors
# re-arranged for plot
tol_muted_adjusted = [
    "#AA4499",
    "#88CCEE",
    "#EE7733",
    "#44AA99",
    "#1f78b4",
    "#CC6677",
    "#117733",
    "#999933",
    "#DDCC77",
    "#CC3311",
    "#882255",
    "#000000",
    "#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)
In [2]:
# this cell is tagged as `parameters` for papermill parameterization
dms_scores = None
alignment = None
metadata = None
out_dir = None
In [3]:
# Parameters
dms_scores = "results/summaries/cell_entry_and_pH_stability.csv"
alignment = "non-pipeline_analyses/Arenavirus_nextstrain/Results/NW_GPC/Alignments/protein_ungapped_no_outgroup.fasta"
metadata = "non-pipeline_analyses/Arenavirus_nextstrain/Results/NW_GPC/metadata.tsv"
out_dir = "results/comparison_to_natural_seqs/"
In [4]:
# # Uncomment for running interactive
# dms_scores = "../results/summaries/cell_entry_and_pH_stability.csv"
# alignment = "../non-pipeline_analyses/Arenavirus_nextstrain/Results/NW_GPC/Alignments/protein_ungapped_no_outgroup.fasta"
# metadata = "../non-pipeline_analyses/Arenavirus_nextstrain/Results/NW_GPC/metadata.tsv"
# out_dir = "../results/comparison_to_natural_seqs/"
In [5]:
def get_natural_sequence_counts(site: int, amino_acid: str, sequences: pd.Series) -> pd.Series:
    """
    Vectorized count of amino acid occurrences at a site across sequences.
    Uses str accessor for fast element-wise comparison instead of Python loop.
    """
    return (sequences.str[site - 1] == amino_acid).sum()

def add_seq_counts(dms_df, seqs_df, query_str, col_name):
    """
    Query sequences and add counts of each mutation
    """

    # Extract subset of sequences sequences
    seqs = seqs_df.query(query_str)["sequence"]

    seq_counts = {
        site: seqs.str[site - 1].value_counts()
        for site in dms_df["site"].unique()         
    }

    dms_df[col_name] = dms_df.apply(
        lambda x: seq_counts[x["site"]].get(x["mutant"], 0), axis=1
    )

    return dms_df.copy()

# Load alignment
natural_seqs_df = pd.DataFrame(
    [{"strain": str(r.id), "sequence": str(r.seq)} for r in AlignIO.read(alignment, "fasta")]
)

# Load metadata and merge with sequences
metadata_df = pd.read_csv(metadata, sep="\t")         
natural_seqs_df = natural_seqs_df.merge(
    metadata_df, how="left", on="strain", validate="one_to_one"
)

# Extract DMS parental sequence
ROMERO_STRAIN = "Romero_P3235_MW227339-1_2017-03-16"
SEQ_LEN = 485
romero_sequence = natural_seqs_df.loc[
    natural_seqs_df["strain"] == ROMERO_STRAIN, "sequence"
].iat[0]

# Create a wildtype df to fill in missing Josiah/wildtype measurements as 0
romero_df = pd.DataFrame({
    "site":             range(1, SEQ_LEN + 1),
    "wildtype":         list(romero_sequence),
    "mutant":           list(romero_sequence),
    "entry in C musculinus TfR1 cells":     0,
    "entry in C callosus TfR1 cells":       0,
    "pH stability":                         0,
    "Romero reference":                  True,
})

# Load DMS scores
scores = pd.read_csv(dms_scores)

# Combine DMS data and reference data
AA_level_df = (
    pd.concat([scores, romero_df], ignore_index=True)
    .sort_values("site")
    .fillna(False)
    .reset_index(drop=True)
)

# Drop gaps and wildtype counts
AA_level_df = (
    AA_level_df
    .query("wildtype != mutant and mutant != '-'")
    .reset_index(drop=True)
)

queries_dict = {
    "candid_related == 'yes' and virus == 'Junin virus'" : "candid_related_counts",
    "candid_related == 'no' and virus == 'Junin virus'" : "Junin_counts",
    "virus == 'Machupo virus'" : "Machupo_counts",
    "virus == 'Guanarito virus'" : "Guanarito_counts",
    "virus == 'Chapare virus'" : "Chapare_counts",
    "virus == 'Sabia virus'" : "Sabia_counts",
}

for query, name in queries_dict.items():
    AA_level_df = add_seq_counts(
        AA_level_df, 
        natural_seqs_df, 
        query, 
        name,
    )

# Drop unseen mutants now, before the bin expansion, to keep frame small.
AA_level_df = (
    AA_level_df
    .query("candid_related_counts > 0 or Junin_counts > 0 or Machupo_counts > 0 or Guanarito_counts > 0 or Chapare_counts > 0 or Sabia_counts > 0")
    .reset_index(drop=True)
)

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

AA_level_df.to_csv(out_dir + "all_arenavirus_GPC_mutations.csv", index=False)

AA_level_df.head(25)

/loc/scratch/10204897/ipykernel_496/780397715.py:63: FutureWarning: Downcasting object dtype arrays on .fillna, .ffill, .bfill is deprecated and will change in a future version. Call result.infer_objects(copy=False) instead. To opt-in to the future behavior, set `pd.set_option('future.no_silent_downcasting', True)`
  .fillna(False)
Out[5]:
site wildtype mutant entry in C callosus TfR1 cells entry in C musculinus TfR1 cells pH stability sequential_site region Romero reference candid_related_counts Junin_counts Machupo_counts Guanarito_counts Chapare_counts Sabia_counts
0 4 F L 0.121 0.1559 -0.004533 4.0 SSP False 0 1 56 16 14 3
1 5 I V -0.7275 -0.5309 -0.01279 5.0 SSP False 0 0 1 0 14 0
2 5 I L -0.121 -0.133 -0.00307 5.0 SSP False 0 0 0 0 0 1
3 5 I F -0.9645 -0.8474 -0.04738 5.0 SSP False 0 0 0 0 0 2
4 8 M F -0.08289 0.03622 0.01414 8.0 SSP False 0 0 56 16 14 3
5 9 Q E -0.4582 -0.2721 -0.002436 9.0 SSP False 0 0 0 0 0 3
6 10 E D -0.4249 -0.3477 -0.01656 10.0 SSP False 0 0 0 16 0 0
7 11 I V -1.107 -0.6762 -0.05926 11.0 SSP False 0 0 0 0 0 3
8 13 T I -0.1068 -0.1987 -0.03291 13.0 SSP False 0 0 0 16 0 0
9 13 T N -0.2939 -0.3443 -0.1045 13.0 SSP False 0 0 0 0 14 3
10 13 T V -0.6867 -0.6824 -0.01471 13.0 SSP False 0 0 56 0 0 0
11 14 F I -0.2889 -0.0889 -0.1203 14.0 SSP False 0 0 0 0 14 3
12 15 L I 0.345 -0.0917 -0.01441 15.0 SSP False 0 0 0 0 14 3
13 15 L F -0.5096 -0.3799 -0.0331 15.0 SSP False 0 0 0 16 0 0
14 16 Q E 0.6121 0.3272 -0.0076 16.0 SSP False 0 0 0 16 0 0
15 16 Q H -1.836 -1.804 0.0296 16.0 SSP False 0 0 0 0 0 3
16 19 L I -0.1448 -0.0845 -0.02416 19.0 SSP False 0 0 0 0 14 3
17 21 I V -0.387 -0.2729 -0.01739 21.0 SSP False 0 0 0 16 0 0
18 24 V I 0.2928 0.462 0.02468 24.0 SSP False 0 0 0 0 14 3
19 24 V A -0.4142 -0.1876 0.001948 24.0 SSP False 0 0 0 16 0 0
20 25 A V -0.8714 -1.08 0.06278 25.0 SSP False 0 0 0 16 0 0
21 27 S T -0.2138 -0.2838 -0.07993 27.0 SSP False 0 0 0 16 0 0
22 29 I L -0.4522 -0.2621 0.01315 29.0 SSP False 0 0 0 16 0 0
23 31 I V -0.07335 -0.1088 -0.006591 31.0 SSP False 0 0 61 1 0 0
24 31 I M -0.5194 -0.3058 0.01158 31.0 SSP False 0 0 0 1 0 0
In [6]:
def plot_distributions(fig, axes, metrics, seq_status, out_dir, fig_name, color, region=None, ylim=None):

    # Plot histplot for each metric
    for index_1, seq in enumerate(seq_status):
        for index_2, metric in enumerate(metrics):
            # Remove NaN/inf values for cleaner plotting
            data_to_plot = None
            if region == None:
                data_to_plot = AA_level_df.query(f"{seq} > 0")[metric]
            else:
                data_to_plot = AA_level_df.query(f"({seq} > 0) and (site >= {region[0]}) and (site <= {region[1]})")[metric]
    
            binrange = None
            if index_2 == 1:
                binrange = (-1,0.25)
            else:
                binrange = (-6,1)# (-3,1)
            
            sns.histplot(
                data=data_to_plot,
                ax=axes[index_2][index_1],
                bins=30,
                binrange=binrange,
                stat="count",
                color=color,
                edgecolor="black",
                alpha=0.5,
                #shrink=0.8,
            )

            if ylim != None:
                axes[index_2][index_1].set_ylim(ylim[index_2][0], ylim[index_2][1])
            
            axes[index_2][index_1].set_title(f"{seq}", fontsize=8)
            axes[index_2][index_1].set_xlabel(f"{metric}", fontsize=8)
            axes[index_2][index_1].set_ylabel("Number of mutations with effect", fontsize=8)
            axes[index_2][index_1].grid(False)
    
    fig.tight_layout()
    
    # Make output dir if doesn't exist
    if not os.path.exists(out_dir):
        os.mkdir(out_dir)
    
    # Save figure
    plt.savefig(out_dir + fig_name)
    plt.show()

Plot full length GPC data

In [7]:
# Define the metrics to facet by
metrics = [
    "entry in C musculinus TfR1 cells",
    # "entry in C callosus TfR1 cells",
    "pH stability",
]
seq_status = [
    "Junin_counts",
    "candid_related_counts",
    # "Machupo_counts",
    # "Guanarito_counts",
    # "Chapare_counts",
    # "Sabia_counts",
    # "other_NW_counts",
]

# Create figure
fig, axes = plt.subplots(2, 2, figsize=(6, 6))
plot_distributions(fig, axes, metrics, seq_status, out_dir, "JUNV_Candid_full_length.svg", color="#000000", region=None)
No description has been provided for this image

NW arenavirus full length

In [8]:
# Define the metrics to facet by
metrics = [
    "entry in C musculinus TfR1 cells",
    # "entry in C callosus TfR1 cells",
    "pH stability",
]
seq_status = [
    # "Junin_counts",
    # "candid_related_counts",
    "Machupo_counts",
    "Guanarito_counts",
    "Chapare_counts",
    "Sabia_counts",
    # "other_NW_counts",
]

# Create figure
fig, axes = plt.subplots(2, 4, figsize=(12, 6), sharey="row")
plot_distributions(fig, axes, metrics, seq_status, out_dir, "NW_GPCs_full_length.svg", color="#000000", region=None, ylim=[[0,65], [0,145]])
No description has been provided for this image

NW SSP, GP1, and GP2

In [9]:
# Create figure
fig, axes = plt.subplots(2, 4, figsize=(12, 6), sharey="row")
plot_distributions(fig, axes, metrics, seq_status, out_dir, "NW_GPCs_SSP.svg", color=tol_muted_adjusted[0], region=[1,58], ylim=[[0,13], [0,21]])
No description has been provided for this image
In [10]:
fig, axes = plt.subplots(2, 4, figsize=(12, 6), sharey="row")
plot_distributions(fig, axes, metrics, seq_status, out_dir, "NW_GPCs_GP1.svg", color=tol_muted_adjusted[1], region=[59,251], ylim=[[0,35], [0,85]])
No description has been provided for this image
In [11]:
fig, axes = plt.subplots(2, 4, figsize=(12, 6), sharey="row")
plot_distributions(fig, axes, metrics, seq_status, out_dir, "NW_GPCs_GP2.svg", color=tol_muted_adjusted[2], region=[252,485], ylim=[[0,21], [0,37]])
No description has been provided for this image

Create combined plot for each region

In [12]:
def plot_distributions_combined(fig, axes, metric, seq_status, out_dir, fig_name, color, region_list, ylim=None):

    # Plot histplot for each metric
    for index_1, seq in enumerate(seq_status):
        for index_2, region in enumerate(region_list):
            # Remove NaN/inf values for cleaner plotting
            data_to_plot = None
            if region == None:
                data_to_plot = AA_level_df.query(f"{seq} > 0")[metric]
            else:
                data_to_plot = AA_level_df.query(f"({seq} > 0) and (site >= {region[0]}) and (site <= {region[1]})")[metric]
    
            binrange = binrange = (-6,1)
            
            sns.histplot(
                data=data_to_plot,
                ax=axes[index_2][index_1],
                bins=30,
                binrange=binrange,
                stat="count",
                color=color[index_2],
                edgecolor="black",
                alpha=0.5,
                #shrink=0.8,
            )

            if ylim != None:
                axes[index_2][index_1].set_ylim(ylim[index_2][0], ylim[index_2][1])
            
            axes[index_2][index_1].set_title(f"{seq}", fontsize=8)
            axes[index_2][index_1].set_xlabel(f"{metric}", fontsize=8)
            axes[index_2][index_1].set_ylabel("Number of mutations with effect", fontsize=8)
            axes[index_2][index_1].grid(False)
    
    fig.tight_layout()
    
    # Make output dir if doesn't exist
    if not os.path.exists(out_dir):
        os.mkdir(out_dir)
    
    # Save figure
    plt.savefig(out_dir + fig_name)
    plt.show()
In [13]:
region_list = [
    [1,58],
    [59,251],
    [252,485],
]

ylim_list = [
    [0,13],
    [0,36],
    [0,21],
]

fig, axes = plt.subplots(3, 4, figsize=(12, 9), sharey="row", sharex="col")
plot_distributions_combined(fig, axes, "entry in C musculinus TfR1 cells", seq_status, out_dir, "NW_GPCs_combined.svg", tol_muted_adjusted[0:3], region_list, ylim_list)
No description has been provided for this image