Techniques for evaluating evolutionary stability of genome segments

Abstract

Evolutionary rate can be estimated for genome segments (or amino acid sequences) of any length, from a single codon (or amino acid position) to an entire protein or an entire genome. Given a reference sequence for a genome of a target at an initial time and one or more strain-specific sequences for respective genomes of one or more strains of the target in a time period after the initial time, a segment of the genome can be selected, and a genetic divergence parameter can be computed for the segment based on genetic distances between the segment in the reference sequence and a corresponding segment in each of the strain-specific sequences. Based on the genetic divergence parameter, an evolutionary rate or evolutionary stability parameter for the segment can be computed and used in design of antiviral therapies or vaccines and / or for estimation of tumor growth.

Description

TECHNIQUES FOR EVALUATING EVOLUTIONARY STABILITY OF GENOME SEGMENTSBACKGROUND

[0001] This disclosure relates to analysis of genome segments, including segments of proteins, and in particular to evaluating evolutionary stability or evolutionary rate of genome segments.

[0002] Many diseases are caused by pathogens such as viruses or other antigens that are made up of proteins. Such pathogens can be rendered ineffective by antibodies that bind to particular sites, referred to as epitopes, along a protein of the pathogen. Epitopes are short chains of amino acids. It is well understood that the immune system, when exposed to a pathogen, can learn to produce appropriate immune response. For example, B-cell epitopes are recognized by secreted antibodies or B-cell receptors, while T-cell epitopes are presented by a major histocompatibility complex (MHC) and recognized by T-cell receptors. This principle underlies much work in vaccine development, as well as development of antibody-based treatments for pathogenic diseases. However, if the pathogen mutates, a particular antibody may lose its effectiveness, leaving people vulnerable to infection by a new strain and / or rendering existing treatments ineffective. For example, mutation of RNA polymerase in an influenza virus can confer resistance to antiviral medications such as favipiravir and oseltamivir. Accordingly, understanding the evolutionary stability of pathogens can help with designing effective vaccinations or treatments.

[0003] In cancer vaccines, selection of appropriate epitopes is crucial for effective treatment. Neoepitopes are more likely to occur on unstable genomic regions, which can be evaluated from the perspective of evolution. Therefore, profiling the evolutionary stability of segments in tumor cells, such as oncogenes and tumor suppressor genes (TSGs) , can contribute to neo-epitope identification for the development of cancer vaccines.

[0004] In terms of diagnosis and prognosis, circulating tumor DNAs (ctDNAs) are considered to be promising biomarkers for cancer diagnosis. Evaluating evolutionary stability of ctDNA has the potential to enable inference of the orders of aberrations and stages of progression as an alternative method for tumor phylogenetics.

[0005] Evolutionary stability can be evaluated by estimating the evolutionary rate based on the rate of observed mutations. Conventional approaches are primarily based on developing a strain phylogeny. Examples of such methods include regression of root-to-tip (RTT) , parameterization of maximum likelihood (ML) trees, and Bayesian variations of these methods. A phylogenetic tree is typically constructed using large genome segments, such as the full sequence of a major surface protein (typically 50-200 amino acids) . The resulting evolutionary rate, therefore, applies to the overall speed of evolution of the protein.SUMMARY

[0006] It has been observed that the evolutionary rate is heterogenous across different locations within a genome. For instance, in the case of pathogens such as viruses, epitopes are often under higher evolutionary selection pressure (since avoiding antibodies confers a survival advantage on the pathogen) and tend to mutate at a higher rate than other regions. Evolutionary rates obtained using conventional techniques fail to account for such differences and may not accurately reflect the rate of evolution of an actual or potential epitope. In addition, better understanding of which portions of a protein are evolutionarily stable may facilitate selection of more stable candidate epitopes as targets for drug or vaccine development. Estimation of evolutionary stability of tumor genomes through analyzing the ctDNA may also provide an indicator of the speed of tumor growth.

[0007] Certain embodiments described herein relate to techniques for estimating evolutionary rate that can be applied to genome segments (or amino acid sequences) of any length, from a single codon (or amino acid position) to an entire protein or an entire genome. In some embodiments, computing an evolutionary rate for a genome segment can include: obtaining a reference sequence for a genome of a target at an initial time; obtaining strain-specific sequences for respective genomes of one or more strains of the target in a time period after the initial time; selecting a segment of the genome to analyze; computing a genetic divergence parameter for the segment based on genetic distances between the segment in the reference sequence and a corresponding segment in each of the strain-specific sequences; and computing, based on the genetic divergence parameter, an evolutionary rate for the segment.

