Traces of read-level recombination in co-infection samples

Supplementary Methods 3. for manuscript ‘Systematic detection of co-infection and intra-host recombination in more than 2 million global SARS-CoV-2 samples’

Author
Affiliation

Orsolya Anna Pipek

Department of Physics of Complex Systems, ELTE Eötvös Loránd University, Budapest, Hungary

Published

July 10, 2023

Contents
In this notebook, we analyse a subset of previously identified co-infection samples for the presence of short reads that overlap variant-defining mutations of both parental strains of the sample and carry both of these mutations simultaneously. We investigate how the genomic distribution of canonical lineage-specific mutations inherently affects the distribution of overlapping reads and thus the detectability of recombination breakpoints. We further examine reads carrying traces of recombination for recombination hotspots, traces of subgenomic RNA, comparability to AF-based analysis results and possible chimeric origin. As an introduction, we discuss the technical and theoretical challenges of subclonal recombinant detection.
Code 
import pandas as pd
import re
import os
import itertools
from collections import Counter
from scipy.stats.stats import pearsonr
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
plt.style.use('seaborn-whitegrid')
import numpy as np
%matplotlib inline

1 Difficulties in detecting recombinants

The detection of subclonal recombinant genomes in co-infection samples is hampered by multiple factors that together make it nigh impossible to reliably distinguish between true evidence of recombination and artefacts. Here we briefly summarize the main causes of this difficulty.

  • Subclonality: Given that co-infection samples usually contain significant amounts of the two original parental strains, recombinant genomes comprise only a small portion of the sequenced viral population. Thus any attempt at identifying recombinants must reckon with decreased coverage and consequently a limited amount of available data.
  • PCR artefacts: PCR amplification is the standard method for the generation of sufficient genetic material prior to sequencing. Most of the sequencing data for SARS-CoV-2 has been produced by a pipeline that incorporates PCR amplification in its initial steps. Admittedly, there are a few experimental setups in which metagenomic sequencing was employed, which eliminates various problematic issues introduced by PCR-based methods, however, the abundance of SARS-CoV-2 genomes in these samples is usually too low to draw meaningful conclusions.
    • Systematic bias in alternate allele frequency distribution: It has been previously shown that alternate allele frequencies measured at defining mutations in artificial samples of mixed variants do not correctly reflect the original mixture proportions in the sample (Bal et al., 2022). Additionally, based on our observations, a systematic bias can be identified in the alternate allele frequency distribution of defining mutations in co-infection samples of specific variant combinations (see Supplementary file 2.). We hypothesize that this effect might be due to the preferential attachment of primers to one of the parental viral strains that carries a set of favourable or lacks a set of disadvantageous mutations. As a result of this, the presence of subclonal recombinant genomes cannot be reliably detected by subtle shifts of alternate allele frequencies along the genome.
    • Chimeras: It has been known for decades that during PCR amplification, PCR-mediated recombination or “chimera formation” systematically occurs (Brakenhoff et al., 1991), generating artificial sequences that are essentially no different from true recombinants. It is virtually impossible to dependably distinguish between these, thus one has to assume that chimera formation is relatively rare, while viral recombination is well-documented in laboratory settings, hence also expected to occur in co-infection samples.
  • Subgenomic RNA: Besides virus genome length RNA, diagnostic samples of SARS-CoV-2 have been shown to carry leader sequence-containing subgenomic RNAs (sgRNAs) as well (Kim et al., 2020, Alexandersen et al., 2020). Recombination occurring in the sgRNAs has no effect on viral evolution and cannot be passed to future viral populations, thus its presence is less relevant than recombination of the genomic RNA. Short reads showing signs of recombination that contain the so-called leader sequence and/or were soft-clipped during alignment can be relatively confidently categorized as originating from sgRNA. However, due to short read lengths and fairly long sequences of the transcriptome, a short read might still be sgRNA-derived, even without the presence of the leader sequence.
  • The low number of defining mutations: Given that, disregarding the relatively low number of defining mutations, parental strains in SARS-CoV-2 co-infection samples are highly similar, recombination events might go completely undetected. The identification of recombination breakpoints is limited to the genomic ranges between defining mutations, thus the uncertainty in their location is extremely high.
  • Uneven distribution of defining mutations: Defining mutations are unevenly distributed across genomic positions, a disproportionately high amount (considering gene lengths) of them are located on genes S and N (Fig. 1), making it very difficult to detect recombination breakpoints occurring in other, less frequently mutated regions of the genome.
  • Short read lengths: Direct evidence of recombination events can come from short reads that simultaneously contain the defining mutations of both parental strains in a co-infection sample. This approach, however, is limited by the relatively short read lengths (100-200 bp) in sequencing data generated by Illumina platforms, as only those defining mutation pairs are overlapped by the same reads that are located close enough on the genome. The recent advances in Nanopore sequencing technologies might provide a solution for this problem, as they usually generate reads ranging from 10 to 100 kbp in length.

In the following analysis, we collect reads overlapping defining mutations of multiple variant stains and show how their detectability is largely influenced by the presence and distribution of variant-defining mutations. We also examine them for traces of recombination and discuss the resulting distribution of recombination breakpoints and how it is biased by the factors described above.

Code 
# gene locations
gene_loc = {
    "ORF1ab": [265, 21555],
    "S": [21562, 25384],
    "ORF3a": [25392, 26220],
    "E": [26244, 26472],
    "M": [26522, 27191],
    "ORF6": [27201, 27387],
    "ORF7a": [27393, 27759],
    "ORF7b": [27755, 27887],
    "ORF8": [27893, 28259],
    "N": [28273, 29533],
    "ORF10": [29557, 29674]}

color_dict = {
    "ORF1ab": "#001219",
    "S": "#10454f",
    "ORF3a": "#005f73",
    "E": "#0a9396",
    "M": "#94d2bd",
    "ORF6": "#e9d8a6",
    "ORF7a": "#ee9b00",
    "ORF7b": "#c96602",
    "ORF8": "#bb3e03",
    "N": "#ae2012",
    "ORF10": "#9b2226" 
}

# mutually exclusive defining mutations in all co-infection samples
df_coinf_allmuts = pd.read_csv("datafile3.csv")

# generate ID for unique mutation
def get_mutation_base(r):
    return str(r["pos"])+"_"+r["ref"]+"_"+r["alt"]

# function to associate mutation with gene
def get_gene_for_pos(p):
    for g in gene_loc.keys():
        if p >= gene_loc[g][0] and p <= gene_loc[g][1]:
            return g

# figure to show defining mutations along the genome        
fig = plt.figure(figsize=(15,3))
gs = fig.add_gridspec(2, hspace=0.5, height_ratios = [0.2, 0.8])
axs = gs.subplots(sharex=True)

yy = 2
for g in gene_loc.keys():
    axs[0].hlines(y = 0, xmin=gene_loc[g][0], xmax=gene_loc[g][1], lw=30, color = color_dict[g])
    if (yy == 2):
        va_c = "bottom"
        ha_c = "left"
    else:
        va_c = "top"
        ha_c = "right"
    axs[0].annotate(g, xy= (np.mean([gene_loc[g][0], gene_loc[g][1]]), yy*0.5), 
                 xytext = (np.mean([gene_loc[g][0], gene_loc[g][1]]), yy), rotation=30, ha="center", va = va_c, fontsize=8,
                arrowprops=dict(arrowstyle="-"))
    yy *= -1
axs[0].set_ylim(-1,1)
axs[0].set_xlim(0,30000)
axs[0].set_axis_off()

mutbase = df_coinf_allmuts.apply(get_mutation_base, axis=1)
axs[1].hist(df_coinf_allmuts.groupby(mutbase).agg({"pos": lambda x: list(x)[0]})["pos"], range=(0,30000), bins=60, 
            color="#023e8a")
axs[1].set_xlabel("Genomic position (binned by 500 bp)", fontsize=10)
axs[1].set_ylabel("Number of unique mutually \nexclusive defining mutations", fontsize=10)
axs[1].set_xlim(0,30000)
axs[1].set_ylim(0,30)
axs[1].grid(visible=False)
axs[1].tick_params(size=5, color="#666666")
axs[1].label_outer()

ymin, ymax = axs[1].get_ybound()
for g in gene_loc.keys():
    axs[1].vlines(x=gene_loc[g][1], ymin=ymin, ymax=ymax, color="black", alpha=0.5, lw=0.5)
axs[1].vlines(x=gene_loc["ORF1ab"][0], ymin=ymin, ymax=ymax, color="black", alpha=0.5, lw=0.5)    
    
plt.show()
plt.close()

# figure to show defining mutation density in genes
gg = Counter(df_coinf_allmuts.groupby(mutbase).agg({"pos": lambda x: list(x)[0]})["pos"].apply(get_gene_for_pos))

c = 0
plt.figure(figsize=(15,1))

yy = 0.05
for g in gene_loc.keys():
    local_dens = gg[g]*1000/(gene_loc[g][1]-gene_loc[g][0])
    plt.hlines(y = 0, xmin=c, xmax=c+local_dens, lw=30, color = color_dict[g])
    if yy > 0:
        ha_c = "left"
        va_c = "bottom"
    else:
        ha_c = "right"
        va_c = "top"
    plt.annotate(g+": " +str(round(local_dens,2)) + "/kbp", 
                 xy=(c+local_dens/2, yy), 
                 ha=ha_c, va = va_c, fontsize=8, color="black", rotation=30)
    yy *= -1
    c += local_dens
    
plt.axis("off")
plt.show()
plt.close()

(a) Distribution of mutually exclusive defining mutations along the SARS-CoV-2 genome, binned by 500 bp. All mutations are included that are considered to be mutually exclusive defining mutations of a variant in any of the variant combinations of the 7,413 co-infection samples. Each specific mutation has been counted once to create this figure. The top panel shows the locations of specific genes along the genome. Vertical lines on the bottom panel indicate gene boundaries.

(b) Mutually exclusive defining mutation density in various genes. Mutations not overlapping any of the canonical genes have been discarded for this figure. The widths of coloured regions are proportional to the density of defining mutations in the genes.

Figure 1: Location and density of mutually exclusive defining mutations in the SARS-CoV-2 genome.

2 Obtaining aligned sequencing data from ENA

The analyses in this notebook use raw, aligned sequencing data (BAM files) to identify reads that carry defining mutations of multiple variants. In order to perform these investigations, the database of the European Nucleotide Archive (ENA) was queried for the location of BAM files of a list of previously selected co-infection samples and data was then locally downloaded for downstream analysis.

Note

