Method for reconstructing tumor clone heterogeneity and evolutionary relationship based on tumor DNA sequencing data

Abstract

The invention belongs to the field of oncology and computational biology, and discloses a method for reconstructing tumor clone heterogeneity and evolutionary relationship based on tumor DNA sequencing data. Comprising the steps of mutation multiplicity and cancer cell fraction calculation, subcloning clustering and recognition, clustering post-processing, tumor cloning evolutionary tree inference, subcloning proportion calculation and tumor sample evolutionary tree inference. According to the method, the mutation multiplicity of SNV is calculated through the Euclidean distance to accurately represent the CCF, and the subcloning clustering result is post-processed to evaluate the credibility and rationality of the subcloning clustering result. Secondly, quickly deducing a tumor cloning evolutionary tree by adopting a greedy strategy, and calculating a subcloning proportion according to an evolutionary relationship;

Description

Technical Field

[0001] The present invention belongs to the field of oncology and computational biology, and specifically relates to a method for reconstructing tumor clonal heterogeneity and evolutionary relationship based on tumor DNA sequencing data. Background Art

[0002] The mainstream theory holds that tumor stem cells have the ability to proliferate indefinitely and differentiate in multiple directions. After a series of somatic mutations, they form the initial tumor clone cell population. These clone populations continue to acquire mutations through natural selection or external pressure, thus forming related but genetically distinguishable offspring subpopulations within a single tumor. Differences in growth rate, invasive ability, and drug sensitivity among these subpopulations lead to intratumor heterogeneity, which in turn leads to tumor treatment failure, drug resistance, and death. Reconstructing tumor clonal heterogeneity is crucial for studying tumor development and designing precise and personalized treatment strategies.

[0003] Existing methods usually use probabilistic mixture models implemented by Markov chain Monte Carlo or Bayesian variational inference to perform unsupervised clustering of variant allele frequencies (VAF) or cancer cell fractions (CCF) of somatic single nucleotide variants (SNVs) to identify clonal structures within tumors. VAF reflects the abundance of mutations in a sample. However, due to the influence of copy number variations (CNVs) and tumor purity, VAF cannot accurately characterize the proportion of cloned cells. CCF, which combines CNV and tumor purity with VAF, can provide more accurate information. CCF represents the proportion of cancer cells carrying mutations at a specific site to all cancer cells. However, the calculation of CCF is highly correlated with the mutation multiplicity of SNVs in cells. Misjudgment of mutation multiplicity can lead to significant differences in CCF, resulting in confusion between clonal mutations and subclonal mutations, thereby inferring incorrect tumor clonal evolutionary relationships. Most existing methods do not consider the impact of CNV on mutation multiplicity, nor do they consider the uncertainty in the reconstruction process in each sample and the rationality of clonal evolutionary relationships. Summary of the invention

[0004] In order to overcome the shortcomings of the prior art, the present invention provides a method for reconstructing tumor clonal heterogeneity and evolutionary relationship based on tumor DNA sequencing data, accurately characterizing CCF by calculating the mutation multiplicity of SNVs through Euclidean distance, and post-processing the subclone clustering results to evaluate their credibility and rationality. Then, a greedy strategy is used to quickly infer the tumor clonal evolutionary tree, and the subclone ratio is calculated based on the evolutionary relationship.

[0005] The above object of the present invention is achieved through the following technical solution: A method for reconstructing tumor clonal heterogeneity and evolutionary relationship based on tumor DNA sequencing data, the steps are as follows:

[0006] 1. Calculate mutation multiplicity and cancer cell fraction: Calculate mutation multiplicity by the Euclidean distance between mutations, and then calculate the cancer cell fraction (CCF) based on VAF, CNV and tumor purity;

[0007] 2. Subclone clustering and identification: Use the Beta mixture model to cluster the CCF of the mutations calculated in step 1 and identify subclone clusters;

[0008] 3. Post-clustering processing: Evaluate the subclone clustering results in step 2, delete outliers, merge clusters with no statistical differences, and handle the situation where clone clusters are misidentified as subclone clusters due to over-analysis of clone clusters;