[0008] The following detailed description, together with the accompanying drawings, will provide a better understanding of the nature and advantages of the claimed invention.BRIEF DESCRIPTION OF THE DRAWINGS

[0009] FIG. 1 shows a flow diagram of a process for determining an evolutionary rate for a genomic segment according to some embodiments.

[0010] FIG. 2 shows a table comparing results of a conventional process for determining an evolutionary rate and processes according to some embodiments.

[0011] FIG. 3 shows a graph of segment-specific evolutionary rate determined according to some embodiments for the spike protein of SARS-CoV-2 as a function of nucleotide position.

[0012] FIG. 4 shows a table listing the fastest-evolving and slowest-evolving (or most conservative) 90-nucleotide segments of the spike protein, as determined using a process according to an embodiment.

[0013] FIG. 5 is a table listing candidate epitopes for the spike protein of SARS-CoV-2, including evolutionary stability quantified according to an embodiment.DETAILED DESCRIPTION

[0014] The following description of exemplary embodiments of the invention is presented for the purpose of illustration and description. It is not intended to be exhaustive or to limit the claimed invention to the precise form described, and persons skilled in the art will appreciate that many modifications and variations are possible. The embodiments have been chosen and described in order to best explain the principles of the invention and its practical applications to thereby enable others skilled in the art to best make and use the invention in various embodiments and with various modifications as are suited to the particular use contemplated.

[0015] Certain embodiments described herein relate to techniques for estimating evolutionary rate that can be applied to genome segments (or amino acid sequences) of any length, including segments corresponding to B-cell epitopes or T-cell epitopes. Some embodiments use a model referred to herein as “genetic divergence through time, ” or “GDT. ”

[0016] As will become apparent, GDT models can avoid unrealistic assumptions that are incorporated in conventional techniques. Examples of such assumptions include a constant evolutionary rate across different genome segments or time periods, or independence of mutations over genomic positions. GDT models can be based on direct calculation of genetic distance and probability functions through time, with no need to fit an entire phylogeny of  viral strains; accordingly, GDT models can be more computationally efficient than conventional techniques.

[0017] According to some embodiments, GDT models can be based on a probabilistic estimation of substitution rate within a genomic segment. In examples herein, the amount of substitution is described by a Poisson random variable. Specifically, suppose that a genetic sequence (e.g., a DNA or RNA sequence of a pathogen) has a number of positions J. For a genomic segment (k) of consecutive nucleotides within the sequence, a substitution rate λkcan be defined as the substitutions per nucleotide site per unit time, shared among the sites in the segment. The segment can start at any location within the sequence and can have any length. Where the segment k has length 1 (a single nucleotide site) , λk is the site-specific substitution rate, and where sequence k has length J, λk is the genomic-level substitution rate. Thus, λk can be understood as an evolutionary rate parameter for segment k.

[0018] In some embodiments, GDT models provide a technique for evaluating λk. For example, if T is the time-to-substitution for any one site within segment k, it can be assumed that T follows an exponential distribution with event rate parameter λk. Accordingly, the probability of observing substitutions at time t at any site in k can be expressed as:

[0019] As long as the substitution rate λk is small (which is the case for known pathogens) , the expected number of substitutions αk in segment k can be estimated using a Poisson approximation of the binomial distribution: αk=lk·pk (t) ,           (2)

[0020] where lk is the number of sites in segment k. Thus, the expected value for the total number of substitutions Yk in segment k from an initial time to time t is a Poisson random variable of parameter αk, with expectation value: E (Yk|αk) =αk .           (3)

[0021] Parameter αk can be determined empirically for a particular pathogen by obtaining a baseline sample sequence at an origin time, obtaining a number of sample sequences at a later time t, and comparing the sequences to determine the number of substitutions in each sample. Defining a “genetic divergence” parameter as the sample mean of the number of substitutions in segment k between the initial time and time t,  can be used as the  maximum likelihood estimator of αk. Accordingly, from Eqs. (1) and (2) , the substitution rate λk can be computed as:

[0022] FIG. 1 shows a flow diagram of a process 100 for computing a substitution rate for a genomic segment according to some embodiments. Portions of process 100 can be implemented using computer systems. At block 102, sequence data is obtained for a “reference” genome representing a target (e.g., a virus or other pathogen or tumor cells) at an initial time. The initial time can be but need not be when the target was first observed or isolated, and any point in time for which a genomic sequence for the target is available can be selected as the initial time. Genetic sequence data can be obtained using appropriate techniques, including known techniques for DNA or RNA sequencing, accessing databases of known sequences, or other techniques.