The original set of previously detected 7,413 co-infection samples was first limited to those 6,999 that contained traces of exactly two variants. (I.e. samples identified as the mixtures of three or more variant strains were discarded.) Then the resulting set of samples was further downsampled to a set of 100 samples for this analysis to decrease computation time and simultaneously preserve the prevalence of specific variant combinations. The 13 co-infection samples identified as putative subclonal recombinants based on alternate allele frequency (AF) shifts along their genome (see Supplementary File 2) were added to this list, along with 5 artificial mixture samples of study PRJNA827817, resulting in altogether 118 samples for this analysis.

Code 
# loading data of samples to analyse
df_samples = pd.read_csv("datafile5.tar.gz")

# function to write ENA query
def get_mixed_samples_ENA_string(list_of_acc):
    return " OR ".join(list(set(['sample_accession="'+k+'"' for k in list_of_acc])))

# ENA query
s = """curl -X POST -H "Content-Type: application/x-www-form-urlencoded" -d 'result=analysis&query={ENA_string}&fields=sample_accession,analysis_alias,analysis_accession,analysis_type,submitted_ftp&format=tsv' "https://www.ebi.ac.uk/ena/portal/api/search" > {output}""".format(ENA_string=get_mixed_samples_ENA_string(list(df_samples["sample_accession"])),
                                                                  output = "download_info_for_read_analysis_20230705.tsv")

# runnig ENA query
!$s

# loading file paths for queried samples
df_download = pd.read_csv("download_info_for_read_analysis_20230705.tsv", sep="\t")

# retrieving file paths for downloading
files_to_download = list(df_download[(df_download["submitted_ftp"].str.contains("output.tar.gz")) | 
                                     (df_download["submitted_ftp"].str.contains(".bam"))]["submitted_ftp"])

# downloading the files
!mkdir -p data_for_recombinant_read_analysis

for fp in files_to_download:
    filelist = fp.split(";")
    ff = [k for k in filelist if ".bam" in k or "output.tar.gz" in k][0]
    !sleep 5
    !wget $ff -P data_for_recombinant_read_analysis/
    
# unzipping tar files
!cd data_for_recombinant_read_analysis/; for ffile in *.tar.gz; do tar xf "$ffile"; done

# sorting, indexing BAM files
!mkdir -p data_for_recombinant_read_analysis/indexed_bams
bam_files = !find "data_for_recombinant_read_analysis/" -type f -name "*.bam"
for bf in bam_files:
    bfs = "data_for_recombinant_read_analysis/indexed_bams/"+bf.split("/")[-1].split(".bam")[0]+"_sorted.bam"
    !samtools sort $bf -o $bfs
    !samtools index $bfs

3 Detection of overlapping and recombinant reads

We principally used the mpileup command of the samtools software to identify reads that overlapped at least one mutually exclusive defining mutation of both the comprising variant strains and carried the alternate allele in both of these genomic positions. The following pipeline was implemented for the analysis:

  1. Get the pileup of the BAM file for all genomic positions in the sample where a mutually exclusive mutation of one of the comprising variant strains is present.
    • filter for only these positions (-l position_list.txt)
    • disable base alignment quality (BAQ) computation (-B)
    • filter for both base (-Q 30) and mapping quality (-q 30)
    • disable sequencing depth cutoff (-d 0)
    • also output read IDs for each detected base (--output-QNAME)
  2. Find all possible pairs of genomic positions (from the above list), where one comprising strain has a mutually exclusive defining mutation in one of the positions and the other strain has one in the other position. Query the pileup for these pairs of genomic positions.
  3. Determine the number of reads overlapping both positions, based on the read IDs provided in the pileup.
  4. Categorize each overlapping read as one of the following:
    • the read carries the mutually exclusive defining mutation of one variant strain,
    • the read carries the mutually exclusive defining mutation of the other variant strain,
    • the read carries the mutually exclusive defining mutation of both variant strains (i.e. supports recombination),
    • the read carries neither mutually exclusive defining mutations.
  5. Calculate the numbers of reads belonging to each group for all relevant mutation pairs in the sample and save the results.
Code 
# functions for processing
def get_sorted_variants(v):
    return ", ".join(sorted(v.split(", ")))

def get_mutation_id(r):
    return str(r["pos"])+"_"+r["ref"]+"_"+r["alt"]+"_"+var_dict[r["variant"]]

def get_short_name_from_fname(fn):
    flist = fn.split(";")
    ff = [k for k in flist if ".bam" in k or "output.tar.gz" in k][0]
    return ff.split("/")[-1].split("_")[0].split(".")[0]

def prepare_sample_details(runid):
    shortname = df_samples[df_samples["runid"] == runid].iloc[0]["shortname"]
    bam_filepath = "data_for_recombinant_read_analysis/indexed_bams/" + shortname + "_sorted.bam"
    var_comb = df_coinf_allmuts[df_coinf_allmuts["runid"] == runid]["variants"].apply(get_sorted_variants).iloc[0]
    # correcting variant combination for single sample with 3 variants (only the two most abundant ones are considered)
    if runid == 4424071: var_comb = get_sorted_variants("Delta_B.1.617.2, Alpha_B.1.1.7")
    position_list_file = "data_for_recombinant_read_analysis/poslists/" + "_".join([var_dict[v] for v in var_comb.split(", ")]) + "_poslist.txt"
    return var_comb, position_list_file, bam_filepath, shortname

def get_mutation_pairs(runid):
    df_tmp = df_coinf_allmuts[df_coinf_allmuts["runid"] == runid]
    all_relevant_muts = list(df_tmp.apply(get_mutation_id, axis=1))
    all_mut_pairs = []
    for p in itertools.combinations(all_relevant_muts, 2):
        if p[0].split("_")[-1] != p[1].split("_")[-1] and p[0].split("_")[0] != p[1].split("_")[0]:
            all_mut_pairs.append(p)
    return all_mut_pairs

def correct_mutation_string(ref, mut):
    # deletion
    if len(ref) > len(mut):
        mut = "-"+str(len(ref)-1)+ref[1:]
    # insertion
    if len(ref) < len(mut):
        mut = "+"+str(len(mut)-1)+mut[1:]
    return mut

def get_base_list(base_string):
    b = re.sub("\^.", "", base_string)
    b = re.sub("\$", "", b)
    if not "+" in b and not "-" in b:
        return list(b)
    else:
        bases = []
        skipchars = 0
        for ci,c in enumerate(b):
            if skipchars == 0:
                if c != "+" and c != "-":
                    bases.append(c)
                # deletion or insertion
                if c == "-" or c == "+":
                    if len(bases) > 0:
                        bases = bases[:-1]
                    delnum = int(re.findall("[0-9]+", b[ci+1:])[0])
                    bases.append(b[ci:ci+delnum+1+len(str(delnum))])
                    skipchars = delnum+len(str(delnum))
            else:
                skipchars -= 1
        return bases

# joining sample IDs with comprising variants
df_samples_to_variants = df_coinf_allmuts.groupby("runid").agg({"variants": lambda x: list(x)[0]}).reset_index()
df_samples = pd.merge(df_samples, df_samples_to_variants, on="runid")
df_samples["variants"] = df_samples["variants"].apply(get_sorted_variants)
# correcting variant combination for single sample with 3 variants (only the two most abundant ones are considered)
df_samples["variants"].mask(df_samples["runid"] == 4424071, "Delta_B.1.617.2, Alpha_B.1.1.7", inplace=True )

# variant coding dictionary
var_dict = {'Alpha_B.1.1.7': "A",
 'Epsilon_B.1.427_429': "E",
 'Gamma_P.1': "G",
 'Iota_B.1.526': "I",
 'Delta_B.1.617.2': "D",
 'Zeta_P.2': "Z",
 'B.1.177': "20E",
 'Omicron_BA.1':"OBA1",
 'Mu_B.1.621': "M", 
 'Lambda_C.37': "L",
 'B.1.617.3': "B16173",
 'Omicron_BA.2.12.1': "OMBA2",
 'Omicron_BA.5': "OMBA5",
 'B.1.1.318': "B11318",
 'Beta_B.1.351': "B",
 'Omicron_BA.3': "OMBA3",
 'B.1.623': "B1623", 
 'Omicron_BA.4': "OMBA4", 
 'Kappa_B.1.617.1': "K", 
 'Eta_B.1.525': "ET",
 'A.23.1': "19A", 
 'Theta_P.3': "T"}

# creating position list files for all types of variant combinations
variant_combinations = list(df_samples["variants"].unique())

!mkdir -p data_for_recombinant_read_analysis/poslists
for vc in variant_combinations:
    var1, var2 = vc.split(", ")
    with open("data_for_recombinant_read_analysis/poslists/" + "_".join([var_dict[v] for v in [var1, var2]]) + "_poslist.txt", "w") as f:
        all_pos = list(df_coinf_allmuts[df_coinf_allmuts["variant"] == var1].sort_values(by="pos")["pos"].unique())
        all_pos += list(df_coinf_allmuts[df_coinf_allmuts["variant"] == var2].sort_values(by="pos")["pos"].unique())
        all_pos = sorted(list(set(all_pos)))
        for p in all_pos:
            f.write("NC_045512.2\t"+str(p)+"\n")
            
# joining sample IDs with sample short names
df_download = df_download[df_download["submitted_ftp"].str.contains(".bam") | 
                          df_download["submitted_ftp"].str.contains("output.tar.gz")]
df_download["shortname"] = df_download["submitted_ftp"].apply(get_short_name_from_fname)
df_samples = pd.merge(df_samples, df_download[["sample_accession", "shortname"]], on = "sample_accession")

# reference fasta
refseq = "/v/projects/ebi-vcf-wfct0p/pipeko/refs/GCF_009858895.2_ASM985889v3_genomic.fna"

# running recombinant read detection for all samples
!mkdir -p data_for_recombinant_read_analysis/recombination_detection_results

