CN110782947A - Identification of cancer drivers based on functional regions of protein sequences - Google Patents

Abstract

The identification of cancer drivers is a key challenge in explaining the mechanisms of cancer occurrence and enabling precision medicine. There are many ways to identify cancer drivers based on individual mutation sites or entire genes. But they ignore a lot of medium-sized functional elements. It is hypothesized that mutations that occur in different regions of the protein sequence have different effects on cancer progression. Here, we developed a novel functional driver region (frDriver) identification method based on Bayesian probability and multiple linear regression models to identify protein regions that can regulate gene expression levels and have high potential for functional impact. Combining gene expression data and somatic mutation data, combined with functional impact scores (SIFT, PROVEAN) as prior knowledge, we identified the most accurate cancer driver regions for predicting gene expression levels. We evaluate the performance of frDriver on the BRCA and GBM datasets of TCGA. The results show that frDriver identifies known cancer drivers and outperforms three other state-of-the-art methods (eDriver, ActiveDriver, and OncodriveCLUST).

Description

Identification of cancer drivers based on functional regions of protein sequences

technical field

The invention relates to data mining in bioinformatics, in particular to a mining of cancer bioinformatics data. Specifically, it relates to methods for identifying novel functionally-driven cancer region associations through Bayesian probability and multiple linear regression models based on protein sequence functional regions.

Background technique

The identification of cancer drivers is a key challenge in explaining the mechanisms of cancer occurrence and enabling precision medicine. There are many ways to identify cancer drivers based on individual mutation sites or entire genes. But they ignore a large number of mid-scale functional elements.

Mutations can affect cellular regulatory processes, signaling processes, and more. Due to the complexity of the mutations, different functional effects result. Compared with germline cells, somatic cells acquire mutations at a faster rate, tens to hundreds of times faster than primary cells. However, only a small subset of mutations provide a selective growth advantage and play a dominant role in the development of tumors, termed tumor drivers; most mutations are neutral and non-functional, termed cancer passengers.

Numerous computational methods and tools have been developed to differentiate cancer drivers from passengers. Algorithms for identifying driver genes are often based on the frequency of genetic mutations. They define river genes by identifying the relationship between genes and background mutations, and how to construct background mutation models has become the focus of research. Background mutation models are used to quantify the accumulation of passenger mutations, and estimates of background mutation rates that are too high or too low can lead to inaccurate results. Synonymous mutations are often used as a fixed background model. However, studies have shown that background mutations are not evenly distributed, so some researchers are considering combining other biometric techniques that may affect mutation frequency. For example, MutSigCV assesses background by considering patient-specific mutational profile signatures and gene-specific signatures (eg, expression, replication). But not all genetic mutation sites are functional, and different sites may have different effects, depending on where they are located. Therefore, some investigators have proposed hotspots, ie, amino acid sites with high mutation frequency, ie, three-base regions of genes that contain many mutations. Chang et al. proposed a binomial statistical model for the detection of repetitively mutated residues, replacing previous methods for detecting repetitively mutated residues. By aligning the morphological protein families of genes across 22 different cancer types, hotspots were identified to identify rare functional mutations.

However, these analyses ignore a large number of medium-sized functional elements, such as protein interfaces or protein subunits. Previous studies have shown that cancer formation is associated with enhanced local function. The study of cancer drivers on a regional basis provides better statistical performance than genes and hotspots. OncodriveCLUST constructs a background model by encoding silent mutations, identifying positions of multiple mutations above a background rate threshold as potentially meaningful cluster seeds, and then including other groups of mutations to form clusters. eDiver identifies regions of the protein rich in somatic missense mutations using the internal distribution of somatic missense mutations across functional regions of the protein and a binomial test to test for bias in the number of mutations observed in that region. DMCM estimates mutation density using data-adaptive bandwidth kernel density estimation (KDE) and identifies variable-length clusters in amino acid sequences. However, these methods only consider local sequence features and ignore gene regulatory effects associated with mutations. Projects such as TCGA provide multi-omics data such as somatic mutation and gene expression, and the combination of different omics data provides better performance than single data. Several studies have begun to systematically predict cancer drivers using gene expression data. For example, xSeq built a statistical model based on a hierarchical Bayesian approach and known interaction networks and applied it to the TCGA pan-cancer dataset to study the impact of somatic mutations on the expression profiles of 12 tumor types. Since most of the currently known networks are error-prone or incomplete, there are still some limitations that tend to bring noise to the results. On the other hand, this method has a certain bias towards undetected networks. At present, according to the degree of protein conservation and structural information, there are many methods to estimate the impact of amino acid changes on protein function by comparing homologous amino acid sequences, and to calculate their functional impact scores, such as SIFT and MutationAssessor. Integrating these functional influence factors into the driver identification process will provide more functional information.

