import psycopg2
import pandas as pd
import os
import matplotlib.pyplot as plt
import numpy as np
import copy
plt.style.use('seaborn-whitegrid')
%matplotlib inline
# for pretty table formatting
from IPython.display import Markdown
import sys
from tabulate import tabulateThe CoVEO database is a PostgreSQL database storing the mutational data (VCF files) of SARS-CoV-2 sequencing samples uploaded to the European COVID-19 Data Portal, maintained partly by the efforts of the Versatile Emerging infectious disease Observatory (VEO) consortium. The dataset is unique in the sense that besides the commonly available consensus sequences of the samples, it also contains low alternate allele-frequency mutations and sequencing depth information in a straightforwardly queryable format, along with sample metadata to allow for simple filtering. Additionally, samples of the database are analysed with a standardized variant calling workflow (available on GitHub) in order to keep technical bioinformatics artefacts at a minimum and to obtain comparable results in spite of multiple sample collectors and various laboratory protocols.
This notebook uses PostgreSQL queries and python code to collect and further analyse data. Upon reasonable request, we provide access to the CoVEO database. Please, e-mail kooplex@elte.hu with any inquiries.
def execute_query(query):
try:
host = "xxx.xxx.xxx.xxx" # for login details and inquiries, e-mail kooplex@elte.hu
port = 1234
dbname = "dbname"
username = "username"
pwd = "pwd"
options = "-c search_path=shema_name"
with psycopg2.connect("host='{}' port={} dbname='{}' user={} password={} options='{}'".format(host,
port,
dbname,
username,
pwd,
options)) as conn:
data = pd.read_sql_query(query, conn)
return data
except Exception as e:
print("Something went wrong.", e)Date of latest modification in the database:
execute_query("select max(ml.stop) from merge_log ml")| max | |
|---|---|
| 0 | 2023-05-19 00:40:13.801316 |
Total number of samples with a human host (with appropriate metadata) in the database:
sql = '''\
SELECT COUNT(DISTINCT(mut.runid)) AS total_number_of_human_samples
FROM vcf_key mut
INNER JOIN n_content nc
ON mut.runid = nc.runid
INNER JOIN metaextension me
ON mut.runid = me.runid
INNER JOIN metadata m
ON mut.runid = m.runid
WHERE m.host_id = 3'''
execute_query(sql)| total_number_of_human_samples | |
|---|---|
| 0 | 3093454 |
Total number of good quality samples with a human host (with appropriate metadata) in the database:
sql = '''\
SELECT COUNT(DISTINCT(mut.runid)) AS total_number_of_goodQ_samples
FROM vcf_key mut
INNER JOIN n_content nc
ON mut.runid = nc.runid
INNER JOIN metaextension me
ON mut.runid = me.runid
INNER JOIN metadata m
ON mut.runid = m.runid
WHERE nc.estimated_n_content <= 0.1 AND me.base_count > 100000 AND m.host_id = 3'''
execute_query(sql)| total_number_of_goodq_samples | |
|---|---|
| 0 | 2172927 |
Time range for the collection date of good quality samples:
sql = '''\
SELECT MIN(m.collection_date) AS earliest, MAX(m.collection_date) AS latest
FROM vcf_key mut
INNER JOIN n_content nc
ON mut.runid = nc.runid
INNER JOIN metaextension me
ON mut.runid = me.runid
INNER JOIN metadata m
ON mut.runid = m.runid
WHERE nc.estimated_n_content <= 0.1 AND me.base_count > 100000 AND m.host_id = 3 AND m.collection_date_valid
AND m.collection_date <= me.first_created AND m.collection_date > '2019-01-01 00:00:00'
'''
execute_query(sql)| earliest | latest | |
|---|---|---|
| 0 | 2019-12-30 | 2022-06-30 |
Instead of using a precompiled list of genomic variations characteristic of each SARS-CoV-2 viral strain, we used the marker table provided by Valieris et al. in which all distinguishing mutations are listed with the number of GISAID samples containing the reference and alternate alleles for each lineage.
As an initial step, we selected the mutations listed in this table that were unique and highly indicative of specific viral strains. More precisely, for each mutation, the lineages with the largest and second-largest prevalence were identified. Genomic variations with a largest prevalence of larger than 80% and a second-largest prevalence of less than 10% were considered as “unique defining mutations” of the lineage with the highest mutational incidence. The number of unique defining mutations of each strain are listed below.
The following code generates Supplementary Table 2. and Supplementary Table 3.
# loading the marker table of Valieris et al.
df_markers_base = pd.read_csv("https://github.com/rvalieris/LCS/raw/master/data/pre-generated-marker-tables/pango-designation-markers-v1.9.tsv.gz", sep="\t")
# functions to calculate largest and second-largest prevalences and their corresponding variants
def get_max_af_variant(r):
return r["sample"].split(", ")[np.argmax(r["af"])]
def get_max_af(r):
return np.max(r["af"])
def get_second_largest_af_variant(r):
return r["sample"].split(", ")[np.argsort(r["af"])[-2]]
def get_second_largest_af(r):
return r["af"][np.argsort(r["af"])[-2]]
# finding unique defining mutations
df_markers_base = df_markers_base[df_markers_base["dp"] > 0]
df_markers_base["af"] = df_markers_base["adalt"]/df_markers_base["dp"]
df_markers_all = df_markers_base.groupby(["chrom", "pos", "ref", "alt"]).agg({"sample": lambda x: ", ".join(x),
"af": lambda x: list(x)}).reset_index()
df_markers_all["largest_af_variant"] = df_markers_all.apply(get_max_af_variant, axis=1)
df_markers_all["largest_af"] = df_markers_all.apply(get_max_af, axis=1)
df_markers_all["second_largest_af_variant"] = df_markers_all.apply(get_second_largest_af_variant, axis=1)
df_markers_all["second_largest_af"] = df_markers_all.apply(get_second_largest_af, axis=1)
df_markers = df_markers_all[(df_markers_all["largest_af"] > 0.8) & (df_markers_all["second_largest_af"] < 0.1)]
# generating and saving Supplementary Table 2.
supptable2 = df_markers[["pos", "ref", "alt",
"largest_af", "largest_af_variant", "second_largest_af", "second_largest_af_variant"]]
supptable2["largest_af"] = supptable2["largest_af"]*100
supptable2["second_largest_af"] = supptable2["second_largest_af"]*100
supptable2.rename(columns = {"pos": "genomic position",
"ref": "reference allele",
"alt": "alternate allele",
"largest_af": "largest prevalence (%)",
"largest_af_variant": "variant with largest prevalence",
"second_largest_af": "second largest prevalence (%)",
"second_largest_af_variant": "variant with second largest prevalence"}, inplace=True)
supptable2.to_csv("SuppTable2.csv", index=False)
# simplifying marker table
df_markers = df_markers[["largest_af_variant", "pos", "ref", "alt"]]
df_markers.rename(columns={"largest_af_variant": "sample"}, inplace=True)
# number of defining mutations for each variant
df_marker_nums = df_markers.groupby("sample", as_index=False).count()[["sample", "pos"]].rename(columns={"sample":"variant",
"pos":"num_unique_defmuts"}).sort_values(by="num_unique_defmuts",
ascending=False).reset_index(drop=True)
# generating and saving Supplementary Table 3.
df_marker_nums.rename(columns={"num_unique_defmuts": "number of unique defining mutations"}, inplace=True)
df_marker_nums.to_csv("SuppTable3.csv", index=False)
# pretty-printing table
Markdown(tabulate(
df_marker_nums.to_numpy(),
headers=["Variant","Number of unique defining mutations"]
))| Variant | Number of unique defining mutations |
|---|---|
| Mu_B.1.621 | 22 |
| Gamma_P.1 | 20 |
| Eta_B.1.525 | 19 |
| Alpha_B.1.1.7 | 19 |
| AV.1 | 18 |
| A.23.1 | 17 |
| Theta_P.3 | 16 |
| B.1.1.318 | 15 |
| Omicron_BA.1 | 14 |
| Lambda_C.37 | 13 |
| B.1.623 | 10 |
| Beta_B.1.351 | 8 |
| B.1.617.3 | 8 |
| B.1.177 | 6 |
| Zeta_P.2 | 6 |
| Kappa_B.1.617.1 | 5 |
| Delta_B.1.617.2 | 5 |
| Iota_B.1.526 | 4 |
| Epsilon_B.1.427_429 | 4 |
| Omicron_BA.2.12.1 | 4 |
| Omicron_BA.5 | 3 |
| Omicron_BA.3 | 2 |
| Omicron_BA.4 | 2 |
The altogether 3,093,454 samples of a human host in the CoVEO database were initially filtered to exclude samples that had a total base count of 100,000 or less to avoid misinterpreting sparse sequencing data. Additionally, to further ensure relatively even coverage of the viral genome, we discarded samples that had a sequencing depth of less than 10 in more than 10% of the 29,903 genomic positions of the reference genome (NCBI ID: NC_045512.2). This filtering step resulted in 2,172,927 remaining samples.
Good-quality samples were considered to be putative co-infection samples if
were present in them. Mutations were not filtered for allele frequency.
The following code generates Additional Datafile 1.
sql = '''\
WITH markers_table AS (SELECT j->>'sample' AS variant,
j->>'pos' AS pos,
j->>'ref' AS ref,
j->>'alt' AS alt
FROM json_array_elements('{markersjson}') AS j),
total_markers AS (SELECT m.variant, COUNT(m.pos) AS num_total_markers FROM markers_table m GROUP BY variant),
all_mutations AS (SELECT mut.runid, mut.pos, mut.ref, mut.alt, mt.variant
FROM vcf_key mut
INNER JOIN markers_table mt
ON mut.pos = CAST(mt.pos AS INTEGER) AND mut.ref = mt.ref AND mut.alt = mt.alt
INNER JOIN n_content nc
ON mut.runid = nc.runid
INNER JOIN metaextension me
ON mut.runid = me.runid
INNER JOIN metadata m
ON mut.runid = m.runid
WHERE nc.estimated_n_content <= 0.1 AND me.base_count > 100000 AND m.host_id = 3),
all_samples AS (SELECT mut.runid, mut.variant, COUNT(mut.pos) AS num_defmuts, tm.num_total_markers
FROM all_mutations mut
INNER JOIN total_markers tm
ON mut.variant = tm.variant
GROUP BY mut.runid, mut.variant, tm.num_total_markers),
samples_to_variant_nums AS (SELECT als.runid, COUNT(als.variant) AS num_of_vars, string_agg(als.variant, ', ') AS variants
FROM all_samples als
WHERE CAST(als.num_defmuts AS FLOAT)/als.num_total_markers > 0.5
GROUP BY als.runid),
mixed_samples AS (SELECT samples_to_variant_nums.runid, samples_to_variant_nums.variants
FROM samples_to_variant_nums
WHERE samples_to_variant_nums.num_of_vars > 1),
mixed_samples_mutations_base AS (SELECT vk.key, vk.runid, vk.pos, vk.ref, vk.alt, v.dp, v.af, m.variant
FROM vcf_key vk
INNER JOIN vcf v
ON vk.key = v.key
INNER JOIN markers_table m
ON vk.pos = CAST(m.pos AS INTEGER) AND vk.ref = m.ref AND vk.alt = m.alt
WHERE vk.runid IN (SELECT runid FROM mixed_samples)),
mixed_samples_mutations AS (SELECT msm.runid, msm.pos, msm.ref, msm.alt, msm.dp, msm.af, msm.variant, mixed_samples.variants
FROM mixed_samples_mutations_base msm
INNER JOIN mixed_samples
ON msm.runid = mixed_samples.runid
GROUP BY msm.runid, msm.pos, msm.ref, msm.alt, msm.dp, msm.af, msm.variant, mixed_samples.variants)
SELECT msm.runid, msm.pos, msm.ref, msm.alt, msm.dp, msm.af, msm.variant, msm.variants
FROM mixed_samples_mutations msm
INNER JOIN metadata
ON msm.runid = metadata.runid
'''.format(markersjson = df_markers[["sample", "pos", "ref", "alt"]].to_json(orient="records"))
putative_coinfection_samples_allmuts = execute_query(sql)
# saving Additional datafile 1.
putative_coinfection_samples_allmuts.to_csv("datafile1.csv", index=False)# sorting list of comprising variants alphabetically
def get_sorted_variants(v):
return ", ".join(sorted(v.split(", ")))
print("Number of candidate co-infection samples:", str(len(putative_coinfection_samples_allmuts["runid"].unique())),
"\nNumber of unique variant combinations:", len(putative_coinfection_samples_allmuts["variants"].apply(get_sorted_variants).unique()))Number of candidate co-infection samples: 29666
Number of unique variant combinations: 1270Given that the number of defining mutations in a given variant greatly affects detectability of co-infection cases, we re-evaluated each sample from the above list of putative co-infection samples based on the following procedure:
For the variant composition of the given sample, we identified “mutually exclusive defining mutations” of the comprising variants. We expanded the list of unique markers to all mutations that had a largest prevalence of at least 80% in GISAID samples in one of the variants of the given composition, but a second-largest prevalence of no more than 10% in all other variants of the given variant composition.
If the candidate sample had less than 50% of the mutually exclusive defining mutations of any of its variants, it was discarded from further analysis.
The following code generates Additional Datafile 2. and Supplementary Table 4.
# function to check if a given marker is a mutually exclusive mutation for any of the variants in the combination
def filt_disjunct_markers(r):
try:
rtmp = np.array(r["af"])
if np.max(rtmp) > 0.8 and np.max(rtmp[rtmp != np.max(rtmp)]) < 0.1:
return r["sample"].split(", ")[np.argmax(rtmp)]
else:
return "not marker"
except:
return "not marker"
# all possible variant combinations of candidate co-infection samples
all_variant_combinations = list(putative_coinfection_samples_allmuts["variants"].apply(get_sorted_variants).unique())
# iterating through all variant combinations
putative_coinfection_samples_allmuts_refined = []
disjunct_markers_all = []
for variants in all_variant_combinations:
# selecting mutually exclusive defining mutations of the given variant combination
variant_list = variants.split(", ")
df_markers_tmp = df_markers_base[df_markers_base["sample"].isin(variant_list)].groupby(["chrom", "pos", "ref", "alt"]).agg({"af": lambda x: list(x),
"sample": lambda x: ", ".join(x)}).reset_index()
df_markers_tmp["marker"] = df_markers_tmp.apply(filt_disjunct_markers, axis=1)
df_markers_tmp = df_markers_tmp[df_markers_tmp["marker"] != "not marker"]
df_markers_tmp = df_markers_tmp[["marker", "pos", "ref", "alt"]]
df_markers_tmp.rename(columns={"marker": "sample"}, inplace = True)
df_markers_tmp["variant_combination"] = variants
disjunct_markers_all.append(df_markers_tmp)
df_markers_tmp = df_markers_tmp[["sample", "pos", "ref", "alt"]]
# collecting samples with the given variant combination
runids_to_variant = list(putative_coinfection_samples_allmuts[putative_coinfection_samples_allmuts["variants"].apply(get_sorted_variants) == variants]["runid"])
# re-evaluating all samples of the given variant combination with mutually exclusive defining mutations
sql = '''\
WITH markers_table AS (SELECT j->>'sample' AS variant,
j->>'pos' AS pos,
j->>'ref' AS ref,
j->>'alt' AS alt
FROM json_array_elements('{markersjson}') AS j),
total_markers AS (SELECT m.variant, COUNT(m.pos) AS num_total_markers FROM markers_table m GROUP BY variant),
mixed_samples_mutations_base AS (SELECT vk.key, vk.runid, vk.pos, vk.ref, vk.alt, m.variant
FROM vcf_key vk
INNER JOIN markers_table m
ON vk.pos = CAST(m.pos AS INTEGER) AND vk.ref = m.ref AND vk.alt = m.alt
WHERE vk.runid IN ({runids_to_collect})),
mixed_samples_mut_counts AS (SELECT mut.runid, mut.variant, COUNT(mut.pos) AS num_defmuts, tm.num_total_markers
FROM mixed_samples_mutations_base mut
INNER JOIN total_markers tm
ON mut.variant = tm.variant
GROUP BY mut.runid, mut.variant, tm.num_total_markers),
final_list AS (SELECT msmc.runid, COUNT(msmc.variant) AS num_of_vars, string_agg(msmc.variant, ', ') AS variants
FROM mixed_samples_mut_counts msmc
WHERE CAST(msmc.num_defmuts AS FLOAT)/msmc.num_total_markers > 0.5
GROUP BY msmc.runid),
final_mut_list AS (SELECT msmb.runid, msmb.pos, msmb.ref, msmb.alt, v.dp, v.af, msmb.variant, msmc.num_defmuts, msmc.num_total_markers, fl.variants
FROM mixed_samples_mutations_base msmb
INNER JOIN vcf v
ON v.key = msmb.key
INNER JOIN final_list fl
ON msmb.runid = fl.runid AND position(msmb.variant in fl.variants) > 0
INNER JOIN mixed_samples_mut_counts msmc
ON msmc.runid = msmb.runid AND msmc.variant = msmb.variant)
SELECT * FROM final_mut_list
'''.format(markersjson = df_markers_tmp[["sample", "pos", "ref", "alt"]].to_json(orient="records"),
runids_to_collect = ", ".join(["'"+str(k)+"'" for k in runids_to_variant]))
df_tmp = execute_query(sql)
if df_tmp.shape[0] > 0:
putative_coinfection_samples_allmuts_refined.append(df_tmp)
# mutually exclusive mutations: merging all tables for different variant combinations
disjunct_markers_all = pd.concat(disjunct_markers_all)
disjunct_markers_all.rename(columns={"variant_combination": "putative combination of variants",
"sample": "mutually exclusive mutation of variant",
"pos": "genomic position",
"ref": "reference allele",
"alt": "alternate allele"}, inplace=True)
disjunct_markers_all = disjunct_markers_all[["putative combination of variants",
"mutually exclusive mutation of variant",
"genomic position",
"reference allele",
"alternate allele"]]
# counting number of samples for each variant combination
var_vc = putative_coinfection_samples_allmuts.groupby("runid",
as_index=False).agg({"variants": lambda x: list(x)[0]})["variants"].apply(get_sorted_variants).value_counts()
# only keeping those for which at least 10 samples were identified
disjunct_markers_all = disjunct_markers_all[disjunct_markers_all["putative combination of variants"].isin(var_vc[var_vc > 10].keys())]
# sorting by variant combination frequency
sorter = list(var_vc[var_vc > 10].keys())
disjunct_markers_all["putative combination of variants"] = disjunct_markers_all["putative combination of variants"].astype("category")
disjunct_markers_all["putative combination of variants"] = disjunct_markers_all["putative combination of variants"].cat.set_categories(sorter)
disjunct_markers_all.sort_values(["putative combination of variants", "mutually exclusive mutation of variant"], inplace=True)
# saving Supplementary Table 4.
disjunct_markers_all.to_csv("SuppTable4.csv", index=False, sep=";")
# samples: merging all tables for different variant combinations
putative_coinfection_samples_allmuts_refined = pd.concat(putative_coinfection_samples_allmuts_refined)
# discarding samples with only one variant remaining with a mutually exclusive defining mutation ratio of 0.5 or larger
putative_coinfection_samples_allmuts_refined = putative_coinfection_samples_allmuts_refined[putative_coinfection_samples_allmuts_refined["variants"].str.contains(",")]
# saving Additional datafile 2.
putative_coinfection_samples_allmuts_refined.to_csv("datafile2.csv", index=False)print("Number of refined candidate co-infection samples:", str(len(putative_coinfection_samples_allmuts_refined["runid"].unique())),
"\nNumber of unique variant combinations for refined samples:", len(putative_coinfection_samples_allmuts_refined["variants"].apply(get_sorted_variants).unique()))Number of refined candidate co-infection samples: 22180
Number of unique variant combinations for refined samples: 711The number of supposed co-infection samples was then determined with multiple thresholds for the ratio of required mutually exclusive defining mutations in the range of 0.5 to 1, and the final filtering limit of 0.8 was chosen for the identification of a total number of 7,700 co-infection samples.
# calculating mutually exclusive defining mutation ratio for all comprising variants for all samples
df_tmp = putative_coinfection_samples_allmuts_refined.groupby(["runid", "variant", "num_defmuts", "num_total_markers", "variants"]).count().reset_index()
df_tmp = df_tmp[["runid", "variant", "num_defmuts", "num_total_markers", "variants"]]
df_tmp["defmut_ratio"] = df_tmp["num_defmuts"]/df_tmp["num_total_markers"]
df_tmp = df_tmp.groupby(["runid", "variants"]).agg({"defmut_ratio": lambda x: list(x),
"variant": lambda x: ", ".join(list(x))}).reset_index()
# function to check if all mutually exclusive defining mutation ratios are larger than tr
def is_all_larger(v, tr):
return all([vv >= tr for vv in v])
# calculating number of co-infection samples for multiple thresholds
thresholds = np.linspace(0.5,1,100)
num_samples = []
for tr in thresholds:
num_samples.append(df_tmp["defmut_ratio"].apply(is_all_larger, args=(tr,)).sum())
# total number of good-quality samples in the database
samples_total_num = 2172927
# calculating ratio of co-infection samples
ratio_samples = [k*100/samples_total_num for k in num_samples]
# plotting the results
fig, ax1 = plt.subplots(figsize=(5, 5))
ax2 = ax1.twinx()
ax1.plot(thresholds, num_samples, color="#03045e", lw=2)
ax1.set_xlim(0.5,1)
ax1.set_ylim(0,25000)
yt = [5000*k for k in range(1,6)]
ax1.set_yticks(yt)
ax1.set_xlabel("Threshold for ratio of mutually exclusive\n defining mutations present", fontsize=12)
ax1.set_ylabel("Number of co-infection samples", fontsize=12)
ax1.annotate(str(df_tmp.shape[0]),
xy = (0.495, df_tmp.shape[0]), fontsize=10,
ha="right", va="center", annotation_clip = False)
ax1.annotate("Number of co-infection samples with \nall mutually exclusive \ndefining mutations present: " + str(num_samples[-1]),
xy = (0.51, 110), fontsize=10, color="black",
ha="left", va="bottom", annotation_clip = False)
ax2.set_yticks([round(k*100/samples_total_num,2) for k in yt])
ax2.set_ylim(0,30000*100/samples_total_num)
ax2.annotate("Percentage of co-infection samples (%)",
xy = (1.12, 0.5), fontsize=12, rotation=-90,
ha="center", va="center", xycoords = "axes fraction")
ax2.fill_between([0.5, 1], 0.2, 0.5, color="#caf0f8", alpha=0.4)
ax2.annotate("Literary range for \nco-infection rate: 0.2-0.5%",
xy = (0.51, 0.22), fontsize=10, color="#00b4d8",
ha="left", va="bottom", annotation_clip = False)
ax2.vlines(x=0.8, ymin=0, ymax=30000*100/samples_total_num, color="#ae2012", linestyle="dashed")
ax2.annotate("Threshold for co-infection\n detection: 0.8",
xy = (0.79, 1.28), fontsize=10, color="#ae2012",
ha="right", va="top", annotation_clip = False)
ax2.annotate("Final nuber of\n co-infection\n samples: " + str(df_tmp["defmut_ratio"].apply(is_all_larger, args=(0.8,)).sum()),
xy = (0.98, 0.48),
fontsize=10, color="black",
ha="right", va="top", annotation_clip = False)
plt.grid()
plt.show()
This corresponds to a co-infection rate of:
print(round(7700*100/samples_total_num,2), "%", sep="")0.35%The following code generates Additional Datafile 3.
# getting list of mutually exclusive defining mutations in co-infection samples
coinf_samples = list(df_tmp[df_tmp["defmut_ratio"].apply(is_all_larger, args=(0.8,))]["runid"])
coinfection_samples_allmuts = putative_coinfection_samples_allmuts_refined[putative_coinfection_samples_allmuts_refined["runid"].isin(coinf_samples)]
# saving Additional datafile 3.
coinfection_samples_allmuts.to_csv("datafile3.csv", index=False)
# variant compositions with at least 50 co-infection samples
k = coinfection_samples_allmuts.groupby(["runid", "variants"]).count().reset_index()["variants"].apply(get_sorted_variants).value_counts()
top_var_comp = k[k>50]
# plotting the results
fig, ax = plt.subplots(figsize=(8, 3))
bars = ax.bar(np.arange(len(top_var_comp)), top_var_comp, width=0.9, color="#caf0f8")
for bar in bars:
height = bar.get_height()
ax.annotate(f'{height}', xy=(bar.get_x() + bar.get_width() / 2, height), xytext=(0, 3),
textcoords="offset points", ha='center', va='bottom')
ax.grid(visible=True)
ax.set_ylabel("Number of \nco-infection samples", fontsize=12)
ax.set_xticks(np.arange(len(top_var_comp)))
ax.set_xticklabels(top_var_comp.keys(), rotation=35, ha="right", va="top", rotation_mode='anchor')
ax.set_ylim(0,3000)
plt.show()
sql = '''\
SELECT l.variant_id, COUNT(DISTINCT(l.runid)) AS number_of_goodQ_samples
FROM vcf_key mut
INNER JOIN n_content nc
ON mut.runid = nc.runid
INNER JOIN metaextension me
ON mut.runid = me.runid
INNER JOIN metadata m
ON mut.runid = m.runid
LEFT OUTER JOIN lineage l
ON mut.runid = l.runid
WHERE nc.estimated_n_content <= 0.1 AND me.base_count > 100000 AND m.host_id = 3
GROUP BY l.variant_id'''
df_allvariants = execute_query(sql)
def get_variant_name_short(r):
if not pd.isnull(r):
return r.split(" (")[0]
else:
return r
df_allvariants["variant_name"] = df_allvariants["variant_id"].apply(get_variant_name_short)
df_allvariants = df_allvariants.groupby("variant_name", as_index=False).agg({"number_of_goodq_samples": sum}).sort_values(by="number_of_goodq_samples", ascending=False)
fig, ax = plt.subplots()
ax.bar(list(df_allvariants[df_allvariants["variant_name"] != "Other variant"]["variant_name"])+["Other variant"],
list(df_allvariants[df_allvariants["variant_name"] != "Other variant"]["number_of_goodq_samples"])+[741728],
color=["#caf0f8"]*(df_allvariants.shape[0]-1)+["lightgray"], width=0.9)
ax.set_yscale("log")
ax.set_xticklabels(list(df_allvariants[df_allvariants["variant_name"] != "Other variant"]["variant_name"])+["Other variant"],
rotation=35, ha="right", va="top", rotation_mode='anchor')
ax.set_ylabel("Number of good-quality samples \nin the CoVEO database assigned to variant", fontsize=12)
plt.show()
plt.close()
The figure below shows the number of co-infection samples that include the given variant in function of the number of good quality samples assigned to the given variant in the database. (Note the logarithmic axes.) It is apparent that the more samples are available from a given variant, the more likely it is to detect co-infection samples that include that variant in their variant composition. The grey line represents a linear relationship between the two values, i.e. a straight line with a slope of 1 in a log-log plot.
df_coinf_variants = coinfection_samples_allmuts.groupby(["runid", "variants"]).count().reset_index()
number_of_coinf_samples = []
for ri,r in df_allvariants.iterrows():
num_coinf = df_coinf_variants[df_coinf_variants["variants"].str.contains(r["variant_name"])].shape[0]
number_of_coinf_samples.append(num_coinf)
df_allvariants["number_of_coinf_samples_with_variant"] = number_of_coinf_samples
avg_factor = (df_allvariants[df_allvariants["variant_name"] != "Other variant"]["number_of_coinf_samples_with_variant"]/df_allvariants[df_allvariants["variant_name"] != "Other variant"]["number_of_goodq_samples"]).mean()
fig, ax = plt.subplots()
ax.scatter(df_allvariants[df_allvariants["variant_name"] != "Other variant"]["number_of_goodq_samples"],
df_allvariants[df_allvariants["variant_name"] != "Other variant"]["number_of_coinf_samples_with_variant"])
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel("Number of good-quality samples \nin the CoVEO database assigned to variant", fontsize=12)
ax.set_ylabel("Number of co-infection samples \nincluding the variant", fontsize=12)
ax.plot(np.linspace(0,df_allvariants["number_of_goodq_samples"].max(),1000),
np.linspace(0,df_allvariants["number_of_goodq_samples"].max(),1000)*avg_factor, color="lightgray")
ax.set_ylim(0,10000);
for ri,r in df_allvariants.iterrows():
if r["variant_name"] == "Other variant":
continue
op = (0,7)
if r["variant_name"] == "Eta":
op = (10,7)
if r["variant_name"] == "Mu":
op = (20,0)
if r["variant_name"] == "Zeta":
op = (-20,3)
if r["variant_name"] == "Gamma":
op = (-20,5)
if r["variant_name"] == "Omicron":
op = (0,-15)
if r["variant_name"] == "Iota":
op = (10,-15)
ax.annotate(r["variant_name"], xy = (r["number_of_goodq_samples"], r["number_of_coinf_samples_with_variant"]),
xytext = op,
textcoords = "offset points", ha="center")
plt.show()
plt.close()
Once the set of co-infection samples have been identified, we further query the database to collect their available metadata.
The following code generates Additional Datafile 4.
sql = """
SELECT m.runid, m.sample_accession, m.experiment_accession, m.study_accession,
m.collection_date, m.collection_date_valid, me.first_created,
c.country_name,
i.instrument_platform, i.instrument_model,
h.host, me.description
FROM metaextension me
LEFT OUTER JOIN metadata m
ON m.runid = me.runid
LEFT OUTER JOIN country c
ON c.id = m.country_id
LEFT OUTER JOIN host h
ON h.id = m.host_id
LEFT OUTER JOIN instrument i
ON i.id = m.instrument_id
WHERE me.runid in ({sample_list_string});
""".format(sample_list_string = ", ".join(["'"+str(k)+"'" for k in coinf_samples]))
coinf_samples_detailed_metadata = execute_query(sql)
# joining with variant info
coinf_samples_detailed_metadata = pd.merge(coinfection_samples_allmuts.groupby("runid", as_index=False).agg({"variants": lambda x: list(x)[0]}),
coinf_samples_detailed_metadata, on="runid", how="outer")
coinf_samples_detailed_metadata["variants"] = coinf_samples_detailed_metadata["variants"].apply(get_sorted_variants)
# saving Additinal Datafile 4.
coinf_samples_detailed_metadata = coinf_samples_detailed_metadata[["runid", "sample_accession", "experiment_accession", "study_accession",
"variants", "collection_date", "collection_date_valid", "first_created", "country_name",
"instrument_platform", "instrument_model", "host", "description"]]
coinf_samples_detailed_metadata.to_csv("datafile4.csv", index=False)To see how co-infection samples are distributed across different studies, we determined the number of good-quality samples of a human host with available mutation information of various study accession IDs in the CoVEO database. Study-specific prevalence rates were calculated as the percentage of identified co-infection samples within a given study.
The following code generates Supplementary Table 1.
sql = '''
SELECT m.study_accession, COUNT(m.runid) AS goodQ_study_samples
FROM metadata m
INNER JOIN n_content nc
ON m.runid = nc.runid
INNER JOIN metaextension me
ON m.runid = me.runid
WHERE nc.estimated_n_content <= 0.1 AND me.base_count >= 100000 AND m.host_id = 3 AND m.study_accession IN ({sample_list_string})
GROUP BY m.study_accession;
'''.format(sample_list_string = ", ".join(["'"+str(k)+"'" for k in coinf_samples_detailed_metadata[~pd.isnull(coinf_samples_detailed_metadata["study_accession"])]["study_accession"].unique()]))
df_study_data = execute_query(sql)
df_study_data = pd.merge(coinf_samples_detailed_metadata.groupby("study_accession", as_index=False).count()[["study_accession",
"runid"]],
df_study_data, on = "study_accession")
df_study_data.rename(columns={"study_accession": "ENA study accession",
"runid": "number of co-infection samples",
"goodq_study_samples": "number of good-quality samples in the study"}, inplace = True)
df_study_data["co-infection prevalence (%)"] = 100*df_study_data["number of co-infection samples"]/df_study_data["number of good-quality samples in the study"]
df_study_data.sort_values(by = "number of good-quality samples in the study", ascending = False, inplace = True)
# saving Supplementary Table 1.
df_study_data.to_csv("SuppTable1.csv", index=False)# pretty-printing table
Markdown(tabulate(
df_study_data[df_study_data["co-infection prevalence (%)"] >= 10].to_numpy(),
headers=["ENA study accession","Number of co-infection samples",
"Number of good quality samples in the study", "Co-infection prevalence (%)"]
))| ENA study accession | Number of co-infection samples | Number of good quality samples in the study | Co-infection prevalence (%) |
|---|---|---|---|
| PRJNA817870 | 117 | 152 | 76.9737 |
| PRJNA817806 | 44 | 101 | 43.5644 |
| PRJNA853723 | 82 | 97 | 84.5361 |
| PRJNA809680 | 5 | 8 | 62.5 |
| PRJNA827817 | 5 | 7 | 71.4286 |
| PRJNA748832 | 1 | 6 | 16.6667 |
| PRJNA698337 | 1 | 6 | 16.6667 |
| PRJNA728440 | 1 | 4 | 25 |
| PRJNA804575 | 1 | 1 | 100 |
The listed studies either contain very few samples or were specifically pre-selected to include large amounts of co-infection cases (studies PRJNA817870, PRJNA827817, PRJNA853723, PRJNA817806, PRJNA809680).
To see how co-infection samples are distributed across different countries, we determined the number of good-quality samples of a human host with available mutation information of various countries in the CoVEO database. Country-specific prevalence rates were calculated as the percentage of identified co-infection samples within a given country.
sql = '''
SELECT c.country_name, COUNT(m.runid) AS goodQ_study_samples
FROM metadata m
INNER JOIN country c
ON c.id = m.country_id
INNER JOIN n_content nc
on nc.runid = m.runid
INNER JOIN metaextension m2
on m2.runid = m.runid
WHERE m2.base_count >= 100000 and nc.estimated_n_content <= 0.1
GROUP BY c.country_name;
'''
df_country = execute_query(sql)
df_country_data = pd.merge(coinf_samples_detailed_metadata.groupby("country_name", as_index=False).count()[["country_name",
"runid"]],
df_country, on = "country_name", how="outer")
df_country_data.replace(np.nan, 0, inplace=True)
df_country_data.rename(columns = {"runid": "number of co-infection samples",
"goodq_study_samples": "number of good-quality samples in country"}, inplace=True)
df_country_data["co-infection prevalence (%)"] = 100*df_country_data["number of co-infection samples"]/df_country_data["number of good-quality samples in country"]