[0009] 4. Inferring the evolutionary tree of tumor clones: Rapidly inferring the evolutionary relationship between subclone clusters based on the infinite site hypothesis theory and the idea of ​​greedy strategy;

[0010] 5. Calculate the subclone ratio: Calculate the subclone ratio of each sample through evolutionary relationships and identify the evolutionary trends of subclone clusters in different samples;

[0011] 6. Infer the evolutionary tree of tumor samples: Find the most recent common ancestor and infer the evolutionary tree of tumor samples.

[0012] Furthermore, the specific calculation method of step 1 is: define P = {p1, p2, ..., p n} represents the CCF of SNV with mutation multiplicity of 1, and defines C = {c1, c2, ..., c m} represents the CCF calculated for different m of a SNV, and the mutation multiplicity of a SNV is calculated by formula (2):

[0013]

[0014] Where M is the major copy number of the SNV and m represents the mutation multiplicity.

[0015] Furthermore, step 2 is specifically as follows: first, pre-filter each sample to exclude SNVs whose CCF is significantly greater than 1.2 due to sequencing noise, and then model the probability that CCF c on N samples belongs to cluster k using formula (3):

[0016]

[0017] where π k is the prior probability that c belongs to k, and Beta is shown in formula (4):

[0018]

[0019] Finally, the Bayesian variational inference function of the R software package SciClone is used to model the Beta mixture model to approximate the posterior probability of formula (3), and then cluster assignment is performed based on the posterior probability.

[0020] Further, step 3 is specifically as follows: calculate the interquartile range of each subclone cluster, regard SNVs exceeding 1.5 times the interquartile range as outliers, remove outliers in each subclone cluster, and clusters with less than 5 mutations. The 95% confidence interval of each cluster in each sample is calculated by bootstrap; a two-sided Wilcoxon test is performed between the clusters of all samples to merge clusters without statistical differences, specifically, first extract the cluster combination of all samples. Secondly, the R package wilcox.test is used to evaluate the statistical significance between different clusters, iterate all cluster combinations, and if there is no statistical difference between the two clusters, one of the cluster IDs is changed to another cluster ID to achieve merging. Finally, the cluster with the highest CCF value is designated as a clonal cluster, and the other clusters are regarded as subclonal clusters. In order to ensure that each sample has a reasonable subclonal structure, it is evaluated whether there is a subclonal cluster with a CCF value exceeding the clonal cluster in each sample. If it exists, the subclonal cluster is merged with the clonal cluster.

[0021] Furthermore, the step 4 is specifically as follows: Assume that n subclone clusters are obtained after clustering post-processing, and there are s samples in total. Define CCF i,j represents the CCF value of the jth cluster in the ith sample. Construct an n×n adjacency matrix A to represent the inclusion relationship between clusters, with all initial values ​​set to 0. For each sample k∈{1,2,...,s}, if CCF k,i >CCF k,j Then A[i,j]+=1. If two clusters violate the crossover rule, that is, CCF k1,i >CCF k1,j , and CCF k2,j >CCF k2,i , then set A[i,j] and A[j,i] to 0. This means that there is no direct parent-child relationship between the two clusters. A directed acyclic ancestral graph is obtained. Then, based on the sum condition (the CCF of the parent cluster must be greater than or equal to the sum of the CCFs of all its direct child clusters) and the idea of ​​the greedy strategy, a clonal evolutionary tree is quickly generated from this directed graph. Define the root node r and find all possible child nodes S of r. r ={j|A[r,j]>0}, and sort the child nodes according to CCF. Then add child nodes one by one until the next added child node will cause a violation of the sum condition:

[0022]

[0023] where ∈ is a customizable CCF error buffer (default is 5% of the CCF of the parent cluster). The above process is then repeated for the child nodes of the root node until all subclone clusters are added to the tree.

[0024] Furthermore, step 5 is specifically as follows: for each cluster of each sample, the average CCF value is calculated, and then the ratio of all clusters in the sample relative to the clone cluster is calculated sample by sample. Combined with the clonal evolution relationship, the CCF value of the parent cluster without the influence of the sub-cluster is calculated, that is, the ratio of subclone clusters with nested relationships is calculated. Finally, the relative ratio of all subclone clusters of each sample can be obtained. Subclone clusters with gradually increasing subclone ratios in different samples may be in the expansion stage and have more growth advantages.