[0023] At block 104, sequence data is obtained for the genome (s) of one or more strains of the target at some later time period, which can cover, e.g., weeks, months, or years after the initial time. There can be, but need not be, a gap between the initial time and the start of the later time period. As at block 102, any period in time for which genomic sequences are available can be selected, and any appropriate technique for obtaining genetic sequence data can be applied.

[0024] At block 106, a segment k within the genome can be selected. For example, a segment that codes the amino acid sequence for a known or candidate epitope for antibody binding can be selected. Such a segment might be, e.g., 12 to 18 nucleotides in length. More generally, the location of a segment within the genome and the length (lk) of the segment can be selected as desired.

[0025] At block 108, genetic divergence parameter can be computed for the segment. In various embodiments,  can be computed by determining, for each strain i′at the later time, a genetic distance Yk (r0, i′) between the segment in strain i′and the initial genome (r0) , then computing a sample mean of the genetic distance according to:

[0026] where nt denotes the number of observed strains at time period t. Various distance metrics can be used to define the genetic distance Yk (r0, i′) ; examples are described below.

[0027] At block 110, an evolutionary rate, e.g., substitution rate λk, can be computed using the genetic divergence parameter. For instance, Eq. (4) above can be applied to compute the evolutionary rate. In some embodiments, evolutionary stability can be computed in addition to or instead of evolutionary rate. Evolutionary stability can be quantified as a parameter having an inverse relationship to evolutionary rate (higher evolutionary rate corresponds to lower evolutionary stability and vice versa) . For example, evolutionary stability of a segment can be quantified as the negative logarithm of the evolutionary rate λk of that segment. (From Eq. (4) , it should be apparent that evolutionary stability can be computed without first computing evolutionary rate λk. )

[0028] For a given sequence and time period t, blocks 106, 108, and 110 of process 100 can be applied to any number of segments of any length, including overlapping segments. Thus, the evolutionary stability of a genome can be characterized based on the evolutionary rates of its segments. In addition or instead, process 100 can be repeated for the same segment with different choices of initial and later time, which may reveal changes in evolutionary rate as a function of time.

[0029] As noted above, at block 108, a genetic distance Yk (r0, i′) between the segment in a later strain i′and the initial genome (r0) is used to compute the genetic divergence parameter  Genetic distance can be defined using a variety of distance metrics. For example, in some embodiments, genetic distance is based on considering only substitutional mutations, in which one base is replaced by another. (This approach is referred to below as “GDT Type S, ” or “GDT-S. ” ) If denotes the genetic sequence of strain i, where j∈ [1, J] denotes a position without gaps in the genomic sequence and x denotes the nucleotide at position j in the genetic sequence of strain i, a genetic distance Yk (i, i′) between strains i and i’ can be defined as:

[0030] where I {·} is a function that returns 1 if the argument is true, 0 otherwise. Eq. (6) corresponds to the Hamming distance calculated on segment k of the genome.

[0031] As another example, in some embodiments, genetic distance is based on considering both substitutional mutations and gaps, i. e., insertions or deletions. (This approach is referred to below as GDT Type SG, ” or “GDT-SG. ” ) As used herein, a “gap length” refers to the number of consecutive positions where insertion or deletion occurs. For a gap length greater than 1, counting mutations using Hamming distance with consideration of each gap may over-estimate the impact of the gap. Accordingly, a penalty term can be introduced to account for length of the gap. For example, genetic distance Yk (i, i′) can be defined as:

[0032] where lg is the length of a continuous gap (i. e., the number of bases inserted or deleted) and gis an extended gap penalty, which can be selected as desired, e.g., g=0.1.

[0033] Eqs. (6) and (7) are examples of computing genetic distance. Other metrics can be substituted.

[0034] To illustrate performance of GDT processes, GDT-S and GDT-SG implementations of process 100 were applied to the SARS-CoV-2 whole genome and compared to results obtained using conventional Bayesian Evolutionary Analysis Sampling Trees (BEAST) software, which is based on phylogenetic trees. Both models were applied to a set of 360 SARS-CoV-2 sequences isolated in California during three time periods: period 1 from May to October of 2020; period 2 from March to August of 2021; and period 3 from November 2021 to April 2022. The reference (ancestral) strain for all time periods was the Wuhan-Hu-1 genome. (All sequences were downloaded from the Global Initiative on Sharing All Influenza Data (GISAID) database. )

