Code
import time
import pandas as pd
import os
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')
import numpy as np
import scipy.stats
from scipy.optimize import curve_fit
%matplotlib inline1 Alternate allele frequency distributions at defining mutations
1.1 Standardized alternate allele frequencies
A naive hypothesis would suggest that the alternate allele frequencies (AFs) measured at mutually exclusive defining mutations in a co-infection sample should directly reflect the variant proportions comprising the sample, i.e. a Delta – Omicron (BA.1) sample with variant ratios 80%-20% (respectively) should have alternate AFs at mutually exclusive Delta-defining mutations of around 0.8 and at mutually exclusive Omicron (BA.1)-defining mutations of around 0.2. However, literary evidence (Bal et al., 2022) proves that this is usually not the case, mostly due to systematic bias introduced by the primers used for sequencing.
Here we investigate the alternate AF distribution in co-infection samples of specific variant compositions. To get comparable results across samples, alternate AFs measured at mutually exclusive variant-defining positions are standardized with the mean and standard deviation of AFs within the given sample for the given variant, so that the mean of standardized values is 0, while their standard deviation is 1. Thus in a sample of variant composition \(A-B\), the \(\overline{AF}_{s,p,A}\) standardized AF in a mutually exclusive defining mutation of variant \(A\), at genomic position \(g\) is:
\[\overline{AF}_{s,g,A} = \frac{AF_{s,g,A} - \mu_{s,A}}{\sigma_{s,A}}\]
where \(\mu_{s,A}\) is the mean and \(\sigma_{s,A}\) is the standard deviation of alternate AFs measured in defining positions of variant \(A\), in sample \(s\).
The figures below show the standardized AFs in samples of the top 4 most frequent variant compositions separately for the mutually exclusive defining mutations of the comprising variants.
Note
Out of the 7,700 previously identified co-infection samples, only those 7,290 are considered here that had exactly two comprising variant strains.
Code
# helper functions for treating variant combinations consistently
def get_sorted_variants(v):
return ", ".join(sorted(v.split(", ")))
def get_num_of_variants(v):
return len(v.split(", "))
# importing all mutually exclusive defining mutations of co-infection samples
df_coinf_allmuts = pd.read_csv("datafile3.csv")
# filtering for only those samples that consist of exactly two parental strains
all_samples_to_check = df_coinf_allmuts[(df_coinf_allmuts["variants"].apply(get_num_of_variants) == 2)]["runid"].unique()
# variant compositions with more than 500 co-infection samples
varcomp_counts = df_coinf_allmuts[df_coinf_allmuts["runid"].isin(all_samples_to_check)].groupby("runid").agg({"variants": lambda x: list(x)[0]}).reset_index()["variants"].apply(get_sorted_variants).value_counts()
varcomp_to_plot = list(varcomp_counts[varcomp_counts>500].keys())
# plotting standardized AFs in defining mutation positions for sample groups of different variant compositions
pos_to_highlight = {"Alpha": [15279], #syn C>T
"Delta": [23604], #C23604G; S:P681R
"Omicron (BA.1)": [11282], #AGTTTGTCTG11282A; ORF1a:L3674-; ORF1a:S3675-; ORF1a:G3676-
"Iota": [9867], # T9867C, ORF1a:L3201P
"Epsilon": [22917] #T>G; S:L452R (increased escape from antibody)
}
for vc in varcomp_to_plot:
variant_list = vc.split(", ")
var1_short = variant_list[0].split("_")[0]
var2_short = variant_list[1].split("_")[0]
if "Omicron" in var1_short:
var1_short = variant_list[0].split("_")[0] + " (" + variant_list[0].split("_")[1] + ")"
if "Omicron" in var2_short:
var2_short = variant_list[1].split("_")[0] + " (" + variant_list[1].split("_")[1] + ")"
figtitle = var1_short + " - " + var2_short + " co-infection samples"
fig = plt.figure(figsize=(15,5))
gs = fig.add_gridspec(2, hspace=0)
axs = gs.subplots(sharex=True, sharey=True)
# fig.suptitle(figtitle, fontsize=16)
pos_list_var1 = []
pos_list_var2 = []
varcomp_samples = list(df_coinf_allmuts[(df_coinf_allmuts["runid"].isin(all_samples_to_check)) &
(df_coinf_allmuts["variants"].apply(get_sorted_variants) == vc)]["runid"].unique())
for s in varcomp_samples:
df_tmp = df_coinf_allmuts[df_coinf_allmuts["runid"] == s]
df_tmp = df_tmp.sort_values(by="pos")
# rescaling AF values
mean_var1 = df_tmp[df_tmp["variant"] == variant_list[0]]["af"].mean()
std_var1 = df_tmp[df_tmp["variant"] == variant_list[0]]["af"].std()
mean_var2 = df_tmp[df_tmp["variant"] == variant_list[1]]["af"].mean()
std_var2 = df_tmp[df_tmp["variant"] == variant_list[1]]["af"].std()
axs[0].plot(df_tmp[df_tmp["variant"] == variant_list[0]]["pos"],
(df_tmp[df_tmp["variant"] == variant_list[0]]["af"]-mean_var1)/std_var1,
color="black", alpha=0.01)
pos_list_var1 += list(df_tmp[df_tmp["variant"] == variant_list[0]]["pos"])
axs[1].plot(df_tmp[df_tmp["variant"] == variant_list[1]]["pos"],
(df_tmp[df_tmp["variant"] == variant_list[1]]["af"]-mean_var2)/std_var2,
color="black", alpha=0.01)
pos_list_var2 += list(df_tmp[df_tmp["variant"] == variant_list[1]]["pos"])
for p in pos_list_var1:
if p in pos_to_highlight[var1_short]:
axs[0].vlines(x = p, ymin=-8, ymax=8, lw=1, color="#c1121f")
else:
axs[0].vlines(x = p, ymin=-8, ymax=8, ls="dashed", lw=1, color="#d6ccc2")
for p in pos_list_var2:
if p in pos_to_highlight[var2_short]:
axs[1].vlines(x = p, ymin=-8, ymax=8, lw=1, color="#c1121f")
else:
axs[1].vlines(x = p, ymin=-8, ymax=8, ls="dashed", lw=1, color="#d6ccc2")
axs[0].annotate(var1_short + " mutually exclusive defining mutations", xy=(500, -7), xycoords='data')
axs[1].annotate(var2_short + " mutually exclusive defining mutations", xy=(500, -7), xycoords='data')
for ax in axs:
ax.set_xlabel("Genomic position", fontsize=14)
ax.set_ylabel("Standardized AF", fontsize=14)
ax.set_xlim(0,30000)
ax.set_ylim(-8,8)
ax.grid(visible=False)
ax.tick_params(size=5, color="#666666")
ax.label_outer()
plt.show()
plt.close() 
1.2 Possible causes of bias in alternate AFs
It is apparent from Fig. 1 that in some genomic positions, standardized AF values tend to markedly differ from zero. We hypothesize that this effect might due to the preferential attachment of primer sequences to genomes containing or lacking specific mutations.
It is also evident that these systematic biases vary based on the specific variant composition under investigation. For example, standardized alternate AFs at the Alpha-defining synonymous 15279:C>T mutation are around zero in Alpha – Iota samples (Fig. 1b), while they can have much lower values in Alpha – Epsilon (Fig. 1c)) and Alpha – Delta (Fig. 1d)) samples. This observation also supports the theory that the primers used during PCR amplification can have a preference for specific genome sequences, which is directly influenced by the set of available genomes (i.e. the variant composition).
To check this assumption, we selected Delta – Omicron (BA.1) co-infection samples for further analysis. We binned the genome into regions of 500 bp and selected those that overlapped at least one Delta- and one Omicron (BA.1)-defining mutation. In these regions, we calculated the mean value of standardized alternate AFs of all mutually exclusive defining mutations that overlapped the given region, separately for Delta- and Omicron (BA.1)-defining mutations and for each sample. These mean values are plotted in Fig. 2a and Fig. 2b, for Delta- and Omicron (BA.1)-defining positions, respectively. In Fig. 2c, these values are displayed simultaneously on the same chart and in Fig. 2d, the sum of these values is shown for each region, for each sample. In regions indicated by the vertical arrows, the mean standardized AFs of Delta- and Omicron (BA.1)-defining mutations are shifted from zero in the opposite direction and their sum shows a narrower distribution around zero. This suggests that in these regions, genomes of one of the two comprising variants are disproportionately amplified, while the relative abundance of the other variant is consequently lowered.
Code
# function for easier standardization of AFs
def rescale_af(r, mean_om, std_om, mean_delta, std_delta):
if r["variant"] == "Delta_B.1.617.2":
return (r["af"]-mean_delta)/std_delta
elif r["variant"] == "Omicron_BA.1":
return (r["af"]-mean_om)/std_om
else:
return np.nan
vc = "Delta_B.1.617.2, Omicron_BA.1"
vc_samples = list(df_coinf_allmuts[df_coinf_allmuts["variants"].apply(get_sorted_variants) == vc]["runid"].unique())
variant_list = vc.split(", ")
pos_list = []
pos_ranges = np.linspace(0,30000,61)
pos_to_plot = [(pos_ranges[i]+pos_ranges[i+1])/2 for i in range(len(pos_ranges)-1)]
xmin_steps = pos_ranges[:-1]
xmax_steps = pos_ranges[1:]
bluecol = "#023e8a"
redcol = "#c1121f"
fig = plt.figure(figsize=(15,8))
gs = fig.add_gridspec(4, hspace=0)
axs = gs.subplots(sharex=True, sharey=True)
for ax in axs:
ax.hlines(y = 0, xmin=0, xmax=30000, linestyle="dashed", lw=1, color="black", alpha=0.1)
for s in vc_samples:
df_tmp = df_coinf_allmuts[df_coinf_allmuts["runid"] == s]
df_tmp = df_tmp.sort_values(by="pos")
# rescaling omicron AF values
mean_om = df_tmp[df_tmp["variant"] == variant_list[1]]["af"].mean()
std_om = df_tmp[df_tmp["variant"] == variant_list[1]]["af"].std()
mean_delta = df_tmp[df_tmp["variant"] == variant_list[0]]["af"].mean()
std_delta = df_tmp[df_tmp["variant"] == variant_list[0]]["af"].std()
df_tmp["stdAF"] = df_tmp.apply(rescale_af, args=(mean_om, std_om, mean_delta, std_delta), axis=1)
sum_to_plot_d = []
sum_to_plot_o = []
sum_to_plot_s = []
for i in range(len(pos_ranges)-1):
df_tmptmp = df_tmp[(df_tmp["pos"] > pos_ranges[i]) & (df_tmp["pos"] <= pos_ranges[i+1])]
if df_tmptmp.shape[0] > 1 and len(df_tmptmp["variant"].unique()) == 2:
sum_to_plot_d.append(df_tmptmp[(df_tmptmp["variant"] == variant_list[0])]["stdAF"].mean())
sum_to_plot_o.append(df_tmptmp[(df_tmptmp["variant"] == variant_list[1])]["stdAF"].mean())
sum_to_plot_s.append(df_tmptmp[(df_tmptmp["variant"] == variant_list[0])]["stdAF"].mean()+df_tmptmp[(df_tmptmp["variant"] == variant_list[1])]["stdAF"].mean())
else:
sum_to_plot_d.append(np.nan)
sum_to_plot_o.append(np.nan)
sum_to_plot_s.append(np.nan)
axs[0].hlines(np.array(sum_to_plot_d)[~np.isnan(np.array(sum_to_plot_d))],
np.array(xmin_steps)[~np.isnan(np.array(sum_to_plot_d))],
np.array(xmax_steps)[~np.isnan(np.array(sum_to_plot_d))],
colors=redcol, alpha = 0.01)
axs[1].hlines(np.array(sum_to_plot_o)[~np.isnan(np.array(sum_to_plot_o))],
np.array(xmin_steps)[~np.isnan(np.array(sum_to_plot_o))],
np.array(xmax_steps)[~np.isnan(np.array(sum_to_plot_o))],
colors=bluecol, alpha = 0.01)
axs[2].hlines(np.array(sum_to_plot_d)[~np.isnan(np.array(sum_to_plot_d))],
np.array(xmin_steps)[~np.isnan(np.array(sum_to_plot_d))],
np.array(xmax_steps)[~np.isnan(np.array(sum_to_plot_d))],
colors=redcol, alpha = 0.01)
axs[2].hlines(np.array(sum_to_plot_o)[~np.isnan(np.array(sum_to_plot_o))],
np.array(xmin_steps)[~np.isnan(np.array(sum_to_plot_o))],
np.array(xmax_steps)[~np.isnan(np.array(sum_to_plot_o))],
colors=bluecol, alpha = 0.01)
axs[3].hlines(np.array(sum_to_plot_s)[~np.isnan(np.array(sum_to_plot_s))],
np.array(xmin_steps)[~np.isnan(np.array(sum_to_plot_s))],
np.array(xmax_steps)[~np.isnan(np.array(sum_to_plot_s))],
colors="black", alpha = 0.01)
axs[0].set_ylabel("Mean\n standardized AF", fontsize=12)
axs[0].annotate("In Delta mutually exclusive defining mutations", xy=(500, -7), xycoords='data')
axs[0].annotate("a)", xy=(500, 5.5), xycoords='data', fontsize = 14)
axs[1].set_ylabel("Mean\n standardized AF", fontsize=12)
axs[1].annotate("b)", xy=(500, 5.5), xycoords='data', fontsize = 14)
axs[1].annotate("In Omicron (BA.1) mutually exclusive defining mutations", xy=(500, -7), xycoords='data')
axs[2].set_ylabel("Mean\n standardized AF", fontsize=12)
axs[2].annotate("c)", xy=(500, 5.5), xycoords='data', fontsize = 14)
axs[2].annotate('In', xy=(500, -7), xycoords='data')
axs[2].annotate('Delta', xy=(900, -7), xycoords='data', color = redcol)
axs[2].annotate('and', xy=(1900, -7), xycoords='data')
axs[2].annotate('Omicron (BA.1)', xy=(2650, -7), xycoords='data', color = bluecol)
axs[2].annotate('mutually exclusive defining mutations', xy=(5500, -7), xycoords='data')
axs[3].set_ylabel("Sum of mean\n standardized AF", fontsize=12)
axs[3].annotate("Sum of mean standardized AFs for the two variants", xy=(500, -7), xycoords='data')
axs[3].annotate("d)", xy=(500, 5.5), xycoords='data', fontsize = 14)
axs[3].annotate("", xy=(24250, -4), xytext=(24250, -7), arrowprops=dict(arrowstyle="->"))
axs[3].annotate("", xy=(26750, -4), xytext=(26750, -7), arrowprops=dict(arrowstyle="->"))
axs[3].annotate("", xy=(28750, -4), xytext=(28750, -7), arrowprops=dict(arrowstyle="->"))
for ax in axs:
ax.set_xlabel("Genomic position (binned by 500 bp)", fontsize=14)
ax.set_xlim(0,30000)
ax.set_ylim(-8,8)
ax.grid(visible=False)
ax.tick_params(size=5, color="#666666")
ax.label_outer()
plt.show()
plt.close() 
1.3 Estimation of standardized AF distributions
Given that our main approach for detecting traces of recombination in co-infection samples is to identify putative breakpoints based on alternate AF shifts along the genome, it is important to distinguish between true signals and ones simply caused by the above detailed inherent variation in AF values. To this end, we estimate the standardized alternate AF distribution at each mutually exclusive defining mutation of each variant for each variant composition separately. In other words, we fit a bimodal probability density distribution to the normalized histogram of the data points at each vertical line in Fig. 1, and do the same for all additionally relevant variant compositions.
The generic form of the function we use for fitting is the following:
def bimodal(x, mu1, s1, mu2, s2):
f = (1/(s1*np.sqrt(2*np.pi)))*np.exp((-1/2)*((x-mu1)/s1)**2)
f += (1/(s2*np.sqrt(2*np.pi)))*np.exp((-1/2)*((x-mu2)/s2)**2)
return f/2The following figures show the results of the fitting procedure for the genomic positions indicated by red vertical lines in Fig. 1.
Code
# bimodal function to fit
def bimodal(x, mu1, s1, mu2, s2):
f = (1/(s1*np.sqrt(2*np.pi)))*np.exp((-1/2)*((x-mu1)/s1)**2)
f += (1/(s2*np.sqrt(2*np.pi)))*np.exp((-1/2)*((x-mu2)/s2)**2)
return f/2
fig = plt.figure(figsize=(12,4))
gs = fig.add_gridspec(2, 4, hspace=0, wspace=0)
axs = gs.subplots(sharex=True, sharey=True)
ax_idx_dict = {'Delta_B.1.617.2, Omicron_BA.1': {'Delta_B.1.617.2': [1,0], 'Omicron_BA.1': [0,0]},
'Alpha_B.1.1.7, Delta_B.1.617.2': {'Delta_B.1.617.2': [1,1], 'Alpha_B.1.1.7': [0,1]},
'Alpha_B.1.1.7, Epsilon_B.1.427_429': {'Epsilon_B.1.427_429': [1,2], 'Alpha_B.1.1.7': [0,2]},
'Alpha_B.1.1.7, Iota_B.1.526': {'Iota_B.1.526': [1,3], 'Alpha_B.1.1.7': [0,3]}}
# min_number_of_datapoints = 100
for vc in varcomp_to_plot:
var1 = vc.split(", ")[0]
var2 = vc.split(", ")[1]
var1_short = var1.split("_")[0]
var2_short = var2.split("_")[0]
if "Omicron" in var1_short:
var1_short = var1.split("_")[0] + " (" + var1.split("_")[1] + ")"
if "Omicron" in var2_short:
var2_short = var2.split("_")[0] + " (" + var2.split("_")[1] + ")"
varcomp_samples = list(df_coinf_allmuts[df_coinf_allmuts["variants"].apply(get_sorted_variants) == vc]["runid"].unique())
AF_list_var1 = []
AF_list_var2 = []
for s in varcomp_samples:
df_tmp = df_coinf_allmuts[df_coinf_allmuts["runid"] == s]
df_tmp = df_tmp.sort_values(by="pos")
df_tmp_var1 = df_tmp[df_tmp["variant"] == var1]
df_tmp_var2 = df_tmp[df_tmp["variant"] == var2]
mean_var1 = df_tmp_var1["af"].mean()
std_var1 = df_tmp_var1["af"].std()
mean_var2 = df_tmp_var2["af"].mean()
std_var2 = df_tmp_var2["af"].std()
p1 = pos_to_highlight[var1_short][0]
if df_tmp_var1[df_tmp_var1["pos"]==p1].shape[0] > 0:
AF_list_var1.append((df_tmp_var1[df_tmp_var1["pos"]==p1].iloc[0]["af"]-mean_var1)/std_var1)
p2 = pos_to_highlight[var2_short][0]
if df_tmp_var2[df_tmp_var2["pos"]==p2].shape[0] > 0:
AF_list_var2.append((df_tmp_var2[df_tmp_var2["pos"]==p2].iloc[0]["af"]-mean_var2)/std_var2)
ax1idx_0, ax1idx_1 = ax_idx_dict[vc][var1]
ax2idx_0, ax2idx_1 = ax_idx_dict[vc][var2]
d = np.histogram(AF_list_var1, bins=100, density=True)
x = 0.5*(d[1][1:] + d[1][:-1])
y = d[0]
popt,pcov=curve_fit(bimodal, x, y, p0=[np.percentile(AF_list_var1, 25), 1, np.percentile(AF_list_var1, 75), 1])
axs[ax1idx_0][ax1idx_1].hist(AF_list_var1, bins=100, density=True, color="#dee2e6");
axs[ax1idx_0][ax1idx_1].plot(x, bimodal(x, *popt), color="#03045e")
d = np.histogram(AF_list_var2, bins=100, density=True)
x = 0.5*(d[1][1:] + d[1][:-1])
y = d[0]
popt,pcov=curve_fit(bimodal, x, y, p0=[np.percentile(AF_list_var2, 25), 1, np.percentile(AF_list_var2, 75), 1])
axs[ax2idx_0][ax2idx_1].hist(AF_list_var2, bins=100, density=True, color="#dee2e6");
axs[ax2idx_0][ax2idx_1].plot(x, bimodal(x, *popt), color="#03045e")
for i in range(2):
for j in range(4):
ax = axs[i][j]
ax.set_xlabel("Standardized AF", fontsize=10)
ax.set_ylabel("Probability density", fontsize=10)
ax.set_xlim(-6,6)
ax.set_ylim(0,1)
ax.grid(visible=False)
ax.tick_params(size=5, color="#666666")
ax.label_outer()
ax.vlines(x=0, ymin=0, ymax=1, color="black", linestyle="dashed", alpha=0.1)
ax.set_yticks([0,0.25,0.5,0.75])
if j != 0:
ax.get_yaxis().set_visible(False)
axs[0][0].annotate('Delta - Omicron (BA.1) samples', xy=(0, 1.1), xycoords='data', annotation_clip=False, ha="center")
axs[0][1].annotate('Alpha - Delta samples', xy=(0, 1.1), xycoords='data', annotation_clip=False, ha="center")
axs[0][2].annotate('Alpha - Epsilon samples', xy=(0, 1.1), xycoords='data', annotation_clip=False, ha="center")
axs[0][3].annotate('Alpha - Iota samples', xy=(0, 1.1), xycoords='data', annotation_clip=False, ha="center")
axs[0][0].annotate('11282: AGTTTGTCTG>A \nORF1a: L3674-; S3675-; G3676-',
xy=(-5.5, 0.78), xycoords='data', ha="left", fontsize=8)
axs[0][1].annotate('15279: C>T \n(synonymous)',
xy=(-5.5, 0.78), xycoords='data', ha="left", fontsize=8)
axs[0][2].annotate('15279: C>T \n(synonymous)',
xy=(-5.5, 0.78), xycoords='data', ha="left", fontsize=8)
axs[0][3].annotate('15279: C>T \n(synonymous)',
xy=(-5.5, 0.78), xycoords='data', ha="left", fontsize=8)
axs[1][0].annotate('23604: C>G \nS: P681R',
xy=(-5.5, 0.78), xycoords='data', ha="left", fontsize=8)
axs[1][1].annotate('23604: C>G \nS: P681R',
xy=(-5.5, 0.78), xycoords='data', ha="left", fontsize=8)
axs[1][2].annotate('22917: T>G \nS: L452R',
xy=(-5.5, 0.78), xycoords='data', ha="left", fontsize=8)
axs[1][3].annotate('9867: T>C \nORF1a: L3201P',
xy=(-5.5, 0.78), xycoords='data', ha="left", fontsize=8)
plt.show()
plt.close() 
2 Identification of putative recombinants
2.1 Hard-filtering of co-infection samples
Clonally recombinant samples would in principle have genomes that were fused together from the appropriate parts of the genomes of parental viral strains at some breakpoint(s), thus would exhibit signs of different sets of mutually exclusive defining mutations of their parental variants before and after the breakpoint(s). In practice, however, the two parental strains usually co-exist with the recombinant strain (with varying ratios) within a single sample. In theory, assuming no bias in AF distributions, putative recombinant breakpoints could be identified from the shifts in AFs observed for the variant-defining mutations of the parents. In an ideal setting, the absolute value of the AF shift corresponds to the ratio of the recombinant genome in the sample (Fig. 4).
To select samples with supposed evidence of the presence of recombinant genomes, we developed a pipeline that detects putative breakpoints in co-infection samples where the mean alternate AF of one set of mutually exclusive defining mutations increases, while the mean alternate AF of the other set of mutually exclusive defining mutations decreases. To filter out presumed artefacts and noise, only those genomic positions were retained as possible breakpoints
- where the absolute AF shift for both variants was 0.05 or larger,
- that separated the genome in a way that resulted in at least four variant-defining mutations of both parental strains on each side of the position,
- where the absolute AF shift for both variants was larger than the standard deviation of AFs on either side of the position,
- which resulted in genomic segments where the standard deviations of AFs in the defining mutations of either variant were lower than 0.1.
Note
In this analysis (as well as previously), all mutually exclusive lineage-defining mutations were considered for the given variant combination to increase statistical power.
Code
# running hard-filtering for all samples
st = time.time()
print("Starting processing...")
i = 1
min_defmut_oneside = 4
min_avg_af_diff = 0.05
max_af_std = 0.1
samples_to_putative_breakpoints = dict()
samples_to_best_breakpoint = dict()
for s in all_samples_to_check:
i += 1
if i/len(all_samples_to_check) >= 0.01:
print(">", end="")
i = 0
df_tmp = df_coinf_allmuts[df_coinf_allmuts["runid"] == s]
pos_all = sorted(list(df_tmp["pos"]))
variant_list = df_tmp.iloc[0]["variants"].split(", ")
putative_breakpoints = []
best_breakpoint = "X"
best_diff = 0
for p in pos_all[min_defmut_oneside+1:-(min_defmut_oneside+1)]:
if df_tmp[df_tmp["pos"]<=p].groupby("variant").count().reset_index().shape[0] < 2:
continue
if df_tmp[df_tmp["pos"]<=p].groupby("variant").count().reset_index()["runid"].min() < min_defmut_oneside:
continue
if df_tmp[df_tmp["pos"]>p].groupby("variant").count().reset_index().shape[0] < 2:
continue
if df_tmp[df_tmp["pos"]>p].groupby("variant").count().reset_index()["runid"].min() < min_defmut_oneside:
continue
df_before = df_tmp[df_tmp["pos"]<=p].groupby("variant").agg({"af": [np.mean, np.std]}).reset_index()
df_after = df_tmp[df_tmp["pos"]>p].groupby("variant").agg({"af": [np.mean, np.std]}).reset_index()
avg_af_diff1 = df_before[df_before["variant"] == variant_list[0]].iloc[0][("af", "mean")]-df_after[df_after["variant"] == variant_list[0]].iloc[0][("af", "mean")]
avg_af_diff2 = df_before[df_before["variant"] == variant_list[1]].iloc[0][("af", "mean")]-df_after[df_after["variant"] == variant_list[1]].iloc[0][("af", "mean")]
max_std1 = np.max([df_before[df_before["variant"] == variant_list[0]].iloc[0][("af", "std")], df_after[df_after["variant"] == variant_list[0]].iloc[0][("af", "std")]])
max_std2 = np.max([df_before[df_before["variant"] == variant_list[1]].iloc[0][("af", "std")], df_after[df_after["variant"] == variant_list[1]].iloc[0][("af", "std")]])
if np.sign(avg_af_diff1) != np.sign(avg_af_diff2) and np.abs(avg_af_diff1) > min_avg_af_diff and np.abs(avg_af_diff2) > min_avg_af_diff and np.abs(avg_af_diff1) > max_std1 and np.abs(avg_af_diff2) > max_std2 and max_std1 < max_af_std and max_std2 < max_af_std:
putative_breakpoints.append(p)
if np.abs(avg_af_diff1) > best_diff and np.abs(avg_af_diff2) > best_diff:
best_diff = np.min([np.abs(avg_af_diff1), np.abs(avg_af_diff2)])
best_breakpoint = p
if len(putative_breakpoints) > 0:
samples_to_putative_breakpoints[s] = putative_breakpoints
samples_to_best_breakpoint[s] = best_breakpoint
print("\nFinished processing.")
print("Number of putative recombinants: " + str(len(samples_to_best_breakpoint)) + "/" + str(len(all_samples_to_check)) + " samples.")
et = time.time()
elapsed_time = et - st
print("-"*100)
print('Execution time: ' + str(round(elapsed_time/60)) + ' minutes') Starting processing...
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Finished processing.
Number of putative recombinants: 46/7290 samples.
----------------------------------------------------------------------------------------------------
Execution time: 87 minutes2.2 Correction for AF-distribution bias
The above initial filtering step resulted in 46 putative subclonal recombinants out of the 7,290 investigated co-infection samples. This drastic decrease in the number of eligible samples underlines the fact that our pipeline uses extremely strict filtering in order to separate true signs of recombination from sequencing noise.
To further refine these findings, we calculated the log-likelihood of observing the measured standardized alternate AFs in the sample given the putative breakpoint and also in a setting with no breakpoints. To determine these quantities, we first fitted bimodal distributions to the standardized alternate AF distribution of each mutually exclusive defining mutation of each variant in each variant combination. Given the \(p(\overline{AF}_{s,g,A}, g, A, A-B)\) probability of observing \(\overline{AF}_{s,g,A}\) standardized alternate AF at a defining mutation of variant \(A\) in a co-infection sample \(s\) with variant composition \(A-B\) at genomic position \(g\), the log-likelihood of measuring the given set of standardized AFs in a sample without a recombination breakpoint is:
\[\log{\mathcal{L}_{nobp; s}} = \sum_{g \in A}\log{p(\overline{AF}_{s,g,A}, g, A, A-B)} + \sum_{g \in B}\log{p(\overline{AF}_{s,g,B}, g, B, A-B)}\]
Assuming a recombination breakpoint at genomic position \(g = b\), the log-likelihood becomes:
\[\log{\mathcal{L}_{bp; s}} = \sum_{g \in A; g \leq b}\log{p(\overline{AF}_{s,g,A}, g, A, A-B)} + \sum_{g \in A; g > b}\log{p(\overline{AF}_{s,g,A}, g, A, A-B)} + \] \[+\sum_{g \in B; g \leq b}\log{p(\overline{AF}_{s,g,B}, g, B, A-B)} + \sum_{g \in B; g > b}\log{p(\overline{AF}_{s,g,B}, g, B, A-B)}\]
Note
Throughout this analysis, we assume that AF distributions at mutually exclusive defining mutations are independent.
Important
Restrictions on the value of \(g\) in the above formulas influence the set of defining mutations considered for the calculation of mean and standard deviation, which in turn affects the value of standardized AFs.
In the following steps, we collect those out of the 46 previously selected putative recombinants, for which the log-likelihood of a no-breakpoint model is lower than that of a breakpoint model.
Code
# fitting bimodal curves to measured AF-distributions at defining mutations
# for each variant combination and each variant separately
variant_combinations = list(df_coinf_allmuts[df_coinf_allmuts["runid"].isin(list(samples_to_best_breakpoint.keys()))]["variants"].apply(get_sorted_variants).unique())
min_number_of_datapoints = 100
variant_distribution_parameters = dict()
i = 1
for vc in variant_combinations:
var1 = vc.split(", ")[0]
var2 = vc.split(", ")[1]
variant_distribution_parameters[vc] = dict()
variant_distribution_parameters[vc][var1] = dict()
variant_distribution_parameters[vc][var2] = dict()
samples_to_check = list(df_coinf_allmuts[df_coinf_allmuts["variants"].apply(get_sorted_variants) == vc]["runid"].unique())
AF_list_var1 = dict()
AF_list_var2 = dict()
for s in samples_to_check:
df_tmp = df_coinf_allmuts[df_coinf_allmuts["runid"] == s]
df_tmp = df_tmp.sort_values(by="pos")
df_tmp_var1 = df_tmp[df_tmp["variant"] == var1]
df_tmp_var2 = df_tmp[df_tmp["variant"] == var2]
mean_var1 = df_tmp_var1["af"].mean()
std_var1 = df_tmp_var1["af"].std()
mean_var2 = df_tmp_var2["af"].mean()
std_var2 = df_tmp_var2["af"].std()
for p in list(df_tmp_var1["pos"]):
if p in AF_list_var1:
AF_list_var1[p].append((df_tmp_var1[df_tmp_var1["pos"]==p].iloc[0]["af"]-mean_var1)/std_var1)
else:
AF_list_var1[p] = [(df_tmp_var1[df_tmp_var1["pos"]==p].iloc[0]["af"]-mean_var1)/std_var1]
for p in list(df_tmp_var2["pos"]):
if p in AF_list_var2:
AF_list_var2[p].append((df_tmp_var2[df_tmp_var2["pos"]==p].iloc[0]["af"]-mean_var2)/std_var2)
else:
AF_list_var2[p] = [(df_tmp_var2[df_tmp_var2["pos"]==p].iloc[0]["af"]-mean_var2)/std_var2]
for p,l in AF_list_var1.items():
if len(l) > min_number_of_datapoints:
d = np.histogram(l, bins=100, density=True)
x = 0.5*(d[1][1:] + d[1][:-1])
y = d[0]
popt,pcov=curve_fit(bimodal, x, y, p0=[np.percentile(l, 25), 1, np.percentile(l, 75), 1])
variant_distribution_parameters[vc][var1][p] = [*popt]
for p,l in AF_list_var2.items():
if len(l) > min_number_of_datapoints:
d = np.histogram(l, bins=100, density=True)
x = 0.5*(d[1][1:] + d[1][:-1])
y = d[0]
popt,pcov=curve_fit(bimodal, x, y, p0=[np.percentile(l, 25), 1, np.percentile(l, 75), 1])
variant_distribution_parameters[vc][var2][p] = [*popt]
# function to calculate log-likelihood for sample
def get_likelihood(pos_list, af_list, variant_comb, variant, bp = None, lowest_number_of_points = 3):
ps = []
params = variant_distribution_parameters[variant_comb][variant]
if bp is None:
if sum([p in params for p in pos_list]) < lowest_number_of_points*2:
return -1000
af_mean = np.mean(af_list)
af_std = np.std(af_list)
for pi, p in enumerate(pos_list):
if p in params:
if (np.abs(af_list[pi]-af_mean) < 1e-10 or af_std < 1e-10):
ps.append(bimodal(0, params[p][0], params[p][1], params[p][2], params[p][3]))
else:
ps.append(bimodal((af_list[pi]-af_mean)/af_std, params[p][0], params[p][1], params[p][2], params[p][3]))
else:
if np.min([sum([p in params and p <= bp for p in pos_list]), sum([p in params and p > bp for p in pos_list])]) < lowest_number_of_points:
return -1000
# before bp
af_mean_before = np.mean(np.array(af_list)[np.array(pos_list) <= bp])
af_std_before = np.std(np.array(af_list)[np.array(pos_list) <= bp])
# after bp
af_mean_after = np.mean(np.array(af_list)[np.array(pos_list) > bp])
af_std_after = np.std(np.array(af_list)[np.array(pos_list) > bp])
for pi, p in enumerate(pos_list):
if p in params:
if p <= bp:
if (np.abs(af_list[pi]-af_mean_before) < 1e-10 or af_std_before < 1e-10):
ps.append(bimodal(0, params[p][0], params[p][1], params[p][2], params[p][3]))
else:
ps.append(bimodal((af_list[pi]-af_mean_before)/af_std_before, params[p][0], params[p][1], params[p][2], params[p][3]))
else:
if (np.abs(af_list[pi]-af_mean_after) < 1e-10 or af_std_after < 1e-10):
ps.append(bimodal(0, params[p][0], params[p][1], params[p][2], params[p][3]))
else:
ps.append(bimodal((af_list[pi]-af_mean_after)/af_std_after, params[p][0], params[p][1], params[p][2], params[p][3]))
return np.log(np.prod(ps))
# finding samples where the breakpoint model has larger log-likelihood than the no-breakpoint model
i = 1
recomb_samples = []
odds_ratios = dict()
recomb_ratios = dict()
for s, bp in samples_to_best_breakpoint.items():
df_tmp = df_coinf_allmuts[df_coinf_allmuts["runid"] == s]
df_tmp = df_tmp.sort_values(by="pos")
variant_comb = df_tmp["variants"].apply(get_sorted_variants).iloc[0]
var1, var2 = variant_comb.split(", ")
# checking for var1
pos_list_c = list(df_tmp[df_tmp["variant"] == var1]["pos"])
af_list_c = list(df_tmp[df_tmp["variant"] == var1]["af"])
var1_ll_nobreakpoint = get_likelihood(pos_list_c, af_list_c, variant_comb, var1)
var1_ll_breakpoint = get_likelihood(pos_list_c, af_list_c, variant_comb, var1, bp)
# checking for var2
pos_list_c = list(df_tmp[df_tmp["variant"] == var2]["pos"])
af_list_c = list(df_tmp[df_tmp["variant"] == var2]["af"])
var2_ll_nobreakpoint = get_likelihood(pos_list_c, af_list_c, variant_comb, var2)
var2_ll_breakpoint = get_likelihood(pos_list_c, af_list_c, variant_comb, var2, bp)
if var1_ll_breakpoint+var2_ll_breakpoint > var1_ll_nobreakpoint+var2_ll_nobreakpoint:
recomb_samples.append(s)
odds_ratios[s] = np.exp((var1_ll_breakpoint+var2_ll_breakpoint)-(var1_ll_nobreakpoint+var2_ll_nobreakpoint))
# calculating recombinant ratio
mean_before_var1 = df_tmp[(df_tmp["variant"] == var1) & (df_tmp["pos"] <= bp)]["af"].mean()
mean_after_var1 = df_tmp[(df_tmp["variant"] == var1) & (df_tmp["pos"] > bp)]["af"].mean()
mean_before_var2 = df_tmp[(df_tmp["variant"] == var2) & (df_tmp["pos"] <= bp)]["af"].mean()
mean_after_var2 = df_tmp[(df_tmp["variant"] == var2) & (df_tmp["pos"] > bp)]["af"].mean()
recomb_ratios[s] = np.mean([np.abs(mean_before_var1-mean_after_var1), np.abs(mean_before_var2-mean_after_var2)])
print("Number of confirmed recombinants: " + str(len(recomb_samples)) + "/" + str(len(samples_to_best_breakpoint)) + " samples.")Number of confirmed recombinants: 13/46 samples.Fig. 5. shows the (non-standardized) AFs measured at defining positions for the remaining 13 samples, along with putative breakpoint positions. Before and after breakpoint means and standard deviations of relevant AFs are marked with horizontal dashed lines and shaded, semi-transparent regions. Odds-ratios of the breakpoint model vs. the no-breakpoint model and estimated ratios of the recombinant genome are also displayed.
Code
# saving list of putative recombinants with their recombination breakpoints
with open("AFbased_recombinants_20230705.csv", "w") as ff:
ff.write("runid,putative_breakpoint\n")
for r in recomb_samples:
ff.write(str(r) + ","+str(samples_to_best_breakpoint[r])+"\n")
# loading SRA IDs of co-infection samples
df_coinf_sra = pd.read_csv("datafile4.csv")
# function to pretty print odds-ratios
def pretty_print(n, digits = 2):
if n < 10:
return str(round(n, digits)), 0
else:
p = 0
while n > 10:
n /= 10
p += 1
return str(round(n, digits)), str(p)
bluecol = "#023e8a"
redcol = "#c1121f"
orangecol = "#F49D37"
greencol = "#0A7055"
variant_color_dict = {"Alpha": orangecol,
"Delta": bluecol,
"Omicron (BA.1)" : redcol,
"Iota": greencol}
fig = plt.figure(figsize=(15,len(recomb_samples)*3.5))
gs = fig.add_gridspec(len(recomb_samples), hspace=0.7)
axs = gs.subplots(sharex=False, sharey=False)
i = 0
for s in recomb_samples:
sra_id = "unknown"
if df_coinf_sra[df_coinf_sra["runid"]==s].shape[0]>0:
sra_id = df_coinf_sra[df_coinf_sra["runid"]==s].iloc[0]["sample_accession"]
best_breakpoint = samples_to_best_breakpoint[s]
df_tmp = df_coinf_allmuts[df_coinf_allmuts["runid"] == s]
variant_list = sorted(df_tmp.iloc[0]["variants"].split(", "))
var1_short = variant_list[0].split("_")[0]
var2_short = variant_list[1].split("_")[0]
if "Omicron" in var1_short:
var1_short = variant_list[0].split("_")[0] + " (" + variant_list[0].split("_")[1] + ")"
if "Omicron" in var2_short:
var2_short = variant_list[1].split("_")[0] + " (" + variant_list[1].split("_")[1] + ")"
axs[i].scatter(df_tmp[df_tmp["variant"] == variant_list[0]]["pos"],
df_tmp[df_tmp["variant"] == variant_list[0]]["af"],
color=variant_color_dict[var1_short], label=var1_short + "-defining mutations")
axs[i].scatter(df_tmp[df_tmp["variant"] == variant_list[1]]["pos"],
df_tmp[df_tmp["variant"] == variant_list[1]]["af"],
color=variant_color_dict[var2_short], label=var2_short + "-defining mutations")
axs[i].legend(bbox_to_anchor=(1, 1.04), loc="lower right",borderaxespad = 0, ncol=2)
t = r'$\mathrm{' + sra_id + '}'
t += '; \mathrm{OR} = '
or_1, or_2 = pretty_print(odds_ratios[s])
if or_2 == 0:
t += or_1
elif or_2 == "1":
t += or_1 + '\cdot 10'
else:
t += or_1 + '\cdot 10^' + or_2
t += '; \mathrm{RR} = ' + str(int(round(recomb_ratios[s]*100, 0))) + '\%$'
axs[i].annotate(t, xy=(0, 1.1), xycoords='data', annotation_clip=False, ha="left", fontsize=12)
axs[i].vlines(x = best_breakpoint, color="black", ymin=0, ymax=1, lw=1)
var1_mean_before = df_tmp[(df_tmp["variant"] == variant_list[0]) & (df_tmp["pos"] <= best_breakpoint)]["af"].mean()
var1_mean_after = df_tmp[(df_tmp["variant"] == variant_list[0]) & (df_tmp["pos"] > best_breakpoint)]["af"].mean()
var2_mean_before = df_tmp[(df_tmp["variant"] == variant_list[1]) & (df_tmp["pos"] <= best_breakpoint)]["af"].mean()
var2_mean_after = df_tmp[(df_tmp["variant"] == variant_list[1]) & (df_tmp["pos"] > best_breakpoint)]["af"].mean()
axs[i].hlines(y=var1_mean_before, xmin=0, xmax = best_breakpoint, color=variant_color_dict[var1_short], ls="dashed", lw=1)
axs[i].hlines(y=var1_mean_after, xmin=best_breakpoint, xmax = 30000, color=variant_color_dict[var1_short], ls="dashed", lw=1)
axs[i].hlines(y=var2_mean_before, xmin=0, xmax = best_breakpoint, color=variant_color_dict[var2_short], ls="dashed", lw=1)
axs[i].hlines(y=var2_mean_after, xmin=best_breakpoint, xmax = 30000, color=variant_color_dict[var2_short], ls="dashed", lw=1)
var1_std_before = df_tmp[(df_tmp["variant"] == variant_list[0]) & (df_tmp["pos"] <= best_breakpoint)]["af"].std()
var1_std_after = df_tmp[(df_tmp["variant"] == variant_list[0]) & (df_tmp["pos"] > best_breakpoint)]["af"].std()
var2_std_before = df_tmp[(df_tmp["variant"] == variant_list[1]) & (df_tmp["pos"] <= best_breakpoint)]["af"].std()
var2_std_after = df_tmp[(df_tmp["variant"] == variant_list[1]) & (df_tmp["pos"] > best_breakpoint)]["af"].std()
axs[i].fill_between(np.linspace(0, best_breakpoint, 100), var1_mean_before-var1_std_before, var1_mean_before+var1_std_before,
alpha=0.1, color=variant_color_dict[var1_short])
axs[i].fill_between(np.linspace(best_breakpoint, 30000, 100), var1_mean_after-var1_std_after, var1_mean_after+var1_std_after,
alpha=0.1, color=variant_color_dict[var1_short])
axs[i].fill_between(np.linspace(0, best_breakpoint, 100), var2_mean_before-var2_std_before, var2_mean_before+var2_std_before,
alpha=0.1, color=variant_color_dict[var2_short])
axs[i].fill_between(np.linspace(best_breakpoint, 30000, 100), var2_mean_after-var2_std_after, var2_mean_after+var2_std_after,
alpha=0.1, color=variant_color_dict[var2_short])
axs[i].grid()
i += 1
for ax in axs:
ax.set_xlabel("", fontsize=0)
ax.set_ylabel("AF", fontsize=14)
ax.set_xlim(0,30000)
ax.set_ylim(0,1)
ax.tick_params(size=5, color="#666666")
axs[-1].set_xlabel("Genomic position", fontsize=14)
plt.show()
plt.close() 