[0025] Further, step 6 is specifically as follows: first, the mutations are grouped in binary form according to the presence of the mutations in the sample, where the presence of the mutation in the sample is 1 and the absence is 0. The group with all 1s is the trunk mutation, and the other groups are branch mutations. In order to avoid too many groups when the sample size is too large, the present invention will delete the group with less than 5 mutations in the group. Then, the number of 1s in the group is calculated to divide the group into different levels, and the high-level group is stipulated to be the parent node of the low-level group. Starting from the leaf node, the respective parent nodes are searched from bottom to top. If the position of 1 of a low-level group is the same as the position of 1 of a high-level group, it means that the high-level group can be used as the parent node of this low-level group. Then, the top-down deletion results in a child node having multiple parent nodes. That is, the definition of a directed edge (V1→V2) indicates that V1 is the ancestor node of V2. If V1→V2→V3 and V1→V3, the edge from the ancestor node V1 to the descendant node V3 will be deleted. Finally, a sample tree with tumor samples as leaf nodes can be obtained.

[0026] The beneficial effects of the present invention compared with the prior art are:

[0027] (1) Correctly identifying mutation multiplicity is the first step in inferring the clonal structure of a tumor. Compared with the traditional method of using a fixed threshold to calculate the mutation multiplicity of SNVs, the present invention uses Euclidean distance to dynamically identify mutation multiplicity, which is more accurate and more interpretable.

[0028] (2) Through pre-filtering before subclone clustering and a series of processing after clustering, the present invention effectively removes the uncertainty in the reconstruction process and ensures the rationality of the clonal evolutionary relationship.

[0029] (3) In terms of inferring the tumor clone evolution tree, the present invention, based on the idea of ​​greedy strategy, can quickly infer the clone evolution tree within a few seconds without traversing all paths from the root node to the leaf node. In addition, when direct inference is not possible, the present invention can automatically delete the subclone cluster with the least number of mutations and then rerun the inference process.

[0030] (4) For multi-sample data, in addition to the clone evolutionary tree, the present invention also provides a sample evolutionary tree to reflect the overall similarity of the samples and the evolutionary path starting from their most recent common ancestor. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] The present invention will be further described below in conjunction with the accompanying drawings and specific embodiments.

[0032] Figure 1 It is the overall flow chart of the present invention;

[0033] Figure 2 Schematic diagram of the method for calculating mutation multiplicity;

[0034] Figure 3 Schematic diagram of the method for inferring clonal evolutionary trees;

[0035] Figure 4 Schematic diagram of the method for inferring the evolutionary tree of samples;

[0036] Figure 5 It is a result schematic diagram of the present invention. DETAILED DESCRIPTION

[0037] The present invention is described in detail below by specific examples, but the protection scope of the present invention is not limited. Unless otherwise specified, the experimental methods used in the present invention are all conventional methods, and the experimental equipment, materials, reagents, etc. used can be obtained from commercial channels.

[0038] Example 1

[0039] A method for reconstructing tumor clonal heterogeneity and evolutionary relationships based on tumor DNA sequencing data, the steps are as follows:

[0040] 1. Calculation of mutation multiplicity and cancer cell fraction: The current standard method for calculating cancer cell fraction (CCF) is to estimate VAF and tumor purity and local copy number. Define ρ to represent tumor purity, T cn represents the local integer copy number in cancer cells, N cn represents the local integer copy number in normal cells, m represents the mutation multiplicity, and the CCF of a SNV is calculated using formula (1):

[0041]

