Disease-specific quantitative trait site recognition method based on multi-omics integration

Abstract

The invention relates to a disease-specific quantitative trait locus identification method based on multi-omics integration. The method comprises the following steps: acquiring variation sites of whole genome sequencing data of a target object, and molecular phenotypes and molecular abundance of molecular phenotype data; determining an association significance probability value of an association pair formed by the variation point and the molecular phenotype based on the variation point and the molecular abundance, and screening a first association pair from the association pair based on the association significance probability value and condition analysis; determining a consistent second association pair in the normal association pair and the disease association pair, and determining a third association pair with a disease interaction effect in the second association pair; calculating a first effect estimation value and a second effect estimation value of each third association pair; and based on the first effect estimation value and the second effect estimation value of the third correlation pair, determining a target correlation pair related to the Parkinson's disease, and taking the target correlation pair as the identified quantitative trait site. By adopting the method, the Parkinson's specific pathogenic heritable variation can be accurately identified.

Description

Technical Field

[0001] This application relates to the field of disease-specific analysis technology, and in particular to a method for identifying disease-specific quantitative trait loci based on multi-omics integration. Background Technology

[0002] Parkinson's disease (PD) is the second most common neurodegenerative disease in the world after Alzheimer's disease. It has extremely high clinical heterogeneity, not only manifesting as typical motor symptoms such as bradykinesia and resting tremor, but also accompanied by a variety of non-motor symptoms such as olfactory dysfunction and cognitive decline. Moreover, the disease progression and pathological characteristics vary significantly among different patients.

[0003] Currently, research on Parkinson's disease typically involves a mixed analysis of genetic variations in normal subjects and those in diseased individuals. Some studies have also conducted quantitative trait locus analyses at the single level of genome and transcriptomics, proteomics, or metabolomics, identifying some genetic variations that influence molecular phenotypes. However, these existing techniques cannot reveal the differences in the regulatory effects of genetic variations between Parkinson's patients and healthy individuals, making it difficult to identify Parkinson's-specific pathogenic genetic variations. Summary of the Invention

[0004] Therefore, it is necessary to provide a disease-specific quantitative trait locus identification method based on multi-omics integration to address the above-mentioned technical problems. This method can accurately identify Parkinson's-specific pathogenic genetic variations.

[0005] A method for identifying disease-specific quantitative trait loci based on multi-omics integration, the method comprising:

[0006] S1. Obtain the variant sites in the whole genome sequencing data of each target object, the molecular phenotypes in the molecular phenotype data, and the molecular abundance of the molecular phenotypes;

[0007] S2. For each target object, based on the molecular abundance of the variant site and the molecular phenotype, determine the association significance probability value of the association pair formed by the variant site and the molecular phenotype, and screen the first association pair from the association pair based on the association significance probability value and conditional analysis; the first association pair includes normal association pairs of healthy objects and disease association pairs of objects with Parkinson's disease.

[0008] S3. Determine a second association pair that is consistent with the normal association pair and the diseased association pair, and determine a third association pair in the second association pair that has a disease interaction effect;

[0009] S4. Calculate the first effect estimate of each of the third association pairs relative to the healthy subjects and the second effect estimate relative to the subjects with Parkinson's disease;

[0010] S5. Based on the first effect estimate and the second effect estimate of each of the third association pairs, determine the target association pairs related to Parkinson's disease in the third association pairs, and use the target association pairs as the identified quantitative trait loci.

[0011] In this application, by identifying a second association pair that is consistent between normal association pairs and disease association pairs, and identifying a third association pair with a disease interaction effect in the second association pair, the first effect estimate of each third association pair relative to healthy subjects and the second effect estimate relative to subjects with Parkinson's disease are calculated. Based on the first and second effect estimates of each third association pair, the target association pair related to Parkinson's disease in the third association pair is identified, and the target association pair is used as the identified quantitative trait locus. This can reveal the difference in the regulatory effect of genetic variations in Parkinson's patients and healthy subjects, eliminate the interference of "false positives", ensure the reliability of the genetic basis, and thus accurately identify the pathogenic genetic variations related to Parkinson's, that is, accurately identify the quantitative trait locus related to Parkinson's.

[0012] In one embodiment, the process of obtaining variant sites in the whole genome sequencing data of the target object in step S1 includes:

[0013] Acquire human reference genome and whole genome sequencing data of various subjects;

[0014] Based on the human reference genome, a genome analysis toolkit was used to detect variants in each of the whole genome sequencing data to obtain the initial variant sites of each of the objects.

[0015] Using any of the initial mutation sites as target mutation sites, determine a first number of objects containing the target mutation sites. When the ratio of the first number to the first total number of objects is less than a first threshold, remove the target mutation sites from the initial mutation sites of the objects until the initial mutation sites have been traversed, and obtain the mutation sites of each object.

[0016] The number of sites in the union set formed by the variant sites of each of the objects is determined. When the ratio of the second number of variant sites of the object to the total number of sites is less than a second threshold, the object to be determined is removed from the objects to obtain the variant sites of the target object. The object to be determined is any one of the objects.

[0017] In this application, by taking any initial mutation site as the target mutation site, a first number of objects with target mutation sites is determined. When the ratio of the first number to the first total number of objects is less than a first threshold, the target mutation sites are removed from the initial mutation sites of the objects until the initial mutation sites are traversed, thus obtaining the mutation sites of each object. The number of sites in the union set formed by the mutation sites of each object is determined. When the ratio of the second number of mutation sites of the undetermined object to the total number of sites is less than a second threshold, the undetermined object is removed from the objects, thus obtaining the mutation sites of the target object. This can remove unqualified mutation sites and objects, thereby making the final determined target association pairs more accurate.

[0018] In one embodiment, the molecular phenotype data is proteomic data, the molecular phenotype is a protein, and the molecular abundance is the protein expression value of the protein. The process of obtaining the molecular phenotype and the molecular abundance of the molecular phenotype in step S1 includes:

[0019] The proteome data of the target object are detected by a protein detection instrument to obtain the initial protein of the target object and the protein expression value of the initial protein;

[0020] The protein expression value of the initial protein is multiplied by the corresponding correction factor to obtain the corrected expression value. The corrected expression value is then logarithmically transformed to obtain the converted expression value of the initial protein. Using any one of the initial proteins as a candidate protein, a third number of target objects containing the candidate protein is determined. When the ratio of this third number to the second total number of target objects is less than a third threshold, the candidate protein is removed from the initial proteins of the target objects. This process continues until all initial proteins have been traversed, resulting in the proteins of each target object. The proteins include a first protein with a converted expression value and a second protein without a converted expression value.

[0021] Based on the conversion expression value of the first protein, the conversion expression value of the second protein is completed using a random tail distribution imputation method to obtain the conversion expression value of each protein.

[0022] Batch correction and anti-logarithmic transformation are performed on the converted expression value of the protein in the target object to obtain the protein expression value of the protein in the target object.

[0023] In this application, the converted expression value of the initial protein is obtained by logarithmic transformation of the corrected expression value. This can reduce data skewness, stabilize variance, and make it suitable for subsequent statistical analysis.

[0024] In one embodiment, the method further includes:

[0025] Principal component analysis was performed based on the protein expression values ​​of the proteins of each target object to obtain the principal components of each target object and the scores and eigenvalues ​​of each principal component;

[0026] Based on the scores and eigenvalues ​​of each principal component of the target object, a first Mahalanobis distance of the target object is calculated; a fourth threshold associated with the number of principal components and a preset confidence level is determined from a preset chi-square distribution table;

[0027] Target objects whose first Mahalanobis distance is greater than the fourth threshold are removed, and the association significance probability value is determined by the mutation sites and protein expression values ​​of the unremoved target objects.

[0028] In this application, principal component analysis is performed based on the protein expression values ​​of each target object to obtain the principal components, scores, and eigenvalues ​​of each principal component. Based on the scores and eigenvalues ​​of each principal component of the target object, the first Mahalanobis distance of the target object is calculated. A fourth threshold, which is associated with the number of principal components and a preset confidence level, is determined from a preset chi-square distribution table. Target objects with a first Mahalanobis distance greater than the fourth threshold are removed. This can further remove outliers and ensure the reliability of subsequent processing results.

[0029] In one embodiment, the molecular phenotype data is metabolomics data, the molecular phenotype is a target metabolite, and the molecular abundance is the concentration of the target metabolite. The process of obtaining the molecular phenotype and the molecular abundance of the molecular phenotype in step S1 includes:

[0030] Samples were taken and mixed from bodily fluid samples of each target object, and the mixed samples were divided into multiple quality control samples. The mass spectrometry response value and concentration of the first candidate metabolite in the multiple quality control samples were detected in batches using a metabolite detection instrument.

[0031] Based on the mass spectrometry response values ​​of the first candidate metabolites in each of the quality control samples, the coefficient of variation of each first candidate metabolite is calculated.

[0032] First candidate metabolites with a coefficient of variation greater than the fifth threshold are removed to obtain second candidate metabolites in each quality control sample. Then, mass spectrometry response values ​​are filled in for the second candidate metabolites in each quality control sample that are missing mass spectrometry response values ​​to obtain the mass spectrometry response value of each second candidate metabolite in each quality control sample.

[0033] Based on the detection order of each quality control sample and the mass spectrometry response value of each second candidate metabolite, the Spearman rank correlation coefficient of each second candidate metabolite is calculated. Second candidate metabolites with Spearman rank correlation coefficients greater than the sixth threshold are eliminated to obtain the target metabolite, and the concentration of the target metabolite of each target object is determined.

[0034] In this application, by eliminating first candidate metabolites with a coefficient of variation greater than the fifth threshold, second candidate metabolites are obtained in each quality control sample. Then, a mass spectrometry response value completion operation is performed on the missing mass spectrometry response values ​​of the second candidate metabolites in each quality control sample to obtain the mass spectrometry response value of each second candidate metabolite in each quality control sample. Based on the detection order of each quality control sample and the mass spectrometry response value of each second candidate metabolite, the Spearman rank correlation coefficient of each second candidate metabolite is calculated. Second candidate metabolites with a Spearman rank correlation coefficient greater than the sixth threshold are eliminated to obtain the target metabolite. This method can eliminate unstable metabolites and metabolites with low accuracy in mass spectrometry response values, thus avoiding interference with subsequent analysis processes.

[0035] In one embodiment, step S2, which involves selecting a first association pair from the association pairs based on the association significance probability value and conditional analysis, includes:

[0036] Principal component analysis was performed based on the molecular abundance of the molecular phenotypes of each target object to obtain the variance contribution rate of each principal component.

[0037] Determine the minimum number of principal components required for the cumulative variance contribution rate to reach the seventh threshold, and determine the screening threshold based on the minimum number of principal components;

[0038] Remove association pairs whose association significance probability value is greater than the screening threshold to obtain initial candidate association pairs;

[0039] From the variant sites of each of the initial candidate association pairs, identify variant sites with abnormal distribution, and remove the initial candidate association pairs with variant sites with abnormal distribution to obtain the initial association pairs;

[0040] Conditional analysis was performed on the initial association pairs, and the initial association pairs with independent genetic effects were selected as the first association pairs.

[0041] In this application, principal component analysis is performed based on the molecular abundance of the molecular phenotypes of each target object to obtain the variance contribution rate of each principal component. The minimum number of principal components required for the cumulative variance contribution rate to reach the seventh threshold is determined, and a screening threshold is determined based on the minimum number of principal components. Association pairs with association significance probability values ​​greater than the screening threshold are eliminated to obtain initial candidate association pairs. From the variant sites of each variant site in the initial candidate association pairs, variant sites with abnormal distribution are identified, and initial candidate association pairs with variant sites with abnormal distribution are eliminated to obtain initial association pairs. Conditional analysis is performed on the initial association pairs to obtain the first association pair. This can avoid unreliable association pairs from affecting the subsequent determination of target association pairs, thereby improving the accuracy and reliability of target association pairs.

[0042] In one embodiment, the process of determining the third association pair includes:

[0043] For each of the second association pairs, a null model is constructed using the molecular abundance of the molecular phenotype in the second association pair as the dependent variable and the attribute information of the target object to which the second association pair belongs as the independent variable; a main effect model is constructed using the molecular abundance of the molecular phenotype as the dependent variable and the attribute information and the variant site as the independent variables; and an interaction effect model is constructed using the molecular abundance of the molecular phenotype as the dependent variable and the attribute information, the variant site, and the interaction term between the variant site and the disease state of the target object to which it belongs as the independent variables.

[0044] Based on the null model, the main effect model, and the interaction effect model of each second association pair, a likelihood ratio test is performed to obtain the interaction P-value for each second association pair; the interaction P-value represents the degree of improvement in the fitting effect of the interaction effect model after the introduction of the interaction term.

[0045] Based on the interaction P-value of each of the second association pairs, the false discovery rate of each of the second association pairs is calculated using the qvalue method, and the second association pairs with a false discovery rate less than the eighth threshold are identified as the third association pairs.

[0046] In this application, for each second association pair, a null model is constructed with the molecular abundance of the molecular phenotype in the second association pair as the dependent variable and the attribute information of the target object to which the second association pair belongs as the independent variable; a main effect model is constructed with the molecular abundance of the molecular phenotype as the dependent variable and the attribute information and variant sites as independent variables; and an interaction effect model is constructed with the molecular abundance of the molecular phenotype as the dependent variable and the attribute information, variant sites, and the interaction term between the variant sites and the disease state of the target object to which they belong as independent variables. Based on the null model, main effect model, and interaction effect model of each second association pair, a likelihood ratio test is performed to obtain the interaction P-value of each second association pair. Based on the interaction P-value of each second association pair, the qvalue method is used to calculate the false discovery rate of each second association pair. In this way, third association pairs with a false discovery rate less than the eighth threshold can be obtained.

[0047] In one embodiment, step S4 includes:

[0048] Based on the maximum likelihood estimation results of the main effect model and interaction effect model of each of the third association pairs, the main effect coefficient of the variant site and the interaction term coefficient between the variant site and the disease state in each of the third association pairs are obtained.

[0049] The main effect coefficient is used as the first effect estimate, and the sum of the main effect coefficient and the interaction term coefficient is used as the second effect estimate.

[0050] In this application, the main effect coefficients of the variant sites and the interaction term coefficients between the variant sites and the disease state in each third association pair are obtained by using the maximum likelihood estimation results based on the main effect model and the interaction effect model of each third association pair. In this way, the main effect coefficients can be used as the first effect estimate, and the sum of the main effect coefficients and the interaction term coefficients can be used as the second effect estimate.

[0051] In one embodiment, the process of determining the target association pair in step S5 includes:

[0052] If the absolute value of the second effect estimate of the third association pair is greater than the absolute value of the first effect estimate, and the absolute value of the first effect estimate is less than the no-effect threshold, then the third association pair is determined to be a target association pair related to Parkinson's disease.

[0053] In this application, by comparing the absolute values ​​of the second effect estimate and the first effect estimate, and by comparing the absolute value of the first effect estimate with the no-effect threshold, target association pairs related to Parkinson's disease can be accurately screened.

[0054] In one embodiment, the method further includes:

[0055] Based on the physical location information of the mutation sites in the chromosome in each of the third association pairs, each of the third association pairs is divided into multiple sets;

[0056] The target set with consistent genetic regulatory effects was screened from the set using colocalization analysis.

[0057] Based on the molecular phenotypes in the target set and the Parkinson's disease status of the target objects to which each of the third association pairs in the target set belongs, the causal effects between molecular phenotypes and the causal effects of molecular phenotypes on disease phenotypes are evaluated by Mendelian randomization analysis, and multiple causal regulatory pathways are obtained.

[0058] Based on the causal regulation path, a regulation network is constructed, and key nodes in the regulation network are selected according to the network centrality index or importance score of each node in the regulation network.

[0059] Functional annotations were performed on the key nodes, and pathway enrichment analysis was conducted on the molecular phenotypes in the key nodes.

[0060] In this application, each third association pair is divided into multiple sets based on the physical location of the variant sites in the chromosome. Colocalization analysis is used to screen target sets with consistent genetic regulatory effects. Based on the molecular abundance of molecular phenotypes in the target sets and the Parkinson's disease status of the target subjects belonging to each third association pair, Mendelian randomization analysis is used to assess the causal effects between molecular phenotypes and the causal effects of molecular phenotypes on disease phenotypes, resulting in multiple causal regulatory pathways. Based on these causal regulatory pathways, a regulatory network is constructed, and key nodes in the regulatory network are screened based on network centrality indicators or importance scores. Functional annotation is performed on the key nodes, and pathway enrichment analysis is conducted on the molecular phenotypes within the key nodes. This allows for the extraction of a causal disease network diagram from massive amounts of genomic data, precisely locating the key molecules and pathways that play a decisive role in the network, thereby revealing the deep molecular mechanisms underlying the development of Parkinson's disease. Attached Figure Description