for rid in list(df_samples["runid"]):
    # getting basic info about sample
    variant_combination, poslist_file, bamfile, sample_short = prepare_sample_details(rid)
        
    # generating pileup
    pileup_large = !samtools mpileup -B -f $refseq -l $poslist_file --output-QNAME -d 0 -q 30 -Q 30 $bamfile
    pileup_dict = {k.split("\t")[1]:k for k in pileup_large[2:]}
    
    # getting list of all relevant mutation pairs
    interesting_pairs_all = get_mutation_pairs(rid)
    
    # iterating through all possible mutation pairs
    var1_only = []
    var2_only = []
    recombinant = []
    nomut = []
    mutpair = []
    total_reads = []
    
    for p in interesting_pairs_all:        
        pos1 = p[0].split("_")[0]
        pos2 = p[1].split("_")[0]
        var1 = p[0].split("_")[-1]
        var2 = p[1].split("_")[-1]
        mut1 = p[0].split("_")[-2].upper()
        mut2 = p[1].split("_")[-2].upper()
        ref1 = p[0].split("_")[1].upper()
        ref2 = p[1].split("_")[1].upper()
        
        # getting their order right (based on genomic position)
        if pos1 > pos2:
            pos1, pos2 = pos2, pos1
            var1, var2 = var2, var1
            mut1, mut2 = mut2, mut1
            ref1, ref2 = ref2, ref1
        
        # correcting alternate allele string
        mut1 = correct_mutation_string(ref1, mut1)
        mut2 = correct_mutation_string(ref2, mut2)

        # parsing pileup line1
        if str(pos1) not in pileup_dict:
            continue
        k = pileup_dict[str(pos1)]
        b1 = k.split("\t")[4].upper()
        bases1 = get_base_list(b1)
        reads1 = k.split("\t")[6].split(",")

        # parsing pileup line2
        if str(pos2) not in pileup_dict:
            continue
        k = pileup_dict[str(pos2)]
        b2 = k.split("\t")[4].upper()
        bases2 = get_base_list(b2)
        reads2 = k.split("\t")[6].split(",")

        bases1 = np.array(bases1)
        reads1 = np.array(reads1)
        bases2 = np.array(bases2)
        reads2 = np.array(reads2)
        
        # getting overlapping reads
        read_list = list(set(reads1).intersection(set(reads2)))
        if len(read_list) == 0:
            continue

        v1_only = 0
        v2_only = 0
        recomb = 0
        no_v = 0
        for r in list(read_list):
            base1 = bases1[reads1 == r][0]
            base2 = bases2[reads2 == r][0]
            if base1.upper() == mut1 and base2.upper() == mut2:
                recomb += 1
            elif base1.upper() in [",", "."] and base2.upper() == mut2:
                v2_only += 1
            elif base2.upper() in [",", "."] and base1.upper() == mut1:
                v1_only += 1
            else:
                no_v += 1
      
        if var1 == var_dict[variant_combination.split(", ")[0]]:
            var1_only.append(v1_only)
            var2_only.append(v2_only)
        else:
            var1_only.append(v2_only)
            var2_only.append(v1_only)
        recombinant.append(recomb)
        nomut.append(no_v)
        mutpair.append(" | ".join(p))
        total_reads.append(len(read_list))

    df_tmp = pd.DataFrame()
    df_tmp["mutation_pair"] = mutpair
    df_tmp[var_dict[variant_combination.split(", ")[0]]+"_reads_only"] = var1_only
    df_tmp[var_dict[variant_combination.split(", ")[1]]+"_reads_only"] = var2_only
    df_tmp["recombinant_reads"] = recombinant
    df_tmp["reads_with_no_mutation"] = nomut
    df_tmp["number_of_overlapping_reads"] = total_reads
    df_tmp["recombinant_ratio"] = df_tmp["recombinant_reads"]/df_tmp["number_of_overlapping_reads"]
    df_tmp.to_csv("data_for_recombinant_read_analysis/recombination_detection_results/" + sample_short + ".csv", index=False)

4 Distribution of overlapping reads along the genome

As previously discussed, the number of reads that overlap mutually exclusive defining mutations of both comprising variants is largely influenced by the density of these mutations along the genome. Hereby we define the location of an overlapping read as the midpoint of the defining mutation pair it overlaps. Fig. 2a depicts the distribution of overlapping read locations with a 500 bp binning with a logarithmic vertical scale. Overlapping reads were counted together for all 118 analysed samples. The left panel of Fig. 2b shows the correlation between the number of overlapping reads at a given genomic range and the number of mutually exclusive defining mutations that fall into that same range. The right panel of Fig. 2b displays the same results for the density of mutually exclusive defining mutations and overlapping reads in specific genes.

