Origins and Effects of Mutations in SARS-CoV-2 Genome

Bioinformatics Program Technical University of Munich Ludwig-Maximilians-Universit¨at M¨unchen Bachelor’s Thesis in Bioinformatics Origins and Effects of Mutations in SARS-CoV-2 Genome Stefan Dvoretskii Bioinformatics Program Technical University of Munich Ludwig-Maximilians-Universit¨at M¨unchen Bachelor’s Thesis in Bioinformatics Entstehen und Auswirkungen der Mutationen im Genom des SARS-CoV-2 Origins and Effects of Mutations in SARS-CoV-2 Genome Author: Stefan Dvoretskii Supervisor: Julien Gagneur, TUM Chair of Computational Molecular Medicine Advisors: Christian Mertes, Julien Gagneur Submitted: 15.09.2020 I confirm that this Bachelor’s thesis is my own work and I have documented all sources and materials used. Abstract Infection waves caused by Coronavirinae break out yearly at seasonal cycles, usually achieving significant number of infected humans, but in most of them only leading to relatively mild clinical symptoms and illness progression. While other viruses of the Coronavirinae subfamily have already been involved in more severe outbreaks in terms of symptoms, SARS-CoV-2 outbreak that has started in November 2019 in Wuhan, China has had an unprecedented outcome both in clinical terms and further societal and everyday life aspects. Alongside with its far-reaching spread, SARS-CoV-2 genome has mutated, acquiring genomic variants that could have effects impacting virus lifestyle, exempli gratia altering its biological machinery processes or even make it resistant to a drug treatment. Variants in the SARS-CoV-2 genome could moreover be indicative of the path, both geographic, temporal and phylogenetic, laid back by the virus bearing the genome. To analyse the variants, as well as their development history and potential impact, I have designed a genomic pipeline operating on SARS-CoV-2 next generation sequencing reads. I show that this pipeline is useful for finding potentially pathogenic variants or spatial and temporal relationships of viral genomes assembled from patient swab probes. This information is first and foremost proposed to be used in clinical settings, for instance as preliminary information for drug treatment adjustments or clinical arrangments of patients. Durch Coronavirinae verursachte Infektionswellen brechen in saisonaler Weise aus und infizieren ¨ublicherweise eine signifikante Anzahl an Menschen. Die meisten Infizierten bekommen jedoch nur milde bis gar keine Symptome. Andere Viren der Subfamilie waren bereits in Erkrankungen im Menschen involviert und haben dabei teilweise st¨arkere Symptome als SARS-CoV-2 verursacht. Die durch SARS-CoV-2 verursachte Pandemie hat aber zu unerwarteten Konsequenzen im klinischen als auch im gesellschaftlichem Bereich und dem allt¨aglichen Leben gef¨uhrt. W¨ahrend der weitreichenden Ausbreitung des Virus ist das Genom von SARS-CoV-2 nebenbei mutiert. Solche genetischen Ver¨anderungen k¨onnen die biologischen Prozesse des Virus ver¨andern, was auch klinisch relevante Folgen haben k¨onnte, zum Beispiel falls der Virus dadurch gegen die Behandlung mit bestimmten Arzneimitteln resistent wird. Varianten im Genom des SARS-CoV-2 k¨onnen außerdem im Sinne der phylogenetischen, zeitlichen bzw. geographischen Entwicklung des Virus informativ sein. Um die Varianten, ihre Entwicklungsgeschichte und m¨ogliche Auswirkungen zu analysieren, habe ich eine genomische Pipeline entwickelt, die mit rohen NGS Reads als Eingabe verschiedene Informationen ¨uber das virale Genom und dessen Varianten ableitet und untersucht. Ich zeige, dass diese Pipeline sowohl potentiell pathogenische Varianten als auch Beziehungen zwischen untersuchten viralen Genomen aufdeckt. Diese Informationen k¨onnen in der ersten Stelle der Verbesserung der Behandlung in der Praxis dienen, zum Beispiel der Behandlung mit Arzneimittteln, beziehungsweise um Kreuzinfektionen fr¨uhzeitig zu erkennen. Acknowledgment I want to kindly thank all the members of Chair of Computational Molecular Medicine at the Technical University of Munich, current and past, for creating a unique working atmosphere and providing a comfortable environment for the research. I would especially like to thank  Julien Gagneur, the Head of the Chair, for his guidance into the scientific working as well as his ideal support on the present project. I am also grateful for believing in my scientific potential and having me as a student apprentice during the last two years in the laboratory. The working experience I have gained during that time has allowed me to be best prepared for this Thesis work.  Christian Mertes, for his guidance in the project, as well as his support all along the way during this project work, let it be scientific, programming or visualization questions.  Alex Karollus, for his help at the last steps of my work as well as his readiness to continue the investigation on the topic of the present Thesis. Apart from this, I want to say thanks to all of the lecturers and methodologists involved in making the B. Sc. Bioinformatics course at the TU Munich/LMU so informative and equipping for the current Thesis. I am also thankful to my classmates for support and cheerful atmosphere. Last but not least, I am grateful to my family members and close friends for supporting me during the work on the current Thesis. Contents List of Tables iii List of Figures iii List of Abbreviations v 1 Introduction 1 1.1 Novel human Coronavirus SARS-CoV-2 . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1 Biological characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.2 Proteinmachinery .................................. 3 1.1.3 Known variants and their impact . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2 NGSRNAsequencing .................................... 5 1.2.1 Twist Fast Hybridization Target Enrichment . . . . . . . . . . . . . . . . . . . 5 1.2.2 Ampliconsequencing................................. 6 1.2.3 Qualitycontrol.................................... 7 1.2.4 Downstreamanalysis................................. 8 1.2.5 Variantcalling .................................... 8 1.3 Viral variants and their correspondence with patient clinical history . . . . . . . . . . 8 2 Material and Methods 10 2.1 Initialqualitycontrol..................................... 10 2.1.1 InitialQCbatchofdata............................... 10 2.1.2 In silico quality assurance and initial variant calling pipeline . . . . . . . . . . . 11 2.2 Mutations and variant effects investigation . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2.1 Variant Effect Predictor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2.2 Mutationentropy................................... 12 2.2.3 Structural features analysis of the protein product . . . . . . . . . . . . . . . . 13 2.3 High-throughputpipeline .................................. 13 3 Results and Discussion 16 3.1 Qualitycontrol ........................................ 16 3.1.1 Pipelineoutput.................................... 16 3.1.2 QC results interpretation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.2 Comparative performance of variant calling . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2.1 bcftoolsvs.lofreq .................................. 19 3.2.2 Our pipeline and CoVpipe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.3 Variantsinvestigation .................................... 24 3.3.1 Overview ....................................... 24 3.3.2 15101C>T variant potentially increases resistance of SARS-CoV-2 to Remdesivirtreatment .................................... 26 3.3.3 24453A>G variant could change the binding of SARS-CoV-2 to ACE2 receptor 27 i Contents Contents 4 Conclusion and Outlook 30 4.1 SummaryandConclusion .................................. 30 4.2 Outlook ............................................ 30 4.2.1 Additionaldata.................................... 30 4.2.2 Quasispeciestheory ................................. 31 Bibliography 32 Appendix i SupplementaryData........................................ i ii List of Tables 2.1 QCbatchsamples.............................. 10 List of Figures 1.1 Coronavirus replication . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Discontinuous transcription . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Twist Hybridization protocol . . . . . . . . . . . . . . . . . . . . . . . . 6 1.4 Amplicon sequencing workflow . . . . . . . . . . . . . . . . . . . . . . . 7 2.1 IGValignmentcheck............................ 12 2.2 V-pipeschema................................ 14 3.1 Read coverage per postition plot . . . . . . . . . . . . . . . . . . . . . . 17 3.2 ReadcoverageECDF............................ 17 3.3 Human reads contamination check . . . . . . . . . . . . . . . . . . . . . 18 3.4 Scatterplot of lofreq and bcftools variants . . . . . . . . . . . . . . . . . 20 3.5 bcftools calls heatmap, V-pipe vs. CoVpipe . . . . . . . . . . . . . . . 21 3.6 Difference between V-pipe-SARS-CoV-2 and CoVpipe calls . . . . . . . 22 3.7 Nextstrain entropy for pipeline calls . . . . . . . . . . . . . . . . . . . . 23 3.8 V-pipe-SARS-CoV-2 calls heatmap . . . . . . . . . . . . . . . . . . . . 25 3.9 Remdesivirbinding............................. 26 3.10 Conserved amino acid 553 with Remdesivir resistance . . . . . . . . . . 28 3.11 Mutation 24453A>G in protein structure . . . . . . . . . . . . . . . . . 29 A.1 Per base Phred score plot . . . . . . . . . . . . . . . . . . . . . . . . . i A.2 Read adapter content, FastQC . . . . . . . . . . . . . . . . . . . . . . . i A.3 Read duplication stats . . . . . . . . . . . . . . . . . . . . . . . . . . . ii A.4 Alignmentstatistics............................. ii A.5 Phred quality per position . . . . . . . . . . . . . . . . . . . . . . . . . iii A.6 Mean Phred quality distribution . . . . . . . . . . . . . . . . . . . . . . iii A.7 Ambiguousbasesplot............................ iv iii 1.2 NGS RNA sequencing Introduction 1.1.3 Known variants and their impact Numerous variants in Non-structural protein 12 (a.k.a. RdRp a.k.a. Pol) were proven to alter Remdesivir treatment outcomes [18], supposedly because of changed folding and, as a result, decreased binding of nsp12 to the RNA template. Above that, some mutations in the nsp14 are found to effectively change this enzyme activity, and retain it from the proofreading. This results in a drastic increase in mutation rate of the virus in-host [19]. There are hints that all of the Non-structural proteins 12 to 16 could be especially susceptible to a significant mutation outcome as members of the viral copying machinery [20]. Furthermore, a mutation in S1 subunit of Spike protein is proven to reduce subunit shedding and increase invectivity [21]. Some screening studies have been performed, assessing effects of multiple mutations in the protein coding regions of SARS-CoV-2 genome [22] [23]. Other less direct ways of missense mutations being effective on the level of protein machinery interactions, like mutations in nsp1 leading to faster cell recovery or mutations in nsp5 making subproteins of pp1ab unseparatable are, of course, also thinkable. 1.2 NGS RNA sequencing The advent of new generation sequencing, exempli gratia Illumina high-throughput sequencing protocols [24], makes the parallel and high-throughput sequencing of probes from clinical patients more efficient and less costly. This fact stands for increased sequencing capabilities in virological studies and hence rapid increase in data that can and needs to be analysed in the field of virology [25]. There are some specific steps in virus NGS protocol like viral enrichment [26], which increases viral signal in the resulting read data. 1.2.1 Twist Fast Hybridization Target Enrichment This Next-generation Sequencing protocol is developed by Twist Biosciences [27]. It is schematically summarized in Figure 1.3. SARS-CoV-2 variant of targeted NGS consists of 6 probes appx. 5kbp each, that enrich for a specific region in the viral genome. This allows to increase viral signal and decrease sequencing costs. 5 1.2 NGS RNA sequencing Introduction Figure 1.3: Targeted NGS as compared to the “normal” NGS. Probe “captures” the zone of interest in the (meta)genome and enriches for its reads. Source: Twist Biosciences. 1.2.2 Amplicon sequencing Amplicon sequencing is a highly targeted approach that enables researchers to analyze genetic variation in specific genomic regions. The ultra-deep sequencing of PCR products (amplicons) allows efficient variant identification and characterization. This method uses oligonucleotide probes designed to target and capture regions of interest, followed by next-generation sequencing (NGS) [28]. See Figure 1.4 for a schematic overview of the protocol workflow. 6 1.2 NGS RNA sequencing Introduction Figure 1.4: Amplicon sequencing workflow. Amplicon sequencing enables a wide range of research applications for the discovery, validation, or screening of genetic variants. Source: Illumina [28]. 1.2.3 Quality control The first step in the down-stream analysis of NGS data is quality control of reads. The steps here could include measuring Phred score across the reads, which gives a hint on the ambiguity of the base calls, GC content of the sequence and percent of reads mapped to the reference genome [29]. Quality control steps are for assessment of 7 1.3 Viral variants and their correspondence with patient clinical history Introduction how good the raw data or the data after initial processing steps is for the downstream analysis. Furthermore, its results also point out the experimental steps that can be improved. This information is used in further experiments. 1.2.4 Downstream analysis After the raw read quality is controlled for, downstream analysis of the reads can begin. It includes, but is not limited to: preprocessing of the reads (including steps like read trimming etc. that suit the sequencing protocol and research questions), de novo genome assembly and alignment to a reference sequence. 1.2.5 Variant calling Once the reads are aligned to a reference genome, one can begin with calling variants. This is done by step-wise analysing the reads covering a same interval and investigating those bearing an alternative sequence (which could be a single base in case of SNV and larger reference/alternative sequence, that would be an InDel) as compared to the reference genome. This information tells to a certain degree of significance how the genetic sequence specified as reference differs from the aligned reads in terms of base content. A variant, described as a sequence substitution at a defined position in the reference genomic sequence, thus acquires a variant frequency - i.e., the fraction of all reads relevant to its motif that supports its alternative sequence. Formally: fvariant =Rvariant R where fvariant is an estimated frequency of variant, Rvariant is the quantity of relevant reads supporting the variant sequence, Ris the quantity of all reads that are relevant to the variant’s position. Further metrics like base read coverage at a given position, quality of variant call, genotype etc. are important in context of a given variant. They can all play role in further quality assessment and interpretation of variants. Good quality of a variant call indicates that it is unlikely to happen under the null model assumption; big read coverage of a variant call makes sequencing errors that bias the alternative allele frequency less likely. Significantly big alternative allele frequency indicates that quantitatively more reads and hence viral copies bear this variant, making it clinically more significant. 1.3 Viral variants and their correspondence with patient clinical history Variant calls are connected with viral evolution both between the hosts (intra-host) as well as in a host (inter-host); both are of importance. We make some hypotheses prior to the variant analysis:  Viral genome variants that increase its viability/resistance to a drug treatment could increase in frequency in one particular host. 8 1.3 Viral variants and their correspondence with patient clinical history Introduction  Multiple infection of one person could be possible.  Two patients, especially those spending significant time in close spatial proximity, could be cross-infecting each other with their own specific strains, which could potentially then be seen based on the genomic measurements and their time points The first deliverable of the variant calling pipeline I construct is the identification of variants in the viral samples. The goal connected to this one is the placement of the identified variants in the context of mutations in other viral samples worldwide [30] [31]. The second goal is the understanding of possible impact of identified variants on clinical outcome of the patients. Last goal is the correlation of genomic variants with patients clinical history. 9 2 Material and Methods 2.1 Initial quality control 2.1.1 Initial QC batch of data The initial batch sent by the collaborators conducting sequencing experiments included 5 patient samples with different viral titers and 3 control samples. It was intended for a primitive QC and determination of the proper protocol for sequencing experiments. This included assessing number of sequencing cycles (for details on sequencing cycles in Illumina NGS workflow see [24]) and viral titer of samples that would be sufficient for the pipeline purposes (i.e., variant calling). Samples in the QC batch contained all log10(viral titer) values in the range from 2 to 6, as well as two dilutions of control viral RNA from the prepping kit. Table 2.1 summarizes information on the samples in the first QC batch. SampleID viral titer S244675 5 ∗102 S244378 5 ∗103 S244564 5 ∗104 S244676 5 ∗105 S244395 5 ∗106 Ctrl4 (greater dilution of control viral RNA) Ctrl2 (lesser dilution of control viral RNA) Ctrl0 (negative control) Table 2.1: Overview of the viral titer in the samples from the first, QC batch. Ctrl4 and Ctrl2 samples contained control viral RNA from Twist prepping kit, whereby concentration of viral RNA in Ctrl4 was planned to be appx. 2 times bigger than in Ctrl2. Ctrl0 sample was a negative control (i.e., not intended to contain any RNA). The prepping kit was from Twist Biosciences with Twist Synthetic SARS-CoV-2 RNA Control 2 (GenBank ID MN908947.3) [27] as the viral reference sequence. Libraries were prepped from all the eight samples and pooled afterwards. For viral enrichment, Twist Fast Hybridization Target Enrichment Protocol was used [26]. Sequencing was performed using Illumina MiSeq system [32] with paired-end reads and 251 sequencing cycles (from each end). Samples were multiplexed. 10 2.1 Initial quality control Material and Methods 2.1.2 In silico quality assurance and initial variant calling pipeline With the purpose of analysing the accordance of experimental protocol with the goal of calling variants, I have assembled an initial variant calling pipeline, that received a codename of ”virus-seq”. It includes following steps: 1. Initial FastQC analysis of raw reads. For that, I have applied fastqc tool, a quality control tool for high throughput sequence data [33]. 2. Read trimming: for that, I have used Trimmomatic suite [34]. Its configuration was based on the manual assessment of the reports from the FastQC analysis. Namely, I have trimmed the first 3 leading bases of each read as well as any sequences longer than 150 basepairs, because the Phred scores [35] after 150 basepairs in read were estimated as ”mediocre” or ”bad” by FastQC. See Figure A.1 for an example per-base Phred quality plot. Another trimming step was to remove Illumina adapters. This was also reasoned by FastQC output for adapter content of reads (see Figure A.2 for an example plot). For adapter removal, I have used Illumina TruSeq3 adapter list [36]. I have applied default Trimmomatic settings for ILLUMINACLIP rule explained in Trimmomatic manual [37] along with the option of keeping both reads after adapter clipping to preserve downstream tools consistency (as they work with paired reads only). Additionally, after the listed trimming steps, I have removed read sequences shorter than 36 basepairs. 3. Reference alignment: for this step, Burows-Wheeler-Aligner against the reference sequence NC 045512v2 [38] [39] was used. The choice of aligner was motivated by the fact that the genome of SARS-CoV-2 is mostly coding and in particular does not contain any spliced-out regions. In this step, I have run bwa mem algorithm [40] with default settings. 4. After the alignment, the pipeline proceeds to apply Picard’s MarkDuplicates tool [41] to remove biological as well as artificial duplicate reads from the alignment. 5. The next step is running samtools statistics tools including depth, flagstats etc. to generate various alignment statistics. 6. In this step, MultiQC is used. It combines output reports from the tools of the previous steps. MultiQC automatically recursively searches the specified directories for output of the other tools using name matching and generates a summarized report over the samples. 7. Variant calling. For this step, I have used a series of [42] calls, that goes as following: bcftools mpileup –max-depth 6000 -f NC 045512.v2.fasta (conversion of reference alignment in the pileup format, maximum of 6000 reads at a given position are taken into account, comparison with NC 045512.v2.fasta [9]), bcftools call - mv -Ov (call variants, use alternative multiallelic and rare variants caller, output in the Variant Calling Format) and bcftools view –include ”INFO/DP≥20” (include only variant calls with the sequencing depth of more than 20 reads). 11 2.2 Mutations and variant effects investigation Material and Methods 2.2 Mutations and variant effects investigation The assesment of the called variants can be separated in two principal steps: 1) assessment of the quality of variant (does its context make sense?) and 2) assessment of the consequences of variant (how far does its effect propagate?). For the first step, I have used IGV browser [43] to check variants ”in the raw data”, i.e. directly in the alignment. Figure 2.1: An example of IGV alignment screenshot for a sample. The tracks are aligned by coordinate, which is specified by the upmost line. The ”Annotations” track contains genomic annotations fetched from a file. The bottom track is the alignment view. The upper part of it (area plot) depicts the per-base read coverage, while the bars below indicate separate reads aligned. The colored lines on the reads illustrate read alternative bases as compared to the reference, each base has its own color. This view lets check variants by assessing how many reads ”have a variant” at a given position (i.e., support the sequence of a variant), and compare it to the called alternative allele frequency of a variant. 2.2.1 Variant Effect Predictor In the second step, I have used Variant Effect Predictor (VEP) [44] to assess the potential outcome of the variants. VEP assessments are sequence-based, and normally also take in account the transcript-exon model; which, however, was not adding new information due to the absence of introns in the SARS-CoV-2 genome. I had to manually generate .gff annotation required by the tool automatically downstream from UCSC browser export [45] - it had to have a proper data format required by VEP. Variant Effect Predictor output includes class of mutation (including synonymous, missense mutations, stop gained mutation classes etc.), amino acid substitution happening as a result of the mutation and position of amino acid in the protein product, which is used in the later steps of analysing protein product structural features. 2.2.2 Mutation entropy I have added Nextstrain browser [31] information on the position entropy in SARSCoV-2 genome. Nextstrain browser syncs data from GISAID sequencing intiative [30], 12 2.3 High-throughput pipeline Material and Methods which is a Global Initiative on Sharing All Influenza Data. It builds a phylogenetic tree from it and further artifacts, one of those being a worldwide entropy track, including information of per-position entropy from the synchronized samples alignment. Higher entropy means that more samples from the observed representative set of the samples worldwide (N=4449 at the accession date) diverge in the nucleotide on the position, whereas lower entropy means quantitatively greater base conservation. This gives a hint on how often a certain position in SARS-CoV-2 genome is mutated worldwide. 2.2.3 Structural features analysis of the protein product I used PyMOL [46] and JyMOL [47] as well as their embedded versions to interactively visualize protein structures. The assessment was done manually for the variants found to most promising based on their frequency, predicted mutation class and entropy. For that, after picking the proper Uniprot ID [48] of the protein product structure (e.g. 6yyt for Non-structural Protein 12) [49]. I analyzed the surface structure of the protein product using the interactive view of protein structures provided by Uniprot [47], especially taking into account whether  the mutation was close or far away from the active center of the protein  the mutation occured in a motif with defined secondary structure (α-helix, βsheet etc.)  the mutation occured in another annotated region of protein - exempli gratia a transmembrane region or a Uniprot domain This way, I was able to find some protein variants potentially having an effect on the clinical outcome of a patient. 2.3 High-throughput pipeline With the knowledge yielded from the first and second QC reports, the experimental collaborators have sent more patient samples (N=563) to a high-throughput sequencing platform, increasing sample quantity sequenced in one batch from tens to hundreds. It was clear that the pipeline had to be designed for a greater flexibility. The most important requirements were segregation and parallel execution of jobs, as well as good capacities for managing dependencies and cluster execution. I have found an actively developed viral NGS processing pipeline called V-pipe [50]. Basically, the pipeline consisted of the steps with the same goals and results as my initial virus-seq pipeline. Figure 2.2 gives an overview of those. 13 2.3 High-throughput pipeline Material and Methods Figure 2.2: Schematic overview of the steps covered in the original V-pipe pipeline. I have left out steps dedicated to quasispecies detection from V-pipe to begin with, as well as changed some minor details of the pipeline, like filtering alignments for samtools flags. I called the new, modified pipeline V-pipe-SARS-CoV-2. It goes as following: Each sample is clearly identified by a subject ID and timestamp, at which it was taken. For each sample, the initial input are 2 read files in FASTQ format [51], one file for forward and reverse reads respectively. This reads are extracted from a gzip archive if provided and put in the read preprocessing software, namely PRINSEQ [52], used for quality control based on sliding window Phred quality approach and Ns threshold, and Trimmomatic [34] restricted to just the read trimming fully identical to the virus-seq algorithm. Once the read preprocessing step is dealt with, the pipeline proceeds to align reads to the reference genome, which is set to be NC 045512.2 (wuhCor1), the same as in virus-seq. We use the same variant calling algorithm as in virus-seq and Variant Effect Predictor with the virus-seq settings. On top of that, the original visualization step was modified in V-pipe-SARS-CoV-2 to obtain an ”all-in-one-place” report, including VEP annotation and Nextstrain entropy of the variants. The architecture of the pipeline, which is built using Snakemake [53], allows for great parallel processing and documentation capabilities. The newest version 14 3.2 Comparative performance of variant calling Results and Discussion Figure 3.5: Heatmap of variant calls made by the V-pipe-SARS-CoV-2 (“our pipeline”). bcftools was used as variant caller. Color intensity stands for variant frequency (from 0.0 to 1.0). Rows are different variants identified by the position and alternative allele sequence. Note that not each variant name/tick is displayed because of the space limitations. Columns stand for samples; their names are hidden here for visualisation purposes. Only variants that are called by both pipelines are displayed. Notably, positive controls like 23403-G, that are known to be highly spread in European patients, are called. 21 3.2 Comparative performance of variant calling Results and Discussion Figure 3.6: Heatmap of variant calls made by our modified version of V-pipe (“our pipeline”) as compared to the CoVpipe calls (“their pipeline”). Color intensity stands for V-pipe-SARS-CoV-2 variant frequency (from 0.0 to 1.0) minus CoVpipe variant frequency. Rows are different variants identified by the position and alternative allele sequence. Note that not each variant name/tick is displayed because of the space limitations. Columns stand for samples; their names are hidden here for visualisation purposes. Only variants that are called by both pipelines are displayed. 22 3.2 Comparative performance of variant calling Results and Discussion Figure 3.7: Nextstrain entropy histograms for variant calls. Upper instagram is for CoVpipe, while bottom is for V-pipe-SARS-CoV-2. Much smaller relative proportion of high-entropy variants that correspond to clades-characteristic variants in SARS-CoV-2 worldwide phylogenetic tree [31] means that CoVpipe calls a lot of “novel” variants not seen in the samples worldwide before. Given the size of the cohort (N=182), such situation is dubious. 23 3.3 Variants investigation Results and Discussion 3.3 Variants investigation Mutations of SARS-CoV-2 can potentially lead to significant clinical outcomes, such as resistance to drug treatment, or at least changes in biological processes that indirectly affect clinical case. Our pipeline allows to call and identify multiple aspects of genetic variants in SARS-CoV-2 samples taken from patients to correlate them with patients’ treatment information and deduce possible coincidences between the two. 3.3.1 Overview In the data from the first two batches and two “high-throughput” batches that have had at least one high-frequency variant (alternative allele frequency in reads >0.1) (N=291), we recover 163 unique variants. From those, above the half is missense, while another quarter is synonymous, and about 10% represent translation stop site gain or loss A.14. The overview of variant call frequencies is given by Figure 3.8. 24 3.3 Variants investigation Results and Discussion Figure 3.8: Heatmap of variant frequencies for all unique variants in V-pipe-SARS-CoV-2 calls. Rows are variant identifiers; some are hidden because of space limitations. Columns represent samples; their names are also hidden because of the space limitations. Rows and columns are clustered using the complete linkage matrix according to Euclidean metric. We clearly see a cluster of “European” characteristic mutations in the most upper side of the heatmap. This variants are absent in some samples because these do not report sufficient coverage at the variant positions. There is also an interesting cluster on the left top of the heatmap just under the “European” rows that includes several variants. It is comprised of samples of the high throughput batches, as well as S14 and S0 from the second QC batch according to Twist protocol, which are known to be similar by our earlier heatmaps (e.g. Figure A.15). Most of the variants are however only seen in single or few samples. This could mean that they disappear with the time, or do not achieve sufficient coverage in the other relevant samples (e.g. samples from the same patient). 25 3.3 Variants investigation Results and Discussion 3.3.2 15101C>T variant potentially increases resistance of SARSCoV-2 to Remdesivir treatment A partial cohort of clinical samples we ran through the pipeline has undergone Remdesivir treatment [61]. S244395 S5 (a.k.a S5 5M) was a sample from one of such patients, together with that having more than satisfying read base coverage and good read quality, making it suitable for highly trustable variant calls. One of such called variants was 15101C>T, which as we could see from Variant Effect Predictor output [44] was in Non-structural protein 12 (further nsp12) of SARS-CoV-2 and very close to the active center of the protein. Non-structural protein Remdesivir ”builds” itself into the active center of nsp12, where the RNA template is normally being “passed through”. Residue 554A, whose second codon has nucleotide coordinate 15101, is positioned very close to the place where Remdesivir binds, as could be seen in Figure 3.9. Figure 3.9: nsp12 molecular structure (UniProt ID: 6yytA) in PyMOL [46] with ligands. Spirals show α-helices, while arrows show β-sheets in the protein structure. Green, pink and cyan colors depict different sub-chains of nsp12. Dark orange spiral with rosa-blue tails is the RNA template and nucleotides on it. Yellow arrow pointing to the yellow compound is an arrow indicating Remdesivir docking site. Red arrow points to the residue 554A, which is affected by the 15101C>T mutation. As could be seen from single disclosed clinical case of the single sample, S244395 S5 26 3.3 Variants investigation Results and Discussion adheres to an individual who has received a 10-day long treatment with Remdesivir, but had a prolonged illness scenario (measurement on day 28 still positive), in contrast to usually shortened cases as a result of Remdesivir treatment. We speculate that this mutation could be changing the shape of the active center of nsp12, which is suggested by the proximity of suscepted amino acid to the active center as well as known conservation and mutational Remdesivir resistance in the neighbor amino acid (553) found by another study [18]. 3.3.3 24453A>G variant could change the binding of SARS-CoV-2 to ACE2 receptor A variant identified only in one sample by V-pipe-SARS-CoV-2, but in more than 90% of the samples (N=182) by CoVpipe (Figure 3.6) involves the amino acid 964K in the Spike protein of SARS-CoV-2. Such divergence in the variant calling results is explained by duplicates removal. It is a missense variant with substitution of Lysine for Arginine. Figure 3.11 sheds a light on the structural position of according amino acid site. This substitution could potentially change the folding of Spike protein, hence changing its functional qualities. Probably the most interesting thing about this variant is that it has a Nextstrain entropy of 0.0. That means that it has not been encountered before in any of the Nextstrain worldwide pick representatives. So if CoVpipe calls are not based on artificial duplicates in this case, this variant could serve as a “barcoding” variant for our samples cohort. As the samples originate from the same hospital, this variant could above all indicate the local origin of the samples. 27 3.3 Variants investigation Results and Discussion Figure 3.10: Amino acid 553 in other members of Coronavirus family is known to make viruses resistant to Remdesivir treatment [18]. This is a neighboring position to amino acid 554. 28 3.3 Variants investigation Results and Discussion Figure 3.11: 24453A>G in the structural view of S protein. Subchains of the protein are depicted in different colors. Red arrow points to the 964K site in the structure, in which codon the variant is found. This amino acid site is contained in an α-helix structure. Pink arrows indicates the binding site of the Spike protein (right side). Given the relative proximity of the amino acid site to the binding site, an amino acid substitution could lead to conformational changes at the S protein binding site. 29 4 Conclusion and Outlook 4.1 Summary and Conclusion In this thesis, I propose a raw reads processing, variant calling and annotation pipeline that allows to make clinical insights into the genetic variation of SARS-CoV-2 between and inside patients. With the information gained, I have uncovered some variants that could have clinical importance, and report the highlight 15101C>T possibly increasing resistance to Remdesivir treatment inside patients. 4.2 Outlook 4.2.1 Additional data In the dataset which was originally used to write this thesis, all samples were anonymized. Clinical metadata would allow us to match virus samples to patients, hence bringing the timescale to our variant frequency observation, which could give us insights on the interpatient evolution of the virus. Timepoints of sequencing would give us the ability to hypothesize on cross-infection as well as to provide more exact estimates of the interpatient mutation propagation speed. Further fields of metadata could be informative, like treatment with drugs, type of clinical syptoms (mild/severe/critical). In particular, the abovementioned annotation would allow to answer e.g. following scientific questions:  How variant frequencies differ based on treatment?  Can we provide evidence for mutations occurring during hospital stay? Can those be caused by the treatment?  Is there evidence for cross-infections?  Are variants in genome connected with clinical case severity? Having the genome sequences of the clinical patients, we could combine our insights from human genomes variation and virus genome variation, potentially uncovering relationships between those that have a definite effect on the clinical case. 30 Appendix Supplementary Data Figure A.1: Phred score quickly detoriates after first 150 basepairs, as well as is relatively lower at the very beginning of reads, which has reasoned our choice of trimming parameters. Figure A.2: We see that adapter content starts to significantly grow after around 70bp in the read. i Appendix Figure A.3: Read duplication stats. QC batch samples are highlighted in green. Figure A.4: Reference genome alignment statistics. QC batch samples are highlighted in green. ii Appendix Figure A.5: Mean per position Phred quality in reads. QC batch samples are highlighted in green. Figure A.6: Mean per position Phred quality in reads. QC batch samples are highlighted in green. iii Appendix Figure A.7: N content per read base position. QC batch samples are highlighted in green. Figure A.8: Adapter content per read base position. QC batch samples are highlighted in green. iv Appendix Figure A.9: Duplicate reads in the samples of first and second data batches. Figure A.10: Duplicate reads in the samples of the second data batch, as provided by FastQC [33]. v Appendix Figure A.11: Reference genome alignment statistics. Second batch samples (17) bars are highlighted. Figure A.12: Viral titer versus reads mapped to SARS-CoV-2 reference genome scatterplot of quality control batch of data. vi Appendix Figure A.13: Called variants amounts of lofreq caller as compared to bcftools. A lot of lofreq calls do not get detected by bcftools at all, although most of them are made on low coverage positions and hence are generally dubious because of possible sequencing errors. vii Appendix Figure A.14: Classes of mutations in the unique variants identified by V-pipe-SARS-CoV-2. “None” class stands for mutations without any fetched annotation. Clear prevalence of missense variants goes alongline with selection hypothesis: variants having effect on clinical outcome persist and achieve high frequency (in this setup frequency >0.1). viii Appendix Figure A.15: Heatmap of variant frequencies for all unique variants in V-pipe-SARS-CoV-2 QC batches calls. Rows are variant identifiers; some are hidden because of space limitations. Columns represent samples. Rows and columns are clustered using the complete linkage matrix according to Euclidean metric. S14 and S0 clearly indicate similar variant frequencies signatures. ix