[0061] Figure 1 This is a diagram illustrating the application environment of a disease-specific quantitative trait locus identification method based on multi-omics integration in one embodiment.

[0062] Figure 2 This is a flowchart illustrating a disease-specific quantitative trait locus identification method based on multi-omics integration in one embodiment.

[0063] Figure 3 This is a technical roadmap for one embodiment;

[0064] Figure 4 This is a schematic diagram of the analysis flow of a disease-specific quantitative trait locus identification method based on multi-omics integration in another embodiment;

[0065] Figure 5 Here is a flowchart of disease-specific analysis in one embodiment;

[0066] Figure 6 This is a schematic diagram illustrating the integration, colocalization analysis, pathway enrichment analysis, and functional annotation of Parkinson's-related target association pairs in one embodiment.

[0067] Figure 7 This is an internal structural diagram of a computer device in one embodiment. Detailed Implementation

[0068] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0069] The disease-specific quantitative trait locus identification method based on multi-omics integration provided in this application can be applied to, for example... Figure 1 In the application environment shown, terminal 102 interacts with server 104 via a wired / wireless channel. A data storage system can store the data that server 104 needs to process. The server acquires variant sites from the whole-genome sequencing data of each target object, molecular phenotypes from the molecular phenotype data, and the molecular abundance of the molecular phenotypes. For each target object, based on the molecular abundance of variant sites and molecular phenotypes, the server determines the association significance probability value of association pairs composed of variant sites and molecular phenotypes, and selects a first association pair from the association pairs based on the association significance probability value and conditional analysis; the first association pair includes normal association pairs for healthy objects and disease association pairs for objects with Parkinson's disease. The server identifies second association pairs that are consistent between the normal and disease association pairs, and identifies third association pairs with disease interaction effects within the second association pairs. The server calculates the first effect estimate of each third association pair relative to healthy objects and the second effect estimate relative to objects with Parkinson's disease. Based on the first and second effect estimates of each third association pair, the server determines the target association pairs related to Parkinson's disease within the third association pairs; the target association pairs are quantitative trait loci. The terminal 102 can be, but is not limited to, various personal computers, laptops, smartphones, tablets, IoT devices, etc. The server 104 can be a single server, a server cluster consisting of multiple servers, or a cloud computing center consisting of multiple servers.

[0070] In one embodiment, such as Figure 2 As shown, a method for identifying disease-specific quantitative trait loci based on multi-omics integration is provided, which can be applied to... Figure 1 Taking server 104 as an example, the following steps are included:

[0071] S1. Obtain the variant sites in the whole genome sequencing data of each target object, the molecular phenotype in the molecular phenotype data, and the molecular abundance of the molecular phenotype;

[0072] Whole-genome sequencing data is extracted from blood samples of the target subjects, while molecular phenotypic data is extracted from bodily fluid samples of the target subjects. Methods for extracting whole-genome sequencing data include, but are not limited to, kit-based methods, magnetic bead methods, and traditional phenol-chloroform methods. Bodily fluid samples include, but are not limited to, blood samples, urine samples, and saliva samples. Whole-genome sequencing data can be sequencing data with a sequencing depth of 30X.

[0073] Molecular phenotypic data includes proteomic and / or metabolomic data. Proteomic data is a digital record of the identity (qualitative), quantity (quantitative), and state (modification, interactions) of proteins. Metabolomic data is a record of the chemical structure identification and concentration determination of metabolites (such as amino acids, lipids, carbohydrates, organic acids, nucleotides, etc.).

[0074] When molecular phenotypic data is proteomic data, the molecular phenotype is protein, and molecular abundance is protein expression value. When molecular phenotypic data is metabolomic data, the molecular phenotype is target metabolite, and molecular abundance is target metabolite concentration.

[0075] Variants can be detected using GATK (Genome Analysis Toolkit). Molecular phenotypes and their molecular abundances can be determined using detection instruments.

[0076] S2. For each target object, based on the molecular abundance of the variant site and the molecular phenotype, determine the association significance probability value of the association pair composed of the variant site and the molecular phenotype, and screen the first association pair from the association pairs based on the association significance probability value and conditional analysis; the first association pair includes normal association pairs of healthy objects and disease association pairs of objects with Parkinson's disease.

[0077] The process involves inputting the molecular abundance of each variant site and molecular phenotype of the target object into association significance probability detection software to obtain the association significance probability value of each association pair. Each variant site of the target object and each molecular phenotype of the target object constitute an association pair. Association significance probability detection software includes, but is not limited to, tensorQTL, MatrixEQTL, and SKAT (Sequence Kernel Association Test).

[0078] Variants are categorized into common and rare variants. Common variants are those with a minimum allele frequency greater than or equal to 1%, while rare variants are those with a minimum allele frequency less than 1%. Since association pairs consist of variants and molecular phenotypes, they can be further classified into pairs corresponding to common variants and pairs corresponding to rare variants. Similarly, when the molecular phenotype data is proteomic or metabolomic, association pairs can be classified into pairs corresponding to proteins and pairs corresponding to target metabolites. The association significance probability of association pairs consisting of common variants and molecular phenotypes can be obtained using linear models such as tensorQTL and MatrixEQTL. The association significance probability of association pairs consisting of rare variants and the molecular abundance of the molecular phenotype can be obtained using methods such as SKAT and tensorQTL.

[0079] S3. Identify the second association pair that is consistent between the normal association pair and the disease association pair, and identify the third association pair that has a disease interaction effect in the second association pair;

[0080] Consistency refers to the consistency of variant sites in normal association pairs and disease association pairs (or the distance between variant sites is less than a preset difference) and the consistency of molecular phenotypes.

[0081] Disease interaction refers to the combined effect of the same gene on the abundance of molecular phenotypes under different disease states.

[0082] S4. Calculate the first effect estimate of each third association pair relative to healthy subjects and the second effect estimate relative to subjects with Parkinson's disease;

[0083] The first effect estimate is the independent effect of the variant sites in the third association pair on the molecular abundance of the molecular phenotype in the healthy population. The second effect estimate is the independent effect of the variant sites in the third association pair on the molecular abundance of the molecular phenotype in the Parkinson's disease population.

[0084] S5. Based on the first and second effect estimates of each third association pair, identify the target association pairs in the third association pairs that are associated with Parkinson's disease, and use the target association pairs as the identified quantitative trait loci.

[0085] Specifically, based on the magnitudes of the first and second effect estimates of the third association pair, the target association pairs related to Parkinson's disease in the third association pair are determined.

[0086] In the aforementioned multi-omics integrated method for identifying disease-specific quantitative trait loci, a second association pair consistent between normal and disease-related association pairs is identified, and a third association pair with disease interaction effects is identified within the second association pair. The first effect estimate of each third association pair relative to healthy subjects and the second effect estimate relative to subjects with Parkinson's disease are calculated. Based on the first and second effect estimates of each third association pair, target association pairs related to Parkinson's disease are identified within the third association pairs. These target association pairs are used as the identified quantitative trait loci. This approach reveals the differences in the regulatory effects of genetic variations between Parkinson's patients and healthy subjects, eliminates false positive interference, ensures the reliability of the genetic basis, and thus accurately identifies pathogenic genetic variations related to Parkinson's disease, i.e., accurately identifies quantitative trait loci related to Parkinson's disease.

[0087] In one embodiment, the process of obtaining variant sites in the whole-genome sequencing data of the target object in step S1 includes:

[0088] Acquire human reference genome and whole genome sequencing data of various subjects;

[0089] Based on the human reference genome, a genome analysis toolkit was used to detect variants in each whole genome sequencing data to obtain the initial variant sites for each object.

[0090] Using any initial mutation site as the target mutation site, determine the first number of objects with the target mutation site. When the ratio of the first number to the first total number of objects is less than a first threshold, remove the target mutation site from the initial mutation sites of the objects until the initial mutation sites have been traversed, and obtain the mutation sites of each object.

[0091] The number of sites in the union set formed by the variant sites of each of the objects is determined. When the ratio of the second number of variant sites of the object to the total number of sites is less than a second threshold, the object to be determined is removed from the objects to obtain the variant sites of the target object. The object to be determined is any one of the objects.

[0092] The human reference genome is a digitized, haploid, linear collection of DNA (Deoxyribonucleic acid) sequences, representing the typical structure and consensus sequence of the human genome. Human reference genomes include, but are not limited to, the hg19 reference genome, the hg38 reference genome, and the T2T reference genome.