[0035] FIG. 2 shows a table 200 comparing evolutionary rates for each time period obtained using a conventional BEAST process and GDT-S and GDT-SG processes according to some embodiments. As shown, in period 1 (during the early stage of SARS-CoV-2 circulation) , all three models yielded similar evolutionary rates. In periods 2 and 3, as more genetic variants appeared, GDT-SG began to provide somewhat higher estimates than BEAST and GDT-S, as a result of the difference in treatment of insertion and deletion events. Table 2 also shows the computing time used by each process (executed using the same hardware) ; it is noted that GDT processes are over 30,000 times faster than BEAST for a sample size of 120 genetic sequences. The computational efficiency of GDT processes opens  up prospects for efficiently and effectively analyzing thousands or even millions of sequences.

[0036] According to some embodiments, segment-based determination of evolutionary rate, as provided by process 100, makes it possible to profile genome-wide evolutionary stability. For instance, process 100 can be applied to multiple segments defined using a sliding window with a segment length lk and step length Q. Blocks 104, 106, and 108 of process 100 can be applied to segments of length lk starting at position 1, position 1+Q, position 1+2Q, and so on. By way of example, a sliding-window implementation of process 100 was applied to the SARS-CoV-2 virus, using a segment length lk=90 nucleotides and step length Q=30 nucleotides. (In this example, the segment length lk is longer than the step length Q, so that segments overlap. However, this is not required. ) To reflect the mutation activities in the subsequent waves of COVID-19 infection generated by the Alpha, Beta, Delta, and Omicron variants, sequence samples isolated in California during the period from March 2020 through February 2022 were used, for a total of 45, 413 sequence samples.

[0037] FIG. 3 shows a graph 300 of evolutionary rate for the spike protein of SARS-CoV-2 as a function of nucleotide position obtained using a sliding-window implementation of process 100. Graph 300 shows that the evolutionary rate is highest in the SD2 (shown at 302) , RBD (shown at 304) , and NTD (shown at 306) regions, while the S2 subunit (shown at 308) is relatively conserved. It is noted that the high-rate regions enclose some nucleotide locations where mutations in the Omicron BA. 5 variant appeared three months after the analysis time period ended. Further quantifying the results, FIG. 4 shows a table 400 listing the ten fastest-evolving and ten slowest-evolving (or most conservative) 90-nucleotide segments of the spike protein. It is noted that profiling of localized evolutionary rates within a genome is not supported by whole-genome based approaches such as BEAST.

[0038] According to some embodiments, segment-based determination of evolutionary rate, as provided by process 100, can be used to evaluate the stability of epitopes of arbitrary length as candidates for developing broadly reactive antibodies and / or vaccines. Evolutionary stability of a segment can be quantified as the negative logarithm of the evolutionary rate of that segment (e.g., λk as computed using process 100) . In various embodiments, candidate epitopes can include known epitopes recorded in existing databases,  such as the Immune Epitope Database and Analysis Resource (IEDB) , and / or computationally predicted epitopes.

[0039] By way of example, IEDB lists 1, 416 candidate epitopes that are linear peptides on SARS-CoV-2 with known T-cell responses in human hosts. Evolutionary rates for these peptides were calculated to identify candidate epitopes that exhibit both high stability (lowest evolutionary rates) and high response frequency (as determined from the known T-cell responses) . A total of 32 candidate epitopes were identified. FIG. 5 is a table 500 listing some of these epitopes, including the protein in which the candidate epitope is located; the amino acid sequence of the epitope; length of the amino acid sequence; start and end positions of the epitope within the protein; and the stability, quantified as the negative logarithm of the evolutionary rate λk as computed using an implementation of process 100. Such information can be used to select a target epitope or a combination of target epitopes for development of a vaccine or antibody therapy.

[0040] While the invention has been described with reference to specific embodiments, those skilled in the art will appreciate that variations and modifications are possible. For instance, the size (number of nucleotides) and location of each segment, and the number of segments, can be selected as desired. The analysis described herein can also be applied to amino acid sequences rather than nucleotide sequences. Techniques described herein can be used to determine evolutionary rate of any genome, protein, or segment thereof. As noted above, the ability to evaluate evolutionary rate or evolutionary stability for small segments of a genome or protein of a pathogen can facilitate identification of target epitopes for vaccines and / or treatments and can also facilitate prediction of future trends in rates of infection and / or disease.

[0041] All processes described herein are also illustrative and can be modified. Operations can be performed in a different order from that described, to the extent that logic permits; operations described above may be omitted or combined; and operations not expressly described above may be added.