Code 
# loading all files with results of recombinant read detection
all_res_files = !ls data_for_recombinant_read_analysis/recombination_detection_results/*.csv

# function to format pvalues
def format_pvalue(pval, digits=3):
    if round(pval, digits) != 0:
        return str(round(pval, digits)), 0
    else:
        i = 0
        while pval < 1:
            pval *= 10
            i += 1
        return str(round(pval, digits)), i
    
all_overlapping_midpoints = []
all_overlapping_genes = []

for fn in all_res_files:
    df_tmp = pd.read_csv(fn)
    for ri, r in df_tmp.iterrows():
        midpoint = np.mean([int(r["mutation_pair"].split(" | ")[0].split("_")[0]), 
                            int(r["mutation_pair"].split(" | ")[1].split("_")[0])])
        all_overlapping_midpoints += [midpoint]*r["number_of_overlapping_reads"]
        all_overlapping_genes += [get_gene_for_pos(midpoint)]*r["number_of_overlapping_reads"]
        
region_overlapping_reads, _ = np.histogram(all_overlapping_midpoints, range=(0,30000), bins=60)
region_defmuts, _ = np.histogram(df_coinf_allmuts.groupby(mutbase).agg({"pos": lambda x: list(x)[0]})["pos"], range=(0,30000), bins=60)

dens_defmuts = [gg[g]*1000/(gene_loc[g][1]-gene_loc[g][0]) for g in gene_loc.keys()]
olg = Counter(all_overlapping_genes)
dens_overlapping_reads = [olg[g]/(gene_loc[g][1]-gene_loc[g][0]) for g in gene_loc.keys()]

# figure to show overlapping reads along the genome        
fig = plt.figure(figsize=(15,3))
gs = fig.add_gridspec(2, hspace=0.5, height_ratios = [0.2, 0.8])
axs = gs.subplots(sharex=True)

yy = 2
for g in gene_loc.keys():
    axs[0].hlines(y = 0, xmin=gene_loc[g][0], xmax=gene_loc[g][1], lw=30, color = color_dict[g])
    if (yy == 2):
        va_c = "bottom"
        ha_c = "left"
    else:
        va_c = "top"
        ha_c = "right"
    axs[0].annotate(g, xy= (np.mean([gene_loc[g][0], gene_loc[g][1]]), yy*0.5), 
                 xytext = (np.mean([gene_loc[g][0], gene_loc[g][1]]), yy), rotation=30, ha="center", va = va_c, fontsize=8,
                arrowprops=dict(arrowstyle="-"))
    yy *= -1
axs[0].set_ylim(-1,1)
axs[0].set_xlim(0,30000)
axs[0].set_axis_off()

axs[1].hist(all_overlapping_midpoints, range=(0,30000), bins=60, 
            color="#2170D6")
axs[1].set_xlabel("Genomic position (binned by 500 bp)", fontsize=10)
axs[1].set_ylabel("Number of\n overlapping reads", fontsize=10)
axs[1].set_xlim(0,30000)
axs[1].set_ylim(1,10**8)
axs[1].grid(visible=False)
axs[1].tick_params(size=5, color="#666666")
axs[1].label_outer()
axs[1].set_yscale('log')

for g in gene_loc.keys():
    axs[1].vlines(x=gene_loc[g][1], ymin=0, ymax=10**8, color="black", alpha=0.5, lw=0.5)
axs[1].vlines(x=gene_loc["ORF1ab"][0], ymin=0, ymax=10**8, color="black", alpha=0.5, lw=0.5)    
    
plt.show()
plt.close()

# figure to show correlations

fig = plt.figure(figsize=(12,4))
gs = fig.add_gridspec(ncols=2, wspace=0.3)
axs = gs.subplots()

axs[0].scatter(np.array(region_overlapping_reads)/(10**6), region_defmuts, alpha=0.5, color="#2170D6")
axs[0].set_xlabel(r'Number of overlapping reads in 500 bp region ($10^6$)')
axs[0].set_ylabel("Number of mutually exclusive \ndefining mutations in 500 bp region")
axs[0].grid(visible=False)
axs[0].tick_params(size=5, color="#666666")

# printing correlation coefficient and p-value
_, xmax = axs[0].get_xbound()
ymin, _ = axs[0].get_ybound()
r,pvalue = pearsonr(region_defmuts, region_overlapping_reads)
pbase, i = format_pvalue(pvalue)
if i == 0:
    t = r'$\mathrm{R} = ' + str(round(r,3)) + '; \mathrm{p-value} = ' + pbase + '$'
else:
    t = r'$\mathrm{R} = ' + str(round(r,3)) + '; \mathrm{p-value} = ' + pbase + '\cdot 10^{-' + str(i) + '}$'

axs[0].annotate(t, xy = (xmax, ymin), ha="right", va="bottom", fontsize=12)

axs[1].scatter(dens_overlapping_reads, dens_defmuts, color="#2170D6")

# axs[1].set_xlim(-100,4000)
axs[1].set_ylim(-1,45)
axs[1].set_xlabel("Density of overlapping reads in gene (1/bp)")
axs[1].set_ylabel("Density of mutually exclusive \ndefining mutations in gene (1/kbp)")
axs[1].grid(visible=False)
axs[1].tick_params(size=5, color="#666666")

xmin, xmax = axs[1].get_xbound()
ymin, ymax = axs[1].get_ybound()

offset_x = (xmax-xmin)/80
offset_y = (ymax-ymin)/80

for i,g in enumerate(gene_loc.keys()):
    if g in ["ORF6", "S", "ORF7b", "ORF8", "N"]:
        axs[1].annotate(g, xy = (dens_overlapping_reads[i], dens_defmuts[i]),
                       xytext = (dens_overlapping_reads[i]+offset_x, dens_defmuts[i]+offset_y))

axs[1].annotate("ORF1ab", xy=([dens_overlapping_reads[i] for i, g in enumerate(gene_loc.keys()) if g == "ORF1ab"][0], 
                           [dens_defmuts[i] for i, g in enumerate(gene_loc.keys()) if g == "ORF1ab"][0]), 
             xytext=(200, 7.5), arrowprops=dict(arrowstyle="-", color="#2170D6"),
            ha = "left", verticalalignment = "center")
axs[1].annotate("ORF3a", xy=([dens_overlapping_reads[i] for i, g in enumerate(gene_loc.keys()) if g == "ORF3a"][0], 
                           [dens_defmuts[i] for i, g in enumerate(gene_loc.keys()) if g == "ORF3a"][0]), 
            xytext=(100, 19), arrowprops=dict(arrowstyle="-", color="#2170D6"),
            ha = "left", verticalalignment = "center")
axs[1].annotate("E", xy=([dens_overlapping_reads[i] for i, g in enumerate(gene_loc.keys()) if g == "E"][0], 
                           [dens_defmuts[i] for i, g in enumerate(gene_loc.keys()) if g == "E"][0]), 
            xytext=(0, 17), arrowprops=dict(arrowstyle="-", color="#2170D6"),
            ha = "left", verticalalignment = "center")
axs[1].annotate("M", xy=([dens_overlapping_reads[i] for i, g in enumerate(gene_loc.keys()) if g == "M"][0], 
                           [dens_defmuts[i] for i, g in enumerate(gene_loc.keys()) if g == "M"][0]), 
            xytext=(300, 15), arrowprops=dict(arrowstyle="-", color="#2170D6"),
            ha = "left", verticalalignment = "center")
axs[1].annotate("ORF7a", xy=([dens_overlapping_reads[i] for i, g in enumerate(gene_loc.keys()) if g == "ORF7a"][0], 
                           [dens_defmuts[i] for i, g in enumerate(gene_loc.keys()) if g == "ORF7a"][0]), 
            xytext=(400, 5), arrowprops=dict(arrowstyle="-", color="#2170D6"),
            ha = "left", verticalalignment = "center")
axs[1].annotate("ORF10", xy=([dens_overlapping_reads[i] for i, g in enumerate(gene_loc.keys()) if g == "ORF10"][0], 
                           [dens_defmuts[i] for i, g in enumerate(gene_loc.keys()) if g == "ORF10"][0]), 
            xytext=(150, 10), arrowprops=dict(arrowstyle="-", color="#2170D6"),
            ha = "left", verticalalignment = "center")

r,pvalue = pearsonr(dens_defmuts, dens_overlapping_reads)
pbase, i = format_pvalue(pvalue)
if i == 0:
    t = r'$\mathrm{R} = ' + str(round(r,3)) + '; \mathrm{p-value} = ' + pbase + '$'
else:
    t = r'$\mathrm{R} = ' + str(round(r,3)) + '; \mathrm{p-value} = ' + pbase + '\cdot 10^{-' + str(i) + '}$'

axs[1].annotate(t, xy = (xmax, ymin), ha="right", va="bottom", fontsize=12)

plt.show()
plt.close()

(a) Distribution of overlapping read locations along the SARS-CoV-2 genome, binned by 500 bp. All overlapping reads of the 118 analysed samples were considered for this figure. Overlapping read location is defined as the midpoint of mutually exclusive defining mutations it overlaps. The top panel shows the locations of specific genes along the genome. Vertical lines on the bottom panel indicate gene boundaries. Note that the scale of the vertical axis is logarithmic.

(b) Left: Correlation of the number of mutually exclusive defining mutations and the number of overlapping reads in 500 bp genomic regions. Right: Correlation of the density of mutually exclusive defining mutations and the density of overlapping reads in specific genes.

Figure 2: Distribution of the locations of reads overlapping mutually exclusive defining mutations of both comprising variants.

5 Traces of recombination in overlapping reads

5.1 Distribution of recombination breakpoints along the genome

We define the location of a recombination breakpoint as the midpoint of the recombination range. The recombination range is the region of the genome bracketed by two mutually exclusive defining mutations (one of each parental strain) that are overlapped by at least one read in which both mutations are present. The distribution of recombination breakpoints is displayed in Fig. 3a for all 118 analysed co-infection samples. In samples where multiple reads support a single breakpoint, the breakpoint is counted only once.

It has recently been shown that recombination breakpoints occur non-uniformly across the viral genome (Turakhia et al., 2022). However, the computation approach (named “RIPPLES”) providing this finding takes consensus sequences as its input. In order to compare its results with our read-level data, we downloaded Supplementary Table 1 of the above manuscript and calculated the correlation between the number of breakpoints detected by RIPPLES and the number of breakpoints supported by overlapping reads with a 500 bp binning along the genome. The results of this comparison are shown in Fig. 3b.

Code 
# collecting read-level data
sample_breakpoint_midpoints = []
sample_breakpoint_genes = []

for fn in all_res_files:
    df_tmp = pd.read_csv(fn)
    for ri, r in df_tmp.iterrows():
        if r["recombinant_reads"] == 0:
            continue
        midpoint = np.mean([int(r["mutation_pair"].split(" | ")[0].split("_")[0]), 
                            int(r["mutation_pair"].split(" | ")[1].split("_")[0])])
        sample_breakpoint_midpoints += [midpoint]
        sample_breakpoint_genes += [get_gene_for_pos(midpoint)]

# collecting data obtained with RIPPLES
# !wget https://static-content.springer.com/esm/art%3A10.1038%2Fs41586-022-05189-9/MediaObjects/41586_2022_5189_MOESM4_ESM.xlsx

df_paper = pd.read_excel("41586_2022_5189_MOESM4_ESM.xlsx")
paper_breakpoints_all = list(df_paper["breakpoint_1"]) + list(df_paper["breakpoint_2"])
paper_breakpoint_midpoints = []
paper_breakpoint_genes = []

for bp in paper_breakpoints_all:
    if not pd.isnull(bp):
        bp1 = int(bp.replace("(","").replace(")","").split(",")[0])
        bp2 = int(bp.replace("(","").replace(")","").split(",")[1])
        paper_breakpoint_midpoints.append((bp1+bp2)/2)
        paper_breakpoint_genes.append(get_gene_for_pos((bp1+bp2)/2))

# calculating correlations
paper_breakpoint_midpoints_h, _ = np.histogram(paper_breakpoint_midpoints, range=(0,30000), bins=60)
sample_breakpoint_midpoints_h, _ = np.histogram(sample_breakpoint_midpoints, range=(0,30000), bins=60)

sample_breakpoint_genes_h = Counter(sample_breakpoint_genes)
paper_breakpoint_genes_h = Counter(paper_breakpoint_genes)

dens_sample_breakpoint = [sample_breakpoint_genes_h[g]/(gene_loc[g][1]-gene_loc[g][0]) for g in gene_loc.keys()]
dens_paper_breakpoint = [paper_breakpoint_genes_h[g]/(gene_loc[g][1]-gene_loc[g][0]) for g in gene_loc.keys()]


# figure to plot both
gs_top = plt.GridSpec(3, 1, height_ratios=[0.1, 0.45, 0.45], top=0.78)
gs_base = plt.GridSpec(3, 1, hspace=0)
fig = plt.figure(figsize=(15,5))

# top panel showing genes
topax = fig.add_subplot(gs_top[0,:])
yy = 2
for g in gene_loc.keys():
    topax.hlines(y = 0, xmin=gene_loc[g][0], xmax=gene_loc[g][1], lw=30, color = color_dict[g])
    if (yy == 2):
        va_c = "bottom"
        ha_c = "left"
    else:
        va_c = "top"
        ha_c = "right"
    topax.annotate(g, xy= (np.mean([gene_loc[g][0], gene_loc[g][1]]), yy*0.5), 
                 xytext = (np.mean([gene_loc[g][0], gene_loc[g][1]]), yy), rotation=30, ha="center", va = va_c, fontsize=8,
                arrowprops=dict(arrowstyle="-"))
    yy *= -1
topax.set_ylim(-1,1)
topax.set_xlim(0,30000)
topax.set_axis_off()

# bottom panel showing overlapping read numbers
ax1 = fig.add_subplot(gs_base[1,:])
ax2 = fig.add_subplot(gs_base[2,:], sharex=ax1)

ax1.hist(sample_breakpoint_midpoints, range=(0,30000), bins=60, 
            color="#5AA1FF")
ax1.set_xlabel("Genomic position (binned by 500 bp)", fontsize=10)
# ax1.set_ylabel("Number of\n recombination breakpoints", fontsize=10)
ax1.set_xlim(0,30000)
ax1.set_ylim(0,400)
ax1.grid(visible=False)
ax1.tick_params(size=5, color="#666666")
ax1.label_outer()

ax1.annotate("Number of recombination breakpoints", xy = (-1500, 0), fontsize=10, rotation=90, ha="left", va="center", 
            annotation_clip=False)

ymin, ymax = ax1.get_ybound()
ax1.annotate("Based on read-level data", xy=(600, ymax*0.9), fontsize=10, ha="left", va="top")
for g in gene_loc.keys():
    ax1.vlines(x=gene_loc[g][1], ymin=ymin, ymax=ymax, color="black", alpha=0.5, lw=0.5)
ax1.vlines(x=gene_loc["ORF1ab"][0], ymin=ymin, ymax=ymax, color="black", alpha=0.5, lw=0.5)    

ax2.hist(paper_breakpoint_midpoints, range=(0,30000), bins=60, 
            color="#98C4FF")
ax2.set_xlabel("Genomic position (binned by 500 bp)", fontsize=10)
# ax2.set_ylabel("Number of\n recombination breakpoints", fontsize=10)
ax2.set_xlim(0,30000)
ax2.set_ylim(0,50)
ax2.grid(visible=False)
ax2.tick_params(size=5, color="#666666")
ax2.label_outer()

ymin, ymax = ax2.get_ybound()
ax2.annotate("Found by RIPPLES (Turakhia et al., 2022)", xy=(600, ymax*0.9), fontsize=10, ha="left", va="top")
for g in gene_loc.keys():
    ax2.vlines(x=gene_loc[g][1], ymin=ymin, ymax=ymax, color="black", alpha=0.5, lw=0.5)
ax2.vlines(x=gene_loc["ORF1ab"][0], ymin=ymin, ymax=ymax, color="black", alpha=0.5, lw=0.5)    

plt.show()
plt.close()

# figure showing correlations
fig = plt.figure(figsize=(12,4))
gs = fig.add_gridspec(ncols=2, wspace=0.3)
axs = gs.subplots()

axs[0].scatter(np.array(sample_breakpoint_midpoints_h), paper_breakpoint_midpoints_h, alpha=0.5, color="#5AA1FF")
axs[0].set_xlabel('Number of recombinant breakpoints\n (read-level analysis)')
axs[0].set_ylabel('Number of recombinant breakpoints\n (obtained by RIPPLES (Turakhia et al., 2022))')
axs[0].grid(visible=False)
axs[0].tick_params(size=5, color="#666666")

# printing correlation coefficient and p-value
_, xmax = axs[0].get_xbound()
ymin, _ = axs[0].get_ybound()
r,pvalue = pearsonr(sample_breakpoint_midpoints_h, paper_breakpoint_midpoints_h)
pbase, i = format_pvalue(pvalue)
if i == 0:
    t = r'$\mathrm{R} = ' + str(round(r,3)) + '; \mathrm{p-value} = ' + pbase + '$'
else:
    t = r'$\mathrm{R} = ' + str(round(r,3)) + '; \mathrm{p-value} = ' + pbase + '\cdot 10^{-' + str(i) + '}$'

axs[0].annotate(t, xy = (xmax, ymin), ha="right", va="bottom", fontsize=12)

axs[1].scatter(dens_sample_breakpoint, dens_paper_breakpoint, color="#5AA1FF")

axs[1].set_xlim(-0.02,0.4)
axs[1].set_ylim(-0.01,0.1)
axs[1].set_xlabel("Density of recombination breakpoints (1/bp)\n (read-level analysis)")
axs[1].set_ylabel("Density of recombination breakpoints (1/bp)\n (obtained by RIPPLES (Turakhia et al., 2022))")
axs[1].grid(visible=False)
axs[1].tick_params(size=5, color="#666666")

xmin, xmax = axs[1].get_xbound()
ymin, ymax = axs[1].get_ybound()

offset_x = (xmax-xmin)/80
offset_y = (ymax-ymin)/80

for i,g in enumerate(gene_loc.keys()):
    if g != "ORF7a":
        axs[1].annotate(g, xy = (dens_sample_breakpoint[i], dens_paper_breakpoint[i]),
                       xytext = (dens_sample_breakpoint[i]+offset_x, dens_paper_breakpoint[i]+offset_y))
    else:
        axs[1].annotate(g, xy = (dens_sample_breakpoint[i], dens_paper_breakpoint[i]),
               xytext = (dens_sample_breakpoint[i]+offset_x, dens_paper_breakpoint[i]-offset_y), va="top")

r,pvalue = pearsonr(dens_sample_breakpoint, dens_paper_breakpoint)
pbase, i = format_pvalue(pvalue)
if i == 0:
    t = r'$\mathrm{R} = ' + str(round(r,3)) + '; \mathrm{p-value} = ' + pbase + '$'
else:
    t = r'$\mathrm{R} = ' + str(round(r,3)) + '; \mathrm{p-value} = ' + pbase + '\cdot 10^{-' + str(i) + '}$'

axs[1].annotate(t, xy = (xmax, ymin), ha="right", va="bottom", fontsize=12)

plt.show()
plt.close()

(a) Upper panel: Number of recombination breakpoints in 500 bp-wide regions of the SARS-CoV-2 genome based on overlapping reads containing multiple defining mutations of the parental strains. Lower panel: Distribution of recombination breakpoints along the genome with a 500 bp binning found by RIPPLES, based on Turakhia et al., 2022.

(b) Left: Relationship between the number of recombination breakpoints indicated by overlapping reads and the number of breakpoints identified by Thurakia et al. 2022 (with the RIPPLES software) in consensus sequences. The genome was binned by 500 bp regions for this analysis. Right: Correlation of the density of recombination breakpoints indicated by overlapping reads and the density of breakpoints identified by Thurakia et al. 2022 in specific genes.

Figure 3: Distribution of recombination breakpoints along the genome.

The above results indicate that the findings of RIPPLES and the results of the read-level analysis are moderately correlated, thus suggesting that both are largely influenced by the inherent distribution of defining mutations along the genome.

5.2 Traces of recombination in artificial mixtures

In theory, artificial mixtures of purified RNA from different viral strains could serve as reliable controls for recombinant detection, as no recombination is expected to occur once viral replication has been terminated. In our set of samples selected for read-level analysis, 7 samples of Bal et al., 2022 (study ID PRJNA817870) and 5 samples of Sovic et al., 2022 (study ID PRJNA827817) are in fact artificial mixtures, thus it is of interest to examine the distribution of recombinant reads along their genomes.

Code 
# getting recombinant read ratios

all_pos = [int(k) for k in np.linspace(1,30000,30000)]
ratios_dict = {"artificial_Delta_Omicron": [],
              "artificial_other": [],
              "non_artificial_Delta_Omicron": [],
              "non_artificial_other": []}

for f in all_res_files:
    sn = f.split("/")[-1].split(".")[0]
    vt = df_samples[df_samples["shortname"] == sn].iloc[0]["variants"]
    si = df_samples[df_samples["shortname"] == sn].iloc[0]["study_accession"]
    
    df_tmp = pd.read_csv(f)
    num_of_overlapping_reads = {k:0 for k in all_pos}
    num_of_recombinant_reads = {k:0 for k in all_pos}

    for ri, r in df_tmp.iterrows():
        posmin = int(r["mutation_pair"].split(" | ")[0].split("_")[0]) 
        posmax = int(r["mutation_pair"].split(" | ")[1].split("_")[0])
        for pi in range(posmin, posmax+1):
            num_of_overlapping_reads[pi] += r["number_of_overlapping_reads"]
            num_of_recombinant_reads[pi] += r["recombinant_reads"]

    counts_overlapping_reads = np.array(list(num_of_overlapping_reads.values()))
    id_nonzero = counts_overlapping_reads > 0
    ratio_recomb_reads = np.zeros_like(all_pos).astype(float)
    ratio_recomb_reads[id_nonzero] += np.array(list(num_of_recombinant_reads.values()))[id_nonzero]/counts_overlapping_reads[id_nonzero]
    
    if vt == "Delta_B.1.617.2, Omicron_BA.1" and si in ["PRJNA817870", "PRJNA827817"]:
        ratios_dict["artificial_Delta_Omicron"] += list(ratio_recomb_reads)
    elif vt != "Delta_B.1.617.2, Omicron_BA.1" and si in ["PRJNA817870", "PRJNA827817"]:
        ratios_dict["artificial_other"] += list(ratio_recomb_reads)       
    elif vt == "Delta_B.1.617.2, Omicron_BA.1" and si not in ["PRJNA817870", "PRJNA827817"]:
        ratios_dict["non_artificial_Delta_Omicron"] += list(ratio_recomb_reads)
    elif vt != "Delta_B.1.617.2, Omicron_BA.1" and si not in ["PRJNA817870", "PRJNA827817"]:
        ratios_dict["non_artificial_other"] += list(ratio_recomb_reads)
        
hist_dict = {}
for k,v in ratios_dict.items():
    if k == "artificial_Delta_Omicron":
        s = df_samples[(df_samples["variants"] == "Delta_B.1.617.2, Omicron_BA.1") & 
                       (df_samples["study_accession"].isin(["PRJNA817870", "PRJNA827817"]))].shape[0]
    elif k == "artificial_other":
        s = df_samples[(df_samples["variants"] != "Delta_B.1.617.2, Omicron_BA.1") & 
                       (df_samples["study_accession"].isin(["PRJNA817870", "PRJNA827817"]))].shape[0]        
    elif k == "non_artificial_Delta_Omicron":
        s = df_samples[(df_samples["variants"] == "Delta_B.1.617.2, Omicron_BA.1") & 
                       (~df_samples["study_accession"].isin(["PRJNA817870", "PRJNA827817"]))].shape[0] 
    elif k == "non_artificial_other":
        s = df_samples[(df_samples["variants"] != "Delta_B.1.617.2, Omicron_BA.1") & 
                       (~df_samples["study_accession"].isin(["PRJNA817870", "PRJNA827817"]))].shape[0]    

    d = np.histogram(v, bins=100, range=(0,1))
    hist_dict[k] = d[0]*100/(s*30000)
    # leaving out zero ratios
    hist_dict[k] = hist_dict[k][1:]
bin_middle = 0.5*(d[1][1:] + d[1][:-1])    
bin_middle = bin_middle[1:]

# plotting figure

fig = plt.figure(figsize=(7,6))
gs = fig.add_gridspec(nrows=4, hspace=0)
axs = gs.subplots()

axs[0].bar(bin_middle, np.sum(hist_dict["artificial_Delta_Omicron"]) - np.cumsum(hist_dict["artificial_Delta_Omicron"]), 
           width=0.007, color="#c1121f", 
            label = "Delta - Omicron (BA.1) artificial samples")


axs[1].bar(bin_middle, np.sum(hist_dict["artificial_other"]) - np.cumsum(hist_dict["artificial_other"]), 
           width=0.007, color="#F3515B", 
            label = "other artificial samples")



axs[2].bar(bin_middle, np.sum(hist_dict["non_artificial_Delta_Omicron"]) - np.cumsum(hist_dict["non_artificial_Delta_Omicron"]), 
           width=0.007, color="#023e8a", 
            label = "Delta - Omicron (BA.1) real samples")


axs[3].bar(bin_middle, np.sum(hist_dict["non_artificial_other"]) - np.cumsum(hist_dict["non_artificial_other"]), 
           width=0.007, color="#2B74D4", 
            label = "other real samples")

for axid, ax in enumerate(axs):
    ax.grid(visible=True)
    ax.tick_params(size=5, color="#666666")
    ax.legend()
    ax.label_outer()
    ymin, ymax = ax.get_ybound()
    ax.vlines(x=0.1, ymin=ymin, ymax=ymax, color="black", lw = 2)  
    vals = [1,2,3]
    if axid == 3: vals = [0.05, 0.10, 0.15]
    ax.yaxis.set_major_locator(mticker.FixedLocator(vals))
    ax.set_yticklabels(['{:.2%}'.format(x/100) for x in vals])

axs[3].set_xlabel("T", fontsize=12)
axs[0].annotate("Mean percentage of genomic positions \nwith recombinant read ratio > T", 
                xy = (-0.16, -1), fontsize=12, rotation=90, 
                      ha="center", va="center", xycoords = "axes fraction")        

plt.show()
plt.close()

Figure 4: The average percentage of genomic positions (per sample) for which the ratio of recombinant reads out of all overlapping ones reaches T. Genomic positions with exactly zero recombinant reads are not shown. Samples were categorized into groups of Delta – Omicron (BA.1) artificial/real and non-Delta – Omicron (BA.1) artificial/real samples. The vertical black line indicates T = 0.1.

Artificial mixture samples do contain a substantial number of recombinant reads, but the prevalence of genomic positions overlapped by a recombinant read ratio of more than 0.1 (vertical black lines) is much lower than in true co-infection samples. This result suggests that breakpoints supported by no more than 10% of the overlapping reads might be considered artefacts due to chimera formation during PCR.

This is further supported by the fact that when calculating the ratio of duplicate reads among recombinant reads in positions with a recombinant read ratio of lower vs. higher than 10%, genomic positions in which recombination was supported by less than 10% of overlapping reads show very high fractions of duplicates, indicating possible evidence of PCR artefacts.

Code 
# getting duplicate read ratios
ratios_dict2 = {"under10percent": [],
              "above10percent": []}

for f in all_res_files:
    sn = f.split("/")[-1].split(".")[0]
    if not os.path.isfile("data_for_recombinant_read_analysis/recombination_sgRNA/"+sn+".csv"):
        continue
        
    num_of_recombinant_reads = {p:0 for p in all_pos}
    num_of_dup_recombinant_reads = {p:0 for p in all_pos}

    df_tmp = pd.read_csv("data_for_recombinant_read_analysis/recombination_sgRNA/"+sn+".csv")
    df_multireads = df_tmp.groupby("mutation_pair", as_index=False).agg({"readID": lambda x: list(x), 
                                                               "is_leader": lambda x: len(list(x))})
    df_multireads = df_multireads[df_multireads["is_leader"]>1]
    all_readids = list(set(", ".join([", ".join(k) for k in df_multireads["readID"]]).split(", ")))
    
    filename = sn+"_readIDs_tocheck.txt"
    bamname = "data_for_recombinant_read_analysis/indexed_bams/" + sn + "_sorted.bam"
    
    with open(filename, "w") as ff:
        for readID in all_readids:
            ff.write(readID+"\n")
            
    k = !samtools view --qname-file $filename $bamname | cut -f1,4
    readID = [ki.split("\t")[0] for ki in k]
    readpos = np.array([int(ki.split("\t")[1]) for ki in k])
    
    is_dup = []
    for rrid in df_tmp["readID"]:
        if rrid in readID and (readpos==readpos[readID.index(rrid)]).sum() > 1:
            is_dup.append(True)
        else:
            is_dup.append(False)

    df_tmp["is_duplicate"] = is_dup
    for rrix, rr in df_tmp.iterrows():
        posmin = int(rr["mutation_pair"].split(" | ")[0].split("_")[0]) 
        posmax = int(rr["mutation_pair"].split(" | ")[1].split("_")[0])
        for pi in range(posmin, posmax+1):
            num_of_recombinant_reads[pi] += 1
            if rr["is_duplicate"]:
                num_of_dup_recombinant_reads[pi] += 1
    
    counts_recombinant_reads = np.array(list(num_of_recombinant_reads.values()))
    id_nonzero = counts_recombinant_reads > 0
    ratio_dup_recomb_reads = np.zeros_like(all_pos).astype(float)
    ratio_dup_recomb_reads[id_nonzero] += np.array(list(num_of_dup_recombinant_reads.values()))[id_nonzero]/counts_recombinant_reads[id_nonzero]

    df_tmp = pd.read_csv(f)
    num_of_overlapping_reads = {k:0 for k in all_pos}
    num_of_recombinant_reads2 = {k:0 for k in all_pos}

    for rrix, rr in df_tmp.iterrows():
        posmin = int(rr["mutation_pair"].split(" | ")[0].split("_")[0]) 
        posmax = int(rr["mutation_pair"].split(" | ")[1].split("_")[0])
        for pi in range(posmin, posmax+1):
            num_of_overlapping_reads[pi] += rr["number_of_overlapping_reads"]
            num_of_recombinant_reads2[pi] += rr["recombinant_reads"]

    counts_overlapping_reads = np.array(list(num_of_overlapping_reads.values()))
    id_nonzero = counts_overlapping_reads > 0
    ratio_recomb_reads = np.zeros_like(all_pos).astype(float)
    ratio_recomb_reads[id_nonzero] += np.array(list(num_of_recombinant_reads2.values()))[id_nonzero]/counts_overlapping_reads[id_nonzero]
    
    ratios_dict2["under10percent"] += list(ratio_dup_recomb_reads[(ratio_recomb_reads < 0.1)*(ratio_recomb_reads > 0)*(counts_recombinant_reads>1)])
    ratios_dict2["above10percent"] += list(ratio_dup_recomb_reads[(ratio_recomb_reads >= 0.1)*(counts_recombinant_reads>1)])

# plotting distributions    
fig, ax = plt.subplots()
parts = ax.violinplot([ratios_dict2["under10percent"], ratios_dict2["above10percent"]], showextrema=False)
parts['bodies'][0].set_facecolor("#c1121f")
parts['bodies'][1].set_facecolor("#023e8a")
parts['bodies'][0].set_alpha(0.7)
parts['bodies'][1].set_alpha(0.7)
ax.boxplot([ratios_dict2["under10percent"], ratios_dict2["above10percent"]], widths = 0.03, 
           patch_artist=True, boxprops=dict(facecolor="white"))
ax.set_ylabel("Ratio of duplicate reads \namong recombinant reads", fontsize=12)
ax.set_xticklabels([r'$0 < T < 0.1$', r'$T \geq 0.1$'], fontsize=12)
plt.show()

Figure 5: The ratio of duplicate reads among recombinant reads in genomic positions with a recombinant read ratio less than vs. at least 10%. Genomic positions with no recombinant reads were excluded from this figure.

5.3 Common recombination breakpoint ranges

To check whether any recombination breakpoints systematically occur in our investigated samples, we plotted the ratio of samples in which sufficient evidence of a breakpoint (supported by at least 10 reads and a recombination read ratio of more than 0.1) can be uncovered for the given genomic position. Co-infection samples of the 4 most common variant combinations were considered for this figure. Artificial samples were discarded during this analysis.

Code 
# count ratio of samples in which the given position is indicated to be in a breakpoint range
all_pos = [int(k) for k in np.linspace(1,30000,30000)]

recomb_ranges_in_samples = dict()
art_samples = list(df_samples[df_samples["study_accession"].isin(["PRJNA817870", "PRJNA827817"])]["runid"])

for fn in all_res_files:
    rid = df_samples[df_samples["shortname"] == fn.split("/")[-1].split(".")[0]].iloc[0]["runid"]
    if rid in art_samples:
        continue
    var_comb = get_sorted_variants(df_coinf_allmuts[df_coinf_allmuts["runid"] == rid].iloc[0]["variants"])
    df_tmp = pd.read_csv(fn)
    
    sample_ratio = 1/((df_samples["variants"] == var_comb) & (~df_samples["runid"].isin(art_samples))).sum()
    
    if var_comb not in recomb_ranges_in_samples:
        recomb_ranges_in_samples[var_comb] = {k: 0 for k in all_pos}
    
    pos_considered = []
    
    for ri, r in df_tmp.iterrows():
        if r["recombinant_reads"] < 10 or r["recombinant_reads"]/r["number_of_overlapping_reads"] < 0.1:
            continue
        posmin = int(r["mutation_pair"].split(" | ")[0].split("_")[0]) 
        posmax = int(r["mutation_pair"].split(" | ")[1].split("_")[0])
        for pi in range(posmin, posmax+1):
            if pi not in pos_considered:
                recomb_ranges_in_samples[var_comb][pi] += sample_ratio
            pos_considered.append(pi)


# figure to show ratio of samples for common variant combinations        
bluecol = "#023e8a"
variants_to_plot = df_samples["variants"].value_counts()[:4].keys().to_list()

# figure to plot both
gs_top = plt.GridSpec(5, 1, height_ratios=[0.1, 0.225, 0.225, 0.225, 0.225], top=0.9)
gs_base = plt.GridSpec(5, 1, hspace=0)
fig = plt.figure(figsize=(15,5))

# top panel showing genes
topax = fig.add_subplot(gs_top[0,:])
yy = 2
for g in gene_loc.keys():
    topax.hlines(y = 0, xmin=gene_loc[g][0], xmax=gene_loc[g][1], lw=30, color = color_dict[g])
    if (yy == 2):
        va_c = "bottom"
        ha_c = "left"
    else:
        va_c = "top"
        ha_c = "right"
    topax.annotate(g, xy= (np.mean([gene_loc[g][0], gene_loc[g][1]]), yy*0.5), 
                 xytext = (np.mean([gene_loc[g][0], gene_loc[g][1]]), yy), rotation=30, ha="center", va = va_c, fontsize=8,
                arrowprops=dict(arrowstyle="-"))
    yy *= -1
topax.set_ylim(-1,1)
topax.set_xlim(0,30000)
topax.set_axis_off()

# bottom panel showing overlapping read numbers
ax1 = fig.add_subplot(gs_base[1,:])
ax2 = fig.add_subplot(gs_base[2,:], sharex=ax1)
ax3 = fig.add_subplot(gs_base[3,:], sharex=ax1)
ax4 = fig.add_subplot(gs_base[4,:], sharex=ax1)


ax4.set_xlabel("Genomic position", fontsize=10)

for ax in (ax1, ax2, ax3, ax4):    
    ax.set_xlim(0,30000)
    ax.set_ylim(0,0.3)
    ax.grid(visible=False)
    ax.label_outer()
    ax.set_yticks([0, 0.1, 0.2])
    ax.tick_params(size=5, color="#666666")
    for g in gene_loc.keys():
        ax.vlines(x=gene_loc[g][1], ymin=0, ymax=10**8, color="black", alpha=0.5, lw=0.5)
    ax.vlines(x=gene_loc["ORF1ab"][0], ymin=0, ymax=10**8, color="black", alpha=0.5, lw=0.5)    

ax1.annotate("Ratio of samples\n with evidence of breakpoint", 
                xy = (-0.05, -1), fontsize=10, rotation=90, 
                      ha="center", va="center", xycoords = "axes fraction")        

snum = ((df_samples["variants"]==variants_to_plot[0]) & (~df_samples["runid"].isin(art_samples))).sum()
ax1.annotate("Delta (B.1.617.2) - Omicron (BA.1) samples (N = "+str(snum)+")", 
             xy = (500, 0.05), fontsize=10, ha="left", va="center", 
            annotation_clip=False)
snum = ((df_samples["variants"]==variants_to_plot[1]) & (~df_samples["runid"].isin(art_samples))).sum()
ax2.annotate("Alpha (B.1.1.7) - Iota (BA.1) samples (N = "+str(snum)+")", 
             xy = (500, 0.05), fontsize=10, ha="left", va="center", 
            annotation_clip=False)
snum = ((df_samples["variants"]==variants_to_plot[2]) & (~df_samples["runid"].isin(art_samples))).sum()
ax3.annotate("Alpha (B.1.1.7) - Epsilon (B.1.427-429) samples (N = "+str(snum)+")", 
             xy = (500, 0.05), fontsize=10, ha="left", va="center", 
            annotation_clip=False)
snum = ((df_samples["variants"]==variants_to_plot[3]) & (~df_samples["runid"].isin(art_samples))).sum()
ax4.annotate("Alpha (B.1.1.7) - Delta (B.1.617.2) samples (N = "+str(snum)+")", 
             xy = (500, 0.05), fontsize=10, ha="left", va="center", 
            annotation_clip=False)

ax4.tick_params(size=5, color="#666666")

all_axs = [ax1, ax2, ax3, ax4]

for varc, dd in recomb_ranges_in_samples.items():
    if varc not in variants_to_plot:
        continue
    axc = all_axs[variants_to_plot.index(varc)]
    axc.plot(dd.keys(), 
             dd.values(),
             color="#2170D6")
    
# adding GISAID deltacron samples breakpoints
redcol = "#9d0208"
ax1.fill_between([210, 6512], 0, 0.6, alpha=0.1, color=redcol)
ax1.fill_between([8393, 10449], 0, 0.6, alpha=0.1, color=redcol)
ax1.fill_between([22028, 22193], 0, 0.6, alpha=0.1, color=redcol)
ax1.fill_between([25000, 25584], 0, 0.6, alpha=0.1, color=redcol)


# breakpoint ranges in GISAID DeltaCron samples (see Supplementary File Z)
ax1.annotate("XF", 
             xy = (350, 0.25), fontsize=10, ha="left", va="center", color = redcol)
ax1.annotate("XS", 
             xy = (8533, 0.25), fontsize=10, ha="left", va="center", color = redcol)
ax1.annotate("XD", 
             xy = (21908, 0.25), fontsize=10, ha="left", va="center", color = redcol)
ax1.annotate("XD", 
             xy = (25040, 0.25), fontsize=10, ha="left", va="center", color = redcol)

ax1.fill_between([22578, 23202], 0, 0.6, alpha=0.1, color=bluecol)
ax1.fill_between([23525, 23854], 0, 0.6, alpha=0.1, color=bluecol)
ax1.fill_between([24130, 24503], 0, 0.6, alpha=0.1, color=bluecol)
ax1.fill_between([26530, 26767], 0, 0.6, alpha=0.1, color=bluecol)

plt.show()
plt.close()

Figure 6: The ratio of various co-infection samples with sufficient evidence of a recombination breakpoint in different genomic positions. For this analysis, co-infection samples were categorized based on their variant composition. Only those variant compositions are shown for which at least 10 samples were available in the raw read analysis pipeline. For each genomic position, the number of samples was calculated in which the given position is part of a recombination breakpoint range based on the evidence of at least 10 short reads and a recombinant read ratio of 0.1 or larger. For each sample, each genomic position was counted a single time only. Regions shaded with light red are recombination breakpoint ranges of the three Pangolin lineages (XF, XS, and XD) in which Delta – Omicron recombination occurs. Areas shaded with light blue indicate intragenic hotspots. Artificial samples were not included in this analysis.

JupyterRequireError: notebook/js/codecell: Timeout. Library 'notebook/js/codecell' is not loaded.

Many of the genomic ranges indicated as putative recombination breakpoint ranges in multiple samples coincide with gene boundaries. Additional intragenic hotspots were regions 22578-23202, 23525-23854, and 24130-24503 in gene S and 26530-26767 in gene M in co-infection samples of Delta – Omicron (BA.1) variants (shaded with light blue above). The majority of recombination hotspots detected from short reads do not correspond to regions of recombination identified from clonal recombinants of the GISAID database (samples assigned to Pango lineages XF, XS and XD, shaded with light red above; for more details, see Supplementary Figure 3).

5.4 Recombination in subgenomic RNA

In order to check if reads carrying traces of recombination originate from subgenomic RNA sequences, we first collected read IDs that were detected as recombinants for each sample, for each relevant pair of mutually exclusive defining mutations and then queried the original BAM files for their details. We checked whether these reads contained the nucleotides of the common 5’-leader sequence attached to sgRNAs during translation and if they were soft-clipped during alignment. If any of these conditions were met, the given read was considered to show signs of sgRNA-origin.

Code 
def get_overlapping_reads_for_sample(s, p, pileup_dict):
    firstone = True
    for pos in p:
        if str(pos) not in pileup_dict:
            overlapping_reads = {}
            continue
        k = pileup_dict[str(pos)]
        reads = k.split("\t")[6].split(",")
        if firstone:
            overlapping_reads = set(reads)
            firstone = False
        else:
            overlapping_reads = overlapping_reads.intersection(set(reads))
    return list(overlapping_reads)

def get_recombinant_reads_for_sample(s, base_comb_muts, overlapping_reads, pileup_dict):
    recombinant_reads = []
    for r in list(overlapping_reads):
        rt = ""
        for m in sorted(base_comb_muts):
            pos = m.split("_")[0]
            var = m.split("_")[-1]
            mut = m.split("_")[-2].upper()
            ref = m.split("_")[1].upper()
            mut = correct_mutation_string(ref, mut)

            k = pileup_dict[str(pos)]
            b = k.split("\t")[4].upper()
            bases = get_base_list(b)
            reads = k.split("\t")[6].split(",")

            bases = np.array(bases)
            reads = np.array(reads)

            base = bases[reads == r][0].upper()
            if base in [",", "."]:
                rt += "R"
            elif base == mut:
                rt += "M"
            else:
                rt += "O"
        if rt == "MM":
            recombinant_reads.append(r)
    return recombinant_reads

!mkdir -p data_for_recombinant_read_analysis/recombination_sgRNA

min_length_of_long_softclip = 10

si = 1

for rid in list(df_samples["runid"]):
    print(si, end=" ")
    si += 1

    # getting basic info about sample
    variant_combination, poslist_file, bamfile, sample_short = prepare_sample_details(rid)
    
    # check if sample has already been processed
    if os.path.isfile("data_for_recombinant_read_analysis/recombination_sgRNA/" + sample_short + ".csv"):
        continue    
    
    # check if recombination detection results are already available
    if not os.path.isfile("data_for_recombinant_read_analysis/recombination_detection_results/" + sample_short + ".csv"):
        continue
    
    # if no recombinant reads were found, also skip
    df_tmp = pd.read_csv("data_for_recombinant_read_analysis/recombination_detection_results/" + sample_short + ".csv")
    df_tmp = df_tmp[df_tmp["recombinant_reads"]>0]
    if df_tmp.shape[0] == 0:
        continue
    
    # generating pileup
    pileup_large = !samtools mpileup -B -f $refseq -l $poslist_file --output-QNAME -d 0 -q 30 -Q 30 $bamfile
    pileup_dict = {k.split("\t")[1]:k for k in pileup_large[2:]}

    # defining empty lists to store data
    df_out = []
    
    # iterating through mutation pairs where recombinant reads were found
    for mp in df_tmp["mutation_pair"]:
        base_comb = " | ".join(sorted(mp.split(" | ")))
        base_comb_muts = base_comb.split(" | ")
        p = [int(k.split("_")[0]) for k in base_comb_muts]
        
        s_overlapping = get_overlapping_reads_for_sample(s, p, pileup_dict)
        s_recombinant = get_recombinant_reads_for_sample(s, base_comb_muts, s_overlapping, pileup_dict)
        
        # saving list of recombinant read IDs to a temporary file
        with open(sample_short + "_recombinant_reads_tmp.txt", "w") as f:
            f.write("\n".join(s_recombinant))

        # collecting read info from BAM file
        k = !samtools view -N $sample_short"_recombinant_reads_tmp.txt" $bamfile
        
        read_info_dict = {s_rid: {"isfound": False, "isleader": False, "issoftclipped": False, "islongsoftclipped": False} 
                          for s_rid in s_recombinant}
#         counted_read_ids = []
        for pp in k:
            if pp.split("\t")[0] not in s_recombinant:
                continue
#             counted_read_ids.append(pp.split("\t")[0])

            read_info_dict[pp.split("\t")[0]]["isfound"] = True

            # checking for leader sequence (it can be in either the read or its mate)
            if "GTAGATCTGTTCTCT" in pp.split("\t")[9]:
                read_info_dict[pp.split("\t")[0]]["isleader"] = True
            
            # checking for softclipped
            read_s = int(pp.split("\t")[3])
            read_e = int(pp.split("\t")[3]) + len(pp.split("\t")[9])
            if read_s <= p[0] and read_e >= p[1]: #if read really overlaps mutation pair (and not its mate)
                res = re.split('(\d+)', pp.split("\t")[5])
                cigar_char = res[2::2]
                cigar_num = [int(k) for k in res[1::2]]
                if "S" in cigar_char:
                    read_info_dict[pp.split("\t")[0]]["issoftclipped"] = True
                    if np.array(cigar_num)[np.array(cigar_char) == "S"].sum() > min_length_of_long_softclip:
                        read_info_dict[pp.split("\t")[0]]["islongsoftclipped"] = True

        s_df = pd.DataFrame()
        s_df["readID"] = [s_rid for s_rid in s_recombinant if read_info_dict[s_rid]["isfound"]]
        s_df["is_leader"] = [read_info_dict[s_rid]["isleader"] for s_rid in s_recombinant if read_info_dict[s_rid]["isfound"]]
        s_df["is_soft_clipped"] = [read_info_dict[s_rid]["issoftclipped"] for s_rid in s_recombinant if read_info_dict[s_rid]["isfound"]]
        s_df["is_soft_clipped_long"] = [read_info_dict[s_rid]["islongsoftclipped"] for s_rid in s_recombinant if read_info_dict[s_rid]["isfound"]]
        s_df["mutation_pair"] = base_comb
        df_out.append(s_df)
    df_out = pd.concat(df_out)
    df_out.to_csv("data_for_recombinant_read_analysis/recombination_sgRNA/" + sample_short + ".csv", index=False)
Code 
all_pos = [int(k) for k in np.linspace(1,30000,30000)]

# loading all files with results of sgRNA detection
all_sgrna_files = !ls data_for_recombinant_read_analysis/recombination_sgRNA/*.csv

sgrna_reads_neither = {k: 0 for k in all_pos}
sgrna_reads_leader = {k: 0 for k in all_pos}
sgrna_reads_sc = {k: 0 for k in all_pos}
sgrna_reads_both = {k: 0 for k in all_pos}

for fn in all_sgrna_files:
    rid = df_samples[df_samples["shortname"] == fn.split("/")[-1].split(".")[0]].iloc[0]["runid"]
    if rid in art_samples:
        continue

    var_comb = get_sorted_variants(df_coinf_allmuts[df_coinf_allmuts["runid"] == rid].iloc[0]["variants"])
    df_tmp = pd.read_csv(fn)
    
    df_tmptmp = df_tmp.groupby(["mutation_pair", "is_leader", "is_soft_clipped"], as_index=False).count()
    
    for ri, r in df_tmptmp.iterrows():
        posmin = int(r["mutation_pair"].split(" | ")[0].split("_")[0])
        posmax = int(r["mutation_pair"].split(" | ")[1].split("_")[0])
        for pi in range(posmin, posmax+1):
            if not r["is_leader"] and not r["is_soft_clipped"]:
                sgrna_reads_neither[pi] += r["readID"]
            elif not r["is_leader"] and r["is_soft_clipped"]:
                sgrna_reads_sc[pi] += r["readID"]
            elif r["is_leader"] and not r["is_soft_clipped"]:
                sgrna_reads_leader[pi] += r["readID"]
            elif r["is_leader"] and r["is_soft_clipped"]:
                sgrna_reads_both[pi] += r["readID"]

counts_neither = np.array(list(sgrna_reads_neither.values()))
counts_leader = np.array(list(sgrna_reads_leader.values()))
counts_both = np.array(list(sgrna_reads_both.values()))
counts_either = np.array(list(sgrna_reads_sc.values())) + counts_leader + counts_both

# figure to ratio of reads containing the leader sequence along the genome     
fig = plt.figure(figsize=(15,3))
gs = fig.add_gridspec(2, hspace=0.5, height_ratios = [0.2, 0.8])
axs = gs.subplots(sharex=True)

yy = 2
for g in gene_loc.keys():
    axs[0].hlines(y = 0, xmin=gene_loc[g][0], xmax=gene_loc[g][1], lw=30, color = color_dict[g])
    if (yy == 2):
        va_c = "bottom"
        ha_c = "left"
    else:
        va_c = "top"
        ha_c = "right"
    axs[0].annotate(g, xy= (np.mean([gene_loc[g][0], gene_loc[g][1]]), yy*0.5), 
                 xytext = (np.mean([gene_loc[g][0], gene_loc[g][1]]), yy), rotation=30, ha="center", va = va_c, fontsize=8,
                arrowprops=dict(arrowstyle="-"))
    yy *= -1
axs[0].set_ylim(-1,1)
axs[0].set_xlim(0,30000)
axs[0].set_axis_off()

id_nonzero = (counts_either+counts_neither) > 0
axs[1].plot(np.array(all_pos)[id_nonzero], (counts_leader+counts_both)[id_nonzero]/(counts_either+counts_neither)[id_nonzero])
axs[1].set_xlabel("Genomic position", fontsize=10)
axs[1].set_ylabel("Percentage of recombinant reads\n with patterns of sgRNA", fontsize=10)
axs[1].set_xlim(0,30000)
axs[1].set_ylim(0,0.3)
axs[1].grid(visible=False)
axs[1].tick_params(size=5, color="#666666")
axs[1].label_outer()

ymin, ymax = axs[1].get_ybound()
for g in gene_loc.keys():
    axs[1].vlines(x=gene_loc[g][1], ymin=ymin, ymax=ymax, color="black", alpha=0.5, lw=0.5)
axs[1].vlines(x=gene_loc["ORF1ab"][0], ymin=ymin, ymax=ymax, color="black", alpha=0.5, lw=0.5)

vals = axs[1].get_yticks()
axs[1].yaxis.set_major_locator(mticker.FixedLocator(vals))
axs[1].set_yticklabels(['{:.0%}'.format(x) for x in vals])
    
plt.show()
plt.close()

Figure 7: The percentage of recombinant reads carrying signs of originating from sgRNA along the genome. Artificial samples were not considered for this analysis.

Leader sequences and soft-clipping in recombinant reads were almost exclusively present when the reads overlapped gene boundaries (S/ORF3a, E/M and ORF8/N). This is in line with our previous knowledge about sgRNA formation and it also suggests that recombinant reads overlapping intragenic regions are unlikely to originate from sgRNA and are rather products of recombination occurring on the genomic RNA.

5.5 Recombinant reads in putative subclonal recombinant samples

Based on measured AF shifts along the genome, 13 samples have been selected as showing putative signs of recombination (see Supplementary File 2). (Note, that 6 of these were artificial mixtures of study PRJNA817870 from Bal et al., 2022.)

The figures below illustrate the distribution of recombinant reads in these 13 samples along with the putative recombination breakpoint that was determined from the location of the AF shift of defining mutations (Supplementary File 2).

Code 
bluecol = "#023e8a"
redcol = "#c1121f"

df_recomb = pd.read_csv("AFbased_recombinants_20230705.csv")
AFrecomb_runids = list(df_recomb["runid"])

gs_top = plt.GridSpec(len(AFrecomb_runids)+1, 1, 
                      height_ratios=[1/(len(AFrecomb_runids)+1)]+[2/(len(AFrecomb_runids)+1)]*len(AFrecomb_runids), 
                      top=0.89)
gs_base = plt.GridSpec(len(AFrecomb_runids)+1, 1, hspace=0)
fig = plt.figure(figsize=(15,15))

# top panel showing genes
topax = fig.add_subplot(gs_top[0,:])
yy = 2
for g in gene_loc.keys():
    gene_min = gene_loc[g][0]
    gene_max = gene_loc[g][1]
    if gene_max < 20000:
        continue
    if gene_min < 20000:
        gene_min = 20000
    topax.hlines(y = 0, xmin=gene_min, xmax=gene_max, lw=30, color = color_dict[g])
    if (yy == 2):
        va_c = "bottom"
        ha_c = "left"
    else:
        va_c = "top"
        ha_c = "right"
    topax.annotate(g, xy= (np.mean([gene_min, gene_max]), yy*0.5), 
                 xytext = (np.mean([gene_min, gene_max]), yy), rotation=30, ha="center", va = va_c, fontsize=8,
                arrowprops=dict(arrowstyle="-"))
    yy *= -1
topax.set_ylim(-1,1)
topax.set_xlim(20000,30000)
topax.set_axis_off()

axes = [fig.add_subplot(gs_base[1,:])]
for i in range(2, len(AFrecomb_runids)+1):
    axes.append(fig.add_subplot(gs_base[i,:], sharex=axes[0]))


axes[-1].set_xlabel("Genomic position", fontsize=10)
axes[-1].tick_params(size=5, color="#666666")

for ax in axes:    
    ax.set_xlim(20000,30000)
    ax.set_yscale("log")
    ax.grid(visible=False)
    ax.label_outer()
    ax.tick_params(size=5, color="#666666")

for rix, rid in enumerate(AFrecomb_runids):
    shortname = df_samples[df_samples["runid"] == rid].iloc[0]["shortname"]
    bp = df_recomb[df_recomb["runid"] == rid].iloc[0]["putative_breakpoint"]
    df_tmp = pd.read_csv('data_for_recombinant_read_analysis/recombination_detection_results/'+shortname+'.csv')
    
    num_of_overlapping_reads = {k:0 for k in all_pos}
    num_of_recombinant_reads = {k:0 for k in all_pos}

    for ri, r in df_tmp.iterrows():
        posmin = int(r["mutation_pair"].split(" | ")[0].split("_")[0]) 
        posmax = int(r["mutation_pair"].split(" | ")[1].split("_")[0])
        for pi in range(posmin, posmax+1):
            num_of_overlapping_reads[pi] += r["number_of_overlapping_reads"]
            num_of_recombinant_reads[pi] += r["recombinant_reads"]

    counts_overlapping_reads = np.array(list(num_of_overlapping_reads.values()))
    id_nonzero = counts_overlapping_reads > 0
    ratio_recomb_reads = np.zeros_like(all_pos).astype(float)
    ratio_recomb_reads[id_nonzero] += np.array(list(num_of_recombinant_reads.values()))[id_nonzero]/counts_overlapping_reads[id_nonzero]
    
    axes[rix].plot(all_pos, counts_overlapping_reads, color=bluecol, linestyle="dotted", lw=1)
    axes[rix].tick_params(axis='y', labelcolor=bluecol)
    
    ax_tmp = axes[rix].twinx()    
    ax_tmp.vlines(x = bp, color="black", ymin=0, ymax=1, lw=2)
    ax_tmp.plot(all_pos, ratio_recomb_reads, color=redcol, lw=1.5)
    ax_tmp.tick_params(axis='y', labelcolor=redcol)
    ax_tmp.set_ylim(0,1)
    ax_tmp.set_yticks([0, 0.4, 0.8])
    ax_tmp.grid(False)
    for g in gene_loc.keys():
        ax_tmp.vlines(x=gene_loc[g][1], ymin=0, ymax=1, color="black", alpha=0.5, lw=0.5)
    ax_tmp.vlines(x=gene_loc["ORF1ab"][0], ymin=0, ymax=1, color="black", alpha=0.5, lw=0.5)    
    
    sname = df_samples[df_samples["runid"]==rid].iloc[0]["sample_accession"]
    if rid in art_samples:
        sname += "\n(artificial mixture)"
    ax_tmp.annotate(sname, 
                    xy = (0.005, 0.94), fontsize=10, ha="left", va="top", color="black", xycoords = "axes fraction")
    
    ax_tmp.annotate(str(num_of_overlapping_reads[bp]), 
                    xy = (bp+50, 0.94), fontsize=10, ha="left", va="top", color=bluecol)
    ax_tmp.annotate(str(round(ratio_recomb_reads[bp-1], 3)), 
                    xy = (bp+50, 0.74), fontsize=10, ha="left", va="top", color=redcol)

    
    
axes[0].annotate("Number of overlapping reads", xy = (-0.05, 1-len(AFrecomb_runids)/2), fontsize=12, rotation=90, 
             ha="center", va="center", color=bluecol,
                 xycoords = "axes fraction")        
axes[0].annotate("Ratio of recombinant reads", xy = (1.05, 1-len(AFrecomb_runids)/2), fontsize=12, rotation=-90, 
             ha="center", va="center", color=redcol,
                 xycoords = "axes fraction")        


plt.show()
plt.close()

Figure 8: The number of overlapping reads (dotted blue lines) and the ratio of ones showing signs of recombination among them (solid red lines) for samples identified as putative subclonal recombinants from shifts in AF distribution. The figure is limited to the 20,000-30,000 genomic position range, as putative recombination breakpoints (shown with black vertical lines) were also confined to this region. Blue and red numbers right after the breakpoints indicate the number of reads overlapping the breakpoint and the ratio of ones carrying signs of recombination among these.

None of the putative breakpoints identified from AF analysis (Supplementary File 2) could be verified from read-level data with a recombinant read ratio of 0.1 or larger (either for artificial or real samples). In many cases, the number of reads overlapping the putative breakpoint was ab ovo very low.