[0093] By inputting the human reference genome and the whole genome sequencing data of each subject into the genome analysis toolkit, the initial variant sites of each subject can be output, which are also the initial variant sites of each whole genome sequencing data.

[0094] There are no kinship relationships among the subjects. The subjects include healthy individuals and those with Parkinson's disease, and both the healthy and those with Parkinson's disease are within a predetermined age range. The onset age of Parkinson's disease in the affected individuals is later than the predetermined age, while the healthy individuals are all individuals without neurological disorders.

[0095] Initial sites in each object that meet the condition "the ratio of the first quantity to the first total quantity of objects is less than the first threshold" will be removed. The removal of initial sites can be implemented using the plink software. The first and second thresholds are preset data.

[0096] Each mutation site corresponds to a mutation location and a mutation type. Different mutation sites have different mutation locations and / or different mutation types. The mutation type refers to the substitution between bases. For example, adenine to guanine is one mutation type, cytosine to guanine is another, guanine to cytosine is yet another, and thymine to cytosine is also a mutation type. Furthermore, the mutation location refers to the chromosome in which the mutation site is located and its position within the chromosome.

[0097] After removing the undetermined objects from the pool of objects, the remaining objects are called target objects. Since the mutation sites of each object have already been obtained, the mutation sites of the target objects can be obtained directly. The removal of undetermined objects that meet the condition "the ratio of the second number of mutation sites to the total number of sites is less than the second threshold" can be achieved using the plink software.

[0098] Furthermore, the human reference genome can be any one of the hg19, hg38, or T2T reference genomes. Based on the human reference genome, a genome analysis toolkit is used to detect variants in each whole-genome sequencing data to obtain the initial variant sites for each object corresponding to the human reference genome. Further, based on the initial variant sites corresponding to the human reference genome, association pairs corresponding to the human reference genome are obtained.

[0099] Furthermore, the method also includes a quality control process, which includes: using bcftools to split multi-allelic genes and delete ambiguous alleles, and using bcftools to identify initial variant sites that meet the criteria. The mutation site was set to deletion, and the bcftools tool was used to remove the heterozygote from the initial mutation site. homozygous After removing the variant sites, the updated initial variant sites for each object are obtained. From these updated initial variant sites, initial variant sites that satisfy the condition "the ratio of the first quantity to the first total quantity of objects is less than a first threshold" are selected. For example, the bcftools tool can be used to... The mutation site was set to deletion, and bcftools was used to convert the heterozygote to a deletion site. homozygous The variant sites were removed. Here, DP is the sequencing depth, GQ is the genotype quality value, QUAL is the variant quality value, and VAF is the variant allele frequency. Setting it to "missing" means that among the initial variant sites of a given object, if a certain initial variant site satisfies... If so, it is determined that the initial mutation site of the object is missing, that is, there is no initial site in the object that satisfies the condition. The variant sites. During quality control, the significance level for the Hardy-Weinberg equilibrium test was uniformly set at 1×10⁻⁶. -6 .

[0100] In this embodiment, by taking any initial mutation site as the target mutation site, a first number of objects with target mutation sites is determined. When the ratio of the first number to the first total number of objects is less than a first threshold, the target mutation sites are removed from the initial mutation sites of the objects until the initial mutation sites are traversed, thus obtaining the mutation sites of each object. The number of sites in the union set formed by the mutation sites of each object is determined. When the ratio of the second number of mutation sites of the undetermined object to the total number of sites is less than a second threshold, the undetermined object is removed from the objects, thus obtaining the mutation sites of the target object. This can remove unqualified mutation sites and objects, thereby making the final determined target association comparison more accurate.

[0101] In one embodiment, the molecular phenotype data is proteomic data, the molecular phenotype is a protein, and the molecular abundance is the protein expression value. The process of obtaining the molecular phenotype and molecular abundance of the molecular phenotype data in step S1 includes:

[0102] The proteome data of the target object is detected by a protein detection instrument to obtain the initial protein and the protein expression value of the initial protein.

[0103] The protein expression value of the initial protein is multiplied by the corresponding correction factor to obtain the corrected expression value. The corrected expression value is then logarithmically transformed to obtain the converted expression value of the initial protein. Using any initial protein as a candidate protein, a third number of target objects containing candidate proteins is determined. When the ratio of this third number to the second total number of target objects is less than a third threshold, candidate proteins are removed from the initial proteins of the target objects. This process continues until all initial proteins have been traversed, resulting in the proteins of each target object. These proteins include a first protein with a converted expression value and a second protein without a converted expression value.

[0104] Based on the conversion expression value of the first protein, the conversion expression value of the second protein is completed using a random tail distribution imputation method to obtain the conversion expression value of each protein;

[0105] Batch correction and anti-logarithmic transformation were performed on the converted expression values ​​of the target protein to obtain the protein expression values ​​of the target protein.

[0106] Protein detection instruments refer to instruments that can detect proteins and protein expression values ​​using proteomics data. For example, a mass spectrometer.

[0107] Initial proteins are all proteins detected from the proteome data. Initial proteins include proteins for which protein expression values ​​were detected and proteins for which protein expression values ​​were not detected.

[0108] The correction factor is a preset value, and each protein corresponds to a specific correction factor. The correction factors for all proteins may be the same or different. When multiplying by the correction factor, only the detected protein expression value is multiplied by the corresponding correction factor.

[0109] The expression for the logarithmic transformation is Y=log b (X), Y is the transformed expression value, X is the corrected expression value, and b is the base.

[0110] Among the proteins of each target, if a protein expression value is detected, then the protein is the first protein with a conversion expression value; if a protein expression value is not detected, then the protein is the second protein without a conversion expression value.

[0111] The random tail distribution imputation method fills in the missing values ​​by simulating the mechanism of "low abundance leading to undetectability" and using low normal distribution random values.

[0112] Since the proteomic data for each target protein were detected in batches, the correction of the converted expression values ​​of the target proteins was also performed in batches. Furthermore, the correction coefficients used for batch correction could be preset, and the correction coefficients for each protein could be the same or different.

[0113] First, batch correction is performed, and then the results of the batch correction are subjected to anti-logarithmic transformation. The expression for anti-logarithmic transformation is Z=b. Q b is the base used in the logarithmic transformation, Z is the protein expression value after antilogarithmic transformation, and Q is the result of batch correction of the transformed expression value.

[0114] In this embodiment, the converted expression value of the initial protein is obtained by logarithmic transformation of the corrected expression value. This can reduce data skewness, stabilize variance, and make it suitable for subsequent statistical analysis.

[0115] In one embodiment, the method further includes:

[0116] Principal component analysis was performed based on the protein expression values ​​of each target object to obtain the principal components of each target object and the scores and eigenvalues ​​of each principal component.

[0117] Based on the scores and eigenvalues ​​of each principal component of the target object, calculate the first Mahalanobis distance of the target object; determine the fourth threshold associated with the number of principal components and the preset confidence level from the preset chi-square distribution table;

[0118] Target objects with a first Mahalanobis distance greater than the fourth threshold are removed, and the association significance probability value is determined by the variant sites and protein expression values ​​of the unremoved target objects.

[0119] Principal component analysis is also known as PCA.

[0120] After performing principal component analysis based on the protein expression values ​​of each target object, the principal components, scores, and eigenvalues ​​of each principal component are obtained.

[0121] The formula for calculating the first Mahalanobis distance is: D is the first Mahalanobis distance, k is the number of principal components, and t is the distance between the principal components and the first principal components. i Let be the score of the i-th principal component. Let be the eigenvalue of the i-th principal component.

[0122] The number of principal components refers to the number of principal components participating in the calculation of the first Mahalanobis distance. The number of principal components can be preset.

[0123] The chi-square distribution table is presented in the form of a two-dimensional matrix. Its core function is to provide critical values ​​at different degrees of freedom and different significance levels. Determining the fourth threshold associated with the number of principal components and the preset confidence level from the pre-defined chi-square distribution table includes: determining the rows representing the number of principal components with degrees of freedom from the chi-square distribution table, the columns representing the significance level as 1 minus the confidence level, and determining the fourth threshold by the value at the intersection of the columns and rows.

[0124] Furthermore, the confidence level is 95%.

[0125] In this embodiment, principal component analysis is performed based on the protein expression values ​​of each target object to obtain the principal components, scores, and eigenvalues ​​of each principal component. Based on the scores and eigenvalues ​​of each principal component of the target object, the first Mahalanobis distance of the target object is calculated. A fourth threshold, which is associated with the number of principal components and a preset confidence level, is determined from a preset chi-square distribution table. Target objects with a first Mahalanobis distance greater than the fourth threshold are removed. This further removes outliers and ensures the reliability of subsequent processing results.

