Science, Faculty of
DSpace
UBCV
Nguyen, Thuy
Master of Science - MSc
University of British Columbia
DNA sequencing gives us insight into how viruses adapt to their host immune systems. Studies of viral populations typically employ deep amplicon sequencing with next-generation reads to capture a detailed sample of genetic variation in a population. The high amount of overlapping sites in a multiple sequence alignment of reads from amplicon sequencing form ideal input for phylogenetic reconstruction, a necessary step for studying evolutionary relations in a population. However, the typical short read lengths of < 600 bp from next generation sequencing technology with the best sequence error rate impose a severe limit on the width of genomic regions for which evolutionary relationships can be analyzed. Shotgun sequencing, in which DNA is fragmented at random positions, is an efficient alternative to amplicon sequencing for covering wider regions of a genome with sufficient depth. Due to the random staggered positions of shotgun reads in a genome, an extremely high percentage of missing data can result in multiple sequence alignment of shotgun sequencing. The absence of sequence homology across the entire set of short reads makes it impossible to reconstruct a phylogenetic tree, limiting the utility of shotgun data for phylogenetic analysis. We developed the Umberjack software pipeline, which employs a 'sliding window' approach to minimize the effect of missing data during phylogenetic reconstruction and obtain evolutionary statistics to detect sites under selection. Using Umberjack to measure a new metric of directional selection I, significant directional selection was detected in treatment-naive HIV populations at sites with previously documented associations with cytotoxic T-lymphocyte (CTL) response. Further, substitutions towards wild-type amino acids were found to occur early within the population's history, but rarely occurred at a site after the appearance of a CTL escape mutation. Measuring the same metric I in drug treated HIV populations, the directional selection due to the constant pressure of drug treatment was much greater than the directional selection from the immune system.
Overcoming Missing Data in Phylogenetic Analysis ofShotgun Sequencing to Detect HIV Adaptation to ImmuneResponsebyThuy NguyenA THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Bioinformatics)The University of British Columbia(Vancouver)August 2016c© Thuy Nguyen, 2016AbstractDNA sequencing gives us insight into how viruses adapt to their host immune sys-tems. Studies of viral populations typically employ deep amplicon sequencing withnext-generation reads to capture a detailed sample of genetic variation in a popu-lation. The high amount of overlapping sites in a multiple sequence alignment ofreads from amplicon sequencing form ideal input for phylogenetic reconstruction,a necessary step for studying evolutionary relations in a population. However, thetypical short read lengths of <600 bp from next generation sequencing technologywith the best sequence error rate impose a severe limit on the width of genomicregions for which evolutionary relationships can be analyzed.Shotgun sequencing, in which DNA is fragmented at random positions, is anefficient alternative to amplicon sequencing for covering wider regions of a genomewith sufficient depth. Due to the random staggered positions of shotgun reads ina genome, an extremely high percentage of missing data can result in multiplesequence alignment of shotgun sequencing. The absence of sequence homologyacross the entire set of short reads makes it impossible to reconstruct a phylogenetictree, limiting the utility of shotgun data for phylogenetic analysis.We developed the Umberjack software pipeline, which employs the ‘slidingwindow’ approach to minimize the effect of missing data during phylogenetic re-construction and obtain evolutionary statistics to detect sites under selection.Using Umberjack to measure a new metric of directional selection I, signifi-cant directional selection was detected in treatment-naive HIV populations at siteswith previously documented associations with cytotoxic T-lymphocyte (CTL) re-sponse. Further, substitutions towards wild-type amino acids were found to occurearly within the population’s history, but rarely occurred at a site after the appear-iiance of a CTL escape mutation. Measuring the same metric I in drug treated HIVpopulations, the directional selection due to the constant pressure of drug treatmentwas much greater than the directional selection from the immune system.iiiPrefaceThe words ‘we’ and ‘our’ throughout the manuscript refers to the work of ThuyNguyen unless otherwise specified.Dr. Art Poon architected the sliding window strategy for phylogenetic analysisof short reads. Dr. Richard Liang wrote the software for generating ancestralselection graphs. Thuy Nguyen implemented and tested Umberjack, and analyzedpatient datasets under the advice of Dr. Art Poon.Aram Karakas, Dr. Rachel McGovern, and other members of the BC Centrefor Excellence in HIV/AIDS prepared and sequenced HIV patient samples.Ethical approval was obtained from Providence Health Care/University of BritishColumbia research ethics board (H04-50276) for the untreated BC patient HOMERcohort and the Netherlands national institutional review board for the Maravirocclinical trial Netherlands cohort.There are no publications at this time.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 HIV Evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1.1 HIV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1.2 HIV Adaptation to Immune System CTL Response . . . . 11.1.3 HIV Adaptation to Drug Treatment . . . . . . . . . . . . 21.2 Phylogenetic Analysis With Short Reads . . . . . . . . . . . . . . 31.2.1 Next Generation Sequencing . . . . . . . . . . . . . . . . 31.2.2 Phylogeny of a Single Population . . . . . . . . . . . . . 41.2.3 Multiple Sequence Alignment . . . . . . . . . . . . . . . 51.2.4 Difficulties in Phylogenetic Reconstruction of HIV . . . . 71.2.5 Phylogenetic Reconstruction Techniques . . . . . . . . . 71.2.6 Tree Error from Missing Data . . . . . . . . . . . . . . . 10v1.2.7 Current Approaches to Phylogenetic Analysis From Miss-ing Data . . . . . . . . . . . . . . . . . . . . . . . . . . . 111.2.8 New Developments to Overcome Missing Data . . . . . . 131.3 Detecting Selection . . . . . . . . . . . . . . . . . . . . . . . . . 141.3.1 Importance of Selection and Evolutionary History . . . . 141.3.2 Evolutionary Statistics . . . . . . . . . . . . . . . . . . . 142 Phylogenetic Analysis With Sliding Windows . . . . . . . . . . . . . 192.1 Proposal of Sliding Windows to Overcome Missing Data . . . . . 192.2 Software Details . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.2.1 Umberjack Use Cases . . . . . . . . . . . . . . . . . . . 202.2.2 Umberjack Computational Requirements . . . . . . . . . 202.2.3 Umberjack Availability . . . . . . . . . . . . . . . . . . . 212.3 Umberjack Implementation . . . . . . . . . . . . . . . . . . . . . 212.3.1 Read Processing . . . . . . . . . . . . . . . . . . . . . . 212.3.2 Window Extraction . . . . . . . . . . . . . . . . . . . . . 232.3.3 Phylogenetic Reconstruction . . . . . . . . . . . . . . . . 242.3.4 Reconstructing Ancestral Sequences . . . . . . . . . . . . 252.3.5 Evolutionary Statistics . . . . . . . . . . . . . . . . . . . 262.4 Method for Validating on Simulated Data . . . . . . . . . . . . . 312.4.1 Selection of Simulation Parameters . . . . . . . . . . . . 312.4.2 Simulated Topology Under Selection . . . . . . . . . . . 322.4.3 Simulated Recombination . . . . . . . . . . . . . . . . . 322.4.4 Simulated Genomes . . . . . . . . . . . . . . . . . . . . 322.4.5 Simulated Shotgun Reads . . . . . . . . . . . . . . . . . 332.4.6 Umberjack settings . . . . . . . . . . . . . . . . . . . . . 332.4.7 Measuring Umberjack dN−dS Accuracy . . . . . . . . . 342.4.8 Measuring Tree Topology Accuracy . . . . . . . . . . . . 342.5 Determining Features that Affect Umberjack Accuracy . . . . . . 352.5.1 Features Considered For Predicting Umberjack Accuracy . 352.5.2 Univariate Correlation Ranking . . . . . . . . . . . . . . 372.5.3 Regression Random Forest . . . . . . . . . . . . . . . . . 372.5.4 Backwards Feature Selection in Regression Random Forest 39vi2.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 392.6.1 Weak Linear Relationship Between Features and ∆ . . . . 392.6.2 Advantages of Regression Random Forests . . . . . . . . 402.6.3 Features Affecting Umberjack Accuracy . . . . . . . . . . 412.6.4 Effect of Recombination . . . . . . . . . . . . . . . . . . 412.6.5 Umberjack Accuracy . . . . . . . . . . . . . . . . . . . . 423 Within-host HIV Evolution . . . . . . . . . . . . . . . . . . . . . . . 473.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.2.1 Dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.2.2 Umberjack Processing . . . . . . . . . . . . . . . . . . . 493.2.3 Umberjack Cleaning . . . . . . . . . . . . . . . . . . . . 503.2.4 Finding CTL Response Associated Amino Acid Polymor-phisms . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.2.5 dN−dS Analysis . . . . . . . . . . . . . . . . . . . . . . 523.2.6 I Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 523.2.7 Timing Infection . . . . . . . . . . . . . . . . . . . . . . 523.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 533.3.1 Genetic Variation Across Sites . . . . . . . . . . . . . . . 533.3.2 Most Subjects in Chronic Stage of Infection . . . . . . . . 543.3.3 Evolving Sites Associated with CTL Response . . . . . . 553.3.4 Host HLA Genotype Drives Diversifying Selection . . . . 563.3.5 Evidence of Selective Sweeps . . . . . . . . . . . . . . . 593.3.6 Reversion to Wild Type . . . . . . . . . . . . . . . . . . . 643.4 Within-host HIV Evolution During Drug Treatment . . . . . . . . 683.4.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . 683.4.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 693.5 Comparing Methods for Detecting Selection . . . . . . . . . . . . 704 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 734.1 Future Directions . . . . . . . . . . . . . . . . . . . . . . . . . . 74Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76viiA Supporting Materials . . . . . . . . . . . . . . . . . . . . . . . . . . 88A.1 Simulated Datasets For Predicting Umberjack Accuracy . . . . . . 88A.2 Untreated Patient Metadata . . . . . . . . . . . . . . . . . . . . . 102A.3 HIV Selective Sweeps in Untreated Patients . . . . . . . . . . . . 108A.4 Effect of Recombination on Umberjack Error . . . . . . . . . . . 108viiiList of TablesTable 2.1 Spearman Correlation Between Feature and ∆ . . . . . . . . . 40Table 2.2 Feature Ranking for Final Random Forest Predicting Umber-jack Accuracy. . . . . . . . . . . . . . . . . . . . . . . . . . . 42Table 3.1 Empirical dN−dS Datasets . . . . . . . . . . . . . . . . . . . 51Table A.1 Umberjack Simulated Dataset Parameter Ranges . . . . . . . . 88Table A.2 Umberjack Simulated Data Parameters . . . . . . . . . . . . . 92Table A.3 Umberjack Recombination Simulated Data Parameters . . . . . 94Table A.4 Untreated Patient Datasets . . . . . . . . . . . . . . . . . . . . 103Table A.5 Untreated Patient Dataset Sequencing Statistics . . . . . . . . 104ixList of FiguresFigure 1.1 Example of Phylogenetic Trees . . . . . . . . . . . . . . . . . 5Figure 1.2 Example of Multiple Sequence Alignment . . . . . . . . . . . 9Figure 1.3 Long Branch Attraction Due to Missing Data . . . . . . . . . 10Figure 1.4 Multiple Sequence Alignment of Typical Shotgun SequencingReads . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13Figure 2.1 Missing Sequence Homology in Shotgun Reads . . . . . . . . 20Figure 2.2 Effect of Recombination on ∆ . . . . . . . . . . . . . . . . . 44Figure 2.3 Umberjack Inferred dN−dS vs Expected dN−dS From Sim-ulated Data . . . . . . . . . . . . . . . . . . . . . . . . . . . 45Figure 2.4 True Umberjack Inaccuracy Versus Random Forest Predictionof Umberjack Inaccuracy . . . . . . . . . . . . . . . . . . . . 46Figure 3.1 Summary of Amino Acid Entropy Along HIV env . . . . . . . 54Figure 3.2 Estimated Durations of Infection at Baseline by Molecular ClockAnalysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55Figure 3.3 Shifts in Matched Non-escape and Escape Mutation Frequen-cies Over Time in HIV nef . . . . . . . . . . . . . . . . . . . 57Figure 3.4 Boxplots of dN− dS at HLA Matched Sites and UnmatchedSites by Gene . . . . . . . . . . . . . . . . . . . . . . . . . . 58Figure 3.5 Site dN − dS per Sample at Sites Associated with CTL Re-sponse for Patient’s HLA . . . . . . . . . . . . . . . . . . . . 60Figure 3.6 Substitutions Tend to Map to Tips of the Within-host Phylogeny 61xFigure 3.7 Nonsynonymous Substitutions Map to Tips of the PhylogenyMore Often Than Synonymous Substitutions . . . . . . . . . 62Figure 3.8 Distributions of Nonsynonymous and Synonymous Substitu-tions in the Tree are More Similar at Sites with Known CTLAssociations . . . . . . . . . . . . . . . . . . . . . . . . . . 63Figure 3.9 Nonsynonymous Substitutions to HLA-matched Amino AcidsTend to Map Deeper in Within-host Phylogenies . . . . . . . 65Figure 3.10 Nonsynonymous Substitutions to Wild-type Residues OccurDeeper in Phylogenies . . . . . . . . . . . . . . . . . . . . . 67Figure 3.11 Violin Plot Distributions of IN and IS for Site Substitutionsfrom Matched Escape Towards Unmatched Wild Type . . . . 68Figure 3.12 Violin and Boxplot of I Statistic of Viral Populations Treatedwith Maraviroc . . . . . . . . . . . . . . . . . . . . . . . . . 70Figure 3.13 Onset of Directional Selection . . . . . . . . . . . . . . . . . 72Figure A.1 ∆ vs Window Total Breakpoint Ratio . . . . . . . . . . . . . . 95Figure A.2 ∆ vs Window-Site Unambiguous Codon Rate . . . . . . . . . 95Figure A.2 ∆ vs Window-Site Amino Acid Depth . . . . . . . . . . . . . 96Figure A.3 ∆ vs Window-Site Substitutions . . . . . . . . . . . . . . . . 96Figure A.3 ∆ vs Window-Site Nonsynonymous Substitutions . . . . . . . 97Figure A.4 ∆ vs Window-Site Synonymous Substitutions . . . . . . . . . 97Figure A.4 ∆ vs Window-Site Expected Nonsynonymous Substitutions PerBranch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98Figure A.5 ∆ vs Window-Site Expected Synonymous Substitutions PerBranch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98Figure A.5 ∆ vs Normalized Window Tree Length . . . . . . . . . . . . . 99Figure A.6 ∆ vs Normalized WRF . . . . . . . . . . . . . . . . . . . . . 99Figure A.6 ∆ vs Normalized Polytomies . . . . . . . . . . . . . . . . . . 100Figure A.7 ∆ vs Codon Distribution P-value . . . . . . . . . . . . . . . . 100Figure A.7 ∆ vs True Site Codon Entropy . . . . . . . . . . . . . . . . . 101Figure A.8 ∆ vs Window-Site Codon Entropy . . . . . . . . . . . . . . . 101Figure A.8 ∆ vs Window-Site Codon Entropy . . . . . . . . . . . . . . . 102xiFigure A.9 Phylogenetic Tree of GAG Consensus Sequence of Samplesfrom Untreated Patients . . . . . . . . . . . . . . . . . . . . . 105Figure A.10 Phylogenetic Tree of NEF Consensus Sequence of Samplesfrom Untreated Patients . . . . . . . . . . . . . . . . . . . . . 106Figure A.11 Phylogenetic Tree of ENV Consensus Sequence of Samplesfrom Untreated Patients . . . . . . . . . . . . . . . . . . . . . 107Figure A.12 Cleaned Umberjack Estimate of Site dN−dS vs True dN−dS 108Figure A.13 Violin Plots Summarizing the Distributions of Site IN− IS . . 109Figure A.14 Effect of Recombination on ∆ . . . . . . . . . . . . . . . . . 110xiiAcknowledgmentsI would like to thank my supervisor Dr. Art Poon for his unfailing patience andtireless efforts in guiding me through this thesis.This work was supported by a CIHR operating grant awarded to Dr. ArtPoon (HOP 111406), by an administrative supplement to NIH NIDA grant R01-DA011591, and by a grant from Genome Canada (Genomics and PersonalizedHealth).xiiiChapter 1Introduction1.1 HIV Evolution1.1.1 HIVHuman Immunodeficiency Virus (HIV) is a double stranded RNA retrovirus thattargets the human immune system. HIV causes the steady decline of CD4 lym-phocytes (white blood cells) through the destruction of infected CD4 cells or stim-ulating cell death pathways for uninfected CD4 cells [21]. Left untreated, HIVinfection will lead to Acquired Immune Deficiency Disorder (AIDS) during whichopportunistic infections can cause mortality [15]. Upon entering a cell, HIV inte-grates its own genome into the genome of the infected cell. Using the cell’s naturalreplication process, it transcribes and translates viral proteins to create progenyviruses that bud out of the infected cell to infect other cells [38]. A full cycle ofreplication occurs approximately once a day for actively replicating HIV.1.1.2 HIV Adaptation to Immune System CTL ResponseThe cytotoxic T-cell (CTL) response is a mechanism from the adaptive immunesystem that destroys infected cells [19, 36]. CTL-mediated vaccines have beenproposed for boosting CTL targeting of HIV infected cells to reduce viral load andthus risk of transmission [36]. In untreated patients, the efficacy of CTL response is1responsible for controlling the viral population after viremia and partially dictatesthe rate of progression to AIDS [70].Viral proteins in an infected cell are degraded by normal cellular processes intopeptides that are specifically bound by human leukocyte antigen (HLA) moleculesand presented on the cell surface, signalling the CTLs to destroy the cell. The HLAmolecules are encoded by alleles at the highly variable Major HistocompatibilityComplex (MHC) class I loci in the human genome. Each HLA allele binds specificviral peptides that are known as CTL epitopes. However, HIV can evade the CTLresponse by developing ‘escape’ mutations within CTL epitopes that prevent itfrom being recognized and bound by the corresponding HLA molecule [36].Following the standard nomenclature [11], we refer to HIV mutations knownto be statistically associated with an HLA allele as ‘matched’ if they occur in a hostcarrying that allele. Matched mutations that are expected to decrease in frequencyin that host are referred to as ‘non-escape’ alleles, and those expected to increaseare ‘escape’ alleles.With a rapid mutation rate of 10−4 mutations/site/day in a genome of size9.5kbp [36], HIV infection rapidly diversifies and accumulates genetic differenceswithin each host [96]. For an effective CTL-based vaccine, the vaccine epitopesused to induce CTL response must be designed such that mutations in the epitoperesult in viruses with poor survival fitness or replicative capacity. Detecting siteswhere CTL imparts a directional selective force on HIV can shed light on how todesign vaccine epitopes to prevent evasion of CTL response.We developed the software package Umberjack to explore HIV adaptation tothe CTL response and how well it can be detected through evolutionary statistics.1.1.3 HIV Adaptation to Drug TreatmentWild-type HIV-1 enters host cells by binding to CCR5 co-receptors [18]. Thedrug Maraviroc inhibits HIV entry into host cells by binding to the host CCR5 co-receptor, which alters the 3D conformational structure such that the HIV envelopeglycoprotein can no longer attach [18, 22]. Viruses with mutations that allow it tobind CXCR4 co-receptors instead of CCR5 proliferate under Maraviroc treatment,rapidly rendering the drug ineffective. In over half of untreated HIV1-B infected2patients, viral populations switch from CCR5 to CXCR4 co-receptor usage, fol-lowed by patient progression towards AIDS [72]. Clinical guidelines typicallyrequire that patient HIV populations be tested for prior CXCR4 co-receptor usagebefore the administration of Maraviroc, and only recommend that the drug be usedas an alternative third agent in combination therapy [74].The portion of the HIV genome responsible for co-receptor binding is predom-inantly encoded in the V3 region of the env gene. There is quite a bit of variancein conditions that lead to co-receptor usage. The few rules-based methods pre-dicting co-receptor usage based on mutations at specific sites have high specificitybut low sensitivity, thus the exact sites under directional selection and sequence ofmutations can not easily be inferred a priori [50].We employed the Umberjack software once again to profile within-host evo-lution in four patients that developed resistance mutations within days of usingMaraviroc. The efficacy of Umberjack estimates as predictors for directional se-lection were compared between the untreated patient dataset and the treated patientdataset.1.2 Phylogenetic Analysis With Short Reads1.2.1 Next Generation SequencingNext-generation sequencing (NGS) is the massive parallelization of nucleic acidsequencing reactions that allows millions of DNA templates to be processed si-multaneously. NGS amplicon sequencing has become a popular method of captur-ing genetic variation in pathogen populations due to its ability to detect variants atminority frequencies by producing a separate sequence for up to thousands of indi-viduals in a population. Amplicons are DNA fragments isolated from a particularregion of the genome and replicated using specific primers. Prior to NGS, viralpopulations were sequenced with Sanger bulk sequencing, which yielded a singleconsensus sequence for an entire population where different polymorphisms fromdifferent viruses would be indicated as a mixture of bases at a site [54]. AlthoughSanger sequencing yields longer reads (∼ 1kbp), the consensus loses the evolution-ary relationship between viruses since mutations cannot be attributed to any spe-3cific virus. Using NGS to obtain deep coverage of an amplicon can yield reads thatcover the isolated region end-to-end, providing ideal input for phylogenetic analy-sis, the study of evolutionary relationships [84]. However, the typical read lengthswith current NGS technologies with acceptable error rates remain short (<600 bp),imposing a severe limit on the width of genomic regions for which variation canbe observed. Pyrosequencing, an NGS technology can produce close to 1kbp readshave error rates of higher than 1% and over 3% error rate in homopolymer regionswith repeated nucleotides [62]. Third generation reads, such as Single-MoleculeSequencing in Real Time (SMRT) sequencing [25] and nanopore sequencing [7]can generate very long reads up to 25kbp long; however, they suffer from sequenc-ing error rates of more than 15% [4].Shotgun sequencing, another NGS method in which reads are generated by ran-domly fragmenting long DNA templates, allows efficient coverage of wide regionsof a genome. However, only a fraction of the total reads will cover any particularregion of the genome. The absence of sequence homology across the entire set ofshort reads makes it extremely difficult to reconstruct a phylogenetic tree, limitingthe utility of shotgun data for phylogenetic analysis.1.2.2 Phylogeny of a Single PopulationA phylogeny of a single population displays the ancestry of individuals sampledfrom the population in an acyclic graph. Phylogenies are also referred to as trees.External nodes represent observed individuals within the population, and internalnodes represent ancestral individuals that must be inferred since they no longerexist. In Figure 1.1, the external nodes are A, B, C, D and internal nodes areE, F, G. An edge connecting two nodes, also known as a branch, indicates anancestor-descendent relationship between individuals. A weighted edge indicatesthe amount of evolution or time between the individuals.In a rooted tree, the ancestor of each individual can be inferred, and the rootdenotes the common ancestor of all individuals in the population. In an unrootedtree, the most common ancestor and thus root node are unspecified (Figure 1.1).4(a) (b)Figure 1.1: Example of phylogenetic trees. Nodes A, B, C, D are individualsobserved within the population. Nodes E, F, G are inferred ancestralindividuals. Node G represents the common ancestor of A, B, C, D. (a)Rooted phylogenetic tree. Common ancestor G is known. (b) Unrootedphylogenetic tree. Common ancestor G is unknown.1.2.3 Multiple Sequence AlignmentAs a precursor to building a tree, the genetic differences between each individ-ual must be known. Genomic sites can be compared between individuals whenthe sequences are stacked in a multiple sequence alignment where each row is thesequence of an individual, and each column is a genomic site (Figure 1.2). Thesequences are vertically aligned by maximizing the similarity between each se-quence at each column. Similar sequences are inferred to be homologous, that is,descended from the same ancestor. It is possible that similar sequences actuallydescended from separate ancestors but share similarity due to convergent evolutionfrom similar evolutionary pressures or due to chance.A multiple sequence alignment requires comparing each sequence against allother sequences resulting in O(N2L2) complexity, where N is the number of se-quences and L is the length of the sequence [57]. One heuristic towards multiplesequence alignment is to align each sequence against a reference. In algorithms5where the reference is hashed, such as BWA-mem, pairwise alignment against areference is typically much faster with complexities of roughly O(NL/k), wherek is the size of a word within the sequence to look up in the reference [52]. Thisheuristic can be a preferred solution for multiple sequence alignment of NGS li-braries, which can contain millions of reads. Using the pairwise alignment po-sitions with respect to the reference, we can obtain ‘pseudo‘ multiple sequencealignments in which each column represents a position in the reference sequence.For example, if a reference sequence is AAA and read sequences are AGG andAGT, the alignments against the reference will be:Reference AAA Reference AAARead1 ACG Read2 AGTUsing the positions with respect to the reference, we obtain the ‘pseudo‘ multiplesequence alignment of the reads:Read1 ACGRead2 AGTHowever, pairwise alignments against a reference lose information regarding themultiple sequence alignment of insertions with respect to the reference. For ex-ample, if a reference sequence is AAA and read sequences are ACGTAA andAGGAA, pairwise alignment against the reference will be:Reference A---AA Reference A--AARead1 ACGTAA Read2 AGGAAThe alignments against the reference do not yield any information regarding thepositions of the inserted bases with respect to other inserted bases. Thus, we cannot discern the best alignment from the potential multiple sequence alignments ofthe inserted bases:Read1 ACGTAA Read1 ACGTAA Read1 ACGTAARead2 AGG-AA Read2 AG-GAA Read2 A-GGAA61.2.4 Difficulties in Phylogenetic Reconstruction of HIVOne of the ways in which HIV populations maintain genetic diversity is through re-combination. HIV experiences recombination when its reverse transcriptase jumpsbetween strands during RNA replication. In order for the recombination to be ge-netically noticeable, either strand must be genetically different from the other inthe sections that swap. This only happens in ‘hybrid’ progeny viruses created bycombining an RNA strand from 2 different strains that infect the same cell [10].In HIV, the rate of genetically noticeable recombinations is reported to be 10−5recombinations/bp/day [77]. Recombination affects reconstruction of evolutionaryreconstruction by altering the true topology on either side of a genomic breakpoint.Thus, a tree reconstructed from an alignment containing a breakpoint will be dif-ferent from the true trees on either side of the breakpoint [94].Further, the HIV effective population size can be quite large. The number ofgenetically unique replicating viruses in an HIV population can reach 103 to 105within a patient [63, 79], requiring the same number of overlapping sequences torepresent each virus in a multiple sequence alignment.1.2.5 Phylogenetic Reconstruction TechniquesThe topology and branch lengths of a phylogenetic tree can be reconstructed byseveral algorithms such as neighbour joining [93] and maximum likelihood [28].Although time consuming, maximum likelihood approaches yield far more accu-rate trees than neighbour joining [17], especially when there is missing sequencedata [104].Neighbour JoiningNeighbour joining is a distance-based approach to phylogenetic reconstruction.It employs a dynamic programming algorithm in which pairs of individuals arerecursively clustered together to form a new subtree until all observed individualsand inferred ancestral individuals are included in a single tree. In each iterationthe pair of unconnected individuals with the minimum pairwise genetic distanceis connected to a new node representing their inferred common ancestor. Branchlengths from the new node to either of the pair is set to a function of the pair’s7genetic distance.Genetic and Evolutionary DistanceThe genetic distance is a function of the probability of mutations at each site. Thereare many mathematical models for calculating genetic distance that take into ac-count various biological biases of mutation. The general time reversible (GTR)nucleotide mutation model estimates a separate parameter for the symmetric mu-tation rates between each pair of nucleotides [102]. The Jukes Cantor model is aspecial subset of the GTR which assumes each nucleotide has the same probabilityand the mutation between each nucleotide occurs at the same rate [44]:P(i|k, t) =14 + 34 e−tµ i = k14 − 14 e−tµ i 6= k(1.1)where P(i|k, t) is the probability that nucleotide k mutates to nucleotide i after timet. µ is the nucleotide mutation rate. The Jukes Cantor model is too simplistic torealistically represent HIV evolution, but it allows for quicker calculations in initialstages of phylogenetic reconstruction algorithms.Maximum LikelihoodThe maximum likelihood approach attempts to selects the tree T = Tmax that maxi-mizes the likelihood of tree T given multiple sequence alignment D. The likelihoodof T given D is equivalent in terminology to the probability of D given T , P(D|T ).Given tree T , the likelihood of T given the sequences at site u (also knownas the site-likelihood) can be calculated recursively by traversing the tree in post-order (tips-first, where children are visited before parents). The site-likelihood ofthe subtree rooted at a node is a function of the site-likelihood of the subtrees rootedby its children.At site u, the probability that node k has sequence character sk is the sum of theprobabilities of mutations from sk to all possible characters at the child nodes.8P(k|sk) = ∑sleft,srightP(sleft|sk, tleft)P(left|sleft)P(sright|sk, tright)P(right|sright) (1.2)where P(k|sk) is the probability that node k has sequence character sk at site u.Node k has children left and right with sequence character states sleft and srightrespectively. P(sleft|sk, tleft) is the probability that character sk mutates to childcharacter sleft after the time elapsed on branch tleft connecting node k and childleft. The mutation probability typically follows a mutational model as described inSection 1.2.5.The site-likelihood of the entire tree is the sum of the probabilities that the rootnode has sequence character sk for every possible sk. Since each site is consideredindependent of other sites, the likelihood of T given D can be expressed as theproduct of the likelihoods of T at each site.The search over all possible trees to find the tree Tmax that maximizes the likeli-hood can be computationally expensive. Search times can be shortened with greedyapproaches that stop searching when newer trees no longer increase likelihood ac-cording to some threshold.Figure 1.2: Example of multiple sequence alignment. Each row representsthe genetic sequence of an individual in the population. Each column isa position in the organism’s genome.91.2.6 Tree Error from Missing DataMissing data in multiple sequence alignments is handled differently by differ-ent phylogenetic reconstruction approaches. Neighbour joining ignores sites withmissing data when calculating pairwise genetic distance. Maximum likelihoodsums over the likelihood of all possible characters for missing data [87].Long branch attraction is an issue common to most phylogenetic reconstruc-tion approaches and occurs when there are wide swaths of missing data for somebut not all individuals in a multiple sequence alignment. If there is little overlapin alignment between two sets of individuals, most phylogenetic reconstructionapproaches would force either sets of individuals into separate elongated lineagesregardless of their true relationship (Figure 1.3) [104].Past studies of phylogenetic reconstruction on missing data have shown that thepercentage of missing data was not as important to tree accuracy as whether therewere sites in the multiple sequence alignment with ‘informative’ sequence from allindividuals. Informative sites contain sufficient sequence diversity to distinguishbetween individuals but evolve slowly enough such that mutations can be capturedwithin the sampling timeframe [105].(a) (b)Figure 1.3: Long branch attraction due to missing data. (a) Multiple sequenc-ing alignment with 2 partitions of missing data. (b) Corresponding treewith 2 long lineages.101.2.7 Current Approaches to Phylogenetic Analysis From MissingDataFilling and Trimming AlignmentsIn multiple sequence alignment of genomes between species, unsequenced or bio-logically missing genes can either be filled in with known sequences from databases,or all genomes trimmed such that only sites common to most species remain [13].When multiple populations are shotgun sequenced with separate libraries, inter-population evolutionary relationships can be found through phylogenetic recon-struction of population consensus sequences [54]; however, within-population evo-lutionary relationships will be lost.Amplicon SequencingStudies of viral populations frequently sequence amplicons that are narrow enoughto be covered almost entirely by reads, which reduces sites in the multiple sequencealignment with missing data across individuals [68, 84]. In order to measure evo-lution across a wider range of the genome, overlapping or adjacent amplicons canbe sequenced [107].PCR primers that amplify a particular portion of the genome require regions oneither side of the amplicon that are conserved across the population. As a result,amplicon primer design can be a limiting step for diverse regions of the genomeor when studying communities that may not share similar genomic features acrossspecies or strains [35].Phylogenetic analysis options for shotgun sequencing are minimal as demon-strated by a recent viral population study in which overlapping amplicons acrossthe genome for an HIV population were shotgun sequenced [107]. Even though se-quencing was available across the entire genome, phylogenetic analysis was onlyavailable for the amplicons narrow enough to fit a full read.Inferring Missing DataForeSeqs is newly released software that infers missing data in a multiple sequencealignment and corrects branch lengths on phylogenies made from incomplete data11[16]. It requires a multiple sequence alignment, the corresponding phylogeny, anda user-supplied list of site boundaries delineating partitions of missing data in themultiple sequence alignment. As such, ForeSeqs is more suited for situations inwhich gaps are at the same site for multiple taxa, such as biologically missinggenes.ForeSeqs uses partitions of the alignment with complete data to scale branchlengths for taxa with missing data. Unknown sequence characters are inferred bysimulating evolution down the tree or through marginal likelihood ancestral recon-struction in which the sequence at a node is selected to maximize the likelihood ofthe connecting node.Since ForeSeqs does not correct the topology of the input tree, any errors fromthe initial phylogenetic reconstruction due to missing data will have knock on ef-fects on the reconstruction of missing data. This proves problematic in its applica-tion to NGS shotgun libraries. An initial phylogenetic reconstruction from shotgunmultiple sequence alignments would tend towards a star topology since each readshares little sequence homology with most other reads (Figure 1.4). For example,the multiple sequence alignment of shotgun libraries analyzed in Chapter 3 rangedbetween 70-85% missing data and none of the sites contained sequence from >50% of reads.Inferring Haplotypes From Single Nucleotide PolymorphismsA haplotype is the sequence along a single DNA strand. By aligning reads to a ref-erence genome, the nucleotide variants with respect to the reference can be used todeconvolute the separate haplotypes within a population. Through reconstructionof the haplotypes, we fill in the missing gaps within a multiple sequence alignmentof shotgun reads. Typically, the variants are modelled as a mixture model, and thevariants belonging to the same haplotype are linked together by clustering. Exist-ing software accomplish haplotyping for cancer genomics in which tumour cellsare made up of multiple subclonal populations, as well as for viral populations[49, 85, 106]. However, haplotype reconstruction software are limited to detecting101 to 102 haplotypes within a population, whereas the number of unique replicat-ing HIV haplotypes can reach 103 to 105 within a patient [63, 79].12Figure 1.4: Subset of a multiple sequence alignment of typical shotgun se-quencing reads from a simulated viral population. Colored blocks rep-resent nucleotide bases within a read. White spaces indicate genomicpositions uncovered by the read. This multiple sequence alignment con-tains over 80% missing bases.1.2.8 New Developments to Overcome Missing DataWe developed the software pipeline Umberjack to address the problem of excessmissing data in shotgun sequencing for phylogenetic analysis. Umberjack uses asliding window technique to obtain windows of alignment with sufficient sequenceoverlap and homology at all sites to reconstruct a tree (Figure 2.1). A windowrefers to a slice of a multiple sequence alignment containing only a subset of se-quences and a subset of contiguous sites. Using the reconstructed evolutionaryhistory from each window, Umberjack detects sites under selection by aggregatingover evolutionary statistics from overlapping windows. Chapter 2 details Umber-jack implementation and its ability to accurately output evolutionary statistics fromshotgun sequencing data.131.3 Detecting Selection1.3.1 Importance of Selection and Evolutionary HistorySelection refers to the deterministic change in frequency of individuals in a popu-lation due to their traits that improve or diminish their capacity to reproduce andsurvive in their environment. At the molecular level, selection operates correspond-ingly to change the frequency of genetic variants that code for the traits. Selectiongenerally acts more strongly at genomic positions that have a high functional im-pact. Thus, detecting genetic signatures of selection can provide the basis for di-agnostic tools and identify functionally important sites. For example, potential ge-netic markers for drug resistance in Mycobacterium tuberculosis have been foundby calculating the strength of selection [27].Leveraging the evolutionary history of pathogens can identify potential targetsfor vaccines and treatments. The immune system and drug treatment act as se-lective forces on pathogens. Antibody-based vaccines induce the immune systemto produce proteins called antibodies that bind and neutralize foreign moleculescalled antigens. In a study of untreated patients that were naturally able to controltheir Human Immunodeficiency Virus (HIV) infections, potential vaccine antigenswere discovered by tracing the series of HIV mutations that subsequently stimu-lated the production of antibodies that neutralized a broad range of HIV variants[55].Although taking frequent longitudinal samples is the most accurate way of de-termining evolutionary history, resource constraints make this infeasible. Phylo-genetic and ancestral reconstruction seek to fill in the gaps in evolutionary historywhen samples are unavailable.1.3.2 Evolutionary StatisticsUmberjack generates several statistics that describe the rate of selection at a codonsite, including dN−dS, I, and Shannon entropy. The first two statistics are phylo-genetically informed and take into account the evolutionary history of a population,leading to less confounding between demographic events and actual selection.Selection statistics that are frequency based, such as Shannon Entropy, are af-14fected by demographic events. Fixation of genetic variants may be due to a popula-tion size bottleneck during which a single founder individual passes on its geneticvariants to its progeny. This founder effect will continue if there is no selectionand there has been insufficient time for random mutations to remove or reduce thegenetic variant from the population. Further, a growing population that is undergo-ing neutral evolution with no selection will experience increasing genetic diversityover time simply due to random mutations [91].Most of the statistics label sites as under diversifying or purifying selection.With regards to within-population evolution at the molecular level, diversifying se-lection indicates a proliferation of polymorphisms. Purifying selection indicatesthe removal of variants unfit for survival in the population’s environment. Direc-tional selection indicates only genotypes coding for a single specific trait survive.dN − dSDiversifying or purifying selection at a site can be measured by the dN−dS statis-tic, where dN is the rate of nonsynonymous substitutions and dS is the rate ofsynonymous substitutions. A site is diversifying if dN − dS > 0, purifying ifdN−dS < 0, and neutral if dN−dS = 0:dN−dS = NE[N]− SE[S]The total nonsynonymous substitutions N and total synonymous substitutions Sare normalized by their respective expected numbers, E[N] and E[S], based on thegenetic code [80]. For example, all possible mutations at each nucleotide positionof codon ATG leads to a different amino acid translation. Thus, each nucleotideposition in ATG is expected to be 100% nonsynonymous. The dN−dS statistic waschosen instead of the typical dN/dS ratio [34], since the latter is more numericallyunstable when dS is small.IWe developed a new metric of site selection inspired by a statistical test of neutral-ity proposed by Fu and Li [32]. Their test is based on the distribution of substitutionevents that have been mapped to a phylogeny. We assume that substitutions on the15tips of the phylogeny are more recent in time. Under purifying selection, substi-tutions should tend to be observed on the tips because they represent lineages thathave not yet been removed by selection. Conversely, if a greater proportion of sub-stitutions are mapped to internal branches of the tree, then the variants from thoseancestral lineages have been allowed to proliferate under diversifying selection.A simple metric of site diversifying and purifying selection can therefore becalculated from the numbers of substitutions at a site that map to internal (SI) andexternal branches (SE), normalized by their respective branch lengths tI , tE . Thisfollows a similar approach as Fu and Li; however, we use branch lengths for nor-malization instead of the Fu and Li mathematical simplifications based on binaryphylogenies where each non-root node has at most two child nodes [32]. The met-ric I ranges in [0, 1], where I = 0.5 indicates neutral evolution, I > 0.5 indicatesdiversifying selection, and I < 0.5 indicates purifying selection.I =SItISItI+ SEtEAncestral lineages that maintain a specific amino acid provide evidence of di-rectional selection at the amino acid level. A nonsynonymous substitution to afavoured amino acid in an ancestral individual followed by synonymous substitu-tions in its descendants could indicate that the lineage is undergoing directionalselection. An existing amino acid preserved through synonymous substitutionsalong a lineage except for nonsynonymous substitutions at the tips indicates thatpurifying selection has not yet had time to excise the spurious mutations.Breaking down substitutions as nonsynonymous and synonymous, we can mea-sure directional selection at the amino acid level by calculating the difference be-tween nonsynonymous substitutions and synonymous substitutions at internal ver-sus external branches, IN−IS. IN−IS ranges in [-1, +1], where IN−IS =±1 impliesstrong directional selection and IN− IS ≈ 0 implies no or weak selection.IN− IS =SINtISINtI+ SENtE−SIStISIStI+ SEStEwhere SIN is the number of nonsynonymous substitutions mapping to internal16branches, SIS is the number of synonymous substitutions mapping to internal branches,SEN is the number of nonsynonymous substitutions mapping to external branches,SES is the number of synonymous substitutions mapping to external branches.It is possible that a site indicated as under selection by IN − IS is actually neu-trally evolving but appears as selected because it is located in a selective sweep inwhich sites adjacent to selected sites are dragged to fixation. Although a favouredmutation imparting a selective advantage will be passed to progeny and increasein frequency, the entire genomic backbone containing the mutation may not bepassed due to recombination. RNA sections are swapped between strands duringrecombination, breaking up the linkage of genomic sequence on either side of thebreakpoint. Since sites close in proximity are less likely to segregate to separatestrands, the genomic sequence adjacent to a favoured mutation will be passed toprogeny and increase in frequency simply due to genomic linkage. Variants at neu-tral sites in a selective sweep will remain fixed in frequency until sufficient timehas passed for mutation and recombination to occur [40].EntropyThe distribution of variant frequencies can be used to infer selection without a phy-logeny. For example, high frequency of a variant indicates purifying selection, anda wide distribution of variant frequencies indicates diversifying selection. In addi-tion to avoiding computational efforts of phylogenetic and ancestral reconstruction,calculating variant frequencies can be performed on any type of sequencing.Shannon entropy is a commonly used summary metric to indicate the distribu-tion of variant frequencies at a genomic site. Shannon entropy ranges from [0, ∞),where 0 entropy indicates all individuals have the same variant at a site:H(X) =n∑i=1P(xi)log2P(xi)where H(X) is the Shannon Entropy of site X . xi is ith variant at site X . P(xi) is theprobability of observing variant i at site X . n is the total number of unique variantsat site X .Metric entropy is a normalized version of Shannon entropy with values in [0,171], which allows comparison of entropy between sites that may be covered by adifferent number of reads:M(X) =∑ni=1 P(xi)log2P(xi)log2nwhere M(X) is the Metric Entropy at site X .18Chapter 2Phylogenetic Analysis WithSliding Windows2.1 Proposal of Sliding Windows to Overcome MissingDataParameters that require a phylogeny for accurate estimation, such as rates of sub-stitution, can be progressively estimated from a series of phylogenies along analignment. This approach is similar to smoothing or ‘sliding window’ methods forestimating a summary statistic over a series of data points, except that data withina window is mapped into tree space from which the statistics are derived. For ex-ample, sliding windows of phylogenies along alignments is a popular approach fordetecting recombination [66, 94].The Umberjack software pipeline applies the sliding window approach towardsphylogenetic analysis of shotgun reads. To reduce missing data used in phyloge-netic reconstruction, windows of multiple sequence alignment are selected suchthat all sequences have sufficient overlap as defined by the user (Figure 2.1). Toevaluate evolutionary statistics at each site along the genome, Umberjack aggre-gates the statistics from overlapping windows. Although the trees from each win-dow will vary depending on the sequence diversity and individuals included withineach window, averaging the site evolutionary statistics over each window replicate19has been tested on simulated data to be robust as long as there is sufficient sequencediversity to provide a phylogenetic signal in the windows (Section 2.6.5).Figure 2.1: Depiction of multiple sequence alignment of shotgun reads. Ref-erence sequence in green. Shotgun reads in grey. Slice of alignmentin orange. White spaces indicate missing sequence homology betweenreads.2.2 Software Details2.2.1 Umberjack Use CasesUmberjack is primarily designed to efficiently process extremely high coveragenext-generation sequencing of a population with a reference genome with minimalrepetition and known open reading frames. As such, Umberjack is ideal for study-ing populations of RNA viruses such as HIV, as they typically have little genomeduplication but display a large amount of diversity that requires deep coverage tocapture [48, 73, 98].2.2.2 Umberjack Computational RequirementsUmberjack is Message Passing Interface (MPI) enabled, allowing it to be run acrossmultiple nodes in a cluster. Although Umberjack is written in Python and R, whichare supported on multiple platforms, it has only been tested on Linux environments.Since Umberjack’s aim is to rapidly process many overlapping windows ofphylogenetic analysis simultaneously, its implementation favours speed over ac-curacy. For example, processing 70 windows for a single sample library for a20genomic region of 3kbp sequenced at 7000x with MiSeq Nextera (shotgun) readstakes roughly 17 minutes on a cluster with 50 CPU cores spread across 4 nodes.2.2.3 Umberjack AvailabilityUmberjack and the third party software it is packaged with are publicly availablefor download under a GNU General Public License (GPL): https://github.com/cfe-lab/Umberjack.2.3 Umberjack ImplementationThe Umberjack pipeline processes reads to form a ‘pseudo’ multiple sequencealignment and slices windows out of the alignment. Each window is sent for phylo-genetic reconstruction, ancestral sequence reconstruction, and calculation of evo-lutionary statistics. The evolutionary statistics from each window are aggregatedinto per-site evolutionary statistics. The pipeline is further detailed in the followingsections.2.3.1 Read ProcessingAs input, Umberjack requires pairwise read alignments against a reference in SAMformat [53] or a multiple sequence alignment in Fasta format. Umberjack uses thealignment information specified in the SAM files (positions with respect to thereference, and aligned characters and qualities) to generate a ‘pseudo’ multiplesequence alignment of the reads (Section 1.2.3).One technique for ensuring that majority consensus insertions with respect toa reference genome are included in the ‘pseudo’ multiple sequence alignment is toemploy an iterative approach to pairwise alignment against a reference, using theprevious iteration’s consensus as the reference in the next iteration (Section 3.2.2).The iterations stop when the consensus no longer changes from the previous iter-ation. We used this iterative approach when applying Umberjack to our datasets.However, Umberjack itself does not perform any alignments.The majority of evolutionary statistics reported by Umberjack are substitutionbased. Insertions that do not occur within the majority of the population will havegaps in the multiple sequence alignment for all the individuals lacking the inser-21tion; that is the majority of the population will have missing data at the insertionsites. If there are no substitutions within the insertions, the substitution-based evo-lutionary statistics will not be able to discern evolution at those sites. However,detecting the total number and locations of insertions can be a way to detect evolu-tion and is worthy of exploring in future releases of the software. Further, we notethat alignment errors from indels or otherwise are a potential bias in Umberjackestimates of evolution in that they would show more diversity where there is none.When paired-end or mate-pair reads are used, Umberjack merges the matesinto a single read using the mate alignment positions reported in the SAM file.Both paired-end and mate-pair reads are reads that are sequenced from both endsof a DNA fragment. In the case of paired-end reads, the DNA is sequenced inwardsstarting from the ends. Mate-pair reads are sequenced outwards starting from a po-sition from within the DNA fragment. SAM alignment positions are taken at facevalue without consideration that mates may align to duplicate portions of the refer-ence. This simple and rapid merging procedure is sufficient for RNA virus popu-lations which experience minimal genome repetition. To handle organisms that doexperience genome repetition, users may use their own multiple sequence align-ments as input. Further, reads that contain the exact same sequence are removed,since they could potentially be PCR (Polymerase Chain Reaction) duplicates inwhich the same DNA fragment from the same individual is amplified, yielding noadditional information regarding the population.All bases that do not meet a minimum Phred-like quality score are masked withN’s. A Phred quality score Q is a log transform of the probability of base errordetermined from internal benchmarking by the sequencing technology [9, 26]:Q =−10log10Perrwhere Perr is the probability that a sequenced base is erroneous. Discordant basesfrom overlapping mates are resolved by accepting the base with the highest qualityscore, and are masked with N’s if both mates do not reach the minimum qualityscore, Qmin. Concordant bases from overlapping mates are allowed a lower Phred22quality score before N-masking, Q′min.P(M1 = M2,M1err,M2err) = P(M1 = M2|M1err,M2err)P(M1err)P(M2err)=39P(M1err)P(M2err)Q′min =−10log10(3)+Qmin <= QM1+QM2where P(M1 = M2,M1err,M2err) indicates the probability that mate 1 and mate 2are concordant and erroneous. P(M1 = M2|M1err,M2err) indicates the probabilitythat mate 1 and mate 2 are concordant given that mate 1 and mate 2 are erroneous.P(M1err) indicates the probability that mate 1 is erroneous. Q′min is the new mini-mum quality threshold for concordant overlapping bases. Qmin is the user-definedminimum quality threshold for a base. QM1 is the quality score for mate 1 whichis dictated by the probability that mate 1 is erroneous. There are 3 permutations ofconcordant erroneous overlapping bases and 9 permutations of erroneous overlap-ping bases.The revised quality score cutoff makes the assumption that the probability oferror at each mate is independent of the other mate within a read and that theprobability of error towards each base is the same, which may not necessarily betrue depending on the sequencing technology [95]. However, it provides a roughapproximation of the increase in confidence of concordant overlapping bases thatdoes not require additional user configuration based on sequencing technology.2.3.2 Window ExtractionThe width of windows and the amount of overlap between windows are user con-figurable. Sequences that do not meet a user defined threshold for minimum per-centage of non-N, non-gap bases within the window are excluded from that win-dow. Windows that do not meet user defined thresholds for minimum number ofsequences are excluded from further phylogenetic analysis. Sequencing error hasthe effect of falsely increasing the rate of nonsynonymous mutations, since themajority of mutations are nonsynonymous. We attempt to reduce the amount of se-quencing error making it to the end analysis by masking bases that have low qualityscores (as configured by the user). Thus, bases that are likely to be erroneous will23be treated as missing data instead of influencing the mutation counts.2.3.3 Phylogenetic ReconstructionUmberjack comes packaged with FastTree2 [88] to perform phylogenetic recon-struction. FastTree2 provides a rapid method of phylogenetic reconstruction thathas been benchmarked favourably in terms of speed with minimal drop in accuracycompared to other maximum likelihood phylogenetic reconstruction software [88].First, FastTree2 calculates the pairwise distance between sequences using asimple Jukes Cantor model that takes into account multiple substitutions per site(Equation 1.1). When FastTree2 encounters an N or a gap, it uses the distanceinformation from the rest of the sequence to infer evolutionary distance [87], andweights sequences by the number of positions without missing data. Using theJukes Cantor pairwise genetic distances, FastTree2 generates an initial tree usingneighbour joining (Section 1.2.5).The initial tree topology and branch lengths are further refined using an approx-imate maximum likelihood phylogenetic reconstruction. Instead of finding the treethat yields the maximum likelihood (see Section 1.2.5), FastTree2 uses heuristicsthat reduce search time by only examining a subset of trees. The tree with thehighest likelihood within the subset is selected, but it may not be the maximumlikelihood tree amongst all possible trees. FastTree2 performs up to two sequentialrounds of optimizations for topology and branch length.FastTree2 refines tree topology using nearest neighbour interchange [20]. Travers-ing through the tree in post-order traversal (tips-first, where children are visitedbefore parents), the subtree rooted at the current node is perturbed by swapping thesibling and child nodes. The trees resulting from the swappings are calculated forlikelihood (see Equation 1.2).After optimization of the tree topology, the branch lengths are refined usingthe GTR nucleotide mutation model (Section 1.2.5). Since different sites mightevolve at different speeds, the set of mutation rates at each site are scaled by asite specific evolution scaling factor. Each site specific evolution scaling factoris selected from the category of a discrete 20-category gamma distribution thatmaximizes the likelihood.242.3.4 Reconstructing Ancestral SequencesIn order to calculate dN−dS and I, we need to reconstruct the sequences of ances-tral individuals to count mutations throughout a population’s evolutionary history.Umberjack uses HyPhy to perform maximum likelihood ancestral reconstructionin which the codon sequence of each ancestral individual is inferred from a giventree and multiple sequence alignment [82].HyPhy first estimates the GTR nucleotide mutation model (Section 1.2.5) thatmaximizes the likelihood of the sequences given the tree. Using the fixed GTRmodel, HyPhy estimates the Muse Gaut codon mutation model [76] that maxi-mizes the likelihood of the sequences given the tree (Section 1.2.5). The MuseGaut codon mutation model (MG94) specifies that the probability of a codon mu-tating to another codon depends on the nucleotide mutation rate, the ratio of non-synonymous to synonymous mutations, and whether the mutation is nonsynony-mous or synonymous. Stop codons indicating the end of transcription (TAG, TAA,TGA) are excluded from the codon mutation model and disallowed during ances-tral reconstruction with the reasoning that a premature stop codon would result in adefective organism. Umberjack has a user configurable option to mask stop codonsin window alignments with NNN to allow ancestral reconstruction to proceed.Since the input phylogenetic tree lengths are in units of nucleotide substitutionsper nucleotide site, the tree needs to be scaled to nucleotide substitutions per codonsite using the branch length scaling factor B. B is estimated at the same time as thecodon model via maximum likelihood.The ancestral sequences are reconstructed using a joint maximum likelihoodapproach based on the Pupko dynamic programming algorithm for ancestral re-construction [89] in which the ancestral sequences are selected to maximize thelikelihood of all nodes for a fixed tree.The ancestral character at each site is calculated independently of other sites.For each site, HyPhy traverses each node starting from the external nodes towardsan arbitrarily chosen root node. At each non-root node Lb, it evaluates the like-lihood of every possible character sa at its parent node La. It keeps track of thesequence character sb for the current node Lb that maximizes the likelihood given25that the parent character is sa.P(Lb|La = sa) = P(sb|sa, tb)∏c∈CP(Lc|Lb = sb)where P(Lb|La = sa) denotes the likelihood of node Lb given that the parent nodeLa has character sa. Node Lb has the set of child nodes C. P(sb|sa, tb) denotes theprobability that character sa will mutate to character sb after time tb. The time tb isobtained from the branch length connecting the node Lb to its parent La.The root character is selected to maximize likelihood of the root:P(Lroot |Lroot = sroot) = P(sroot)∏c∈CP(Lc|Lroot = sroot)where P(Lroot |Lroot = sroot) is the likelihood of the root given that it has charactersroot . The root has the set of child nodes C. P(sroot) is the prior probability of thecharacter sroot .The final set of ancestral characters at each site are found by traversing thetree from root towards the tips, and selecting the node character that yields themaximum likelihood for the parent node.2.3.5 Evolutionary StatisticsUmberjack uses a custom HyPhy Batch Language implementation of the SingleLikelihood Ancestor Counting (SLAC) [78, 81] method for calculating dN− dS(Section 1.3.2). HyPhy allows programmers to utilize its functionality via its APIs(application programming interfaces) written in HyPhy Batch Language. AlthoughHyPhy contains an implementation of SLAC, we created a custom SLAC im-plementation for cleaner integration with Umberjack’s workflows, including con-straining tree topologies to those obtained from FastTree2 and ensuring consis-tent clade names. Most importantly, Umberjack’s implementation of SLAC avoidsHyPhy’s built-in resolution of ambiguous codons, which improved accuracy indN−dS estimations from NGS shotgun sequencing data which has more N’s andgaps than the population mixture Sanger sequencing data that HyPhy was origi-nally designed to handle.For each codon site, SLAC calculates dN−dS by counting the nonsynonymous26and synonymous substitutions along branches in the phylogeny, and normalizes thecounts by the expected nonsynonymous and expected synonymous substitutions.dN is the rate of nonsynonymous substitutions per expected nonsynonymoussubstitutions at a codon site. dS is the rate of synonymous substitutions per ex-pected synonymous substitutions at a codon site.dN =NTET [N]dS =STET [S]where NT is the total nonsynonymous substitutions in tree T , and ET [N] is the totalexpected nonsynonymous substitutions T . ST is the total synonymous substitutionsT , and ET [S] is the total expected synonymous substitutions T .The method to calculate dS is analogous to the method for dN. We will focusonly on nonsynonymous substitutions for simplicity.Counting Nonsynonymous SubstitutionsThe nonsynonymous substitutions are counted along the branches in the tree. Ifthe path between the parent and child codon of a branch requires more than onenucleotide mutation, the branch nonsynonymous substitutions are calculated as theaverage nonsynonymous substitutions across all possible paths between the parentand child codon with 1-nucleotide mutations at each step.NT = ∑b∈T1|M| ∑m∈MNm (2.1)where branch b is in tree T . M is the set of paths from the parent to child codon ofbranch b. Each path m consists of steps of 1-nucleotide mutations, where mutationsto stop codons are disallowed. Nm is the total nonsynonymous mutations duringpath m.27For example, if branch b had parent codon AAG and child codon GAA:M ={AAG 1 synonymous−−−−−−−→ AAA 1 nonsynonymous−−−−−−−−−→ GAA,AAG1 nonsynonymous−−−−−−−−−→ GAG 1 synonymous−−−−−−−→ GAA}Nb =1|M|(NAAG→AAA→GAA+NAAG→GAG→GAA)=12(1+1)=1Counting Expected Nonsynonymous Substitutions in a TreeThe expected nonsynonymous substitutions along a tree is weighted by each branchlength:ET [N] =1|T | ∑b∈TtbEb[N] (2.2)where branch b has branch length tb in tree T . Eb[N] is the expected nonsynony-mous substitutions in branch b. |T | is the total length of tree T .Counting Expected Nonsynonymous Substitutions in a BranchThe expected nonsynonymous substitutions along a branch is an average of theexpected nonsynonymous mutations across all paths from the parent to child codon.Each path from the parent to child codon consists of steps of 1-nucleotide mutationswhere stop codons are disallowed.Eb[N] =1|M| ∑m∈MEm[N]where Eb[N] is the expected nonsynonymous substitutions in branch b. M is theset of paths from the parent to child codon. Em[N] is the expected nonsynonymousmutations in path m from the parent codon to child codon.The expected nonsynonymous substitutions in a path, Em[N], is an average of28the expected nonsynonymous substitutions across each codon in the path:Em[N] =1|C| ∑c∈mEc[N]where path m is comprised of the set of codons resulting from 1-nucleotide muta-tion steps from the parent to the child codon, including the parent and child codonthemselves.Counting Expected Nonsynonymous Substitutions in a CodonThe expected nonsynonymous mutations in a codon is the sum of the expectednonsynonymous mutations at each codon position. The expected nonsynonymousmutations at each position ranges in [0, 1]. Thus the total expected nonsynonymousmutations in a codon ranges in [0, 3].Ec[N] =∑i∈cEi[N] =∑i∈c∑j 6=iwi jNi jwhere i is a nucleotide at one of the three nucleotide positions in the codon. wi j isa weight assigned to the mutation i→ j based on the mutation probability obtainedfrom the MG94 codon mutation model. Ni j is the total nonsynonymous mutationsin codon c resulting from mutation i→ j, excluding stop codons.For example, let’s say the parent codon is AAG. Mutation AAG→ TAG is notallowed since TAG is a stop codon. Thus, the expected nonsynonymous mutationsat the first codon position is:Ei[N] =∑j 6=iwi jNi jEfirst A[N] = ∑j 6=AwA jNA j= wACNCAG+wACNGAG= wAC(1)+wAC(1)29The total expected nonsynonymous mutations in codon AAG is:EAAG[N] =∑i∈c∑j 6=iwi jNi j=[wACNCAG+wAGNGAG]+[wACNACG+wATNATG+wAGNAGG]+[wGANAAA+wGCNAAC+wGTNAAT]=[wAC(1)+wAG(1)]+[wAC(1)+wAT(1)+wAG(1)]+[wGA(0)+wGC(1)+wGT(1)]SLAC PowerAs a counting method, SLAC has low power at sites with few substitutions result-ing from low diversity [81]. Although maximum likelihood methods exist that havebeen documented to have more power for populations with fewer sequences or lit-tle diversity, such as Fixed Effects Likelihood (FEL) where dN/dS is a parameterto be estimated via maximum likelihood at each site, those methods were bench-marked to be at least 6 times slower than SLAC for alignments with as little as 100sequences, and at least 10 times slower for alignments with 1000 sequences. Con-sidering viral populations are deep sequenced closer to 10000X and that selectionat each site would be recalculated for every overlapping window, likelihood meth-ods were abandoned for the initial release of the Umberjack pipeline. However,incorporating more powerful selection detection methods may be incorporated infuture releases.Internal Versus External SubstitutionsI (Section 1.3.2) is calculated similarly to dN−dS in that substitutions are countedalong the phylogeny, but substitutions are counted separately along internal branchesand external branches. Instead of normalizing by the expected nonsynonymoussubstitutions and synonymous substitutions, the substitution counts are normalizedby the total internal and external branch lengths.30Average Over WindowsAs each site can be covered by multiple overlapping windows, site dN − dS iscalculated as the average site dN − dS across all windows, weighted by the sitedepth of unambiguous codons in each window.Site I is calculated as the average site I across all windows, weighted by thesite substitutions (along the phylogeny) in each window. Site IN and site IS arecalculated analogously to I except that windows are weighted by the window-site’snonsynonymous substitutions and synonymous substitutions.2.4 Method for Validating on Simulated DataIn order to test the accuracy of substitution counts across the phylogeny, Umber-jack’s estimates of dN−dS were compared to known values from simulated data.50 populations were simulated under parameter settings within ranges informed byempirical HIV measurements. Each population contained 1000 observed viruseswith genomes of 9kbp. The observed viruses in the simulated populations weresequenced in silico using parameter ranges measured from MiSeq Nextera (shot-gun) libraries produced by the BC Centre for Excellence in HIV/AIDS. Umberjackwas run on each read library resulting in 82 502 Umberjack window-site dN−dSestimates aggregated into 14 640 site dN−dS estimates.2.4.1 Selection of Simulation ParametersIn order to generate simulated data that covered the most representative combina-tion of parameters, we sampled the parameters using a Latin hypercube sampling[101]. Latin hypercube sampling creates an N-dimensional grid, where each di-mension represents a parameter. For each sample, a point is selected within thegrid to maximize the distance between each point. Parameter values for each sim-ulated dataset were taken from the coordinates of the selected points in the Latinhypercube. Parameter ranges can be found in Table A.1. The parameters for simu-lated datasets used to test Umberjack accuracy can be found in Table A.2.312.4.2 Simulated Topology Under SelectionThe topology of a population under selection was generated using an ancestralselection graph [45]. Under neutral evolution with no selection, the time from 2individuals to their common ancestor follows an exponential distribution wherethe time between branching events on the tree increases closer to the root [75]. Forlineages under selection, the ancestral selection graph corrects the distribution oftime between branching events using the selection rate [45]. Selection rates wereset to 0.01 per individual per generation in all simulated datasets. Average selectioncoefficients in envelope proteins of HIV populations in untreated HIV patients havebeen estimated to be 0.008-0.02 per individual per generation [77].2.4.3 Simulated RecombinationRecombination was simulated by altering the tree topology at each recombinationbreakpoint within the genome by randomly pruning and regrafting a subtree to an-other location within the tree. In addition to the 50 simulated populations generatedto predict Umberjack accuracy, we generated an extra 27 simulated populations tofocus solely on the effects of recombination by fixing the parameters for high levelsof diversity and high sequence quality. Although the effective HIV recombinationrate is 10−5 recombinations/bp/generation [77], we wanted to determine the effectsof light to heavy recombination using recombination rates in [10−4, 0 recombina-tions/bp/generation], equivalent to a breakpoint occurring up to every 2 codons inthe genome within the time course of the infection. Parameters for the recombina-tion datasets are found in Table A.3.2.4.4 Simulated GenomesThe branch lengths of the trees created by the ancestral selection graphs were inunits of generations. Keeping the topologies, the entire tree lengths were scaled tounits of substitutions/site according to the simulated population parameters: treelength = mutation rate x generations.To simulate the different mutation rates of genes, the genome was randomlysectioned in up to 3 sections. The mutation rate in each gene was selected fromthe range [2× 10−5, 8× 10−4] mutations/bp/generation. Mutation rates for HIV32sequenced in plasma have been reported at 10−5 - 10−4 mutations/bp/generation[14, 65].Using the trees as input, we simulated genomes for every observed and an-cestral individual in a population using INDELible [30]. Codon site dN/dS val-ues followed a discrete gamma distribution with 60 categories (shape=1.5, rate=3,scale= 13 ). Site dN/dS values in inter-patient HIV datasets [71] have been reportedto follow a smooth gamma distribution between [0, 3] with shape parameters in[0.2, 0.6] and rate parameters in [0.27, 1.1] .2.4.5 Simulated Shotgun ReadsMiSeq shotgun paired-end 2x250bp reads were simulated for the observed indi-viduals using a customized version of ART read simulation software that featuredadapter contamination [42]. When read lengths are longer than the DNA fragmentbeing sequenced, the sequencer will also sequence the primers attached to the DNAfragment, leading to adapter contamination in the read.Empirical MiSeq Nextera (shotgun) sequence quality distributions measuredfrom libraries generated by the BC Centre for Excellence in HIV/AIDS laboratorywere used as quality guides. These quality distributions resulted in an average se-quence error of 1% and <0.01% indel error across all simulated reads. Simulatedfragment lengths followed a Gaussian distribution with average fragment lengthbetween 104 - 500 bp and standard deviation fixed at 100bp. The libraries gener-ated at the BC Centre for Excellence in HIV/AIDS had mean fragment sizes in 104-147 bp and standard deviations in 64 -100 bp.The genome of each observed individual was sequenced between 0x-8x depthcoverage. Shotgun libraries in our real datasets were estimated as having 2x -25xdepth coverage per virus.Each library was aligned against the population reference genome sequenceusing BWA-mem [52].2.4.6 Umberjack settingsDepending on the Latin hypercube sampling, Umberjack was run on the simulatedpopulations using one of the possible settings:33Setting Config a Config bWindow Size 150bp 300bpMin Window Width Threshold for Read Inclusion in Window 0.7 0.875Min Window Depth Threshold for Window Inclusion 10 10Min Phred Score Threshold for N-Masking 15 202.4.7 Measuring Umberjack dN−dS AccuracyThe accuracy in Umberjack’s estimation of dN− dS at the window-site level wasquantified by the squared error:∆= (xˆ− x)2 (2.3)where xˆ and x were the window-site predicted values of dN−dS and known valuesof site dN−dS, respectively.2.4.8 Measuring Tree Topology AccuracyEvery tree inferred from a window in the Umberjack pipeline was compared againstthe true tree(s) of the simulated population. If there were no recombination break-points within the window, there was only one true tree to compare against theinferred window tree. If there were recombination breakpoints within the window,the true tree for each recombinant section within the window was compared againstthe inferred window tree. Even without recombination, topology and branch lengthdifferences between the true and inferred tree could occur due to sampling, se-quencing error, or heuristics in the phylogenetic reconstruction algorithm.The similarity between an inferred window tree and a true tree was measuredusing their weighted Robinson Foulds distance WRF [92]. Lower WRF values in-dicate higher similarity, with WRF = 0 indicating equivalent trees. WRF measuresthe distance between tree A and tree B by summing up the lengths of branches in Anot found in B, the lengths of branches in B not found in A, and the absolute differ-ence in lengths of branches found in both trees. Branch tA in tree A is equivalent tobranch tB in tree B if the 2 disjoint sets of tip nodes located on either side of branchtA are the same as the 2 disjoint sets of tip nodes located on either side of branch34tB.The overall accuracy of an inferred window tree was calculated as WRF , theaverage WRF between the inferred window tree and each true tree, weighted bythe width of the window represented by each true tree.To fairly compare a true tree against an inferred window tree, the true tree wasdownsampled to include only the individuals captured within the window reads. Acopy of each individual was made in the true tree every time that individual wassequenced in the window. The branch lengths in the true tree were scaled accordingto the diversity found in the slice of true sequences corresponding to the true treeand window. Using FastTree2, the branch lengths of the true tree were chosen tomaximize the likelihood of the slice of true sequences given the constrained truetree topology.2.5 Determining Features that Affect UmberjackAccuracyIn order to determine which features (also known as covariates) affected Umber-jack accuracy and rank their effects, we employed feature selection techniques:ranking univariate correlations and backwards feature selection in a regression ran-dom forest. Initially, backwards feature selection in general linear regression mod-els was also employed but abandoned due to poor model fits. Further details aregiven in Section 2.6.1.We simulated 50 populations with parameters described in Table A.2. Usingthe known dN− dS values from the simulated data, we measured Umberjack ac-curacy using the ∆ metric (Equation 2.3) and ranked features on their ability topredict ∆.2.5.1 Features Considered For Predicting Umberjack Accuracy(A) Window Total Breakpoint Ratio: We define the breakpoint ratio as theratio of bases on either side of the breakpoint within a window, with thelower number as the numerator. The total breakpoint ratio is the sum ofbreakpoint ratios from all breakpoints within a window. This metric is asummary statistic indicating the number of breakpoints within a window and35how close they are to the middle of the window.(B) Window-Site Unambiguous Codon Rate: Rate of unambiguous codonsper read at a codon site in a window.(C) Window-Site Amino Acid Depth: Depth of unambiguous amino acids at acodon site in a window.(D) Window-Site Substitutions: Total nucleotide substitutions at a codon sitein a window. Substitutions are counted between ancestral sequences alongthe phylogeny inferred by the Umberjack pipeline.(E) Window-Site Nonsynonymous Substitutions: Total nonsynonymous sub-stitutions at a codon site in a window. Substitutions are counted between an-cestral sequences along the phylogeny inferred by the Umberjack pipeline.This metric is calculating using Equation 2.1 from SLAC Section 2.3.5.(F) Window-Site Synonymous Substitutions: Total synonymous substitutionsat a codon site in a window, counted similarly to ‘Window-Site Nonsynony-mous Substitutions’.(G) Normalized Window-Site Expected Nonsynonymous Substitutions: Ex-pected nonsynonymous substitutions at a codon site in a window, normalizedby the total tree length, where the tree length is in units of nucleotide substi-tutions per codon site. Expected substitutions are counted between ancestralsequences along the phylogeny inferred by the Umberjack pipeline. Thismetric is equivalent to ET [N]Tree Length (see Equation 2.2). The values rangein [0, 3] and can be considered similar to the average number of nucleotidepositions within a codon that could cause a mutation to a different aminoacid.(H) Normalized Window-Site Expected Synonymous Substitutions: This met-ric is similar to ‘Window-Site Expected Nonsynonymous Substitutions PerBranch’, except it counts synonymous substitutions instead of nonsynony-mous substitutions.36(I) Normalized Window Tree Length: Length of the window tree inferred bythe Umberjack pipeline, normalized by the total reads in the window. Unitsare in nucleotide substitutions/nucleotide site.(J) Normalized WRF : Average distance between window tree inferred by Um-berjack pipeline and true trees corresponding to the window, normalizedby total reads within the window. Units are in nucleotide substitutions/nu-cleotide site. See Section 2.4.8.(K) Normalized Polytomies: Total polytomies in the window tree inferred bythe Umberjack pipeline, normalized by total reads in the window.(L) Codon Distribution P-value: Probability that the codon frequency distri-bution at a codon site within a window is the same as the codon frequencydistribution at the same codon site within the genomes of all extent individu-als in the simulated population. The probability is calculated using the G-testof independence between a set of distributions [67].(M) True Site Codon Entropy: Shannon entropy of codons at a codon site ofthe genomes of all extent individuals in the simulated population.(N) Window-Site Codon Entropy: Shannon entropy of unambiguous codons ata codon site in a window.(O) Window-Site Sequence Error Rate: Rate of erroneously sequenced basesper read at a codon site in a window.2.5.2 Univariate Correlation RankingThe effect of each feature on ∆ was ranked by their Spearman’s correlation [100], anon-parametric measure of the monotonicity of the relationship between two vari-ables. Since each feature was considered independently of other features, correla-tion only provided a rough picture of their effect on ∆.2.5.3 Regression Random ForestSince features had poor correlations and poor linear model fits to ∆, we used aregression random forest model to obtain more accurate predictions of ∆. A back-37wards feature selection algorithm produced a regression random forest using thebest set of features to predict ∆. The features can be found in Table 2.2. A re-gression random forest is a collection of regression trees. A regression tree decideswhich linear regression model should be used to predict the response of a particularinput datapoint. Each tip of a regression tree represents a separate linear regressionmodel. When predicting the response for an input datapoint, the regression tree istraversed from root to tip. Each node is assigned a feature, and the connected edgesrepresent conditions based on a feature threshold value. The path to the next nodeis determined by whether the input datapoint satisfies an edge’s condition [8].To train (i.e. generate) a regression tree, training data with known response val-ues are recursively partitioned to form each node. The training data are partitionedbased on a feature threshold value, and a linear regression model is fit to the datain each partition. The feature and its threshold value are selected to minimize theresiduals of the linear regression models [8].Random forests avoid overfitting by selecting a feature at each node from arandom subset of features, and training each regression tree with a random subsetof training data. Further, the random forest uses the average prediction amongst itscollection of regression trees as the final prediction [8].To gauge the importance of a feature, the values of that feature were permutedacross observations in each training dataset. For each tree, the mean squared er-ror between the true ∆ and the predicted ∆ˆ, MSE(∆, ∆ˆ), was calculated for eachregression tree using permuted and original training data as input for predictions.The feature importance score was calculated as the average increase in MSE(∆, ∆ˆ)due to permuting the feature, normalized by the standard deviation of MSE(∆, ∆ˆ).We used the Breiman (1999) test for feature importance significance, which treatsfeature importance scores as z-scores. The predictive performance of a randomforest was measured as the total MSE(∆, ∆ˆ) from each regression tree.For all random forest models, we used the R randomForest package [56] withdefault parameters, with the exception of ntree, which we chose as 501 to breakties:• mtry = the total random features to select from at each node = floor(one-thirdof total features)38• ntree = total regression trees = 501• permutations = 12.5.4 Backwards Feature Selection in Regression Random ForestThe backwards feature selection algorithm employed for the regression randomforests used cross-validation to decide upon the ideal model size, M, then fit a finalrandom forest model using the highest M ranking features.The ideal model size was chosen as the model size that yielded the best pre-dictive performance across 5 cross-validations. The training datasets were resam-pled between cross-validations. In each cross-validation, an initial random for-est was trained using all available features. A new random forest was trained foreach model size X, using only the X most important features. The feature impor-tance ranks were calculated from the initial random forest and fixed for subsequentmodels. The predictive performance of the model size was calculated as the to-tal MSE(∆, ∆ˆ) using the test datasets from the corresponding random forest in allcross-validations.The final set of features were selected according to their revised reranking fromall the cross-validations. A final random forest was trained using all the data.Model selection was performed using R caret package [31] and R randomFor-est package [56].2.6 Results2.6.1 Weak Linear Relationship Between Features and ∆The features considered for predicting the Umberjack accuracy (Section 2.5.1)demonstrated poor linear relationship with the ∆ metric, causing issues when fit-ting linear models. The Spearman’s correlation between each feature and ∆ wasfairly low, with absolute values ≤ 0.3 (Table 2.1).Certain default R implementations for generalized linear models (speedglmpackage, Gamma distributed response GLM) were unable to fit models with col-inear features, requiring manual exclusion of colinear features during backwards39feature selection. Some nested models were unable to be fit using default R im-plementations due to inability to find convergence of fitted values. These issuesmay have been avoided by using specialized linear regression models that penalizecollinear or low importance features such as Lasso or Ridge regression [41, 103].However, we decided to employ random forests to handle the possible non-linearrelations between each feature with ∆.Visual inspection of scatterplots of ∆ versus each feature in A.1 to A.8 indicatefairly noisy regressions with a wide spread of error.Feature Corr Lower UpperDiversityTrue Site Codon Entropy 0.16 0.15 0.16Window-Site Codon Entropy 0.15 0.15 0.16Normalized Window Tree Length -0.063 -0.07 -0.056Window-Site Nonsynonymous Substitutions 0.059 0.052 0.065Normalized Window-Site Expected Synonymous Substitutions -0.057 -0.064 -0.05Window-Site Synonymous Substitutions 0.035 0.028 0.042Normalized Polytomies 0.021 0.014 0.028Normalized Window-Site Expected Nonsynonymous Substitutions -0.013 -0.019 -0.0059SamplingCodon Distribution P-value -0.34 -0.34 -0.33Window-Site Amino Acid Depth -0.038 -0.045 -0.031Window-Site Unambiguous Codon Rate -0.0096 -0.016 -0.0028Reconstruction & Sequence ErrorWindow-Site Sequence Error Rate 0.16 0.16 0.17Normalized WRF -0.03 -0.037 -0.023Window Total Breakpoint Ratio 0.013 0.0067 0.02Table 2.1: Spearman Correlation Between Feature and ∆; lower and upper95% confidence interval.2.6.2 Advantages of Regression Random ForestsWe focused on random forests due to their built-in mitigation of feature collinear-ity. Since each node split is chosen from a random subset of features, random40forests are less likely to train models comprised mostly of collinear features. Re-gression random forests are also inherently able to model the non-linear relation-ships between the features through their piecewise linear modeling on partitions oftraining data. In addition, the built-in ensemble of regression trees reduces overfit-ting. Given that there were 82502 datapoints in the simulated datasets created totest Umberjack accuracy, the ease of parallelization of each regression tree in therandom forest helped reduce the runtime for model fitting.2.6.3 Features Affecting Umberjack AccuracyDue to the lack of strong linear relations between features and ∆, the feature rank-ing results from the univariate correlations approach and linear regression approachwere discarded from further analysis. The final regression random forest producedafter backwards feature selection was able to explain 77% of the variation in ∆(Figure 2.4).The model contained 11 features as shown in Table 2.2. Most notably, theranking of the feature importance tells us that Umberjack is most sensitive to lowgenetic diversity since it employs a substitution counting technique. Although anapproximate p-value was provided for every feature importance, they should onlybe considered a rough guide. The p-values are based on an approximation of theincrease in model MSE due to feature permutation across regression trees with anormal distribution.2.6.4 Effect of RecombinationA tree reconstructed from a multiple sequence alignment containing a recombina-tion breakpoint will be a combination of the true trees on either side of the break-point. Figure 2.2 and Figure A.14 demonstrate how an increase of recombinationled to an increase in error of inferred trees (as measured by WRF) and consequentlyan increase in error of Umberjack dN−dS estimates (as measured by ∆). However,at the recombination rates typical of HIV (1.4±0.6×10−5 recombinations/bp/gen-eration) [77] or lower, Umberjack estimates with high error were dominated by lowdiversity (as measured by tree length).41Feature Feature Importance Approx p-valueWindow-Site Substitutions 7.7 1.59e-14Window-Site Synonymous Substitu-tions6.3 4.08e-10Normalized Window Tree Length 3.8 1.31e-04Window-Site Sequence Error Rate 3.7 2.24e-04Codon Distribution P-value 3.7 2.56e-04Window-Site Codon Entropy 3.4 6.29e-04Normalized WRF 3.4 7.17e-04Window-Site Unambiguous CodonRate3.3 1.02e-03Window-Site Nonsynonymous Substi-tutions3.3 1.08e-03Window-Site Amino Acid Depth 2.8 4.90e-03True Site Codon Entropy 2.8 4.96e-03Table 2.2: Feature Ranking for Final Random Forest Predicting UmberjackAccuracy.The above results were taken from the 27 simulated population dataset con-centrating on recombination (Table A.3). However, similar results occured in the50 simulated population dataset created with varied parameters chosen from Latinhypercube sampling (Table A.2). Although the rate of recombination simulated inthe latter datasets was smaller, low site diversity (as measured by site substitutions)was the most importantly ranked feature affecting Umberjack accuracy (Table 2.2).Since Umberjack breaks up the genome into shorter windows, Umberjack canreduce the effect of recombination on estimates of site evolutionary statistics ifthe average distance between recombination breakpoints is larger than the windowsize.2.6.5 Umberjack AccuracyMany non-ideal scenarios were simulated to test the limits of Umberjack’s accu-racy. The Lin’s concordance coefficient between Umberjack’s estimate of dN−dSat each window site versus the true site dN−dS for all simulations was fairly lowat 0.59 (95%CI = [0.58-0.60]). Lin’s concordance coefficient ranges in [-1, +1],where 1 indicates perfect concordance, -1 indicates perfect discordance, and 0 in-42dicates no concordance.Umberjack concordance was > 0.83 for ideal scenarios representative of goodlibrary preparation and sufficient population diversity (Figure 2.3). These idealscenarios had an average 15% adapter contamination, 7 substitutions per site, 0.7xdepth coverage per individual, <1% sequencing error, and branch length of 5.0×10−3 substitutions per site.43(a)(b)Figure 2.2: Effect of recombination on ∆. Each point represents a window ofUmberjack estimates on simulated datasets focused on recombination(Table A.3). Black line indicates linear regression fit. Y-axis limitedto Window Average ∆ ≤ 1.0. (a) Points are shaded according the thelength of the window tree in units of nucleotide substitutions/site. (b)Points are shaded according to the recombination rate ρ of population inunits of recombinations/site/generation. Interval denotes average codondistance between breakpoints across the genome.44Figure 2.3: Umberjack inferred dN − dS vs expected dN − dS from simu-lated data. Each point represents a window-site extracted from a sim-ulated dataset. Black line indicates y=x. Red line indicates correlationof Umberjack inferred dN−dS to expected dN−dS for datasets repre-senting non-ideal scenarios such as poor library preparation, insufficientsampling, or low population diversity. Blue line indicates correlation ofUmberjack inferred dN−dS to expected dN−dS for datasets represent-ing ideal scenarios.45Figure 2.4: True Umberjack inaccuracy versus random forest prediction ofUmberjack inaccuracy. Umberjack inaccuracy was measured as thesquared error (∆, Equation 2.3) between Umberjack’s estimate of dN−dS and expected dN − dS. Each point represents a single site from asimulated dataset. Red line indicates y=x. Blue line indicates fitted lineof ∆ vs random forest predicted ∆. Inset: expansion of plot in the regionof ∆= (0,3) and random forest predicted ∆= (0,3).46Chapter 3Within-host HIV Evolution3.1 BackgroundThe host-specific adaptive immune response is a primary source of selection thatshapes the genetic variation of HIV [86]. For example, a vigorous cytotoxic Tlymphocyte (CTL) response (Section 1.1.2) within the first few weeks of infectionis associated with a substantial decline in viral load [6]. During the acute stage ofHIV infection, which lasts approximately for the first 6 months of infection, viralloads rapidly increase to a peak of around 107 copies/mL [36]. As a host CTLcell population targets certain HIV epitopes, new HIV escape mutations emergeand sweep. The CTL population size decreases when the number of infected cellsit binds decreases. Another CTL population targeting different sets of epitopestakes its place, resulting in another set of HIV escape mutations that sweep [33].Transition to a chronic stage of infection is associated with a diminishing rate ofCTL escape mutations, whereupon the viral loads decline to a ‘set point’ levelaveraging around 3×105 copies/mL, with substantial variation among individuals[33, 36]. To evaluate the role of the CTL response in shaping the evolution of HIVwithin hosts, we processed and analyzed shotgun sequencing of viral RNA isolatedfrom plasma of 31 treatment-naive HIV patients with high viral loads.473.2 Methods3.2.1 DatasetThe study used frozen plasma samples from BC HOMER cohort individuals whohad provided informed consent to donate and store the remainder of their routinelycollected plasma samples for the purposes of research relevant to the pathogenesisand treatment of HIV under the Experimental HIV Monitoring Program. Sam-ples were chosen from individuals on the basis of the following criteria: first, theindividuals had at least two samples available before the first therapy start date;second, the samples were collected at least 90 days apart (median 377 days) toincrease the chances of measurable evolution occurring between observations; andthird, the samples were associated with a minimum viral load of 5×104 copies/mL.The last criterion was employed to reduce the probability of template resampling,in which a single template is represented by multiple sequences. Not only doestemplate resampling skew minor variant frequency counts, it also introduces falsediversity during phylogenetic reconstruction if there is sequencing error amongstthe resampled templates [60]. We note that templating resampling is difficult todiscern bioinformatically from true biological proliferation of a variant. There arewet lab methods such as primer IDs that can tag templates to determine how manyunique templates there are in a population; however, these methods were not em-ployed for this dataset [108]. Although patient CTL responses were not explicitlymeasured through ‘wet lab’ assays, the high viral loads implied adaptation to thehost-specific immune responses by the respective HIV populations. Most patientseither demonstrated an increase in viral load or a change within assay measure-ment error. The exceptions were patients BC11 and BC22, for whom a log10-folddecrease in viral load of <-0.67 and <-0.68 was observed, respectively. Furtherdetails on patient measurements can be found in Table A.4.In addition, the genotypes for HLA genes HLA-A, HLA-B, and HLA-C forthese individuals had been previously determined using the procedure documentedin [11]. Briefly, the HLA genes were sequenced with Sanger sequencing, andthe HLA genotypes classified at the resolution of allele group (coarser) or protein(finer). Allele groups were further imputed down to the protein level using the48software HLA Completion [58] for individuals for whom allele groups from allgenes were available. Five patients were missing sequences for HLA-A and HLA-C and thus only their HLA-B allele groups were known. HLA Completion wasunable to compute suitable protein resolutions for one patient due to maximumlikelihood convergence issues; consequently, the HLA genotypes for this patientcould only be resolved at the allele group level.HIV RNA was extracted from the plasma samples using a NucliSENS easy-MAG instrument (bioMe´rieux, St. Laurent, Quebec, Canada). Nested RT-PCRamplification reactions in which a shorter region was isolated from a wider regionin two rounds of PCR were performed as described in previous work [1] usingprimers specific to regions of the HIV-1 genome encoding gag p17, nef, and thethird variable loop of env gp120. The amplified samples were used to generateIllumina Nextera XT libraries and sequenced on an Illumina MiSeq instrument us-ing a paired-end 2 x 250bp reagent kit (version 2), resulting in a random ‘shotgun’distribution of short reads across each target region. Median depth of coveragereached 8200-, 10000-, and 65000-fold for genes gag, nef, and env, respectively.Further details on the sequencing libraries can be found in Table A.5.A phylogenetic tree was constructed using the sample consensus sequences foreach patient at each time point in the original dataset. All patients were infectedwith HIV1B subtype of HIV. Samples that did not phylogenetically cluster togetherby patient indicated potential cross-contamination and were discarded, forming thefinal 31 patient dataset ( Figure A.9, Figure A.10, Figure A.11).3.2.2 Umberjack ProcessingReads from each patient sample were iteratively aligned using BWA-mem [52]against the HXB2 reference (GenBank accession K03455.1), where the previousiteration’s consensus was used as the reference in the next iteration. Re-alignmentwas halted after 10 iterations or if the consensus sequence did not change fromthe previous iteration. When forming the consensus sequence, gaps that inducedframeshifts (gaps not in groups of three) were shifted together into groups of threeas long as they did not introduce premature stop codons. The SAM files from everypatient sample alignment were used as inputs for Umberjack to estimate dN− dS49and I per codon site, Section 1.3.2, (window size = 300bp, window overlap = 30bp,read bases with quality score <20 masked with N, reads with >12.5% missingbases within a window excluded, windows <50 sequences excluded).We note that PCR amplification of populations has been known to introduce ar-tificial recombination that is indistinguishable from true biological recombination[99]. Due to the short average fragment size of the dataset (<200bp), we did notemploy bioinformatic tools to detect the rate of recombination, as there would beinsufficient sequence length to obtain confident estimates. We relied on excludingsites in which random forest predictors deemed that Umberjack produced inaccu-rate estimates.3.2.3 Umberjack CleaningRandom forest predictors were tested for their ability to predict Umberjack errorsolely using features that could be extracted from real datasets. Features whichwould most likely be unknown in an experimental context, such as ‘Window-SiteSequence Error Rate’ and ‘Normalized WRF’, were excluded from backwardsfeature selection. The final random forest trained on the best performing real-dataset features (‘Normalized Window-Site Expected Synonymous Substitutions’,‘Window-Site Synonymous Substitutions’, ‘Normalized Window-Site ExpectedNonsynonymous Substitutions’, ‘Window Unambiguous Codon Rate’, ‘Normal-ized Window Tree Length’) explained 19% of the variance in error within Umber-jack estimates of window-site dN−dS.Although the random forest did not predict error well for dN−dS calculationsat the window-site level, it was still useful in removing unreliable window-sitedN − dS estimates to improve the accuracy of site dN − dS estimates averagedacross windows. Using an empirically derived threshold, Umberjack estimates ofwindow-site dN−dS with an excessive predicted error by the random forest modelwere excluded from further processing. Umberjack site dN−dS was calculated byaveraging estimates across the cleaned windows. Removing dubious window-siteestimations reduced the total number of site dN−dS estimates from the simulateddataset by 1%. Based on a simulation analysis in the previous chapter, the R2for cleaned estimates of site dN − dS compared to true site dN − dS was 80%50Dataset Patients Samples PerPatientTime Span(post infec-tion)Total SeqShankarappaet al. [97]11 10-14 6-12 years 1300 envLiao et al. [55] 1 10 1-4 months 295 envLiu et al. [61] 2 2, 23 8 days-4 years >155 gag,>194 nef,>165 envTable 3.1: Composition of published HIV-1 sequence data sets used to derivean empirical distribution of site-specific dN−dS estimates.(Figure A.12).The error threshold (∆ < 10.7) was based on the IQR (Interquartile Range) ofempirical dN−dS values taken from other empirical datasets of longitudinal sam-ples of untreated patients where several HIV viruses were isolated per timepointand sequenced using conventional capillary-based methods (Table 3.1).Further, the indel-rich V4 region of the env gene (amino acid coordinates 395-413 with respect to ENV protein in HXB2 strain) was excluded from the analysisto avoid erroneous nonsynonymous substitutions that could arise from the spuriousalignment of non-homologous insertions.3.2.4 Finding CTL Response Associated Amino Acid PolymorphismsHIV amino acid polymorphisms associated with CTL response by Carlson et al. ina prior study [11] were located within the multiple sequence alignments and recon-structed ancestral sequences were generated by Umberjack from the patient sam-ple libraries. The CTL-associated HIV polymorphisms were classified as ‘MatchedEscape’ if the polymorphism associated with evading a host-specific CTL responseand a selective HLA allele was present in that patient, ‘Matched Nonescape’ if itwas associated with susceptibility to the CTL response and a selective HLA allelewas present, or ‘Unmatched’ if no selective HLA alleles were present.513.2.5 dN−dS AnalysisSite dN−dS was estimated for each sample separately and across samples (pooleddata) for each patient. Sites under significant selection according to dN−dS weredetermined by performing binomial tests on substitutions counts using ‘nonsyn-onymous’ and ‘synonymous’ as the response categories and the expected synony-mous substitutions as the probability of success [80].3.2.6 I AnalysisSite I and IN− IS were calculated for each sample separately and across both time-point samples for each patient. Sites under significant selection according to I weredetermined by performing permutation tests in which substitutions were randomlyreassigned as either ‘internal’ or ‘external’ using the total internal and externalbranch lengths as weights. The p-value was calculated as the fraction of permuta-tion trials in which the Iˆ calculated from permuted substitutions was further from0.5 (the neutral threshold value) than the original I calculated from non-permuteddata. Site substitutions were permuted for 104 independent trials in each windowcovering the site. The total site substitutions in each window were kept constantduring each trial.In order to perform multivariate regression on IN − IS statistics, IN − IS valueswere transformed from [−1,+1] to [0,1] so that it could be modeled as a zero-oneinflated beta distribution. Further, all p-values from each separate zero-one inflatedbeta regression were Benjamini-Hochberg corrected.3.2.7 Timing InfectionSince clinical information for the dates of HIV-1 infection was unavailable, weestimated these dates using the software package BEAST (Bayesian EvolutionaryAnalysis by Sampling Trees) [23]. BEAST employs a Bayesian Markov chainMonte Carlo (MCMC) method to generate a random sample from the joint pos-terior distribution of rooted trees and rates of evolution in units of time, givensequences labeled with dates of sample collection. Using the most recent commonancestor as a proxy for the infecting virus, the date of infection was set to the dateat the root of the tree. Since the Bayesian MCMC procedure was computationally52intensive, the date of infection was only inferred for one window of alignment perpatient. Both baseline and followup sequences were pooled into the same align-ment and windows were extracted by Umberjack (window size = 300bp, windowoverlap = 30bp, read bases with Phred quality score < 20 masked with N, readswith > 12.5% missing bases within a window excluded, windows < 50 sequencesexcluded). Windows to analyze with BEAST were selected on a score based on acombination of the highest average site entropy across the window, highest aver-age site depth coverage across the window, and lowest percentage of missing bases(N’s, gaps). 85% of the selected windows encompassed the V3 variable region inenv. After running BEAST for two chains of 3× 108 MCMC samples for eachpatient, mutual convergence of both chains on the estimated date of infection wasevaluated through visual inspection of trace plots displaying estimates. The first107 MCMC samples from each chain were discarded as a ‘burn-in’ period.3.3 Results3.3.1 Genetic Variation Across SitesAs expected, sequences covering HIV env gp120, which encodes the surface en-velope glycoprotein, tended to be more variable. Using amino acid entropy toquantify variation, there was a significant difference in variability between genes,with entropy in both nef and env greater than gag (t = 11.9 and t = 40.6, respec-tively; likelihood ratio test, P < 10−15). Entropy in the follow-up samples alsovaried over the genome in a consistent pattern amongst subjects. For example, theamino acid entropy along the portion of HIV env covered by the NGS data tendedto be greatest near the disulfide loop regions (V3, V4 and V5), and lowest in re-gions associated with conserved functional motifs such as the GPGR motif withinV3 and the CD4 binding loop (Figure 3.1). Note that we excluded the V4 regionfrom any further analysis regarding selection due to issues with multiple sequencealignment of excessive indels.53V3GPGR V4 V50.00.51.01.52.0200 300 400 500HXB2 reference codon coordinateAmino Acid EntropyFigure 3.1: Summary of amino acid entropy along HIV env. Each trace cor-responds to a different individual’s follow-up sample.3.3.2 Most Subjects in Chronic Stage of InfectionBased on BEAST timing estimates of the most recent common ancestor, 23/31subjects were in the chronic stage of HIV infection by the baseline sampling date,with a median infection date of 2.6 years prior to baseline (Figure 3.2, Table A.4).The timing of infection for 7/31 subjects was <9 months, which we classified asclose to the boundary of acute/chronic stages of infection. The viral loads and CD4counts from these subjects were also ambiguously acute/chronic. Their minimumviral load was 5×106. The limit of quantification on the viral load assay was 106copies/mL, which several of the patients hit; however, this study population hadbeen selected on the basis of having samples with high viral loads. This can biasresults in that patients with high viral loads indicate proliferation of viral popula-tions and adaptation of the viral population to the host immune system. CD4 countsassociated with these samples were within 330-541 cells/mm3. The BEAST runsfor one of the patients did not converge on a timing estimate, and their infectiontiming information was excluded from further analysis. The median estimated rateof evolution across subjects was 8.3×10−5 mutations per base per day, which was54consistent with published estimates of the HIV mutation rate (3.5×10−5 [64]) andthe hypervariability of the V3 region.0 2 4 6 8 10 12051015202530Years before baselineSubjectFigure 3.2: Estimated durations of infection at baseline by molecular clockanalysis. Each point represents the median estimate for a given subject.The thick line segment indicates the interquartile range, and the thinline indicates the 95% credible interval. A dashed line is drawn at sixmonths to indicate which subjects may have been at an acute or earlystage of infection at baseline.3.3.3 Evolving Sites Associated with CTL ResponseTo study evolution at the protein level in these data, we extracted the amino acidfrequency distributions per site from windows of the NGS alignments for each geneand subject. We used a G-test statistic to evaluate changes in these frequency distri-butions per site between the baseline and follow-up samples. To adjust for multiplecomparisons, we used a simple Bonferroni correction such that α = 1.9×10−6 toclassify sites with a significant change in the amino acid frequency distributionover time, which will be referred to as ‘evolving sites’ for brevity. Despite thisconservative procedure, 5354 (20%) out of 26186 tests were considered signifi-cant. Evolving sites were more likely than not to contain an HIV polymorphismmatched with at least one of the subject’s HLA alleles, insinuating that CTL medi-55ated response drove selection at those sites in those individuals (Odds Ratio=1.3,P = 9.7×10−7).Matched escape mutations were significantly more likely to increase in fre-quency than matched non-escape mutations, and this trend was the most pro-nounced in HIV nef (P = 0.015). Additionally, the increase of matched escapemutations over time was significantly greater for mutations restricted by HLA-Balleles than alleles at the other HLA loci (P = 0.023; Figure 3.3). This is con-sistent with previous studies that determined that HLA-B alleles were the mosteffective at viral load suppression [36].To investigate these associations in greater detail, we fit a generalized linearmixed model on the classification of evolving sites using a binomial link function,with subject as a random effect on model intercepts, rejecting simpler models onthe basis of the Akaike information criterion. In addition to terms correspondingto known HLA associations and matched status within subjects, we observed sig-nificant fixed effects of genes relative to HIV gag (nef : z = 2.2, P = 0.02; env:z = 26.4, P < 10−12); in other words, evolving sites were observed significantlymore often in nef and env. Based on the model and residual deviances, however,the mixed model explained only about 8% of variation in the binomial outcomeof being classified as an evolving site. Considering env is also known as a highlyvariable gene in general, and adaptive to the human immune response, it is possiblethat the increased significance of env compared to nef is due to the larger numberof substitutions found in env. Further, it is possible that the majority of the sitesfound under significant selection were hitchhiking variants caught in a selectivesweep; however, it is difficult to disentangle true selection from a mere sweep.3.3.4 Host HLA Genotype Drives Diversifying SelectionThe next objective was to investigate whether the application of phylogenetic meth-ods could improve or supplement the study of HIV evolution within hosts beyondwhat was attained using frequency distributions alone (section 3.3.3).Site specific patterns of diversifying selection were analyzed using the dN−dSstatistic (Figure 3.4) across samples for each patient using Umberjack. On average,most sites tended to experience purifying selection (median dN − dS = −0.13).56Matched Nonescape Matched Escapelllllllllllllllllllllllllllllllllllll llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll llllllllllllllllllll lllllllllllllllllllllllllllllllllll ll llllllllllllllllllllllllllllllllllllllllllllllllllllllll llllllllllllllllllllllllllllllllll lllllllllllllllllll0.000.250.500.751.000.000.250.500.751.000.000.250.500.751.00HLA−AHLA−BHLA−CBaseline Followup Baseline FollowupSampleFrequencyFigure 3.3: Shifts in matched non-escape and escape mutation frequenciesover time in HIV nef across all subjects. Black lines indicate the changein the frequency of the amino acid polymorphisms from baseline tofollow-up. There was a small but significant tendency for non-escapemutations to decline, and for escape mutations to increase, which wasthe most apparent for mutations restricted by HLA-B alleles. Red dotsindicate the overall mean frequency in each group.Roughly 36% of all sites across genes and subjects had values of dN − dS ex-ceeding zero (diversifying selection), of which only about 0.97% were consid-ered significant based on a q-value cutoff of 0.05 (Benjamini-Hochberg correc-tion). dN− dS varied amongst genes, with significantly greater values observedon average for HIV nef than either env or gag (t = 6.0, P = 1.9× 10−9). Therewas also a significant association between dN − dS and whether that site had aCTL-associated polymorphism matched to the subject’s HLA genotype (t = 5.8,P = 7.0× 10−9), which is consistent with the CTL mediated response driving thediversification of HIV within hosts. Finally, dN− dS was marginally associated57GAG NEF ENVlllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll−4−2024MatchedUnmatchedMatchedUnmatchedMatchedUnmatchedHLA MatcheddN − dSFigure 3.4: Boxplots of dN−dS at HLA Matched Sites and Unmatched Sitesby Gene. Red dots indicate means.with an indicator variable of whether the subject was classified with an acute orearly stage of HIV infection (defined here as being within the first year of infec-tion) on the basis of the molecular clock analysis (t = 2.2, P = 0.026; section3.3.2). The association between dN−dS and HLA matched sites still existed afterstratifying by gene (Mann-Whitney U , gag, P = 1.1×10−8; nef, P = 3.7×10−11;env, P = 4.3×10−3; Figure 3.4).This last result is consistent with the bulk of CTL mediated selection occurringin the earlier period of an HIV infection. At the same time, it may be an artifact ofthe transient nature of diversifying selection in the relatively constant environmentpresented within a single host. Once the host-specific immune response has beenstimulated by the infection, the virus population undergoes a number of ‘selectivesweeps’ (Figure 3.3). During the sweep, an elevated amino acid substitution rate ismanifested by an elevated dN−dS rate indicating selection. After the conclusionof the selective sweep, however, we do not expect to observe any further evolutionat the amino acid level at the sites targeted by the immune response and dN −dS will revert to lower levels consistent with purifying selection. If a sample istaken at the point of a sweep when synonymous substitutions begin to overtakenonsynonymous substitutions to maintain favoured amino acids, dN−dS may even58indicate neutral evolution. Consequently, dN− dS alone is not a sufficient metricfor quantifying the response to selection by HIV within hosts [46].The temporal nature of this statistic is exemplified when we examined the per-sample dN− dS values in Figure 3.5 as opposed to the across-samples dN− dSvalues in Figure 3.4. We compared sites with matched polymorphisms to sites withno matched polymporphisms (Figure 3.5). Here, dN− dS hovers around neutralor purifying at ‘Matched’ sites and dN−dS is purifying at ‘Matched’ for gag andnef, and a slight increase of Matched dN−dS upon followup in env. However, thedifferences between Matched and Unmatched groups is insignificant (P>0.55, lin-ear regression), though we expect that CTL response will imbue selective pressureon Matched sites.3.3.5 Evidence of Selective SweepsThe I statistic suffered from similar issues as dN−dS with respect to aggregatingover the substitutions. Similar to the original Fu & Li statistic [32], I does notdifferentiate between nonsynonymous and synonymous mutations. Consequently,directional selection at the amino acid level may appear neutral. After Benjamini-Hochberg multiple test correction, only 3% of sites were under significant diversi-fying selection and only 19% of sites were under significant purifying selection. Inthis section, we demonstrate that IN − IS is a better measure of relative directionalselection at the amino acid level. Estimates of IN − IS (section 1.3.2) were looselycorrelated with measures of diversifying selection dN− dS (Kendall’s ρ = 0.42,P < 2.2×10−16).Across all patients and sites, nonsynonymous substitutions were significantlymore likely than synonymous substitutions to be mapped to tips of the phylogeny(Figure 3.6; IN < IS, Mann-Whitney U = 2.37×108, P < 2.2×10−16). This resultis consistent with purifying selection being the dominant mode of selection withinhosts. In other words, most nonsynonymous mutations are deleterious. This trendvaried significantly among genes. For instance, from Figure 3.7 and Figure A.13,median IN− IS ≈ 0 in env, indicating that nonsynonymous and synonymous muta-tions were similarly distributed between internal and external branches. However,in gag, median IN − IS = −0.11 indicating that gag nonsynonymous mutations59(a) gag (b) nef(c) envFigure 3.5: Site dN− dS per sample at sites associated with CTL responsefor patient’s HLA. Site dN− dS calculated separately for baseline andfollowup samples. Violin plots show density of sample-sites at eachvalue of dN − dS. Inner boxplots show IQR and median. ‘Matched’groups refer to any site associated with CTL response in study cohort([11]) sharing patient’s HLA. ‘Unmatched’ groups refer to any site withno CTL association for study cohorts sharing patient’s HLA.60Figure 3.6: Substitutions tend to map to tips of the within-host phylogeny.Density plots comparing how often site nonsynonymous and synony-mous substitutions map to internal branches of phylogeny. Red linesindicate site IN , blue lines indicate site IS. IN and IS range in [0, 1],where 0 means all substitutions were found in the tips, and 1 means allsubstitutions were found in internal branches.tended to occur more recently than synonymous substitutions, implying that gagwas under stronger purifying selection than env (t = 19.6, P = 1.5×10−83, linearregression). gag epitopes have been previously reported to have a stronger im-munogenic in CTL response than env [37]. Alternatively, the higher indel rate inenv could have induced more alignment errors, driving up artificial mutations seenat the tips of the env phylogeny.To determine the effect of HLA pressure on HIV evolution, we stratified sitesby whether they contained a Matched polymorphism at anytime during the popula-tion history, as observed at sequencing or reconstructed by the Umberjack pipeline.An entire amino acid site was labeled ‘Matched’ as long as it contained at least oneMatched polymorphism, and ‘Unmatched’ if it never contained a Matched poly-morphism. IN − IS at Matched sites increased significantly with median ≈ 0 in all61Figure 3.7: Nonsynonymous substitutions map to tips of the phylogeny moreoften than synonymous substitutions. The violin plots summarize theoverall distributions of IN and IS statistics across subjects, broken downby gene. A lower I statistic indicates that the substitutions tend to maponto the tips of the phylogeny, which is consistent with purifying se-lection. Black point ranges indicate median and IQR. White squaresindicate mean.genes, and mean equal to −0.03, −0.02, and 0.08 for gag, nef, and env respec-tively (t = 5.4, P = 6.0× 10−8). These greater values of IN − IS indicate that thenonsynonymous substitutions at these HLA-restricted sites mapped earlier in thephylogenies, which is consistent with recent selective sweeps (Figure 3.8). Thedifference between site IN and IS was significant for env, the most variable genestudied in this thesis (env: t = 2.9, P = 0.054). Further, IN for Matched sites wereconsistently higher than IN for Unmatched sites across each gene (gag: t = 9.7,P = 2.6×10−21, nef : t = 6.4, P = 6.2×10−10, env: t = 2.6, P = 0.013), leadingto more evidence that CTL pressure was driving directional selection and causingfavoured amino acids to overtake the population earlier on.Since overall site IN − IS does not differentiate between amino acids, it is bestused as a relative measure of directional selection, as IN− IS = 0 does not necessar-ily mean that a site is neutrally evolving. When sites experience early mutations tofavoured amino acids at sites that may continue to diversify through synonymous62Figure 3.8: Distributions of nonsynonymous and synonymous substitutionsin the tree are more similar at sites with known CTL associations. Theviolin plots summarize the distributions of site IN − IS per patient forsites with known CTL associations. Sites are stratified by whether thesite contained a Matched amino acid at any point in the reconstructedwithin-host phylogeny. Each constituent site IN − IS represents the dif-ference between proportion of site nonsynonymous and synonymoussubstitutions mapped to internal branches of the phylogeny containingboth baseline and followup sequences from a single patient. As longas the site in the patient samples contained at least 1 Matched aminoacid in any of the observed sequences or inferred ancestral sequences,the site was labeled ”Matched”. A lower I statistic indicates that thesubstitutions tend to map onto the tips of the phylogeny, which is con-sistent with purifying selection. The filled circle and whiskers indicatethe medians and IQR respectively. White squares indicate means.63substitutions, recent mutations away from the favoured amino acids at the tips willdecrease site IN , driving IN− IS towards parity. If the site is truly under directionalselection, the unfavoured variants will be purged from the population. Therefore,in order to truly see directional selection by comparing internal vs external sub-stitutions, we need to explicitly label the amino acid polymorphisms at each site.Instead of calculating overall site IN− IS and then stratifying sites as Matched/Un-matched, we stratified substitutions at each site as Matched/Unmatched accordingto the specific amino acid resulting from the substitution. From Figure 3.9, it isevident that nonsynonymous substitutions to Matched amino acids tended to occurearlier than synonymous substitutions, implying that lineages carrying the Matchedamino acids proliferated in the ancestral population under positive selection. Con-versely, nonsynonymous substitutions to Unmatched amino acids occurred moreoften at the tips than the synonymous substitutions, implying that existing aminoacids conferring a fitness advantage were maintained and new deleterious aminoacids were purged before they could proliferate (gag: t = 9.6, P= 7.6×10−2, nef :t = 10.7, P = 1.4×10−26, env: t = 3.6, P = 2.6×10−3).3.3.6 Reversion to Wild TypeThere were no significant differences when the distributions of IN and IS shown inFigure 3.9 were further stratified by substitutions to new Matched Escape versusMatched NonEscape amino acids (P > 0.16). The Carlson et al. study, from whichwe extracted the HLA allele - HIV polymorphism associations for this analysis,reported that Matched NonEscape amino acids were typically (>80%) the same asthe wild-type amino acid, which was defined as the HIV-1 subtype B consensus se-quence as described by the Los Alamos HIV Sequence Database [47]. Reversion towild type epitopes has been noted in patients infected with strains containing mis-matched escape mutations [51, 107]. These reversions took years to emerge, whichwould be consistent with our dataset of chronic patients. Breaking down IN and ISby substitutions towards a wild-type or mutant amino acid at each site (Figure 3.10)revealed that reversions to wild-type occurred significantly earlier in the phyloge-nies (gag: t = 14.6, P= 1.5×10−83, nef : t = 20.4, P= 1.2×10−90, env: t = 15.1,P = 2.6× 10−51), implying a selective advantage relative to the mutant residues.64(a) gag (b) nef(c) envFigure 3.9: Nonsynonymous substitutions to HLA-matched amino acids tendto map deeper in within-host phylogenies. Violin plots summarize thedistributions of IN and IS, stratified by substitutions to Matched or Un-matched amino acids at each site. Each constituent site IN and IS respec-tively represents the proportion of site nonsynonymous and synonymoussubstitutions (to the specific residues) mapped to internal branches ofthe phylogeny containing both baseline and followup sequences froma single patient. A substitution was labeled ‘Matched’ if it resulted inan amino acid Matched with the patient HLA, and ‘Unmatched’ oth-erwise. Black point ranges indicate median and IQR. White squaresindicate mean. Width of violins represents total matched or unmatchedsites.65Although unmatched amino acids typically appeared only at the tips, mutationsto unmatched wild-type appeared to be strongly selected in the absence of CTLpressure very early on in the evolutionary history of the infection (Mann-WhitneyU , P < 1.4× 10−168). Wild-type amino acids always appeared earlier in the phy-logeny than mutant amino acids, but amongst Matched polymorphisms, only nefwild-types were significantly earlier than mutants (P= 6.1×10−39). There was nosignificant difference between the appearance of Matched or Unmatched wild-typeamino acids (Mann-Whitney U , P > 0.29).Only 0.038% of amino acid substitutions exhibited a Matched Escape aminoacid mutating towards an Unmatched wild-type amino acid, indicating that CTLselective pressure continued and was usually stronger than the wild-type fitnessadvantage. At sites in which the Matched Escape reverted to an Unmatched wild-type (Figure 3.11), nef imparted the strongest drive towards wild-type reversion.Amino acid substitutions to wild-types were significantly maintained by synony-mous substitutions (IN > IS, Mann-Whitney U , Benjamini-Hochberg corrected P=5.4× 10−5). In env, nonsynonymous substitutions to Unmatched wild type weregenerally earlier than synonymous substitutions, but there were enough nonsyn-onymous substitutions at the tips to deem Unmatched wild-type amino acids neu-trally selected with a median IN = 0.49 (Mann-Whitney U , Benjamini-Hochbergcorrected , IN < IS, P = 0.0034). gag had the highest baseline mean frequency ofwild-type amino acids (94%) which were maintained by synonymous substitutionsdeep in the tree, but there were insufficient nonsynonymous substitutions to deter-mine if there was significant active amino acid-level selection towards wild-type(Mann-Whitney U , Benjamini-Hochberg corrected , IN < IS, gag P = 0.89). Al-though the precise timing of each substitution could not be reconstructed due tocomputational complexity of Bayesian inference, the wild-type variant to whichthe Matched Escape variant mutated towards already existed either as a major orminor variant at baseline sampling for all genes.66(a) gag (b) nef(c) envFigure 3.10: Nonsynonymous substitutions to wild-type residues occurdeeper in the phylogenies. The violin plot distributions of IN and IS arestratified at each site by substitutions to wildtype amino acid (HIV1-B subtype consensus). Black point ranges indicate median and IQR.White squares indicate mean.67Figure 3.11: Violin plot distributions of IN and IS for site substitutions to-wards Unmatched wild type amino acid from Matched Escape aminoacid. HIV1-B subtype consensus used as wild type. Black point rangesindicate median and IQR. White squares indicate mean.3.4 Within-host HIV Evolution During Drug Treatment3.4.1 BackgroundIn addition to characterizing putative selection from immune response, we soughtto validate how well Umberjack could detect selection in a dataset exhibiting veryclear and strong directional selection. As such, we compared the treatment-naivedataset to a previously published Maraviroc clinical drug trial dataset [68].We employed Umberjack to profile within-host evolution of the 4 patients un-dergoing HIV drug treatment in this trial. Patients were sampled before and af-ter the addition of Maraviroc to an existing drug regimen. Over a median of 5368days, 4-5 longitudinal samples were taken per patient. Samples were amplicon se-quenced with single 454 reads. Although the amplicon sequencing resulted in farmore sequence overlap than would be present in shotgun sequencing, low-qualityclipping lead to variable length reads with jagged alignments, requiring the use ofUmberjack to infer evolutionary statistics at end positions of the alignment. En-hanced Sensitivity Trofile Assays (ESTA) for three out of four patients indicatedthat their viral populations used both CCR5 and CXCR4 co-receptors before start-ing Maraviroc treatment. The co-receptor usage for each virus was also inferredwith Geno2Pheno [3, 68], a machine learning algorithm that uses viral sequenceto build its predictions. Each patient displayed a marked increase in viruses usingCXCR4 between the baseline and final sampling.3.4.2 ResultsOverall, sites from the Maraviroc treated HIV populations were under slightly pu-rifying selection (median I = 0.47). Even though there were several longitudinalsamples per patient, I calculated from the substitutions accumulated from base-line at each sample was not a predictive measure of whether the following samplewould experience a sweep, since sweeps occurred at different rates amongst sitesand patients. That is, I >0 at a site in the previous sample would not mean that thenext sample would not incur a sweep. However, when examining I across all sam-ples, sites in which a variant swept throughout the population in the drug treateddataset were associated with higher I-statistics (P = 7.8× 10−16, Mann WhitneyU) (Figure 3.12).Further, the stronger directional selection due to drug pressure is evident inthe higher I of the drug treated dataset (median I=0.47) compared to the untreateddataset (median I = 0.22), indicating that drug resistance mutations proliferated inthe treated population more so than CTL escape mutations in the untreated pop-ulation. Drug treatment imparts a constant selection pressure, whereas the CTLimmune response varies depending on the CTL cell population size within the pa-tient, leading to variable and weaker drive towards CTL escape [5].69Figure 3.12: Violin and boxplot of I statistic for sites in ENV gene (V3 re-gion) of viral populations treated with Maraviroc.‘Sweep’ and ‘No sweep’ refer to sites in which a variant proliferated throughoutthe entire population by the final sample and sites in which a variant did notproliferate throughout the entire population by the final sample, respectively.3.5 Comparing Methods for Detecting SelectionFrequency based methods of detecting selection can confound selection with de-mographic events such as contraction and expansion of population size [91]. How-ever, site frequency is easily calculated and does not require long reads to obtainevolutionary relationships. Phylogenetic methods will take evolutionary relation-ships into account to avoid demographic confounding [2, 91]. However, whencounting methods are used to evaluate the distribution of mutations in the phy-logeny, there is low power to detect selection at sites with low diversity [80]. More-over, in the untreated patient dataset examined in this study, only 7% of reads werelong enough (>263bp) to sufficiently fill a window for phylogenetic analysis due70to a shorter median fragment size of 109bp (IQR 63-189bp). The median depth ofcoverage of sites that passed window-inclusion thresholds was 362-fold whereasthe median depth of coverage across all sites was 7306-fold. This resulted in dras-tically lower percentages of sites considered under significant selection in compar-ison to a frequency-based method. Thus, while it seems intuitively true that mostsites would not be under directional selection, it is possible that a reduced power ofthe method, combined with fewer observations per site, prevented the phylogeneticmethod from detecting sites under weak directional selection.A disadvantage of aggregating site substitutions across all internal and externalbranches in I, IN , and IS statistics is that we lose the connection of substitutionsalong lineages. If our assumption that selection is similar across all lineages ata site is inaccurate, then it is possible for lineages exposed to purifying selectionor directional selection in different parts of the phylogeny will cancel each otherout. However, the lower resolution of the I selection statistic is more robust tosampling and topological errors during window phylogenetic construction, as wellas limited amino acid diversity. For example, the IN − IS statistic computed forsubject BC27, window 271-570bp in HIV-1 gag, codon site 147 was−0.01, whichis basically indistinguishable from a neutrally-evolving site. Thus, we fail to de-tect strong directional selection towards the HLA-matched amino acid (leucine),to which substitutions were mapped onto multiple lineages in the within-host phy-logeny (Figure 3.13). Breaking down the window-site by substitutions to Matchedamino acids shows that the Matched mutations appear early during the popula-tion history and are maintained deep within the tree by synonymous substitutions(Matched IN = 0.71, IS = 0.66). Further, all substitutions to unmatched aminoacids only occur at the tips with IN,Unmatched = IS,Unmatched = 0.Retaining our focus on this patient-window-site combination, we note the util-ity of deep sequencing compared to conventional bulk Sanger sequencing that pre-dominates the study of HIV adaptation to the CTL-mediated immune response.With deep sequencing, we were able to pick up the change in amino acid con-sensus sequence from baseline to followup as isoleucene (I)→ leucine (L). In thebaseline sample, 79% of sequences carried the codon ATA encoding I at this posi-tion. Subsequently in the follow-up sample, only 44% of sequences carried ATA,while 56% encoded L (TTA 37%, CTA 18%, CTG & CTT & TTG <1%). Depend-71Figure 3.13: Onset of directional selection. All colors not listed in legendare Unmatched amino acid lineages. Phylogeny created from windowof alignment containing baseline and followup sequences from patientBC27. Window coordinates 271-570bp, 1-based with respect to HXB2gag. Majority baseline sampling tips code for isoleucine (I) and getreplaced by blue or purple tips at follow sampling, which code forleucine (L).ing on the mixture calling algorithm, bulk Sanger sequencing may not have beenable to differentiate between I and L amino acids at followup.From the various statistics implemented and examined, we recommend the useof IN and IS for detection of directional selection for within-host populations overdN− dS and frequency based methods. Given a longitudinal sequencing datasetwith sufficient fragment lengths and genetic diversity, it will be superior to dN−dSin detecting temporal changes in nonsynonymous and synonymous substitutions,and will be able to determine if high prevalence variants were simply inheritedduring genetic bottlenecks.72Chapter 4ConclusionsUmberjack (Chapter 2) overcomes the phylogenetic challenges of missing homol-ogy in shotgun sequencing and has been validated on simulated datasets represent-ing a wide range of sequence quality parameters and sampling scenarios of HIVinfections.After successfully applying Umberjack to quantifying HIV evolution withinhosts in the context of the CTL immune response (Chapter 3), significant asso-ciations with patient HLA using the phylogenetic approaches were found in ad-dition to a frequency based approach. Although gag and nef were >95% con-served at most sites, Umberjack was able to detect directional selection in their sitesknown to be under CTL pressure. Further, we found evidence that the HIV popu-lations from all patients experienced reversions to the wild type HIV1-B consensussequence close to the most recent ancestor of the sampled populations. Thesewild type amino acids were further maintained through synonymous substitutions.However, very little reversion to unmatched wild-type occurred from matched es-cape variants, indicating that sites seldom experienced switches in selective pres-sures between general fitness requirements and evasion of the CTL response. Indrug treatment datasets (Section 3.4), the phylogenetic statistics generated by Um-berjack were able to highlight sites undergoing directional selection, as well asshowcase the stronger selective pressures of drug treatment compared to the im-mune response at the chronic stage of infection.734.1 Future DirectionsAlthough second generation high throughput short reads can reveal the effect ofevolutionary pressures on viral populations, their shortcomings in read length re-quire additional processing such as Umberjack in order to be useful. With thirdgeneration reads coming into play, such as 25kbp Single Molecule Real Timereads [25] and Oxford Nanopore reads [90], we have the ability to cover the en-tire genomes of entire viruses such as HIV [29]. Although their high error rates ofover 15% for Single Molecule Real Time reads and over 30% for Oxford Nanoporereads makes them currently prohibitive for phylogenetic analysis [43], they pose apromising method to connect distant regions of the genome, which we can leverageto examine genetic properties such as linkage disequilibrium. New techniques toreduce the error rate of Single Molecular Real Time reads have come forward, suchas circularized consensus sequencing in which DNA is circularized and sequencedmultiple times in a loop. Using this technique, error rates have been reported to bereduced to 2.5%, albeit with a reduction in read length to 2.5kbp [43].Using phylogenetic profiling, we can examine regions of viral genomes thatremain conserved through various evolutionary pressures, allowing us to designbetter vaccines. At present, there have been no successful antibody-based or CTL-based vaccines. The most successful antibody-based vaccine thus far, the ThaiRV144 vaccine, demonstrated 31% effectiveness in prevention of HIV infection[39]. The CTL-based Step vaccine demonstrated no overall protection against HIVinfection, and only small reductions in viral load during acute phase within certainpatient HLA groups. Comparison of genetic distance between the Step vaccineepitope inserts and the HIV populations of the vaccinated and placebo patientsrevealed that vaccinated populations were more distant from the vaccine insert,suggesting a higher rate of escape mutations amongst vaccinated patients [24].CTL based vaccine designs follow either a mosaic design which uses multi-ple common epitope variants as vaccine inserts to induce a broad CTL response,or conserved design which uses conserved epitopes to reduce mutations for CTLescape [69]. A successful vaccine design should not only take into account evo-lution within the HIV population, but evolution within the host immune system aswell. Different lines of T-cells recognize different epitopes, and competition ex-74ists amongst these lineages such that T-cells that successfully target more epitopesproliferate to become the dominant T-cells within the immune system. DominantT-cell lines have been reported to target highly mutating epitopes more than con-served epitopes, especially during acute infection, allowing escape mutations torapidly expand within the HIV population. Supporters for conserved design vac-cines argue that dominance of T-cell lines that prefer highly mutating epitopesbecomes less of an issue if vaccine epitopes are all conserved, since all T-cell re-sponses would focus on conserved regions, reducing early escape [59]. Further,mosaic designs would be ineffective against infections in which matched CTL es-cape mutations are transmitted between patients [12]. Using phylogenetic analysisto confirm that epitopes are conserved due to fitness as opposed to founder effectshelps ensure the best epitopes are selected for a successful vaccine.75Bibliography[1] C. S. Alexander, W. Dong, K. Chan, N. Jahnke, M. V. O’Shaughnessy,T. Mo, M. A. Piaseczny, J. S. Montaner, and P. R. Harrigan. Hiv proteaseand reverse transcriptase variation and therapy outcome inantiretroviral-naive individuals from a large north american cohort. AIDS,15(5):601–7, Mar 2001. → pages 49[2] N. H. Barton. Genetic hitchhiking. Philosophical Transactions of theRoyal Society of London B: Biological Sciences, 355(1403):1553–1562,2000. → pages 70[3] N. Beerenwinkel, M. Da¨umer, M. Oette, K. Korn, D. Hoffmann, R. Kaiser,T. Lengauer, J. Selbig, and H. Walter. Geno2pheno: estimating phenotypicdrug resistance from hiv-1 genotypes. Nucleic acids research, 31(13):3850–3855, 2003. → pages 69[4] C. Bleidorn. Third generation sequencing: technology and its potentialimpact on evolutionary biodiversity research. Systematics and Biodiversity,14(1):1–8, 2016. → pages 4[5] S. Bonhoeffer, R. M. May, G. M. Shaw, and M. A. Nowak. Virus dynamicsand drug therapy. Proceedings of the National Academy of Sciences, 94(13):6971–6976, 1997. → pages 69[6] P. Borrow, H. Lewicki, B. H. Hahn, G. M. Shaw, and M. B. Oldstone.Virus-specific CD8+ cytotoxic T-lymphocyte activity associated withcontrol of viremia in primary human immunodeficiency virus type 1infection. J Virol, 68(9):6103–10, Sep 1994. → pages 47[7] D. Branton, D. W. Deamer, A. Marziali, H. Bayley, S. A. Benner, T. Butler,M. Di Ventra, S. Garaj, A. Hibbs, X. Huang, et al. The potential andchallenges of nanopore sequencing. Nature biotechnology, 26(10):1146–1153, 2008. → pages 476[8] L. Breiman. Random forests. Machine learning, 45(1):5–32, 2001. →pages 38[9] W. Brockman, P. Alvarez, S. Young, M. Garber, G. Giannoukos, W. L. Lee,C. Russ, E. S. Lander, C. Nusbaum, and D. B. Jaffe. Quality scores and snpdetection in sequencing-by-synthesis systems. Genome research, 18(5):763–770, 2008. → pages 22[10] D. S. Burke. Recombination in hiv: an important viral evolutionarystrategy. Emerging infectious diseases, 3(3):253, 1997. → pages 7[11] J. M. Carlson, C. J. Brumme, E. Martin, J. Listgarten, M. A. Brockman,A. Q. Le, C. Chui, L. A. Cotton, D. J. Knapp, S. A. Riddler, et al.Correlates of protective cellular immunity revealed by analysis ofpopulation-level immune escape pathways in HIV-1. Journal of virology,pages JVI–01998, 2012. → pages 2, 48, 51, 60[12] J. M. Carlson, V. Y. Du, N. Pfeifer, A. Bansal, V. Y. Tan, K. Power, C. J.Brumme, A. Kreimer, C. E. DeZiel, N. Fusi, et al. Impact of pre-adaptedhiv transmission. Nature medicine, 22(6):606–613, 2016. → pages 75[13] K. Chen and L. Pachter. Bioinformatics for whole-genome shotgunsequencing of microbial communities. PLoS Comput Biol, 1(2):e24, 2005.→ pages 11[14] J. M. Cuevas, R. Geller, R. Garijo, J. Lo´pez-Aldeguer, and R. Sanjua´n.Extremely high mutation rate of hiv-1 in vivo. PLoS Biol, 13(9):e1002251,2015. → pages 33[15] A. G. Dalgleish, P. Beverley, P. R. Clapham, D. H. Crawford, M. F.Greaves, and R. A. Weiss. The cd4 (t4) antigen is an essential componentof the receptor for the aids retrovirus. Nature, 312(5996):763–767, 1983.→ pages 1[16] D. Darriba, M. Weiß, and A. Stamatakis. Prediction of missing sequencesand branch lengths in phylogenomic data. Bioinformatics, page btv768,2016. → pages 12[17] A. De Bruyn, D. P. Martin, and P. Lefeuvre. Phylogenetic reconstructionmethods: an overview. Molecular Plant Taxonomy: Methods andProtocols, pages 257–277, 2014. → pages 777[18] H. Deng, R. Liu, W. Ellmeier, S. Choe, D. Unutmaz, M. Burkhart, P. D.Marzio, S. Marmon, R. E. Sutton, C. M. Hill, et al. Identification of amajor co-receptor for primary isolates of hiv-1. 1996. → pages 2[19] K. Deng, M. Pertea, A. Rongvaux, L. Wang, C. M. Durand, G. Ghiaur,J. Lai, H. L. McHugh, H. Hao, H. Zhang, et al. Broad ctl response isrequired to clear latent hiv-1 due to dominance of escape mutations.Nature, 517(7534):381–385, 2015. → pages 1[20] R. Desper and O. Gascuel. Fast and accurate phylogeny reconstructionalgorithms based on the minimum-evolution principle. Journal ofcomputational biology, 9(5):687–705, 2002. → pages 24[21] G. Doitsh, N. L. Galloway, X. Geng, Z. Yang, K. M. Monroe, O. Zepeda,P. W. Hunt, H. Hatano, S. Sowinski, I. Mun˜oz-Arias, et al. Cell death bypyroptosis drives cd4 t-cell depletion in hiv-1 infection. Nature, 505(7484):509–514, 2014. → pages 1[22] P. Dorr, M. Westby, S. Dobbs, P. Griffin, B. Irvine, M. Macartney, J. Mori,G. Rickett, C. Smith-Burchnell, C. Napier, et al. Maraviroc (uk-427,857), apotent, orally bioavailable, and selective small-molecule inhibitor ofchemokine receptor ccr5 with broad-spectrum anti-humanimmunodeficiency virus type 1 activity. Antimicrobial agents andchemotherapy, 49(11):4721–4732, 2005. → pages 2[23] A. J. Drummond and A. Rambaut. BEAST: Bayesian evolutionary analysisby sampling trees. BMC Evol Biol, 7:214, 2007.doi:10.1186/1471-2148-7-214. → pages 52[24] P. T. Edlefsen, P. B. Gilbert, and M. Rolland. Sieve analysis in hiv-1vaccine efficacy trials. Current opinion in HIV and AIDS, 8(5), 2013. →pages 74[25] J. Eid, A. Fehr, J. Gray, K. Luong, J. Lyle, G. Otto, P. Peluso, D. Rank,P. Baybayan, B. Bettman, et al. Real-time dna sequencing from singlepolymerase molecules. Science, 323(5910):133–138, 2009. → pages 4, 74[26] B. Ewing, L. Hillier, M. C. Wendl, and P. Green. Base-calling of automatedsequencer traces usingphred. i. accuracy assessment. Genome research, 8(3):175–185, 1998. → pages 22[27] M. R. Farhat, B. J. Shapiro, K. J. Kieser, R. Sultana, K. R. Jacobson, T. C.Victor, R. M. Warren, E. M. Streicher, A. Calver, A. Sloutsky, D. Kaur,78J. E. Posey, B. Plikaytis, M. R. Oggioni, J. L. Gardy, J. C. Johnston,M. Rodrigues, P. K. C. Tang, M. Kato-Maeda, M. L. Borowsky,B. Muddukrishna, B. N. Kreiswirth, N. Kurepina, J. Galagan, S. Gagneux,B. Birren, E. J. Rubin, E. S. Lander, P. C. Sabeti, and M. Murray. Genomicanalysis identifies targets of convergent positive selection in drug-resistantmycobacterium tuberculosis. Nat Genet, 45(10):1183–9, Oct 2013.doi:10.1038/ng.2747. → pages 14[28] J. Felsenstein. Evolutionary trees from dna sequences: a maximumlikelihood approach. Journal of molecular evolution, 17(6):368–376, 1981.→ pages 7[29] E. B. Fichot and R. S. Norman. Microbial phylogenetic profiling with thepacific biosciences sequencing platform. Microbiome, 1(1):1, 2013. →pages 74[30] W. Fletcher and Z. Yang. Indelible: a flexible simulator of biologicalsequence evolution. Molecular biology and evolution, 26(8):1879–1888,2009. → pages 33[31] M. K. C. from Jed Wing, S. Weston, A. Williams, C. Keefer, A. Engelhardt,T. Cooper, Z. Mayer, B. Kenkel, the R Core Team, M. Benesty,R. Lescarbeau, A. Ziem, L. Scrucca, Y. Tang, and C. Candan. caret:Classification and Regression Training, 2016. URLhttp://CRAN.R-project.org/package=caret. R package version 6.0-64. →pages 39[32] Y.-X. Fu and W.-H. Li. Statistical tests of neutrality of mutations. Genetics,133(3):693–709, 1993. → pages 15, 16, 59[33] V. V. Ganusov, N. Goonetilleke, M. K. Liu, G. Ferrari, G. M. Shaw, A. J.McMichael, P. Borrow, B. T. Korber, and A. S. Perelson. Fitness costs anddiversity of the cytotoxic t lymphocyte (ctl) response determine the rate ofctl escape during acute and chronic phases of hiv infection. Journal ofvirology, 85(20):10518–10528, 2011. → pages 47[34] N. Goldman and Z. Yang. A codon-based model of nucleotide substitutionfor protein-coding dna sequences. Molecular biology and evolution, 11(5):725–736, 1994. → pages 15[35] J. M. Gonzalez, M. C. Portillo, P. Belda-Ferre, and A. Mira. Amplificationby pcr artificially reduces the proportion of the rare biosphere in microbialcommunities. PLoS One, 7(1):e29973, 2012. → pages 1179[36] P. J. Goulder and D. I. Watkins. HIV and SIV CTL escape: implications forvaccine design. Nature Reviews Immunology, 4(8):630–640, 2004. →pages 1, 2, 47, 56[37] P. J. Goulder and D. I. Watkins. Impact of mhc class i diversity on immunecontrol of immunodeficiency virus replication. Nature ReviewsImmunology, 8(8):619–630, 2008. → pages 61[38] Y. Han, M. Wind-Rotolo, H.-C. Yang, J. D. Siliciano, and R. F. Siliciano.Experimental approaches to the study of hiv-1 latency. Nature ReviewsMicrobiology, 5(2):95–106, 2007. → pages 1[39] B. F. Haynes, P. B. Gilbert, M. J. McElrath, S. Zolla-Pazner, G. D. Tomaras,S. M. Alam, D. T. Evans, D. C. Montefiori, C. Karnasuta, R. Sutthent, et al.Immune-correlates analysis of an hiv-1 vaccine efficacy trial. New EnglandJournal of Medicine, 366(14):1275–1286, 2012. → pages 74[40] J. Hermisson and P. S. Pennings. Soft sweeps: molecular populationgenetics of adaptation from standing genetic variation. Genetics, 169(4):2335–52, Apr 2005. doi:10.1534/genetics.104.036947. → pages 17[41] A. E. Hoerl and R. W. Kennard. Ridge regression: Biased estimation fornonorthogonal problems. Technometrics, 12(1):55–67, 1970. → pages 40[42] W. Huang, L. Li, J. R. Myers, and G. T. Marth. Art: a next-generationsequencing read simulator. Bioinformatics, 28(4):593–594, 2012. → pages33[43] X. Jiao, X. Zheng, L. Ma, G. Kutty, E. Gogineni, Q. Sun, B. T. Sherman,X. Hu, K. Jones, C. Raley, et al. A benchmark study on error assessmentand quality control of ccs reads derived from the pacbio rs. Journal of datamining in genomics & proteomics, 4(3), 2013. → pages 74[44] T. H. Jukes and C. R. Cantor. Evolution of protein molecules. Mammalianprotein metabolism, 3(21):132, 1969. → pages 8[45] S. M. Krone and C. Neuhauser. Ancestral processes with selection.Theoretical population biology, 51(3):210–237, 1997. → pages 32[46] S. Kryazhimskiy and J. B. Plotkin. The population genetics of dN/dS.PLoS Genet, 4(12):e1000304, Dec 2008.doi:10.1371/journal.pgen.1000304. → pages 5980[47] C. Kuiken, T. Leitner, B. Foley, B. Hahn, P. Marx, F. McCutchan,S. Wolinsky, B. Korber, G. Bansal, and W. Abfalterer. Hiv sequencecompendium 2009. Los Alamos, New Mexico: Los Alamos NationalLaboratory, Theoretical Biology and Biophysics, 2009. → pages 64[48] A. S. Lauring, J. Frydman, and R. Andino. The role of mutationalrobustness in rna virus evolution. Nature reviews Microbiology, 11(5):327–336, 2013. → pages 20[49] J. Lee, P. Mu¨ller, S. Sengupta, K. Gulukota, and Y. Ji. Bayesian featureallocation models for tumor heterogeneity. In Statistical Analysis forHigh-Dimensional Data, pages 211–232. Springer, 2016. → pages 12[50] T. Lengauer, O. Sander, S. Sierra, A. Thielen, R. Kaiser, et al.Bioinformatics prediction of hiv coreceptor usage. Nature biotechnology,25(12):1407–1410, 2007. → pages 3[51] A. Leslie, K. Pfafferott, P. Chetty, R. Draenert, M. Addo, M. Feeney,Y. Tang, E. Holmes, T. Allen, J. Prado, et al. Hiv evolution: Ctl escapemutation and reversion after transmission. Nature medicine, 10(3):282–289, 2004. → pages 64[52] H. Li. Aligning sequence reads, clone sequences and assembly contigswith bwa-mem. arXiv preprint arXiv:1303.3997, 2013. → pages 6, 33, 49[53] H. Li, B. Handsaker, A. Wysoker, T. Fennell, J. Ruan, N. Homer, G. Marth,G. Abecasis, R. Durbin, et al. The sequence alignment/map format andsamtools. Bioinformatics, 25(16):2078–2079, 2009. → pages 21[54] B. Liang, M. Luo, J. Scott-Herridge, C. Semeniuk, M. Mendoza,R. Capina, B. Sheardown, H. Ji, J. Kimani, B. T. Ball, et al. A comparisonof parallel pyrosequencing and sanger clone-based sequencing and itsimpact on the characterization of the genetic diversity of hiv-1. PloS one, 6(10):e26745, 2011. → pages 3, 11[55] H.-X. Liao, R. Lynch, T. Zhou, F. Gao, S. M. Alam, S. D. Boyd, A. Z. Fire,K. M. Roskin, C. A. Schramm, Z. Zhang, et al. Co-evolution of a broadlyneutralizing hiv-1 antibody and founder virus. Nature, 496(7446):469–476,2013. → pages 14, 51[56] A. Liaw and M. Wiener. Classification and regression by randomforest. Rnews, 2(3):18–22, 2002. → pages 38, 3981[57] D. J. Lipman, S. F. Altschul, and J. D. Kececioglu. A tool for multiplesequence alignment. Proceedings of the National Academy of Sciences, 86(12):4412–4415, 1989. → pages 5[58] J. Listgarten, Z. Brumme, C. Kadie, G. Xiaojiang, B. Walker,M. Carrington, P. Goulder, and D. Heckerman. Statistical resolution ofambiguous hla typing data. PLoS Comput Biol, 4(2):e1000016, 2008. →pages 49[59] M. K. Liu, N. Hawkins, A. J. Ritchie, V. V. Ganusov, V. Whale,S. Brackenridge, H. Li, J. W. Pavlicek, F. Cai, M. Rose-Abrahams, et al.Vertical t cell immunodominance and epitope entropy determine hiv-1escape. The Journal of clinical investigation, 123(1):380–393, 2013. →pages 75[60] S. L. Liu, A. G. Rodrigo, R. Shankarappa, G. H. Learn, L. Hsu,O. Davidov, L. P. Zhao, and J. I. Mullins. HIV quasispecies andresampling. Science, 273(5274):415–6, Jul 1996. → pages 48[61] Y. Liu, J. McNevin, J. Cao, H. Zhao, I. Genowati, K. Wong,S. McLaughlin, M. D. McSweyn, K. Diem, C. E. Stevens, et al. Selectionon the human immunodeficiency virus type 1 proteome following primaryinfection. Journal of virology, 80(19):9519–9529, 2006. → pages 51[62] C. Luo, D. Tsementzi, N. Kyrpides, T. Read, and K. T. Konstantinidis.Direct comparisons of illumina vs. roche 454 sequencing technologies onthe same microbial community dna sample. PloS one, 7(2):e30087, 2012.→ pages 4[63] F. Maldarelli, M. Kearney, S. Palmer, R. Stephens, J. Mican, M. A. Polis,R. T. Davey, J. Kovacs, W. Shao, D. Rock-Kress, et al. Hiv populations arelarge and accumulate high genetic diversity in a nonlinear fashion. Journalof virology, 87(18):10313–10323, 2013. → pages 7, 12[64] L. M. Mansky and H. M. Temin. Lower in vivo mutation rate of humanimmunodeficiency virus type 1 than that predicted from the fidelity ofpurified reverse transcriptase. J Virol, 69(8):5087–94, Aug 1995. → pages55[65] L. M. Mansky and H. M. Temin. Lower in vivo mutation rate of humanimmunodeficiency virus type 1 than that predicted from the fidelity ofpurified reverse transcriptase. Journal of virology, 69(8):5087–5094, 1995.→ pages 3382[66] D. Martin, D. Posada, K. Crandall, and C. Williamson. A modifiedbootscan algorithm for automated identification of recombinant sequencesand recombination breakpoints. AIDS Research & Human Retroviruses, 21(1):98–102, 2005. → pages 19[67] J. H. McDonald. Handbook of biological statistics, volume 3. SparkyHouse Publishing Baltimore, MD, 2014. → pages 37[68] R. A. McGovern, J. Symons, A. F. Poon, P. R. Harrigan, S. F. van Lelyveld,A. I. Hoepelman, P. M. van Ham, W. Dong, A. M. Wensing, andM. Nijhuis. Maraviroc treatment in non-r5-hiv-1-infected patients results inthe selection of extreme cxcr4-using variants with limited effect on the totalviral setpoint. Journal of Antimicrobial Chemotherapy, page dkt153, 2013.→ pages 11, 68, 69[69] A. J. McMichael and W. C. Koff. Vaccines that stimulate t cell immunity tohiv-1: the next step. Nature immunology, 15(4):319, 2014. → pages 74[70] A. J. McMichael and S. L. Rowland-Jones. Cellular immune responses toHIV. Nature, 410(6831):980–987, 2001. → pages 2[71] A. G. Meyer and C. O. Wilke. The utility of protein structure as a predictorof site-wise dn/ds varies widely among hiv-1 proteins. Journal of TheRoyal Society Interface, 12(111):20150579, 2015. → pages 33[72] M. Mild, R. R. Gray, A. Kvist, P. Lemey, M. M. Goodenow, E. M. Fenyo¨,J. Albert, M. Salemi, J. Esbjo¨rnsson, and P. Medstrand. High intrapatienthiv-1 evolutionary rate is associated with ccr5-to-cxcr4 coreceptor switch.Infection, Genetics and Evolution, 19:369–377, 2013. → pages 3[73] R. Mills, M. Rozanov, A. Lomsadze, T. Tatusova, and M. Borodovsky.Improving gene annotation of complete viral genomes. Nucleic acidsresearch, 31(23):7041–7055, 2003. → pages 20[74] J. Montaner. Therapeutic guidelines for antiretroviral (arv) treatment ofadult hiv infection, 2015. → pages 3[75] P. A. P. Moran et al. The statistical processes of evolutionary theory. Thestatistical processes of evolutionary theory., 1962. → pages 32[76] S. V. Muse and B. S. Gaut. A likelihood approach for comparingsynonymous and nonsynonymous nucleotide substitution rates, withapplication to the chloroplast genome. Molecular Biology and Evolution,11(5):715–724, 1994. → pages 2583[77] R. A. Neher and T. Leitner. Recombination rate and selection strength inhiv intra-patient evolution. PLoS Comput Biol, 6(1):e1000660, 2010. →pages 7, 32, 41[78] M. Nei and T. Gojobori. Simple methods for estimating the numbers ofsynonymous and nonsynonymous nucleotide substitutions. Molecularbiology and evolution, 3(5):418–426, 1986. → pages 26[79] P. S. Pennings, S. Kryazhimskiy, and J. Wakeley. Loss and recovery ofgenetic diversity in adapting populations of hiv. PLoS Genet, 10(1):e1004000, 2014. → pages 7, 12[80] S. L. K. Pond and S. D. Frost. Not so different after all: a comparison ofmethods for detecting amino acid sites under selection. Molecular biologyand evolution, 22(5):1208–1222, 2005. → pages 15, 52, 70[81] S. L. K. Pond and S. D. Frost. Not so different after all: a comparison ofmethods for detecting amino acid sites under selection. Molecular biologyand evolution, 22(5):1208–1222, 2005. → pages 26, 30[82] S. L. K. Pond and S. V. Muse. Hyphy: hypothesis testing usingphylogenies. In Statistical methods in molecular evolution, pages 125–181.Springer, 2005. → pages 25[83] A. F. Poon, L. C. Swenson, E. M. Bunnik, D. Edo-Matas, H. Schuitemaker,A. B. van’t Wout, and P. R. Harrigan. Reconstructing the dynamics of hivevolution within hosts from serial deep sequence data. PLoS Comput Biol,8(11):e1002753, 2012. → pages 102[84] A. F. Y. Poon, L. C. Swenson, W. W. Y. Dong, W. Deng, S. L.Kosakovsky Pond, Z. L. Brumme, J. I. Mullins, D. D. Richman, P. R.Harrigan, and S. D. W. Frost. Phylogenetic analysis of population-basedand deep sequencing data to identify coevolving sites in the nef gene ofHIV-1. Mol Biol Evol, 27(4):819–32, Apr 2010.doi:10.1093/molbev/msp289. → pages 4, 11[85] S. Prabhakaran, M. Rey, O. Zagordi, N. Beerenwinkel, and V. Roth. Hivhaplotype inference using a propagating dirichlet process mixture model.IEEE/ACM Transactions on Computational Biology and Bioinformatics(TCBB), 11(1):182–191, 2014. → pages 12[86] D. A. Price, P. J. Goulder, P. Klenerman, A. K. Sewell, P. J. Easterbrook,M. Troop, C. R. Bangham, and R. E. Phillips. Positive selection of HIV-184cytotoxic T lymphocyte escape variants during primary infection.Proceedings of the National Academy of Sciences, 94(5):1890–1895, 1997.→ pages 47[87] M. N. Price, P. S. Dehal, and A. P. Arkin. Fasttree: computing largeminimum evolution trees with profiles instead of a distance matrix.Molecular biology and evolution, 26(7):1641–1650, 2009. → pages 10, 24[88] M. N. Price, P. S. Dehal, and A. P. Arkin. Fasttree 2–approximatelymaximum-likelihood trees for large alignments. PloS one, 5(3):e9490,2010. → pages 24[89] T. Pupko, I. Pe’er, R. Shamir, and D. Graur. A fast algorithm for jointreconstruction of ancestral amino acid sequences. Mol Biol Evol, 17(6):890–6, Jun 2000. → pages 25[90] J. Quick, A. R. Quinlan, and N. J. Loman. A reference bacterial genomedataset generated on the minion portable single-molecule nanoporesequencer. Gigascience, 3(1):1, 2014. → pages 74[91] A. Ramı´rez-Soriano, S. E. Ramos-Onsins, J. Rozas, F. Calafell, andA. Navarro. Statistical power analysis of neutrality tests underdemographic expansions, contractions and bottlenecks with recombination.Genetics, 179(1):555–567, 2008. → pages 15, 70[92] D. Robinson and L. R. Foulds. Comparison of phylogenetic trees.Mathematical biosciences, 53(1):131–147, 1981. → pages 34[93] N. Saitou and M. Nei. The neighbor-joining method: a new method forreconstructing phylogenetic trees. Molecular biology and evolution, 4(4):406–425, 1987. → pages 7[94] M. Salminen, J. Carr, D. Burke, and F. McCutchan. Identification ofbreakpoints in intergenotypic recombinants of HIV type 1 by bootscanning.AIDS research and human retroviruses, 11(11):1423, 1995. → pages 7, 19[95] M. Schirmer, R. D’Amore, U. Z. Ijaz, N. Hall, and C. Quince. Illuminaerror profiles: resolving fine-scale variation in metagenomic sequencingdata. BMC Bioinformatics, 17(1):1, 2016. → pages 23[96] R. Shankarappa, J. B. Margolick, S. J. Gange, A. G. Rodrigo, D. Upchurch,H. Farzadegan, P. Gupta, C. R. Rinaldo, G. H. Learn, X. He, X. L. Huang,and J. I. Mullins. Consistent viral evolutionary changes associated with the85progression of human immunodeficiency virus type 1 infection. J Virol, 73(12):10489–502, Dec 1999. → pages 2[97] R. Shankarappa, J. B. Margolick, S. J. Gange, A. G. Rodrigo, D. Upchurch,H. Farzadegan, P. Gupta, C. R. Rinaldo, G. H. Learn, X. He, et al.Consistent viral evolutionary changes associated with the progression ofhuman immunodeficiency virus type 1 infection. Journal of virology, 73(12):10489–10502, 1999. → pages 51[98] E. Simon-Loriere and E. C. Holmes. Gene duplication is infrequent in therecent evolutionary history of rna viruses. Molecular biology andevolution, page mst044, 2013. → pages 20[99] R. Smyth, T. Schlub, A. Grimm, V. Venturi, A. Chopra, S. Mallal,M. Davenport, and J. Mak. Reducing chimera formation during pcramplification to ensure accurate genotyping. Gene, 469(1):45–51, 2010. →pages 50[100] C. Spearman. ” general intelligence,” objectively determined and measured.The American Journal of Psychology, 15(2):201–292, 1904. → pages 37[101] M. Stein. Large sample properties of simulations using latin hypercubesampling. Technometrics, 29(2):143–151, 1987. → pages 31[102] S. Tavare´. Some probabilistic and statistical problems in the analysis of dnasequences. Lectures on mathematics in the life sciences, 17:57–86, 1986.→ pages 8[103] R. Tibshirani. Regression shrinkage and selection via the lasso. Journal ofthe Royal Statistical Society. Series B (Methodological), pages 267–288,1996. → pages 40[104] J. J. Wiens. Missing data and the design of phylogenetic analyses. Journalof biomedical informatics, 39(1):34–42, 2006. → pages 7, 10[105] J. J. Wiens and M. C. Morrill. Missing data in phylogenetic analysis:reconciling results from simulations and empirical data. SystematicBiology, page syr025, 2011. → pages 10[106] O. Zagordi, A. Bhattacharya, N. Eriksson, and N. Beerenwinkel. Shorah:estimating the genetic diversity of a mixed sample from next-generationsequencing data. BMC bioinformatics, 12(1):1, 2011. → pages 1286[107] F. Zanini, J. Brodin, L. Thebo, C. Lanz, G. Bratt, J. Albert, and R. A.Neher. Population genomics of intrapatient hiv-1 evolution. eLife, 4:e11282, 2016. → pages 11, 64[108] S. Zhou, C. Jones, P. Mieczkowski, and R. Swanstrom. Primer id validatestemplate sampling depth and greatly reduces the error rate ofnext-generation sequencing of hiv-1 genomic rna populations. Journal ofvirology, 89(16):8540–8555, 2015. → pages 4887Appendix ASupporting MaterialsA.1 Simulated Datasets For Predicting UmberjackAccuracyFeature Min MaxGenerations 128 8192Recombination Rate (recombinations/bp/generation) 0 5.1e-5Mutation Rate 1 (mutations/bp/generation) 5.6e-7 8.4e-4Mutation Rate 2 (mutations/bp/generation) 5.6e-7 8.4e-4Mutation Rate 3 (mutations/bp/generation) 5.6e-7 8.4e-4Umberjack Config a bRead Coverage Per Individual 0 3Mean Sequence Fragment Size (bp) 104 500Table A.1: Parameter ranges for Latin hypercube sampling to generate sim-ulated datasets to test Umberjack accuracy. Three randomly allocatedcontiguous sections of the genome were assigned a different mutationrate selected from the Latin hypercube sampling to simulate differentmutation rates per gene. Umberjack Configuration a = [Window size =150bp, Min window width coverage = 0.7, Min window read depth = 10,Min phred quality score = 15]. Umberjack Configuration b = [Windowsize = 300bp, Min window width coverage = 0.875, Min window readdepth = 10, Min phred quality score = 20].88Table A.2Generations MutationRate 1MutationRate 2MutationRate 3Recomb.RateBreakpoints Coverage MeanFragmentSize283 1.03e-05 3.63e-05 1.87e-05 3.53e-05 3 5 1351230 6.27e-05 5.54e-06 1.44e-05 1.08e-05 4 0 250307 4.62e-05 6.33e-05 6.39e-05 2.17e-05 2 2 353412 2.94e-05 5.37e-05 6.53e-05 2.43e-05 3 5 220784 2.73e-05 4.26e-05 5.61e-05 1.70e-05 4 0 2774491 5.68e-05 3.15e-06 7.43e-05 2.23e-06 3 1 120474 7.05e-06 3.04e-05 4.82e-05 0.00e+00 0 2 4252343 4.92e-05 5.68e-05 1.31e-05 5.69e-06 4 2 1153614 2.13e-05 1.42e-05 3.91e-05 2.77e-06 3 3 245675 3.32e-05 7.56e-05 8.17e-05 9.88e-06 2 2 487502 3.43e-05 1.67e-05 5.03e-05 2.66e-05 4 6 192881 8.34e-05 2.43e-05 1.12e-05 7.57e-06 2 0 1701060 5.35e-05 8.26e-06 7.34e-06 3.14e-06 1 0 4031174 4.20e-05 5.15e-05 4.43e-05 1.14e-05 4 2 2054817 3.81e-05 6.25e-05 5.75e-05 1.38e-06 2 3 1646284 3.45e-06 3.17e-05 3.78e-05 5.30e-07 1 0 1842620 5.27e-05 1.06e-05 3.41e-05 1.27e-06 1 3 141Continued on next page89Table A.2 – Continued from previous pageGenerations MutationRate 1MutationRate 2MutationRate 3Recomb.RateBreakpoints Coverage MeanFragmentSize775 7.14e-05 3.76e-05 4.54e-05 1.29e-05 3 4 4341813 1.19e-05 2.92e-05 7.80e-05 5.52e-06 3 1 289714 6.57e-05 5.12e-06 1.02e-05 1.40e-05 3 2 2141308 6.39e-05 6.76e-05 3.20e-05 5.10e-06 2 0 3581000 8.05e-06 4.55e-05 5.90e-05 6.67e-06 2 4 377367 6.69e-05 4.47e-05 6.96e-05 0.00e+00 0 3 4092788 7.09e-05 7.42e-05 6.91e-06 3.59e-06 3 1 285339 4.14e-05 5.84e-05 2.74e-05 1.97e-05 2 6 3821597 7.64e-05 1.33e-05 6.10e-05 6.26e-06 3 1 1581525 3.92e-05 2.23e-05 5.35e-05 4.37e-06 2 7 3915044 5.92e-05 8.30e-05 7.28e-05 1.32e-06 2 3 317317 2.36e-05 7.84e-05 2.12e-05 2.10e-05 2 0 3433424 7.32e-05 6.94e-05 7.95e-05 1.95e-06 2 2 370547 7.76e-05 5.56e-05 7.70e-05 1.83e-05 3 6 4935833 8.19e-05 9.96e-06 2.65e-05 5.71e-07 1 7 1787772 3.58e-05 1.82e-05 1.84e-05 4.29e-07 1 1 483268 4.80e-05 2.75e-05 8.41e-05 1.24e-05 1 0 335Continued on next page90Table A.2 – Continued from previous pageGenerations MutationRate 1MutationRate 2MutationRate 3Recomb.RateBreakpoints Coverage MeanFragmentSize5668 2.61e-05 6.49e-05 1.54e-05 5.88e-07 1 1 261417 1.79e-05 7.20e-05 4.34e-06 0.00e+00 0 4 3027164 7.92e-05 1.96e-05 7.12e-05 0.00e+00 0 1 4401396 2.31e-05 5.07e-05 2.71e-06 0.00e+00 0 0 3293116 5.04e-05 3.91e-05 2.20e-05 3.21e-06 3 1 4566981 6.78e-05 8.22e-05 6.16e-05 0.00e+00 0 1 2311963 5.05e-06 4.69e-05 2.51e-05 5.09e-06 3 0 1052076 6.09e-05 3.40e-05 4.28e-05 1.61e-06 1 4 3232441 1.94e-05 2.59e-05 4.10e-05 1.37e-06 1 0 4711753 1.47e-05 2.05e-05 3.55e-05 1.90e-06 1 5 4154360 3.15e-05 6.05e-05 4.74e-05 7.65e-07 1 1 1503024 5.55e-05 4.00e-05 6.63e-05 2.20e-06 2 1 2663865 1.26e-05 7.65e-05 3.09e-05 8.62e-07 1 1 298561 1.60e-05 6.93e-05 2.95e-05 1.19e-05 2 0 228916 7.54e-05 8.08e-05 6.78e-05 0.00e+00 0 5 461624 4.40e-05 4.88e-05 5.17e-05 1.07e-05 2 1 450Continued on next page91Table A.2 – Continued from previous pageGenerations MutationRate 1MutationRate 2MutationRate 3Recomb.RateBreakpoints Coverage MeanFragmentSizeTable A.2: Parameters for generating simulated datasets to predict Umberjack accuracy. Each row represents a sim-ulated population and its simulated paired-end MiSeq sequence library. Recombination rate units in recombina-tions/bp/generation. Mutation rate units in mutations/bp/generation. Read coverage is per extent individual in thepopulation. Sequencing fragment size in bp. Genome size = 900bp. Selection rate = 0.01/generation. Extentpopulation size = 1000. 2x250bp MiSeq paired-end reads. Sequencing fragment size standard deviation = 100bp.Umberjack Configuration a = Window size = 150bp, Min window width coverage = 0.7, Min window read depth= 10, Min phred quality score = 15, Umberjack Configuration b = Window size = 300bp, Min window widthcoverage = 0.875, Min window read depth = 10, Min phred quality score = 20.92Table A.3Name Recombination RateRecombo1 1.29e-06Recombo2 1.29e-06Recombo3 1.29e-06Recombo4 1.29e-06Recombo5 1.29e-06Recombo6 1.29e-06Recombo7 0.00e+00Recombo8 0.00e+00Recombo9 0.00e+00Recombo10 6.45e-07Recombo11 6.45e-07Recombo12 6.45e-07Recombo13 3.23e-06Recombo14 3.23e-06Recombo15 3.23e-06Recombo16 6.45e-06Recombo17 6.45e-06Recombo18 6.45e-06Recombo19 6.45e-05Continued on next page93Table A.3 – Continued from previous pageName Recombination RateRecombo20 6.45e-05Recombo21 6.45e-05Recombo22 1.29e-04Recombo23 1.29e-04Recombo24 1.29e-04Recombo25 1.29e-05Recombo26 1.29e-05Recombo27 1.29e-05Table A.3: Parameters for generating simulated datasets to predict Umberjack accuracy under varying recombination.Each row represents a simulated population and its simulated paired-end MiSeq sequence library. Recombinationrate units in recombinations/bp/generation. Mutation rate units in mutations/bp/generation. Read coverage is perextent individual in the population. Sequencing fragment size in bp. Genome size = 930bp. Selection rate =0.01/generation. Mutation Rate = 4e-5 mutations/bp/generation. Generations = 5000. Extent population size =100. 2x250bp MiSeq paired-end reads. Read coverage per individual = 2x. Sequencing fragment size mean =375bp. Sequencing fragment size standard deviation = 75bp. Umberjack Configuration b = Window size = 300bp,Min window width coverage = 0.875, Min window read depth = 10, Min phred quality score = 20.94(a) Zoomed in to 1.5 × IQR of ∆ (b) All dataFigure A.1: ∆ vs Window Total Breakpoint Ratio. Refer to Feature (A) forfeature description.(a) Zoomed in to 1.5 × IQR of ∆ (b) All dataFigure A.2: ∆ vs Window-Site Unambiguous Codon Rate. Refer to Feature(B) for feature description.95(c) Zoomed in to 1.5 × IQR of ∆ (d) All dataFigure A.2: ∆ vs Window-Site Amino Acid Depth. Refer to Feature (C) forfeature description.(a) Zoomed in to 1.5 × IQR of ∆ (b) All dataFigure A.3: ∆ vs Window-Site Substitutions. Refer to Feature (D) for featuredescription.96(c) Zoomed in to 1.5 × IQR of ∆ (d) All dataFigure A.3: ∆ vs Window-Site Nonsynonymous Substitutions. Refer to Fea-ture (E) for feature description.(a) Zoomed in to 1.5 × IQR of ∆ (b) All dataFigure A.4: ∆ vs Window-Site Synonymous Substitutions. Refer to Feature(F) for feature description.97(c) Zoomed in to 1.5 × IQR of ∆ (d) All dataFigure A.4: ∆ vs Window-Site Expected Nonsynonymous Substitutions PerBranch. Refer to Feature (G) for feature description.(a) Zoomed in to 1.5 × IQR of ∆ (b) All dataFigure A.5: ∆ vs Window-Site Expected Synonymous Substitutions PerBranch. Refer to Feature (H) for feature description.98(c) Zoomed in to 1.5 × IQR of ∆ (d) All dataFigure A.5: ∆ vs Normalized Window Tree Length. Refer to Feature (I) forfeature description.(a) Zoomed in to 1.5 × IQR of ∆ (b) All dataFigure A.6: ∆ vs Normalized WRF . Refer to Feature (J) for feature descrip-tion.99(c) Zoomed in to 1.5 × IQR of ∆ (d) All dataFigure A.6: ∆ vs Normalized Polytomies. Refer to Feature (K) for featuredescription.(a) Zoomed in to 1.5 × IQR of ∆ (b) All dataFigure A.7: ∆ vs Codon Distribution P-value. P-values on x-axis are log10transformed. Refer to Feature (L) for feature description.100(c) Zoomed in to 1.5 × IQR of ∆ (d) All dataFigure A.7: ∆ vs True Site Codon Entropy. Refer to Feature (M) for featuredescription.(a) Zoomed in to 1.5 × IQR of ∆ (b) All dataFigure A.8: ∆ vs Window-Site Codon Entropy. Refer to Feature (N) for fea-ture description.101(c) Zoomed in to 1.5 × IQR of ∆ (d) All dataFigure A.8: ∆ vs Window-Site Codon Entropy. Refer to Feature (O) for fea-ture description.A.2 Untreated Patient MetadataWe estimated that 1.5% of the RNA templates contained in a patient plasma samplemade it to sequencing for NEF and GAG, and 4.4% for ENV. The estimates arebased on the fraction of volume of patient plasma used in viral RNA extraction andPCR amplification using a similar procedure specified in [83]. The percentage ofNEF and GAG templates sent to PCR amplification and sequencing was one thirdof the templates sent for ENV.102Table A.4Min Mean Max 25% 50% 75%Baseline Viral Load (copies/mL) 2.1×104 1.2×105 7.5×105 9.0×104 1.0×105 1.0×105Followup Viral Load (copies/mL) 2.1×104 1.1×104 4.2×105 7.8×104 1.0×105 1.0×105Baseline CD4 Counts (cells/mm3) 190 388 560 325 405 473Followup CD4 Counts (cells/mm3) 20.0 251 460 185 300 345Years Infected at Baseline 0.35 3.8 13 1.2 2.6 5.7Months Between Samples 3.0 16 51 8.4 12 17Table A.4: Untreated Patient Measurement Quantiles. 32 of the patient sample viral loads hit a viral load assay mea-surement upper limit of 105 copies/mL, and 1 patient sample viral load hit another assay limit of 7.5× 105copies/mL. None were remeasured using assays with a higher limits. Viral loads that hit upper measurementlimits were set to the limit in quantile calculations. CD4 samples were taken within 30 days of viral load samples.Baseline CD4 counts were missing for 14 patients, and followup CD4 counts were missing for 12 patients. ‘YearsInfected at Baseline’ are estimated from BEAST timing of the most recent common ancestor of sequencing reads.Patient BC11 is missing an estimate for ‘Years Infected at Baseline’ since the BEAST runs did not converge.103Table A.5gag nef envPrimer Positions 734 - 1833 bp 8776 - 9593 bp 6945 - 7373 bp, 8325 -8775bpMin Templates 310 310 920Mean Templates 1700 1700 5200Max Templates 11000 11000 3300025% Gene Depth Coverage 4500X 6200X 3400X50% Gene Depth Coverage 8200X 10000X 6500X75% Gene Depth Coverage 13000X 14000X 13000X25% Template DepthCoverage2.8X 4.1X 0.71X50% Template DepthCoverage5.6X 7.4X 1.6X75% Template DepthCoverage10X 11X 3.3XMean Fragment Size 137bp 147bp 104bpStd Dev Fragment Size 94.5bp 100bp 64.7bpTable A.5: Sequencing Statistics for Untreated Patient Samples. Primer positions are with respect to HXB2 referencestrain (Accession K03455.1)104Figure A.9: Phylogenetic tree of GAG consensus sequence of samples from untreated patients. Node values indicatebootstrap support. Branch length units in nucleotide substitutions per site.105Figure A.10: Phylogenetic tree of NEF consensus sequence of samples from untreated patients. Node values indicatebootstrap support. Branch length units in nucleotide substitutions per site.106Figure A.11: Phylogenetic tree of ENV consensus sequence of samples from untreated patients. Node values indicatebootstrap support. Branch length units in nucleotide substitutions per site.107Figure A.12: Cleaned Umberjack estimate of site dN− dS vs true dN− dS.Each point represents a site estimate of dN− dS. Sites for which theRandom Forest predictor deemed the Umberjack estimate inaccurateare excluded. Red line is y=x. Blue line is fitted line.A.3 HIV Selective Sweeps in Untreated PatientsA.4 Effect of Recombination on Umberjack Error108Figure A.13: Violin plots summarizing the distributions of site IN− IS statis-tic across subjects, broken down by gene. A IN− IS < 0 indicates thatnonsynonymous substitutions occur later than synonymous substitu-tions, which is consistent with purifying selection. Black point rangesindicate median and IQR. White squares indicate mean.109(a)(b)Figure A.14: Effect of recombination on ∆. Each point represents a windowof Umberjack estimates on simulated datasets focused on recombina-tion (Table A.3). Black line indicates linear regression fit. (a) Pointsare shaded according the the length of the window tree in units of nu-cleotide substitutions/site. (b) Points are shaded according to the re-combination rate ρ of population in units of recombinations/site/gen-eration. Interval denotes average codon distance between breakpointsacross the genome.110
Thesis/Dissertation
eng
Bioinformatics
Vancouver : University of British Columbia Library
Attribution-NonCommercial 4.0 International
Graduate
Overcoming missing data in phylogenetic analysis of shotgun sequencing to detect HIV adaptation to immune response
Text