[0042] Among them, VAF, ρ and local copy number can all be quantitatively determined by sequencing data. However, m cannot be directly determined from these data. The present invention introduces a new method for calculating mutation multiplicity. According to the theory of infinite site hypothesis (each site in the genome will only mutate once at most during evolution), the mutation multiplicity of normal diploid cells and cells with an integer copy number of 1 is always equal to 1. Therefore, the CCF of all SNVs with a mutation multiplicity of 1 is first calculated. For m, there are multiple choices. Select m that can make the Euclidean distance between the CCF of the SNV and the CCF of the n closest SNVs the shortest. Specifically, define P = {p1, p2, ..., p n} represents the CCF of SNV with mutation multiplicity of 1. Define C = {c1, c2, ..., c m} represents the CCF calculated for different m of a SNV, and the mutation multiplicity of a SNV is calculated by formula (2):

[0043]

[0044] Where M is the major copy number of the SNV.

[0045] 2. Subclone clustering: First, each sample is pre-filtered to exclude SNVs whose CCF is significantly greater than 1.2 due to sequencing noise. Then, the present invention models the probability that CCF c on N samples belongs to cluster k using formula (3):

[0046]

[0047] where π k is the prior probability that c belongs to k, and Beta is shown in formula (4):

[0048]

[0049] Finally, the Bayesian variational inference function of the R software package SciClone is used to model the Beta mixture model to approximate the posterior probability of formula (3), and then cluster assignment is performed based on the posterior probability.

[0050] 3. Post-clustering processing: After the initial subclone structure identification, post-processing was performed. First, the interquartile range of each subclone cluster was calculated, and SNVs exceeding 1.5 times the interquartile range were considered outliers. Outliers in each subclone cluster and clusters with less than 5 mutations were removed. The 95% confidence interval for each cluster in each sample was calculated by bootstrap. Since SNVs are shared among all samples, it is impossible to change the cluster assignment of some SNVs in one sample without affecting their assignment in other samples. Therefore, a two-sided Wilcoxon test was performed between the clustering clusters of all samples to merge clusters that were not statistically different. Specifically, the clustering cluster combinations of all samples were first extracted. Secondly, the R package wilcox.test was used to evaluate the statistical significance between different clustering clusters, and all clustering cluster combinations were iterated. If there was no statistical difference between the two clusters, one of the cluster IDs was changed to another cluster ID to achieve the merger. Ultimately, the cluster with the highest CCF value was designated as a clonal cluster, and the other clusters were regarded as subclonal clusters. To ensure that each sample had a reasonable subclonal structure, each sample was evaluated for subclonal clusters with CCF values ​​exceeding those of the clonal clusters. If present, the subclonal clusters were merged with the clonal clusters.

[0051] 4. Inferring the tumor clone evolution tree: Assume that after clustering, n subclone clusters are obtained, with a total of s samples. Define CCF i,j represents the CCF value of the jth cluster in the ith sample. Construct an n×n adjacency matrix A to represent the inclusion relationship between clusters, with all initial values ​​set to 0. For each sample k∈{1,2,...,s}, if CCF k,i >CCF k,j Then A[i,j]+=1. If two clusters violate the crossover rule, that is, CCF k1,i >CCF k1,j , and CCF k2,j >CCF k2,i , then set A[i,j] and A[j,i] to 0. This means that there is no direct parent-child relationship between the two clusters. This step can obtain a directed acyclic ancestral graph. Then, based on the sum condition (the CCF of the parent cluster must be greater than or equal to the sum of the CCFs of all its direct child clusters) and the idea of ​​the greedy strategy, a clonal evolutionary tree is quickly generated from this directed graph. Define the root node r and find all possible child nodes S of r. r ={j|A[r,j]>0}, and sort the child nodes according to CCF. Then add child nodes one by one until the next added child node will cause a violation of the sum condition:

[0052]

[0053] where ∈ is a customizable CCF error buffer (default is 5% of the CCF of the parent cluster). The above process is then repeated for the child nodes of the root node until all subclone clusters are added to the tree.

[0054] 5. Calculate the subclone ratio: For each cluster of each sample, calculate its average CCF value, and then calculate the ratio of all clusters in the sample relative to the clone cluster for each sample. Combined with the clonal evolution relationship, calculate the CCF value of the parent cluster without the influence of the sub-cluster, that is, calculate the ratio of subclone clusters with nested relationships. Finally, the relative ratio of all subclone clusters of each sample can be obtained. Subclone clusters with gradually increasing subclone ratios in different samples may be in the expansion stage and have more growth advantages.