[0042] Data analysis and computational operations of the kind described herein can be implemented in computer systems that may be of generally conventional design, such as a desktop computer, laptop computer, tablet computer, mobile device (e.g., smart phone) , or the like. Such systems may include one or more processors to execute program code (e.g., general-purpose microprocessors usable as a central processing unit (CPU) and / or  special-purpose processors such as graphics processors (GPUs) that may provide enhanced parallel-processing capability) ; memory and other storage devices to store program code and data; user input devices (e.g., keyboards, pointing devices such as a mouse or touchpad, microphones) ; user output devices (e.g., display devices, speakers, printers) ; combined input / output devices (e.g., touchscreen displays) ; signal input / output ports; network communication interfaces (e.g., wired network interfaces such as Ethernet interfaces and / or wireless network communication interfaces such as Wi-Fi) ; and so on.

[0043] Computer programs incorporating features of the present invention that can be implemented using program code may be encoded and stored on various computer readable storage media; suitable media include magnetic disk or tape, optical storage media such as compact disk (CD) or DVD (digital versatile disk) , flash memory, and other non-transitory media. (It is understood that “storage” of data is distinct from propagation of data using transitory media such as carrier waves. ) Computer readable media encoded with the program code may include an internal storage medium of a compatible electronic device and / or external storage media readable by the electronic device that can execute the code. In some instances, program code can be supplied to the electronic device via Internet download or other transmission paths.

[0044] While various circuits and components are described herein with reference to particular blocks, it is to be understood that these blocks are defined for convenience of description and are not intended to imply a particular physical arrangement of component parts. The blocks need not correspond to physically distinct components, and the same physical components can be used to implement aspects of multiple blocks. Components described as dedicated or fixed-function circuits can be configured to perform operations by providing a suitable arrangement of circuit components (e.g., logic gates, registers, switches, etc. ) ; automated design tools can be used to generate appropriate arrangements of circuit components implementing operations described herein. Components described as processors or microprocessors can be configured to perform operations described herein by providing suitable program code. Various blocks might or might not be reconfigurable depending on how the initial configuration is obtained. Embodiments of the present invention can be realized in a variety of apparatus including electronic devices implemented using a combination of circuitry and software.

[0045] Accordingly, although the invention has been described with respect to specific embodiments, it will be appreciated that the invention is intended to cover all modifications and equivalents within the scope of the following claims.

Claims

1.A computer-implemented method comprising:obtaining a reference sequence for a genome of a target at an initial time;obtaining strain-specific sequences for respective genomes of one or more strains of the target in a time period after the initial time;selecting a segment of the genome to analyze;computing a genetic divergence parameter for the segment based on genetic distances between the segment in the reference sequence and a corresponding segment in each of the strain-specific sequences; andcomputing, based on the genetic divergence parameter, an evolutionary rate or evolutionary stability parameter for the segment.2.The computer-implemented method of claim 1 wherein the target is a pathogen.3.The computer-implemented method of claim 2 wherein the segment is selected to correspond to a candidate epitope of the pathogen.4.The computer-implemented method of claim 1 further comprising:performing the acts of selecting a segment of the genome, computing the genetic divergence parameter, and computing a respective evolutionary rate or evolutionary stability parameter for each of a plurality of different segments of the genome.5.The computer-implemented method of claim 4 wherein the different segments are selected according to a sliding window.6.The computer-implemented method of claim 4 wherein the target is a pathogen and the different segments are selected to include candidate epitopes of the pathogen.7.The computer-implemented method of claim 6 further comprising:ranking the candidate epitopes based at least in part on comparing the respective evolutionary rates or evolutionary stability parameters for the different segments.8.The computer-implemented method of claim 4 further comprising:selecting a target epitope or a combination of target epitopes for development of a vaccine or antibody therapy based at least in part on the respective evolutionary rates or evolutionary stability parameters for the different segments.9.The computer-implemented method of claim 1 wherein computing the genetic divergence parameter includes computing an average of a distance metric between the segment in the reference sequence and the corresponding segment in each of the strain-specific sequences.10.The computer-implemented method of claim 9 wherein the distance metric is a Hamming distance based on a number of substitutions.11.The computer-implemented method of claim 9 wherein the distance metric includes a Hamming distance based on a number of substitutions modified by a penalty term that is computed based on a number of consecutive sites of nucleotide deletions or insertions.12.The computer-implemented method of claim 1 wherein the target is a genomic segment of tumor cells.13.The computer-implemented method of claim 12 wherein the segment is selected to correspond to a candidate neoepitype of the tumor cells.14.The computer-implemented method of claim 1 wherein the target is a ctDNA derived from tumor cells.15.The computer-implemented method of claim 14 wherein the segment is selected to estimate the stages of tumor growth or cancer progression.16.A system comprising:a memory; anda processor coupled to the memory and configured to perform a method according to any one of claims 1 to 15.17.A computer-readable storage medium having stored therein program code instructions that, when executed by a processor in a computer system, cause the processor to perform a method according to any one of claims 1 to 15.

Images