[0126] In one embodiment, the molecular phenotype data is metabolomics data, the molecular phenotype is the target metabolite, and the molecular abundance is the concentration of the metabolite. The process of obtaining the molecular phenotype and the molecular abundance of the molecular phenotype in step S1 includes:

[0127] Samples were taken and mixed from bodily fluid samples of each target subject, and the mixed samples were divided into multiple quality control samples. The mass spectrometry response value and concentration of the first candidate metabolite in the multiple quality control samples were detected in batches using a metabolite detection instrument.

[0128] Based on the mass spectrometry response values ​​of the first candidate metabolites in each quality control sample, the coefficient of variation of each first candidate metabolite is calculated.

[0129] First candidate metabolites with a coefficient of variation greater than the fifth threshold are removed to obtain second candidate metabolites in each quality control sample. Then, mass spectrometry response values ​​are filled in for metabolites with missing mass spectrometry response values ​​in the second candidate metabolites of each quality control sample to obtain the mass spectrometry response value of each second candidate metabolite in each quality control sample.

[0130] Based on the detection order of each quality control sample and the mass spectrometry response value of each second candidate metabolite, the Spearman rank correlation coefficient of each second candidate metabolite is calculated. Second candidate metabolites with Spearman rank correlation coefficients greater than the sixth threshold are eliminated to obtain the target metabolite, and the concentration of the target metabolite for each target object is determined.

[0131] The metabolite detection instrument measures the metabolomics data of quality control samples and the metabolomics data of body fluid samples from each subject. Since the instrument cannot detect too much data at once, batch testing is necessary. In batch testing, each batch includes metabolomics data from both the quality control samples and the subject's body fluid samples.

[0132] Among the first candidate metabolites of each quality control sample, there were metabolites with mass spectrometry response values ​​and metabolites without mass spectrometry response values.

[0133] The formula for calculating the coefficient of variation is: CV is the coefficient of variation. The standard deviation of the mass spectrometry response values ​​of the first candidate metabolite in each quality control sample is given. This represents the mean mass spectrometry response value of the first candidate metabolite in each quality control sample. A smaller coefficient of variation indicates that the metabolite is stable and the detection results of the metabolite detection instrument are reliable.

[0134] The fifth and sixth thresholds are pre-set data.

[0135] The mass spectrometry response value completion process includes: identifying missing metabolites with missing mass spectrometry response values ​​from the first candidate metabolites of each quality control sample; identifying target quality control samples from the quality control samples where mass spectrometry response values ​​of missing metabolites are detected; determining the minimum mass spectrometry response value from the missing metabolites in each target quality control sample; and taking half of the minimum value as the mass spectrometry response value of the missing metabolite with the missing mass spectrometry response value. For example, if no mass spectrometry response value of metabolite 1 is detected in quality control sample A, but the minimum mass spectrometry response value of metabolite 1 detected in other quality control samples is X, then X / 2 is taken as the mass spectrometry response value of metabolite 1 in quality control sample A.

[0136] The formula for calculating the Spearman rank correlation coefficient of the second candidate metabolite is as follows: S represents the testing order of each quality control sample, and X represents... j Let R() represent the mass spectrometry response value of the second candidate metabolite j in each quality control sample, where R() represents the rank vector and cov() represents the covariance of the rank vector. Let S be the standard deviation of the rank vector R(S). R(X) is a rank vector j The standard deviation of ).

[0137] If the Spearman rank correlation coefficient of the second candidate metabolite is greater than the sixth threshold, it indicates that the mass spectrometry response value of the second candidate metabolite is not determined by the biological differences of the sample itself, but by the detection order. Therefore, the mass spectrometry response value of the second candidate metabolite is unreliable and needs to be removed.

[0138] PCA analysis was performed on the quality control samples and body fluid samples from each target subject to examine the degree of aggregation of the quality control samples. Since the quality control samples were collected before, during, and after the detection, their aggregation degree can reflect the performance stability of the metabolite detection instrument during the detection period. Secondly, if the detection stability of the metabolite detection instrument is poor, there will be the lowest correlation between the first and last quality control samples. Therefore, the correlation of the mass spectrometry response values ​​of the first and last quality control samples can also reflect the detection stability of the metabolite detection instrument.

[0139] The concentration of the target metabolite for each target object can be obtained using a concentration detection instrument.

[0140] Furthermore, the fifth threshold is 0.2.

[0141] Furthermore, the sixth threshold is 0.3.

[0142] In this embodiment, by eliminating first candidate metabolites with a coefficient of variation greater than the fifth threshold, second candidate metabolites are obtained in each quality control sample. Then, a mass spectrometry response value completion operation is performed on the missing mass spectrometry response values ​​of the second candidate metabolites in each quality control sample to obtain the mass spectrometry response value of each second candidate metabolite in each quality control sample. Based on the detection order of each quality control sample and the mass spectrometry response value of each second candidate metabolite, the Spearman rank correlation coefficient of each second candidate metabolite is calculated. Second candidate metabolites with a Spearman rank correlation coefficient greater than the sixth threshold are eliminated to obtain the target metabolite. This method can eliminate unstable metabolites and metabolites with low accuracy in mass spectrometry response values ​​to avoid interfering with subsequent analysis processes.

[0143] In one embodiment, the method further includes: performing principal component analysis based on the concentration of target metabolites of each target object to obtain target principal components and their scores and eigenvalues ​​for each target object; calculating the second Mahalanobis distance of the target objects based on the scores and eigenvalues ​​of each target principal component; determining an object screening threshold associated with the number of target principal components and a preset confidence level from a preset chi-square distribution table; and removing target objects whose second Mahalanobis distance is greater than the object screening threshold, so as to determine the significance probability value of the association based on the concentration of target metabolites of the unremoved target objects. Further, the confidence level is 95%.

[0144] In one embodiment, step S2, which involves selecting a first association pair from the association pairs based on the association significance probability value and conditional analysis, includes:

[0145] Principal component analysis was performed based on the molecular abundance of the molecular phenotypes of each target object to obtain the variance contribution rate of each principal component.

[0146] Determine the minimum number of principal components required for the cumulative variance contribution rate to reach the seventh threshold, and determine the screening threshold based on the minimum number of principal components;

[0147] Remove association pairs whose association significance probability value is greater than the screening threshold to obtain initial candidate association pairs;

[0148] Identify the variant sites with anomalous distributions from the variant sites of the initial candidate association pairs, and remove the initial candidate association pairs with variant sites with anomalous distributions to obtain the initial association pairs.

[0149] Conditional analysis was performed on the initial association pairs, and the initial association pairs with independent genetic effects were selected as the first association pairs.

[0150] The minimum number of principal components required for the cumulative variance contribution rate to reach the seventh threshold refers to the sum of the variance contribution rates of each principal component in sequence. If the cumulative variance contribution rate after the sum reaches the seventh threshold, the summation stops, and the number of principal components currently participating in the summation is the minimum number of principal components.

[0151] The formula for the filtering threshold is: N is the number of minimum principal components.

[0152] Distribution anomalies refer to situations where the distribution pattern of variant sites deviates from expected statistical regularities or biological common sense. For example, if a variant site involves more than 7 chromosomes, then the association pairs containing that variant site are all initial candidate association pairs that need to be removed.

[0153] The conditional analysis of the initial association pairs, and the designation of the initial association pairs with independent genetic effects as the first association pairs, includes: performing conditional analysis on the initial association pairs to obtain the conditional analysis results for each initial association pair; based on the conditional analysis results of each initial association pair, determining the initial association pairs containing variant sites with independent genetic effects, and designating the initial association pairs containing variant sites with independent genetic effects as the first association pairs. The conditional analysis results are used to characterize whether the variant sites in the initial association pairs are variant sites with independent genetic effects. An independent genetic effect refers to the influence of a variant site on a disease trait or disease risk, which exists independently of other variant sites and is not a spurious association caused by a "linkage" relationship with other variant sites. The conditional analysis can be performed using the GCTA-COJO (Genome-wide Complex Trait Analysis - Conditional & Jointassociation analysis) module in the GCTA (Genome-wide Complex Trait Analysis) software toolkit. By screening the initial association pairs through conditional analysis, the first association pair is obtained. This can eliminate the interference of linkage disequilibrium on the association signal, thereby identifying the first association pair of independent genetic effects.