[0055] 6. Inferring the evolutionary tree of tumor samples: First, the mutations are grouped in binary form according to the presence of mutations in the sample, with 1 for the presence of mutations in the sample and 0 for the absence of mutations. The group with all 1s is the trunk mutation, and the other groups are branch mutations. In order to avoid too many groups when the sample size is too large, the present invention will delete the group with less than 5 mutations in the group. Then the number of 1s in the group is calculated to divide the group into different levels, and the high-level group is stipulated to be the parent node of the low-level group. Starting from the leaf node, the respective parent nodes are searched from bottom to top. If the position of 1 in a low-level group is the same as the position of 1 in a high-level group, it means that the high-level group can be used as the parent node of this low-level group. Then the top-down deletion results in a child node having multiple parent nodes. That is, the definition of a directed edge (V1→V2) indicates that V1 is the ancestor node of V2. If V1→V2→V3 and V1→V3, the edge from the ancestor node V1 to the descendant node V3 will be deleted. Finally, a sample tree with tumor samples as leaf nodes can be obtained.

[0056] The above-described embodiments are only preferred embodiments of the present invention, but not all feasible embodiments of the present invention. For those skilled in the art, any obvious changes made thereto without departing from the principles and spirit of the present invention should be considered to be included in the scope of protection of the claims of the present invention.

Claims

1. A method for reconstructing tumor clonal heterogeneity and evolutionary relationships based on tumor DNA sequencing data, characterized in that: Here are the steps: S1. Calculation of mutation multiplicity and cancer cell fraction: Mutation multiplicity was calculated by the Euclidean distance between mutations, and CCF was calculated based on VAF, CNV and tumor purity; S2. Subclone clustering and identification: Use the Beta mixture model to cluster the CCF of the mutations calculated in step S1 and identify subclone clusters; S3. Post-clustering processing: evaluate the subclone clustering results in step S2, delete outliers, and merge clusters with no statistical differences; Dealing with situations where over-profiling of clonal clusters results in clonal clusters being misidentified as subclonal clusters; S4. Inferring the evolutionary tree of tumor clones: Rapidly inferring the evolutionary relationship between subclone clusters based on the infinite site hypothesis theory and the idea of ​​greedy strategy; S5. Calculate subclone ratio: Calculate the subclone ratio of each sample through evolutionary relationships and identify the evolutionary trends of subclone clusters in different samples; S6. Infer the evolutionary tree of tumor samples: Find the most recent common ancestor and infer the evolutionary tree of tumor samples.

2. The method for reconstructing tumor clonal heterogeneity and evolutionary relationship based on tumor DNA sequencing data according to claim 1, characterized in that: The specific calculation method of step S1 is: define P = {p1, p2, ..., p n } represents the CCF of SNV with mutation multiplicity of 1, and defines C = {c1, c2, ..., c m } represents the CCF calculated for different m of a SNV, and the mutation multiplicity of a SNV is calculated by formula (1): Where M is the major copy number of the SNV and m represents the mutation multiplicity.

3. The method for reconstructing tumor clonal heterogeneity and evolutionary relationship based on tumor DNA sequencing data according to claim 1, characterized in that: The step S2 is specifically as follows: first, each sample is pre-filtered to exclude SNVs whose CCF is significantly greater than 1.2 due to sequencing noise. Then, the probability that CCF c belongs to cluster k on N samples is modeled using formula (2): where π k is the prior probability that c belongs to k, Beta is shown in formula (3): Finally, the Bayesian variational inference function of the R software package SciClone is used to build the Beta mixture model to approximate the posterior probability of formula (2), and then cluster assignment is performed based on the posterior probability.