In summary, existing methods do not adequately consider the role of a large number of medium-sized functional elements in identifying cancer drivers, and rarely compare individual amino acids to estimate the impact of amino acid sequence changes on protein function and calculate their functional impact scores. This orientation identifies cancer driver regions.

SUMMARY OF THE INVENTION

In view of the problems existing in the above methods and the role of regional protein sequences in identifying cancer driving factors, the present invention proposes a method for identifying new functionally driven cancer regions based on protein sequences to identify genes that can regulate gene expression levels and have high functional impact. potential protein regions. Protein domains and intrinsically disordered regions were used as candidate regions. Combining gene expression data and somatic mutation data, combined with functional impact scores (SIFT, PROVEAN) as prior knowledge, we identified the most accurate cancer driver regions for predicting gene expression levels. We evaluate the performance of frDriver on the BRCA and GBM datasets of TCGA. The described method steps include:

1. Gene expression for differential comparison

Use gene expression data and somatic mutation data in the TCGA database to screen out genes with key functions;

2. Data preprocessing stage

About 3000 genes that have been confirmed to have key functions, such as signal transduction genes, cell cycle control genes, and regulatory process-related genes, were screened; in order to make the data evenly distributed, the represented data were quantified and normalized. The logarithm of all gene expression values was taken and then all samples were quantified or normalized with a preprocessing language (Rlanguage); only samples co-existing in both datasets were kept by matching samples in mutation data and gene expression data.

3. Candidate driver regions

Two types of functional domains were selected as candidates: protein domains and intrinsic disorder domains (IDRs); a protein domain is a conserved portion of a specific protein sequence that is capable of functioning and surviving independently of the rest of the protein chain ; Intrinsically disordered domains are important functional regions capable of adapting to cellular conditions and exhibiting many different structures. Using these regions as candidate regions provides an adequate explanation for the effect of functional mutations.

4. Predicting functional characteristics of single amino acid variants

SIFT evaluates the degree of hazard of amino acid changes through homologous sequence alignment according to the positions of amino acid changes in the amino acid sequence and their physicochemical properties. Each position is scanned in an alignment constructed by homologous sequences, and each amino acid will be evaluated Probabilities of occurrence are recorded in a proportional probability matrix; PROVEAN evaluates the functional impact of variants based on sequence alignment scores. This method assesses changes in pairwise sequence similarity due to amino acid variation through sequence alignment scores. The delta alignment score is defined as follows:

Δ(P, v, S) = R(P', S) - A(P, S)

The average value within the cluster and within the set of supporting sequences is obtained, and the resulting unbiased average is used as the score of PROVEAN, as shown in the following formula:

5. Evaluate regional functional potential

To estimate the probability that a single amino acid variation affects the level of gene expression, the functional potential of the mutation is defined as:

Automatically learn feature weight parameters from data. Add a sigmoid function to the result to prevent infinite expansion of FP;

The probability that a mutation in region r affects the gene expression profile is defined as:

FP represents the probability that the jth mutation affects the expression of gene g, that is, the functional potential of the mutation calculated by the equation;

6. Build a regional regulation model

In N samples, based on candidate regions and gene expression values, we model as follows:

Calculate the correlation coefficient w, which can best reflect the relationship between candidate regions and gene expression levels. Calculate w by taking the maximum value of p(w|y).

p(w|y)=Kp(y|w)p(w)

Add L1 regularization parameter. When the RFP value increases, D tends to decrease, so when learning the correlation coefficient w, w tends to a non-zero value;

7. Identify tumor-specific driver regions

The global minima of the above optimization problem are solved by solving the local minima. There is an L1 regularization term, which is not differentiable at zero, so we choose the iterative descent algorithm. The algorithm updates one dimension at a time, and iteratively optimizes the weight β and the correlation coefficient w. Based on the obtained correlation coefficient w, a correlation matrix between specific genes and candidate regions is constructed to represent the degree of correlation between each pair of genes and regions. The score of the candidate region is calculated based on the number of related genes, i.e. the count of non-zero weight coefficients. If the score of a candidate region is greater than a given threshold, the region is predicted to be a cancer driver region;

Description of drawings

Figure 1: Construction of Multiple Linear Regression Models for Regions and Genes

Figure 2: Correlation of Regional Scores with Functional Scores

Figure 3: The results of the four methods on the BRCA and GBM datasets are overlaid in a Venn diagram

Figure 4: Top 25 GO biological process annotations in BRCA and top 30 KEGG pathway annotations

Figure 5: GBM top 25 GO biological process annotations and KEGG pathway annotation top 30

Figure 6: Boxplot of relative expression of MYH7B in normal and BRCA samples

Detailed ways

In order to make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below in conjunction with experiments. It should be understood that the specific embodiments described herein are only used to explain the present invention, but not to limit the present invention.

The hardware environment is mainly a PC host. Among them, the CPU of the PC host is Intel(R) Core(TM) i7-6700, 3.40GHz, the memory is 32GB RAM, and the 64-bit operating system is used. The software takes Windows 7 as the platform and is implemented in R language in the RStudio environment. The RStudio version is 1.1.142, and the R language version is 3.5.0.

The data used are two relatively well-established cancer data sets downloaded from TCGA, Breast Invasive Carcinoma (BRCA) and Glioblastoma Multiforme (GBM), which are relatively related to gene research. There are more data available to verify the results. Information on the number of samples is shown in the first column of the table below. It should be noted here that the two types of cancer data are calculated by performing all steps separately, but for simplicity, they are explained together here.

1. Gene expression for differential comparison

Use gene expression data and somatic mutation data in the TCGA database to screen out genes with key functions;

2. Data preprocessing stage

About 3000 genes that have been confirmed to have key functions, such as signal transduction genes, cell cycle control genes, and regulatory process-related genes, were screened; in order to make the data evenly distributed, the represented data were quantified and normalized. The logarithm of all gene expression values was taken and then all samples were quantified or normalized with a preprocessing language (Rlanguage); only samples co-existing in both datasets were kept by matching samples in mutation data and gene expression data.

3. Candidate driver regions

Two types of functional domains were selected as candidates: protein domains and intrinsic disorder domains (IDRs); a protein domain is a conserved portion of a specific protein sequence that is capable of functioning and surviving independently of the rest of the protein chain ; Intrinsically disordered domains are important functional regions capable of adapting to cellular conditions and exhibiting many different structures. Using these regions as candidate regions provides an adequate explanation for the effect of functional mutations.

The obtained results are shown in Table 1:

4. Predicting functional characteristics of single amino acid variants

SIFT evaluates the degree of hazard of amino acid changes through homologous sequence alignment according to the positions of amino acid changes in the amino acid sequence and their physicochemical properties. Each position is scanned in an alignment constructed by homologous sequences, and each amino acid will be evaluated Probabilities of occurrence are recorded in a proportional probability matrix; PROVEAN evaluates the functional impact of variants based on sequence alignment scores. This method assesses changes in pairwise sequence similarity due to amino acid variation through sequence alignment scores. The delta alignment score is defined as follows:

Δ(P, v, S) = R(P', S) - A(P, S)

The average value within the cluster and within the set of supporting sequences is obtained, and the resulting unbiased average is used as the score of PROVEAN, as shown in the following formula:

5. Evaluate regional functional potential

To estimate the probability that a single amino acid variation affects the level of gene expression, the functional potential of the mutation is defined as:

Automatically learn feature weight parameters from data. Add a sigmoid function to the result to prevent infinite expansion of FP;

The probability that a mutation in region r affects the gene expression profile is defined as:

FP represents the probability that the jth mutation affects the expression of gene g, that is, the functional potential of the mutation calculated by the equation;

6. Build a regional regulation model

In N samples, based on candidate regions and gene expression values, we model as follows:

Calculate the correlation coefficient w, which can best reflect the relationship between candidate regions and gene expression levels. Calculate w by taking the maximum value of p(w|y).

p(w|y)=Kp(y|w)p(w)

Add L1 regularization parameter. When the RFP value increases, D tends to decrease, so when learning the correlation coefficient w, w tends to a non-zero value;

7. Identify tumor-specific driver regions

The global minima of the above optimization problem are solved by solving the local minima. There is an L1 regularization term, which is not differentiable at zero, so we choose the iterative descent algorithm. The algorithm updates one dimension at a time, and iteratively optimizes the weight β and the correlation coefficient w. Based on the obtained correlation coefficient w, a correlation matrix between specific genes and candidate regions is constructed to represent the degree of correlation between each pair of genes and regions. The score of the candidate region is calculated based on the number of related genes, i.e. the count of non-zero weight coefficients. If the score of a candidate region is greater than a given threshold, the region is predicted to be a cancer driver region;

Table 1: Number of mutations, patients and candidate regions in BRCA and GBM datasets

1. Result analysis and verification

Table 1 shows the number of mutations, patients, and candidate regions included in the two cancer data finally obtained. Compared with the BRCA dataset (753 cases), the GBM dataset has fewer samples with only 150 cases, so there are fewer missense mutations and candidate regions obtained according to the above procedure, 7023 and 67, respectively.

Applying the frDriver method to the BRCA and GBM datasets, with both SIFT and PROVEAN functional scores as prior knowledge, predicts regions that can accurately predict gene expression levels as cancer drivers. Figure 2 shows the correspondence between the scores of the 532 candidate regions and the prior knowledge on the BRCA dataset. The regional score (blue bar) is the number of related genes, and the score of the corresponding region (i.e. the score located in the region) represents the score of the perceptual information located in the region (i.e. the score located in the region), so that they can be clearly observed. correlation between. We found that locations with high functional scores corresponded to regions with high scores; whereas those regions with low functional scores were associated with relatively fewer genes. This fully demonstrates that the functional priors we introduce bring important information to our cancer driver predictions. Among them, the p53 DNA-biding domain located on the TP53 gene has a higher function impact score, because the missense mutation located in this region itself is likely to affect the protein function, and the number of missense mutations in this region is relatively large, with a total of 141 mutations. Therefore, the candidate regions are more likely to affect gene expression levels.

Table 2: List of driver regions used for the BRCA dataset

We predicted 20 cancer drivers across 532 regions in the BRCA dataset. Table 2 lists the genes, region names, locations and brief descriptions of their functions for these 20 regions. The PF number is the protein functional region in the Pfam database, and the iur number is the internal disorder region (IDR). The location refers to the start and end location of the area. The regions we identified were distributed between 22 amino acids and 3912 amino acids. Among them, the domain length is between 122-911 amino acids, which is much smaller than the overall IDR. Most of the IDRs are several thousand in length, and only the regions located on the SEMFG2 and PHKA2 genes are less than 1000 in length. Among the genes in the 20 regions, 9 genes, TP53, PIK3CA, CBFB, MAP2K4, MAP3K1, FOXA1, PTEN, ERBB2, CTCF, are included in the Cancer Gene Census (CGC) list and are known as cancer driver genes, or genes There is a strong role in cancer, but extensive evidence is lacking. Only one of these 9 regions is an IDR, the others are domains, annotated by the Pfam database. This may be because the IDR is an indeterminate region in the three-dimensional structure of the protein with multiple as yet unproven functions. Thus, our method is able to identify known cancer driver genes and their specific internal regions. For example, the DNA-binding domain of the TP53 gene P53 located in cancer is a driver, and studies of p53 have shown that the inhibitory domain leads to activation of the cell migration by the chemokine CXCL13 and increased expression of the CXCR5 chemokine receptor gene.

Table 3: List of driver regions used for GBM dataset

We performed the same experiments on the GBM dataset and identified 17 regions out of 67 candidates and considered them to be cancer drivers. As with the BRCA data, we list 17 regions, including their gene names, region numbers, positions on the amino acid sequence, and brief descriptions of their functions (Table 3). The length of this region is a minimum of 76 amino acids and a maximum of 3912 amino acids, with a total of 14 domain names and 3 domain names. .By comparing with the CGC list, the genes TP53, PTEN, EGFR, IDH1, GRM3, PIK3CA, SMC1A in this region are known cancer driver genes or identified potential driver genes. The three regions PF00757, PF01030 and PF00069 are located on the same gene EGFR. GBM is associated with mutation-induced overexpression of EGFR, and EGFR-specific mutations are frequently observed. Somatic mutations, including EGFR, lead to persistent cell activation, which can lead to uncontrolled cell division. In addition, the other three methods used for comparison (ActiveDriver, OncodriveCLUST and eDriver) can also identify the EGFR gene. In addition, studies have shown that the genes we know, such as TCHH, FLG, KLK15, SEMA3C, GABRA6, are associated with cancer.

We compared the frDriver approach with three state-of-the-art approaches, eDriver, ActiveDriver, and OncodriveCLUST, from the perspective of the genes in which this region resides. Using the database DriverDB, a list of results from other methods on the BRCA and GBM datasets was obtained. On the BRCA dataset, eDriver identified 273 driver genes, ActiveDriver identified 220 driver genes, and OncodriveCLUST identified 26 driver genes. frDriver found 20 driver genes. Ten of the 20 genes were identified by at least one other method, making them more illustrative as driver genes (Fig. 3a). In the GBM dataset, eDriver, ActiveDriver and OncodriveCLUST identified 69, 60 and 8 genes, respectively. Using Venn diagram analysis, 7 of the 15 genes identified by frDriver were identified by at least another method (Fig. 3b). In addition, EGFR and TP53 genes were simultaneously identified by four methods and included in the CGC list. Thus, frDriver provides complementary results for other approaches to more accurately identify cancer drivers.

Table 4: Number of cancer samples and mutations

The CGC list contains 723 known cancer driver genes, or a high probability of being a cancer driver gene, but not enough evidence yet. We compared the recognition accuracy of the four methods. Table 4 shows the number of driver genes identified by the four methods, and the number of genes that overlapped with the CGC lists on the BRCA dataset and GBM dataset. On the BRCA dataset and GBM dataset, the accuracy of frDriver is 45% and 46.67%, which are the highest among the four methods. On the BRCA dataset, ActiveDriver and eDriver achieve 13.64% and 12.82% accuracy, respectively. The number of CGC genes they found was 30 and 35, respectively, more than the other two methods, but due to the large number of candidate rivergenes, activedriver and dedriver had high false positives. The same phenomenon also occurs on the GBM dataset. OncodriveCLUST has low accuracy in identifying a small number of driver genes. Therefore, frDriver is superior to other sintered manufacturing methods and the identified driver genes are more reliable.

Claims (8)

1. A protein sequence-based functional region driver identification, characterized in that the implementation steps are: (1) Differential comparison of gene expression, using gene expression data and somatic mutation data in the TCGA database to screen out genes with key functions; (2) Preprocess the sample data, which includes three steps of sample matching, missing value processing, and data standardization, to obtain samples that coexist in the GBM and BRCA data sets, and list them as sample names; (3) Candidate driver regions, two types of functional regions are selected as candidate regions: protein domains and intrinsically disordered domains (IDRs); (4) To predict the functional characteristics of single amino acid variation, the functional impact score of single amino acid variation can be evaluated and calculated according to sequence evolution conservation, protein structure characteristics and amino acid physicochemical properties. SIFT and PROVEAN methods were selected to predict functional effects of single amino acid mutations to measure functional potential; (5) Evaluate the functional potential of the region. The single amino acid variation in the candidate region has a functional effect on the protein due to the functional effect of the single point mutation in the region. The weight of the single amino acid can be calculated to solve the functional potential of the region; (6) Build a regional regulation model, establish a multiple linear regression model between gene expression levels and candidate regions, show that there is a linear relationship between them, and identify regions that can more accurately predict gene expression levels as cancer drivers; (7) To clarify the tumor-specific driver regions, and based on the obtained correlation coefficient w, a correlation matrix between specific genes and candidate regions was constructed to represent the degree of correlation between each pair of genes and regions. The score of the candidate region is calculated based on the number of related genes, i.e. the count of non-zero weight coefficients. If the score of a candidate region is greater than a given threshold, the region is predicted to be a cancer driver region. 2. The cancer-driven identification based on protein sequence functional regions according to claim 1, characterized in that the method is in the stage of differential comparison of gene expression: (1) Select gene expression data and somatic mutation data, and perform consistent processing of expression data and mutation data; (2) Screen about 3,000 genes that have been confirmed to have key functions, such as signal transduction genes, cell cycle control genes, and regulatory process-related genes. 3. The cancer driver identification based on protein sequence functional regions according to claim 1, characterized in that the method is in the data preprocessing stage: (1) Screening about 3,000 genes that have been confirmed to have key functions, such as signal transduction genes, cell cycle control genes, and regulatory process-related genes; (2) In order to make the data distribution uniform, quantify and normalize the represented data. Take the logarithm of all gene expression values, and then quantify or normalize all samples in R; (3) By matching the samples in the mutation data and gene expression data, only the samples that coexist in the two datasets are kept. 4. The protein sequence functional region-based cancer driver identification according to claim 1, characterized in that the method candidate driver region stage: (1) Two types of functional regions were selected as candidate regions: protein domains and intrinsically disordered domains (IDRs); (2) A protein domain is a conserved part of a specific protein sequence that can function and survive independently of the rest of the protein chain; (3) Intrinsically disordered domains are important functional regions that can adapt to cell conditions and exhibit many different structures. Using these regions as candidate regions provides an adequate explanation for the effect of functional mutations. 5. The identification of cancer drivers based on functional regions of protein sequences according to claim 1, characterized in that the method predicts the functional characteristic stage of single amino acid variation: (1) SIFT evaluates the degree of hazard of amino acid changes through homologous sequence alignment according to the positions of amino acid changes in the amino acid sequence and their physicochemical properties. The probability of occurrence of each amino acid is recorded in a proportional probability matrix; (2) PROVEAN evaluates the functional impact of variants based on sequence alignment scores. This method assesses changes in pairwise sequence similarity due to amino acid variation through sequence alignment scores. The delta alignment score is defined as follows: △(P, v, S)=R(P′,S)-A(P,S) The average value within the cluster and within the set of supporting sequences is obtained, and the resulting unbiased average is used as the score of PROVEAN, as shown in the following formula: . 6. The cancer-driven identification based on protein sequence functional regions according to claim 1, characterized in that the method evaluates the functional potential stage of the region: (1) In order to estimate the probability that a single amino acid variation affects the gene expression level, the functional potential of the mutation is defined as:

Automatically learn feature weight parameters from data. Add a sigmoid function to the result to prevent infinite expansion of FP; (2) The probability that a mutation in region r affects the gene expression profile is defined as follows: . FP represents the probability that the jth mutation affects the expression of gene g, that is, the functional potential of the mutation calculated by the equation. 7. The cancer driver recognition based on protein sequence functional regions according to claim 1, characterized in that the method constructs a region regulation model stage: (1) In N samples, according to candidate regions and gene expression values, we model as follows:

(2) Calculate the correlation coefficient w, which can best reflect the relationship between candidate regions and gene expression levels. Calculate w by taking the maximum value of p(w|y). p(w|y)=Kp(y|w)p(w)

Add L1 regularization parameter. When the RFP value increases, D tends to decrease, so when learning the correlation coefficient w, w tends to a non-zero value. 8 . The identification of cancer drivers based on protein sequence functional regions according to claim 1 , wherein the method defines the stage of tumor-specific driver regions, and solves the global minima of the optimization problem by solving local minima. 9 . There is an L1 regularization term, which is not differentiable at zero, so we choose iterative descent algorithm, a non-gradient optimization method [35]. The algorithm updates one dimension at a time, that is, iteratively optimizes the weights β and correlation coefficients w. Regularization parameters D0, D1, E are solved by cross-checking.

Images