[0154] In this embodiment, principal component analysis is performed based on the molecular abundance of the molecular phenotype of each target object to obtain the variance contribution rate of each principal component. The minimum number of principal components required for the cumulative variance contribution rate to reach the seventh threshold is determined, and a screening threshold is determined based on the minimum number of principal components. Association pairs with association significance probability values ​​greater than the screening threshold are eliminated to obtain initial candidate association pairs. From the variant sites of each variant site in the initial candidate association pairs, variant sites with abnormal distribution are identified, and initial candidate association pairs with variant sites with abnormal distribution are eliminated to obtain the first association pair. This can avoid unreliable association pairs from affecting the subsequent determination of target association pairs, thereby improving the accuracy and reliability of target association pairs.

[0155] In one embodiment, the process of determining the third association pair includes:

[0156] For each second association pair, a null model is constructed with the molecular abundance of the molecular phenotype in the second association pair as the dependent variable and the attribute information of the target object to which the second association pair belongs as the independent variable; a main effect model is constructed with the molecular abundance of the molecular phenotype as the dependent variable and the attribute information and variant sites as the independent variables; and an interaction effect model is constructed with the molecular abundance of the molecular phenotype as the dependent variable and the attribute information, variant sites, and the interaction term between the variant sites and the disease state of the target object to which they belong as the independent variables.

[0157] Based on the null model, main effect model and interaction effect model of each second association pair, the likelihood ratio test is performed to obtain the interaction P value of each second association pair; the interaction P value represents the degree of improvement of the fit of the interaction effect model after the introduction of the interaction term.

[0158] Based on the interaction P-value of each second association pair, the false discovery rate of each second association pair is calculated using the qvalue method, and the second association pair with a false discovery rate less than the eighth threshold is identified as the third association pair.

[0159] The constructed null model, main effect model, and interaction effect model are as follows: Zero model: ; Main effects model: ; Interaction effect model: ;

[0160] Y represents the molecular abundance of the molecular phenotype, SNP represents the variant site, D represents the disease state, A represents age, S represents sex, and G represents the molecular abundance of the phenotype. ba Indicates the carrier status of GBA (Glucocere brosidase gene), PC1 to PC 10 The first 10 principal components, This represents the interaction coefficient between the variant site and the disease state. For residuals, Main effect coefficient This represents the main effect coefficient of the disease. The intercept is... , , and All are covariate coefficients. Attribute information includes age, gender, and GBA portability.

[0161] The false discovery rate (FDR) for each second association pair is calculated by performing multiple hypothesis testing on the interaction p-value using the qvalue method.

[0162] Furthermore, the eighth threshold is 0.05.

[0163] In this embodiment, for each second association pair, a null model is constructed with the molecular abundance of the molecular phenotype in the second association pair as the dependent variable and the attribute information of the target object to which the second association pair belongs as the independent variable; a main effect model is constructed with the molecular abundance of the molecular phenotype as the dependent variable and the attribute information and variant sites as the independent variables; and an interaction effect model is constructed with the molecular abundance of the molecular phenotype as the dependent variable and the attribute information, variant sites, and the interaction term between the variant sites and the disease state of the target object to which they belong as the independent variables. Based on the null model, main effect model, and interaction effect model of each second association pair, a likelihood ratio test is performed to obtain the interaction P-value of each second association pair. Based on the interaction P-value of each second association pair, the qvalue method is used to calculate the false discovery rate of each second association pair. In this way, third association pairs with a false discovery rate less than the eighth threshold can be obtained.

[0164] In one embodiment, step S4 includes:

[0165] Based on the maximum likelihood estimation results of the main effect model and interaction effect model of each third association pair, the main effect coefficient of the variant site and the interaction term coefficient between the variant site and the disease state in each third association pair are obtained.

[0166] The main effect coefficient is used as the estimate of the first effect, and the sum of the main effect coefficient and the interaction term coefficient is used as the estimate of the second effect.

[0167] By performing maximum likelihood estimation on the main effect model and interaction effect model of each third association pair, the main effect coefficient of the variant site and the interaction term coefficient between the variant site and the disease state in each third association pair can be obtained.

[0168] In this embodiment, the main effect coefficients of the variant sites and the interaction term coefficients between the variant sites and the disease state in each third association pair are obtained by using the maximum likelihood estimation results based on the main effect model and the interaction effect model of each third association pair. In this way, the main effect coefficients can be used as the first effect estimate, and the sum of the main effect coefficients and the interaction term coefficients can be used as the second effect estimate.

[0169] In one embodiment, the process of determining the target association pair in step S5 includes:

[0170] If the absolute value of the second effect estimate of the third association pair is greater than the absolute value of the first effect estimate, and the absolute value of the first effect estimate is less than the no-effect threshold, then the third association pair is determined to be a target association pair related to Parkinson's disease.

[0171] Among them, the variant sites in the target association pairs only have a significant impact on the molecular phenotype in the Parkinson's disease state.

[0172] Furthermore, if the absolute value of the first effect estimate of the third association pair is greater than the absolute value of the second effect estimate, and the absolute value of the second effect estimate is less than the no-effect threshold, the third association pair is determined to be a health-specific regulatory association pair. The variant sites in the health-specific regulatory association pair have a significant impact on the molecular phenotype in a healthy state, but the regulatory effect is lost or weakened in the Parkinson's disease state.

[0173] Furthermore, if the second effect estimate of the first association pair has the opposite sign to the first effect estimate, then the direction of the influence of the variant site on the molecular phenotype in the third association pair is reversed in the disease state.

[0174] In this embodiment, by comparing the absolute values ​​of the second effect estimate and the first effect estimate, and by comparing the absolute value of the first effect estimate with the no-effect threshold, target association pairs related to Parkinson's disease can be accurately screened.

[0175] In one embodiment, the method further includes:

[0176] Based on the physical location information of the variant sites in the chromosome in each third association pair, each third association pair is divided into multiple sets;

[0177] Colocation analysis was used to screen out target sets with consistent genetic regulatory effects from the set;

[0178] Based on the molecular phenotypes in the target set and the Parkinson's disease status of each third association pair in the target set, the causal effects between molecular phenotypes and the causal effects of molecular phenotypes on disease phenotypes were evaluated using Mendelian randomization analysis, resulting in multiple causal regulatory pathways.

[0179] Based on the causal regulation path, a regulation network is constructed, and key nodes in the regulation network are selected according to the network centrality index or importance score of each node in the regulation network.

[0180] Functional annotations were performed on key nodes, and pathway enrichment analysis was conducted on molecular phenotypes within the key nodes.

[0181] In each set, the variant sites of the third association pairs are in the same or adjacent positions.

[0182] Colocalization analysis can assess whether variant sites in different target sets are driven by the same underlying causal variant, thereby obtaining a target set with consistent genetic regulatory effects. A target set with consistent genetic regulatory effects means that the variant sites in the target set point to the same causal mechanism at two different biological levels (such as gene expression and disease risk).

[0183] Based on the molecular phenotypes in the target set and the Parkinson's disease status of each third association in the target set for their respective target subjects, this study uses Mendelian randomization analysis to assess the causal effects between molecular phenotypes and the causal effects of molecular phenotypes on disease phenotypes, deriving multiple causal regulatory pathways. The specific process includes: constructing causal relationship models between different omics molecular phenotypes and between them and disease phenotypes based on the target set; using the target set as an instrumental variable, evaluating the causal effects of proteins on metabolites and the causal effects of proteins or metabolites on disease phenotypes using Mendelian randomization analysis. During the analysis, various robustness testing methods can be used to verify the causal inference results, including the weighted median method, weighted pattern method, MR-Egger regression, and Steiger directionality test, to ensure that the causal effect estimation results are not affected by multi-level effects or single instrumental variables, thereby improving the reliability and robustness of the inference results. Based on the above causal effect analysis, a multi-level causal regulatory pathway (variant site → protein → metabolite → disease phenotype) is established, which acts on the disease phenotype step by step through the variant site via protein and metabolite, thereby clarifying the direction and intensity of the role of different omics molecules in the occurrence and development of the disease.