4. The method for reconstructing tumor clonal heterogeneity and evolutionary relationship based on tumor DNA sequencing data according to claim 1, characterized in that: The step S3 is specifically as follows: the interquartile range of each subclone cluster is calculated, SNVs exceeding 1.5 times the interquartile range are regarded as outliers, outliers in each subclone cluster and clusters with less than 5 mutations are removed; the 95% confidence interval of each cluster in each sample is calculated by bootstrap; a two-sided Wilcoxon test is performed between the clusters of all samples to merge clusters without statistical differences, specifically, firstly, the cluster combinations of all samples are extracted; secondly, the statistical significance between different clusters is evaluated using the R software package wilcox.test, all cluster combinations are iterated, if there is no statistical difference between two clusters, one of the cluster IDs is changed to another cluster ID to achieve the merger, finally, the cluster with the highest CCF value is designated as the clonal cluster, and the other clusters are regarded as subclonal clusters, in order to ensure that each sample has a reasonable subclonal structure, it is evaluated whether there is a subclonal cluster with a CCF value exceeding the clonal cluster in each sample; if so, the subclonal cluster is merged with the clonal cluster.

5. The method for reconstructing tumor clonal heterogeneity and evolutionary relationship based on tumor DNA sequencing data according to claim 1, characterized in that: The step S4 is specifically as follows: assuming that n subclone clusters are obtained after clustering post-processing, and there are s samples in total; defining CCF i,j represents the CCF value of the jth cluster in the ith sample; construct an n×n adjacency matrix A to represent the inclusion relationship between clusters, with all initial values ​​being 0; for each sample k∈{1,2,...,s}, if CCF k,i >CCF k,j Then A[i,j]+=1; if two clusters violate the crossover rule, that is, CCF k1,i >CCF k1,j , and CCF k2,j >CCF k2,i , then set A[i,j] and A[j,i] to 0; indicating that the two clusters have no direct parent-child relationship; obtain a directed acyclic ancestral graph; then based on the sum condition that the CCF of the parent cluster must be greater than or equal to the sum of the CCFs of all its direct child clusters and the idea of ​​the greedy strategy, quickly generate a clonal evolutionary tree from this directed graph; define the root node r, and find all possible child nodes S of r r ={j|A[r,j]>0}, and sort the child nodes according to CCF; then add the child nodes one by one until the next added child node will cause a violation of the sum condition: where ∈ is a customizable CCF error buffer, which defaults to 5% of the CCF of the parent cluster; the above process is then repeated for the child nodes of the root node until all subclone clusters are added to the tree.

6. The method for reconstructing tumor clonal heterogeneity and evolutionary relationship based on tumor DNA sequencing data according to claim 1, characterized in that: The step S5 is specifically as follows: for each cluster of each sample, the average CCF value is calculated, and then the ratio of all clusters in the sample relative to the clone cluster is calculated sample by sample. Combined with the clonal evolution relationship, the CCF value of the parent cluster without the influence of the sub-cluster is calculated, that is, the ratio of sub-clone clusters with nested relationships is calculated. Finally, the relative ratio of all sub-clone clusters of each sample can be obtained. Sub-clone clusters with gradually increasing sub-clone ratios in different samples may be in the expansion stage and have more growth advantages.

7. The method for reconstructing tumor clonal heterogeneity and evolutionary relationship based on tumor DNA sequencing data according to claim 1, characterized in that: The step S6 is specifically as follows: first, the mutations are grouped in binary form according to the presence of the mutations in the sample, where the presence of the mutation in the sample is 1, and the absence of the mutation is 0; the group with all 1s is the trunk mutation, and the other groups are branch mutations; to avoid too many groups when the sample size is too large, the present invention deletes the group with less than 5 mutations in the group; then, the number of 1s in the group is calculated to divide the group into different levels, and it is stipulated that the high-level group is the parent node of the low-level group; Start from the leaf nodes and search for their respective parent nodes from bottom to top; If the position of 1 of a low-level group is the same as the position of 1 of a high-level group, it means that the high-level group can be the parent node of the low-level group; then top-down deletion results in a child node having multiple parent nodes; that is, the definition of a directed edge (V1→V2) indicates that V1 is the ancestor node of V2. If V1→V2→V3 and V1→V3, the edge from the ancestor node V1 to the descendant node V3 will be deleted; finally, a sample tree with tumor samples as leaf nodes can be obtained.

Images