[0184] Based on causal regulatory pathways, a regulatory network was constructed. Key nodes in the network were selected based on network centrality or importance scores. This involved: systematically integrating all causal regulatory pathways to construct a cross-omics regulatory network; using proteins, metabolites, and their corresponding variant sites in each causal regulatory pathway as network nodes, and identifying regulatory relationships confirmed by causal inference and co-localization analysis as network edges, forming a multi-level network structure with variant sites as starting points, molecular phenotypes as intermediate nodes, and disease phenotypes as endpoints; identifying key nodes and regulatory modules through network topology analysis, including proteins and metabolites that regulate multiple downstream molecules or disease phenotypes, as well as nodes that repeatedly appear in different causal regulatory pathways; and weighting each node based on its genetic regulatory strength, causal effect magnitude, and frequency of repetition in multiple causal regulatory pathways to screen for key nodes that significantly and robustly influence disease phenotypes. This multi-omics regulatory network can be used to systematically analyze the overall regulatory mechanisms of multi-omics molecules in disease development and provides an intuitive reference for key regulatory molecules and potential intervention targets.

[0185] Functional annotation of key nodes and pathway enrichment analysis of molecular phenotypes at these nodes include: functional annotation and biological interpretation of key proteins and metabolites based on regulatory networks; gene function annotation and regulatory region annotation of key variant sites; and pathway enrichment analysis of key proteins and metabolites to determine their involvement in biological processes and metabolic pathways. Based on this, and combining multi-omics causal regulatory pathways and networks, a comprehensive analysis is conducted to elucidate the potential mechanisms by which genetic variations influence disease development through multi-level molecular regulation, thereby achieving a systematic mechanistic explanation from genetic variation to molecular phenotype to disease phenotype.

[0186] In this application, each third association pair is divided into multiple sets based on the physical location of the variant sites in the chromosome. Colocalization analysis is used to screen target sets with consistent genetic regulatory effects from these sets. Based on the molecular phenotypes in the target sets and the Parkinson's disease status of the target individuals belonging to each third association pair, Mendelian randomization analysis is used to assess the causal effects between molecular phenotypes and the causal effects of molecular phenotypes on disease phenotypes, resulting in multiple causal regulatory pathways. Based on these causal regulatory pathways, a regulatory network is constructed, and key nodes in the network are screened based on network centrality indicators or importance scores. Functional annotation is performed on these key nodes, and pathway enrichment analysis is conducted on the molecular phenotypes within them. This allows for the extraction of a causal disease network diagram from massive amounts of genomic data, precisely locating the key molecules and pathways that play a decisive role in the network, thereby revealing the deep molecular mechanisms underlying the development of Parkinson's disease.

[0187] In one embodiment, each target association pair can be divided into multiple sets based on the physical location information of the variant sites in the chromosome, so as to perform colocalization analysis, pathway enrichment analysis and functional annotation through the multiple sets divided by each target association pair.

[0188] The technical roadmap of this application is as follows: Figure 3 As shown. Specifically, the first association pair is determined, and differential analysis is performed on the normal association and disease association pairs in the first association pair to obtain the second association pair. The third association pair is then screened, the target association pair is screened, and pathway enrichment analysis and functional annotation are performed on key nodes.

[0189] The analysis flowchart of this application is as follows: Figure 4 As shown in the figure. Significance analysis refers to identifying target association pairs related to Parkinson's disease. Colocalization analysis can be performed using COLOC, and pathway enrichment analysis can be performed using the hypergeometric distribution test.

[0190] The flowchart of the disease-specific analysis in this application is as follows: Figure 5 As shown in the diagram, the integration, colocalization analysis, pathway enrichment analysis, and functional annotation of target association pairs related to Parkinson's disease are illustrated below. Figure 6 As shown.

[0191] It should be understood that although the steps in the flowcharts of the embodiments described above are shown sequentially according to the arrows, these steps are not necessarily executed in the order indicated by the arrows. Unless explicitly stated herein, there is no strict order restriction on the execution of these steps, and they can be executed in other orders. Moreover, at least some steps in the flowcharts of the embodiments described above may include multiple steps or multiple stages. These steps or stages are not necessarily completed at the same time, but can be executed at different times. The execution order of these steps or stages is not necessarily sequential, but can be performed alternately or in turn with other steps or at least some of the steps or stages of other steps.

[0192] Based on the same inventive concept, this application also provides a device for identifying disease-specific quantitative trait loci based on multi-omics integration to implement the aforementioned method for identifying disease-specific quantitative trait loci based on multi-omics integration. The solution provided by this device is similar to the solution described in the above method. Therefore, the specific limitations of one or more embodiments of the device for identifying disease-specific quantitative trait loci based on multi-omics integration provided below can be found in the limitations of the method for identifying disease-specific quantitative trait loci based on multi-omics integration described above, and will not be repeated here.

[0193] In one embodiment, a disease-specific quantitative trait locus identification device based on multi-omics integration is provided, comprising:

[0194] The information acquisition unit is used to acquire variant sites in the whole genome sequencing data of each target object, molecular phenotypes in the molecular phenotype data, and molecular abundance of molecular phenotypes.

[0195] The first screening unit is used to determine the association significance probability value of the association pair consisting of the variant site and the molecular phenotype for each target object based on the molecular abundance of the variant site and the molecular phenotype, and to screen the first association pair from the association pair based on the association significance probability value and conditional analysis; the first association pair includes normal association pairs of healthy objects and disease association pairs of objects with Parkinson's disease.

[0196] The second screening unit is used to identify the second association pairs that are consistent between the normal association pairs and the disease association pairs, and to identify the third association pairs that have disease interaction effects among the second association pairs.

[0197] The effect value calculation unit is used to calculate the first effect estimate relative to healthy subjects and the second effect estimate relative to subjects with Parkinson's disease for each third association pair;

[0198] The third screening unit is used to determine the target association pairs related to Parkinson's disease in the third association pairs based on the first effect estimate and the second effect estimate of each third association pair. The target association pairs are quantitative trait loci.

[0199] Each module in the aforementioned disease-specific quantitative trait locus identification device based on multi-omics integration can be implemented entirely or partially through software, hardware, or a combination thereof. These modules can be embedded in or independent of the processor in a computer device in hardware form, or stored in the memory of a computer device in software form, so that the processor can call and execute the corresponding operations of each module.

[0200] In one embodiment, a computer device is provided, which may be a server, and its internal structure diagram may be as follows: Figure 7As shown, the computer device includes a processor, memory, and a network interface connected via a system bus. The processor provides computational and control capabilities. The memory includes non-volatile storage media and internal memory. The non-volatile storage media stores the operating system, computer programs, and a database. The internal memory provides an environment for the operation of the operating system and computer programs in the non-volatile storage media. The database stores various types of data. The network interface communicates with external terminals via a network connection. When executed by the processor, the computer program implements a disease-specific quantitative trait locus identification method based on multi-omics integration.

[0201] Those skilled in the art will understand that Figure 7 The structure shown is merely a block diagram of a portion of the structure related to the present application and does not constitute a limitation on the computer device to which the present application is applied. Specific computer devices may include more or fewer components than those shown in the figure, or combine certain components, or have different component arrangements.

[0202] In one embodiment, a computer device is also provided, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps in the above method embodiments.

[0203] In one embodiment, a computer-readable storage medium is provided having a computer program stored thereon that, when executed by a processor, implements the steps in the above method embodiments.

[0204] In one embodiment, a computer program product is provided, including a computer program that, when executed by a processor, implements the steps in the above method embodiments.

[0205] It should be noted that the user information (including but not limited to user device information, user personal information, etc.) and data (including but not limited to data used for analysis, data stored, data displayed, etc.) involved in this application are all information and data authorized by the user or fully authorized by all parties.

[0206] Those skilled in the art will understand that all or part of the processes in the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium. When executed, the computer program can include the processes of the embodiments described above. Any references to memory, databases, or other media used in the embodiments provided in this application can include at least one of non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical memory, high-density embedded non-volatile memory, resistive random access memory (ReRAM), magnetic random access memory (MRAM), ferroelectric random access memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory can include random access memory (RAM) or external cache memory, etc. By way of illustration and not limitation, RAM can take many forms, such as Static Random Access Memory (SRAM) or Dynamic Random Access Memory (DRAM). The databases involved in the embodiments provided in this application may include at least one type of relational database and non-relational database. Non-relational databases may include, but are not limited to, blockchain-based distributed databases. The processors involved in the embodiments provided in this application may be general-purpose processors, central processing units, graphics processing units, digital signal processors, programmable logic devices, quantum computing-based data processing logic devices, etc., and are not limited to these.

[0207] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0208] The embodiments described above are merely illustrative of several implementation methods of this application, and while the descriptions are specific and detailed, they should not be construed as limiting the scope of this patent application. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of this application, and these all fall within the protection scope of this application. Therefore, the protection scope of this application should be determined by the appended claims.

Claims

1. A method for identifying disease-specific quantitative trait loci based on multi-omics integration, characterized in that, The method includes: S1. Obtain the variant sites in the whole genome sequencing data of each target object, the molecular phenotypes in the molecular phenotype data, and the molecular abundance of the molecular phenotypes; S2. For each target object, based on the molecular abundance of the variant site and the molecular phenotype, determine the association significance probability value of the association pair formed by the variant site and the molecular phenotype, and screen the first association pair from the association pair based on the association significance probability value and conditional analysis; the first association pair includes normal association pairs of healthy objects and disease association pairs of objects with Parkinson's disease. S3. Determine a second association pair that is consistent with the normal association pair and the diseased association pair, and determine a third association pair in the second association pair that has a disease interaction effect; S4. Calculate the first effect estimate of each of the third association pairs relative to the healthy subjects and the second effect estimate relative to the subjects with Parkinson's disease; S5. Based on the first effect estimate and the second effect estimate of each of the third association pairs, determine the target association pairs related to Parkinson's disease in the third association pairs, and use the target association pairs as the identified quantitative trait loci.

2. The method according to claim 1, characterized in that, The process of obtaining variant sites in the whole genome sequencing data of the target object in step S1 includes: Acquire human reference genome and whole genome sequencing data of various subjects; Based on the human reference genome, a genome analysis toolkit was used to detect variants in each of the whole genome sequencing data to obtain the initial variant sites of each of the objects. Using any of the initial mutation sites as target mutation sites, determine a first number of objects containing the target mutation sites. When the ratio of the first number to the first total number of objects is less than a first threshold, remove the target mutation sites from the initial mutation sites of the objects until the initial mutation sites have been traversed, and obtain the mutation sites of each object. The number of sites in the union set formed by the variant sites of each of the objects is determined. When the ratio of the second number of variant sites of the object to the total number of sites is less than a second threshold, the object to be determined is removed from the objects to obtain the variant sites of the target object. The object to be determined is any one of the objects.

3. The method according to claim 1, characterized in that, The molecular phenotypic data is proteomic data, the molecular phenotype is a protein, and the molecular abundance is the protein expression value of the protein. The process of obtaining the molecular phenotype and the molecular abundance of the molecular phenotype in step S1 includes: The proteome data of the target object are detected by a protein detection instrument to obtain the initial protein of the target object and the protein expression value of the initial protein; The protein expression value of the initial protein is multiplied by the corresponding correction factor to obtain the corrected expression value. The corrected expression value is then logarithmically transformed to obtain the converted expression value of the initial protein. Using any one of the initial proteins as a candidate protein, a third number of target objects containing the candidate protein is determined. When the ratio of this third number to the second total number of target objects is less than a third threshold, the candidate protein is removed from the initial proteins of the target objects. This process continues until all initial proteins have been traversed, resulting in the proteins of each target object. The proteins include a first protein with a converted expression value and a second protein without a converted expression value. Based on the conversion expression value of the first protein, the conversion expression value of the second protein is completed using a random tail distribution imputation method to obtain the conversion expression value of each protein. Batch correction and anti-logarithmic transformation are performed on the converted expression value of the protein in the target object to obtain the protein expression value of the protein in the target object.

4. The method according to claim 3, characterized in that, The method further includes: Principal component analysis was performed based on the protein expression values ​​of the proteins of each target object to obtain the principal components of each target object and the scores and eigenvalues ​​of each principal component; Based on the scores and eigenvalues ​​of each principal component of the target object, a first Mahalanobis distance of the target object is calculated; a fourth threshold associated with the number of principal components and a preset confidence level is determined from a preset chi-square distribution table; Target objects whose first Mahalanobis distance is greater than the fourth threshold are removed, and the association significance probability value is determined by the mutation sites and protein expression values ​​of the unremoved target objects.

5. The method according to claim 1, characterized in that, The molecular phenotype data is metabolomics data, the molecular phenotype is the target metabolite, and the molecular abundance is the concentration of the target metabolite. The process of obtaining the molecular phenotype and the molecular abundance of the molecular phenotype in step S1 includes: Samples were taken and mixed from bodily fluid samples of each target object, and the mixed samples were divided into multiple quality control samples. The mass spectrometry response value and concentration of the first candidate metabolite in the multiple quality control samples were detected in batches using a metabolite detection instrument. Based on the mass spectrometry response values ​​of the first candidate metabolites in each of the quality control samples, the coefficient of variation of each first candidate metabolite is calculated. First candidate metabolites with a coefficient of variation greater than the fifth threshold are removed to obtain second candidate metabolites in each quality control sample. Then, mass spectrometry response values ​​are filled in for the second candidate metabolites in each quality control sample that are missing mass spectrometry response values ​​to obtain the mass spectrometry response value of each second candidate metabolite in each quality control sample. Based on the detection order of each quality control sample and the mass spectrometry response value of each second candidate metabolite, the Spearman rank correlation coefficient of each second candidate metabolite is calculated. Second candidate metabolites with Spearman rank correlation coefficients greater than the sixth threshold are eliminated to obtain the target metabolite, and the concentration of the target metabolite of each target object is determined.

6. The method according to claim 1, characterized in that, Step S2, which involves selecting the first association pair from the association pairs based on the association significance probability value and conditional analysis, includes: Principal component analysis was performed based on the molecular abundance of the molecular phenotypes of each target object to obtain the variance contribution rate of each principal component. Determine the minimum number of principal components required for the cumulative variance contribution rate to reach the seventh threshold, and determine the screening threshold based on the minimum number of principal components; Remove association pairs whose association significance probability value is greater than the screening threshold to obtain initial candidate association pairs; From the variant sites of each of the initial candidate association pairs, identify variant sites with abnormal distribution, and remove the initial candidate association pairs with variant sites with abnormal distribution to obtain the initial association pairs; Conditional analysis was performed on the initial association pairs, and the initial association pairs with independent genetic effects were selected as the first association pairs.

7. The method according to claim 1, characterized in that, The process of determining the third association pair includes: For each of the second association pairs, a null model is constructed using the molecular abundance of the molecular phenotype in the second association pair as the dependent variable and the attribute information of the target object to which the second association pair belongs as the independent variable; a main effect model is constructed using the molecular abundance of the molecular phenotype as the dependent variable and the attribute information and the variant site as the independent variables; and an interaction effect model is constructed using the molecular abundance of the molecular phenotype as the dependent variable and the attribute information, the variant site, and the interaction term between the variant site and the disease state of the target object to which it belongs as the independent variables. Based on the null model, the main effect model, and the interaction effect model of each second association pair, a likelihood ratio test is performed to obtain the interaction P-value for each second association pair; the interaction P-value represents the degree of improvement in the fitting effect of the interaction effect model after the introduction of the interaction term. Based on the interaction P-value of each of the second association pairs, the false discovery rate of each of the second association pairs is calculated using the qvalue method, and the second association pairs with a false discovery rate less than the eighth threshold are identified as the third association pairs.

8. The method according to claim 7, characterized in that, Step S4 includes: Based on the maximum likelihood estimation results of the main effect model and interaction effect model of each of the third association pairs, the main effect coefficient of the variant site and the interaction term coefficient between the variant site and the disease state in each of the third association pairs are obtained. The main effect coefficient is used as the first effect estimate, and the sum of the main effect coefficient and the interaction term coefficient is used as the second effect estimate.

9. The method according to claim 8, characterized in that, The process of determining the target association pair in step S5 includes: If the absolute value of the second effect estimate of the third association pair is greater than the absolute value of the first effect estimate, and the absolute value of the first effect estimate is less than the no-effect threshold, then the third association pair is determined to be a target association pair related to Parkinson's disease.

10. The method according to claim 1, characterized in that, The method further includes: Based on the physical location information of the mutation sites in the chromosome in each of the third association pairs, each of the third association pairs is divided into multiple sets; The target set with consistent genetic regulatory effects was screened from the set using colocalization analysis. Based on the molecular phenotypes in the target set and the Parkinson's disease status of the target objects to which each of the third association pairs in the target set belongs, the causal effects between molecular phenotypes and the causal effects of molecular phenotypes on disease phenotypes are evaluated by Mendelian randomization analysis, and multiple causal regulatory pathways are obtained. Based on the causal regulation path, a regulation network is constructed, and key nodes in the regulation network are selected according to the network centrality index or importance score of each node in the regulation network. Functional annotations were performed on the key nodes, and pathway enrichment analysis was conducted on the molecular phenotypes in the key nodes.

Images