Method, system, equipment and medium for analyzing biological characteristic correlation in biomics data

Abstract

The invention discloses a method, a system, equipment and a medium for analyzing biological characteristic correlation in biomics data, and belongs to the technical field of bioinformatics. In particular to a method, a system, equipment and a medium for analyzing correlation of biomics data. According to the method, the system, the equipment and the medium, a high-performance and parallel block matrix algorithm is used for processing large-scale correlation and clustering tasks, excellent performance is shown in the aspect of calculating a correlation matrix, and downstream analysis can exceed computer memory resources; according to the method, the approximate P value can be quickly and effectively estimated without repeated complex cumulative distribution function calculation, and the corresponding approximate correction P value can be quickly and effectively estimated based on the P value sampling subset without being based on all the P values; and a top-k algorithm is realized, and faster execution and lower memory consumption are realized by utilizing a partial sorting method based on heaps.

Description

Technical Field

[0001] The present application relates to the technical field of bioinformatics, and in particular, to a method, system, device and medium for analyzing the correlation of biological features in bio-omics data. Background Art

[0002] Correlation analysis is widely used in bio-omics data to infer potential interaction networks within biological systems. In transcriptomics, researchers construct gene regulatory networks by analyzing the coordinated expression of different genes, which enables the identification of key transcriptional regulatory mechanisms and functional modules. In proteomics, correlation quantifies the strength of interactions between proteins and promotes the discovery of new protein complexes and their functions by analyzing protein activity or expression levels. In microbiomics, microbial co-occurrence relationships can be detected based on relative abundance, revealing complex interactions between microbial communities. In addition to its application in single-omics data, correlation analysis can also integrate multi-omics data for comprehensive analysis of complex biological systems.

[0003] In the past decade, the development of multi-omics technologies has generated a large amount of biological data, making it possible to comprehensively analyze biological processes and diseases. Many large studies have collected multi-omics datasets containing millions of tests. The high dimensionality and complexity of such data pose great computational challenges for correlation analysis.

[0004] There are many software packages for correlation analysis, including Numpy, Pandas, WGCNA, Cupy, RAPIDS, and CorALS. Among them, Numpy and Pandas are standard software in data analysis. Numpy uses C language to implement efficient multidimensional array operations and provides the numpy.corrcoef function to calculate the Pearson correlation coefficient matrix. Pandas provides a series of user-friendly statistical functions for tabular data structures. WGCNA focuses on the biomedical field and is specifically optimized for Pearson correlation coefficient calculation. Cupy, RAPIDS, and CorALS are generally used to analyze large-scale data. Cupy and RAPIDS use GPUs to accelerate correlation calculations, but require dedicated hardware support and are subject to memory constraints. CorALS uses a Numpy-based framework that combines dedicated vector projection with space segmentation technology to achieve faster correlation coefficient matrix calculations and top-k correlation approximations. Although these software packages perform well in specific scenarios, they still have limitations, such as high memory consumption and dependence on hardware support. Summary of the invention

[0005] In order to solve at least one of the above technical problems, the present application aims to provide a technical solution for correlation analysis of high-dimensional data, which can use high-performance, parallel algorithms to handle large-scale correlation and clustering tasks, so that downstream analysis can exceed computer memory resources. In order to achieve this purpose, the technical solution adopted by the present application is as follows: The first aspect of the present application provides a method for analyzing the correlation of biological features in biomics data, comprising the following steps: S1, obtaining one or more bio-omics data of multiple samples, and obtaining one or more matrices including biological feature quantitative data; S2, for any two matrices, if the product of the number of biometric features of the two matrices exceeds a first preset threshold, at least one of the two matrices is divided into a plurality of sub-matrices along the biometric feature dimension, and one of the matrices or each of its sub-matrices is subjected to matrix calculation respectively with the other matrix or each of its sub-matrices; S3, repeat step S2 until all matrices are calculated to obtain a correlation matrix, wherein the correlation matrix includes correlation coefficient vectors between each pair of the biometric features.

[0006] In the present application, the bio-omics includes but is not limited to one or more of genomics, transcriptomics, proteomics, epigenomics, metagenomics, and viromics (ie, multi-omics).

[0007] In this application, a bio-omics can obtain at least one matrix including biological features and sample numbers / serial numbers. The biological features include but are not limited to the biological information of the source of the biological sample (such as human, including age, height, weight, etc.), DNA (including cfDNA), and its corresponding modifications such as methylation, or corresponding variants; RNA, protein, virus, bacteria, etc. The format is as follows:

[0008] Different rows represent different biological features, and different columns represent different samples.

[0009] In this application, the correlation coefficient between two biometric features refers to the correlation coefficient of multiple biometric feature pairs, which represents the correlation between multiple pairs of two features. The correlation coefficients of different biometric feature pairs may be the same.

[0010] It should be noted that, in the present application, the purpose of dividing the matrix into blocks is to be able to complete the correlation analysis of a large number of (e.g., more than 1.6 million) biological features on a computer with a small memory (e.g., within 30GB). Based on this, if the number of biological features is small, different bio-omics matrices can also be merged. In some embodiments of the present application, the preset threshold depends on the size of the computing memory, that is, the result of multiplying the two matrices cannot exceed the memory upper limit (first preset threshold) of the computer on which the matrix calculation is to be performed.

[0011] In some embodiments of the present application, the any two matrices include any one matrix and itself, that is, the two matrices are the same matrix, so as to ensure that the correlation between any two biometric features can be obtained. Of course, it also means that the correlation between any one biometric feature and itself can be obtained.

[0012] In some embodiments of the present application, the matrix calculation is to convert the redundant calculation of the correlation coefficients between two biological features into matrix multiplication by using an appropriate vector transformation method. This conversion method not only simplifies the calculation process, but also can use efficient linear algebra routines, such as OpenBLAS, to effectively take advantage of the multi-core processors and SIMD instructions of modern CPUs, and further improve the calculation efficiency.

[0013] In some embodiments of the present application, combined with steps S2 and S3, the problem of computational resource limitations faced in large-scale data correlation analysis can be solved. By dividing the matrix into blocks, only a small portion of data is processed each time. In some specific embodiments of the present application, combined with lazy evaluation technology, the result of calculating and returning the submatrix according to the request can more effectively manage the use of memory. The correlation calculation of the submatrix utilizes vector transformations and highly parallelized modern linear algebra routines, which can significantly improve computing performance.

[0014] In some preferred embodiments of the present application, only one of the matrices may be partitioned (divided into multiple sub-matrices) to improve the versatility and scalability of the above method, and further calculate the correlation between each sub-matrix and the other matrix. Each row in the correlation matrix represents the correlation coefficient between a biological feature of the matrix and all biological features of the other matrix, which is convenient for subsequent analysis and application.

[0015] In some embodiments of the present application, in step S3, the correlation coefficient is selected from one of the Pearson correlation coefficient, the Spearman correlation coefficient and the Kendall correlation coefficient.

[0016] In some embodiments of the present application, the correlation matrix also includes P Value vector and / or correction P The method further comprises the following steps: S4, according to P Value or correction P The correlation matrix is ​​screened by the value to obtain the P Value or correction P The number of biometric pairs whose values ​​are less than the second preset threshold is denoted as d; S5, obtaining the first k biometric pairs with the largest correlation coefficient or obtaining the biometric pairs with a correlation coefficient greater than a third preset threshold, wherein 0 <k<d。

[0017] In some embodiments of the present application, the P The value is approximated based on linear interpolation. The specific steps include: Pre-acquisition L indivual P The values ​​are sorted from small to large. , and obtain its corresponding square correlation coefficient: ; For any correlation coefficient r , which is approximately P value Calculated using the following formula:

[0018] in,

[0019] In some specific embodiments of the present application, the L indivual P The values ​​from small to large are: , ,…,1.

[0020] In some specific embodiments of the present application, the following formula is used to obtain P The corresponding correlation coefficient r :

[0021]

[0022] Further obtain the square correlation coefficient Therefore, the present application can prepare a P Value and correlation coefficient r or statistical variables t The mapping table can be directly based on the calculated correlation coefficient r or statistical variables t Get approximation P value, thereby increasing the calculation speed.

[0023] In some embodiments of the present application,b Represents an extremely small probability value, much smaller than , which can be selected by those skilled in the art. In some embodiments of the present application, b Select 2.2e-16.

[0024] In some embodiments of the present application, the correction P The value is also approximated based on linear interpolation. The specific steps include: M The correlation coefficient of biometric pairs P The values ​​are sorted in ascending order as follows: ,…, ; use P The value correction method obtains M indivual P Correction of value P value : ,…, ; For one P Value, its approximate correction P value Calculated using the following formula:

[0025] in, M >10000,

[0026] In one embodiment of the present application, M The correlation coefficient of biometric pairs P The values ​​are sampled from the bio-omics data to be computed.

[0027] In another embodiment of the present application, the correlation coefficient can be directly established based on the same method as above. r or statistical variables t and mapping table to further speed up.

[0028] In some embodiments of the present application, the P The value correction method is selected from at least one of Bonferroni, Holm, Hochberg, Benjamini-Hochberg and Benjamini-Yekutieli. These correction methods are all based on small samples.

[0029] It should be noted that the Bonferroni, Holm and Hochberg methods require that P Multiply the value by , while the Benjamini-Hochberg method requires P Multiply the value by The Benjamini-Yekutieli method requires P Multiply the value by , and N It is considered as the inspection quantity for calculation.

[0030] In some embodiments of the present application, the second threshold value can be selected as needed, and is generally set to 0.05.

[0031] In some implementation schemes of the present application, in step S5, the top k biometric pairs with the largest correlation coefficients are obtained using a small top heap. In the present application, the k value is set as needed, and any integer greater than 0 can be selected theoretically. However, based on the purpose of biometric pair screening, a k value much smaller than d is generally set (such as 10, 100, 1000, etc., while the d value may still be in the tens of thousands). On the one hand, it is due to memory limitations, and on the other hand, it is meaningless to set the k value too large. For example, if k is set to be larger than the size of the correlation matrix, the result can be returned without any sorting algorithm, and the purpose of screening biometric pairs cannot be achieved.

[0032] In other embodiments of the present application, a third threshold may be set to screen out biometric pairs with a correlation coefficient greater than the third threshold, which are generally considered to be biometric pairs with a relatively strong correlation.

[0033] The second aspect of the present application provides a system for analyzing the correlation of biological features in bio-omics data, comprising the following modules: An input module, for obtaining bio-omics data of a plurality of samples, and for further obtaining one or more matrices including biometric data; A block matrix calculation module is connected to the input module and is used to implement the following steps: For any two matrices, if the product of the number of biometric features of the two matrices exceeds a preset threshold, at least one of the two matrices is divided into multiple sub-matrices along the biometric dimension, and one of the matrices or each of its sub-matrices is subjected to matrix calculation with the other matrix or each of its sub-matrices; Repeat the previous step until all matrices are calculated to obtain a correlation matrix, which includes correlation coefficient vectors between the biological features.

[0034] In some embodiments of the present application, in the block matrix calculation module, the correlation coefficient is selected from one of the Pearson correlation coefficient, the Spearman correlation coefficient and the Kendall correlation coefficient.

[0035] In some embodiments of the present application, the correlation matrix also includesP Value vector and / or correction P Value vector, the system further comprises the following modules: A screening module is connected to the block matrix calculation module and is used to P Value or correction P The correlation matrix is ​​screened by the value to obtain the P Value or correction P The number of biometric pairs whose values ​​are less than the second preset threshold is denoted as d; A top-k selection module is connected to the screening module and is used to obtain the first k biometric pairs with the largest correlation coefficient or to obtain the biometric pairs with a correlation coefficient greater than a third preset threshold, wherein 0 <k<d。

[0036] In some embodiments of the present application, the P The value is estimated based on linear interpolation. The specific steps include: obtaining L indivual P The values ​​are sorted from small to large. , and obtain its corresponding square correlation coefficient: ; For any correlation coefficient r , which is approximately P value Calculated using the following formula:

[0037] in,

[0038] In some embodiments of the present application, the correction P The value is also approximated based on linear interpolation. The specific steps include: M The correlation coefficient of biometric pairs P The values ​​are sorted in ascending order as follows: ,…, ; use P The value correction method obtains M indivual P Correction of value P value : ,…, ; For one P Value, its approximate correction P value Calculated using the following formula:

[0039] in, M >10000,

[0040] In some embodiments of the present application, the P The value correction method is selected from at least one of Bonferroni, Holm, Hochberg, Benjamini-Hochberg and Benjamini-Yekutieli.

[0041] The third aspect of the present application provides a computer device, comprising: a memory for storing a computer program; and a processor for implementing the steps of any method described in the first aspect of the present application when executing the computer program.

[0042] The fourth aspect of the present application provides a computer-readable storage medium, on which a computer program is stored. When the computer program is executed by a processor, the steps of any method described in the first aspect of the present application are implemented.

[0043] Beneficial effects of this application Compared with the prior art, this application has the following beneficial effects: Correlation analysis is a ubiquitous statistical method in multi-omics comprehensive analysis. In recent years, the emergence of multi-omics technologies, including high-throughput sequencing and mass spectrometry, has led to the accumulation of a large amount of biomedical data. However, as the number of features and samples in the dataset increases, the computational complexity and memory consumption of correlation analysis grow rapidly. Existing tools often cannot handle the computational requirements of building correlation networks for high-dimensional data.

[0044] The method, system, device and medium of the present application utilize high-performance, parallel block matrix algorithms to process large-scale correlation and clustering tasks. Compared with the prior art, the method, system, device and medium of the present application demonstrate superior performance in calculating correlation matrices, enabling downstream analysis to exceed computer memory resources.

[0045] The method, system, device and medium of the present application can effectively estimate the approximate P Value or approximate correction P The method can evaluate the significance of correlation coefficients 66 times faster than the existing technology without the need for complex cumulative distribution function calculations, which is a significant improvement.

[0046] The method, system, device and medium of the present application implement a top-k algorithm, which uses a partial sorting method based on a small top heap to achieve faster execution and lower memory consumption. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1A schematic diagram of the process of top-k selection using the PyNetCor method is shown. Matrices A and B represent the input data, with rows corresponding to biological features (e.g., genes, species) and columns corresponding to samples. Matrix B can be the same as matrix A. (a) The input matrix is ​​preprocessed and transformed, dividing matrix A into multiple smaller blocks along the feature dimension while keeping matrix B unchanged. (b) Correlations are calculated using one submatrix of A and the full matrix B at a time, and the remaining submatrices are processed sequentially to optimize memory usage. (c) Finally, a small top heap structure is used to store the highest correlation coefficients, and a heap sort is used to obtain the top k correlation coefficients.

[0048] Figure 2 The running time of correlation calculation using the block matrix algorithm is shown. The running time of the unblocked matrix algorithm and the block matrix algorithm for calculating the Pearson correlation matrix for a metagenomic dataset with 50,000 features and 500 samples is compared using 1-16 cores.

[0049] Figure 3 The figure shows the running time comparison of PyNetCor and other tools for calculating the Spearman correlation matrix on the metagenomic dataset. (a) and (b) are on the MHGF dataset, and (c) and (d) are on the EMP dataset; (a) and (c) represent different numbers of features, and (b) and (d) represent different numbers of samples. The dotted line represents the parallel calculation using 32 threads, and the solid line represents the calculation using a single thread.

[0050] Figure 4 The figure shows the running time and memory consumption comparison of PyNetCor and other tools for top-k (top 1%) correlation selection on metagenomic datasets with different feature numbers. (a) and (b) are on the MHGF dataset, and (c) and (d) are on the EMP dataset; (a) and (c) represent the running time, and (b) and (d) represent the memory consumption. The missing lines indicate program crashes due to exceeding the system memory limit. The dotted line represents parallel calculation using 32 threads, and the solid line represents calculation using a single thread.

[0051] Figure 5 shows the estimated approximation in pyNetCor P Value and approximate correction P The performance of pyNetCore. P Comparison of the running time of the value approximation method and the classic method of CorALS on correlation coefficient datasets of different sizes. The parallel computation using 32 threads is shown. pyNetCor is the only one that can parallelize the computation. P (b) Pearson, Spearman, and Kendall correlation calculations by pyNetCorP (c) Absolute error distribution of multiple hypothesis testing correction using approximate methods provided by pyNetCor, including Bonferroni, Benjamini-Hochberg (BH), Benjamini-Yekutieli (BY), Holm and Hochberg. No significant error was observed in the approximate corrected P value obtained by Bonferroni correction. The absolute errors of other methods were all less than 5e-4.

[0052] Figure 6 Network analysis of microbe-metabolite interactions in IBD is shown. (a) Heat map of the difference in Spearman correlation between disease and control samples of omics data (microbiome, metabolome) from the IBDMDB database. The top 1% of microbe-metabolite pairs with the largest correlation difference are shown, and 10 different microbes are identified. Red indicates increased differential correlation under disease conditions compared to the control group, while blue indicates decreased differential correlation. (b) Significant association network of 10 different microbes and metabolites (FDR < 0.05). Rectangular nodes represent microbes and circular nodes represent metabolites. Red edges indicate positive correlations and blue edges indicate negative correlations, with darker colors representing stronger associations. DETAILED DESCRIPTION

[0053] Unless otherwise indicated, implied from the context, or customary in the prior art, all parts and percentages in this application are based on weight, and the tests and characterization methods used are all synchronized with the filing date of this application. Where applicable, the contents of any patent, patent application or disclosure involved in this application are fully incorporated herein by reference, and their equivalent patent families are also introduced as references, especially the definitions of relevant terms in the art disclosed in these documents. If the definition of a specific term disclosed in the prior art is inconsistent with any definition provided in this application, the definition of the term provided in this application shall prevail.

[0054] In order to make the technical problems, technical solutions and beneficial effects solved by the present application more clearly understood, the present application is further described in detail below in conjunction with embodiments.

[0055] The following examples are used to demonstrate the preferred embodiments of the present application. It will be understood by those skilled in the art that the techniques disclosed in the following examples represent techniques discovered by the inventors that can be used to implement the present application, and therefore can be considered as preferred embodiments of the present application. However, it should be understood by those skilled in the art based on this specification that many modifications can be made to the specific embodiments disclosed herein, and still the same or similar results can be obtained without departing from the spirit or scope of the present application.

[0056] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs, and the disclosure and materials cited therein are hereby incorporated by reference.

[0057] Those skilled in the art will recognize, or be able to ascertain using no more than routine experimentation, many technical equivalents to the specific embodiments of the invention described herein. Such equivalents are intended to be encompassed by the claims.

[0058] The experimental methods in the following examples are conventional methods unless otherwise specified. The instruments and equipment used in the following examples are conventional laboratory instruments and equipment unless otherwise specified; the experimental materials used in the following examples are purchased from conventional biochemical reagent stores unless otherwise specified.

[0059] Example 1 Correlation analysis of macrogenomics data When analyzing metagenomics datasets with a large number of features, the traditional methods described in the background technology often require a lot of memory and rely on high-performance computing clusters.

[0060] For example, when calculating the correlation between two features of metagenomic sequencing data, the number of features may reach more than one million, and the calculated correlation matrix is ​​very large, which can easily lead to memory overflow. In theory, calculating the correlation matrix (1.6M×1.6M) of a metagenomic data set with 1.6 million features requires more than 20TB of memory.

[0061] To this end, in this embodiment, the inventors propose a method (PyNetCor) suitable for processing large-scale data (such as metagenomics data) for correlation analysis, which combines a block matrix algorithm, a heap-based partial sorting algorithm, and a linear interpolation-based P The value approximation method achieves efficient correlation analysis while significantly reducing memory consumption.

[0062] To evaluate the performance of PyNetCor, the inventors used two datasets: A metagenomic HUMAnN3 gene family (MHGF) dataset was generated using HUMAnN3 from more than 500 metagenomic samples for the analysis of microbial community functions. The process involved mapping DNA reads to the UniProtReference Clusters database (UniRef90). Due to the large number of zeros in metagenomic data, the inventors excluded features with zero variance or zero values ​​exceeding 95% to mitigate the impact of zero inflation on the data.

[0063] The other is 16S rRNA gene sequencing data from the Earth Microbiome Project (EMP). The EMP dataset includes 96 studies from different environments around the world, including 23,828 samples. Through 16S rRNA gene sequencing and amplicon data analysis, an operational taxonomic unit (OTU) table containing 307,572 features was generated.

[0064] Combination Figure 1 , PyNetCor is described in detail as follows: 1. Block Matrix Algorithm For two matrices A (m a ×n)、B (m b ×n), where m a and m b They refer to the number of biometric features, and n is the number of samples.

[0065] S1, preprocessing and vector transformation of matrices A and B. After transformation, matrix multiplication can be used directly to calculate the correlation; S2, divide the matrix A into multiple smaller blocks along the feature dimension, and the matrix B remains unchanged; S3, applies lazy evaluation technology, multiplying a sub-block of matrix A with the full matrix B each time it is requested; S4, screening significantly correlated features based on correlation results; S5, repeat S1-S4 until all sub-blocks of matrix A are calculated.

[0066] Among them, in step S1, the preprocessing and vector transformation refers to feature centering processing of matrices A and B so that their mean is 0 and their variance is 1 / n. The correlation calculation can be expressed as the scalar product of two preprocessed biometric features. After the conversion, matrix multiplication can be directly used to improve the calculation efficiency.

[0067] In step S2, the size of the smaller block is 500 by default, and can be adjusted to other values ​​by those skilled in the art.

[0068] Using the block matrix approach, correlation calculations for a metagenomic dataset with 1.6 million biological features can be completed in 3.5 hours using only 8 CPUs and 23GB of memory.

[0069] To verify the advantages of the block matrix method, the inventors used 6 high-dimensional data sets downsampled from the MHGF data set, with the number of features ranging from 50,000 to 1.6 million. The Pearson correlation matrix was calculated using the block matrix algorithm on a personal computer equipped with 8 cores and 32GB RAM (as shown in Table 1).

[0070] Table 1 Time and memory consumption for calculating the Pearson correlation matrix using the block matrix algorithm on a personal computer

[0071] As can be seen from Table 1, the block matrix algorithm can build large-scale correlation networks without relying on expensive high-performance computing resources. Specifically, when processing a metagenomics dataset of 1.6 million features, theoretically more than 20TB of memory is required to calculate the complete correlation matrix. However, using a single-core block matrix algorithm, only 22GB of memory is required, and the calculation can be completed within 9 hours; in multi-threaded mode (parallel), the block matrix algorithm can further improve the computational efficiency and shorten the running time to 3.5 hours.

[0072] Since the block matrix algorithm divides a large matrix into multiple small matrices of a specified size, the correlation calculation results of each small matrix are returned only when requested, thereby achieving controllable memory usage. Therefore, in the block matrix algorithm, the running time may increase. In order to evaluate the performance of the block matrix algorithm, the inventors measured its running time under different numbers of threads and compared it with the traditional method without the block matrix algorithm. Figure 2 As shown in Figure 2, the running time of the block matrix algorithm is about 0.6-24.2% higher than that of the traditional method, mainly due to the additional overhead of multiple memory allocations and deallocations. However, considering that the block matrix algorithm can significantly reduce memory consumption, these minor running time losses are negligible.

[0073] To verify the performance of PyNetCor in correlation analysis, the inventors compared PyNetCor with existing tools, including CorALS, Numpy, Pandas, and WGCNA. Since the single-core implementations of PynetCor and CorALS showed better performance than other tools, the inventors further extended the comparison to a multi-threaded environment using 32 cores. Figure 3 Shown are the running times of various tools for computing the Spearman correlation matrix when processing metagenomic datasets with varying numbers of features and samples.

[0074] For the MHGF dataset (70,000 features), only pyNetCor and CorALS can complete the computation within 10 minutes, with PynetCor consistently showing superior performance. In multithreaded mode, pyNetCor only takes 21.1 seconds, 3 times faster than CorALS and 260 times faster than Pandas. It is worth noting that these performance gaps are expected to widen further as the number of samples increases. For example, on the EMP dataset (70,000 features, >20,000 samples), pyNetCor is 4.3 times faster than CorALS and 1816 times faster than Pandas.

[0075] In addition, PyNetCors also supports the calculation of the Pearson correlation matrix. The running time analysis is shown in Table 2.

[0076] Table 2 Comparison of running time of different methods for calculating Pearson correlation matrix on the EMP dataset

[0077] While pyNetCor and NumPy show comparable runtimes on large-sample (>1000 samples) datasets, pyNetCor shows better performance on small-sample (<500 samples) datasets.

[0078] 2. Heap-based partial sorting algorithm Differential correlation analysis is to obtain the changes in the correlation between features in different states (diseases or different time periods, etc.), among which the most significantly correlated features may involve key biological processes. Before calculating their differences, the correlation matrix of the two groups of samples is calculated first.

[0079] Also refer to Figure 1 , based on the block matrix algorithm, further includes searching for the k (top-k) feature pairs with the highest correlation between the features of matrix A and matrix B, using the top-k correlation (differential correlation) search algorithm. Specifically, the above block matrix algorithm further includes: S6, use the small top heap to obtain the top-k feature pairs.

[0080] In high-dimensional data sets, many weak correlations may occur between features due to random noise. Therefore, the inventors prioritize the core top-k most strongly correlated feature pairs. The inventors iterate all correlations based on the aforementioned high-performance, memory-efficient block matrix algorithm, and retain the top k feature pairs with the highest correlation at each round of calculation. Subsequently, the inventors use the heap sort algorithm to sort the top-k correlations. The average time complexity of the top-k search process is ( n logk ) and space complexity (k).

[0081] In order to verify the performance of the heap-based partial sorting algorithm, the inventors first directly used Numpy and WGCNA to implement top-k selection, that is, first calculate the complete correlation matrix, and then obtain the top-k correlations through an efficient sorting algorithm. However, due to the huge memory consumption of this method, the inventors were unable to run the direct implementation based on Numpy and WGCNA on a data set with more than 70,000 features ( Figure 4 ).

[0082] Therefore, pyNetCor provides a method for top-k correlation coefficient selection, which uses block matrix calculation and heap-based partial sorting algorithm, while CorALS extracts approximate top-k correlations through eigenvector projection and ball tree. On the MHGF dataset, pyNetCor only needs 150s running time and 2.3GB of memory to use a single thread ( Figure 4 (a and b) Searching for the top 1% of relevance among 70,000 features is 34 times faster than CorALS and uses only 39% of the maximum memory of CorALS. The same trend is observed on the EMP dataset ( Figure 4 It is worth noting that pyNetCor maintains very small memory consumption even with multi-core parallel acceleration.

[0083] The following results fully demonstrate that pyNetCor is capable of performing efficient correlation analysis on high-dimensional data such as metagenomics.

[0084] 3. Based on linear interpolation P Value Estimation Building large correlation networks usually requires hypothesis testing of correlation coefficients to determine the reliability of the correlation. Existing tools mainly focus on the computational performance of correlation matrix calculations, resulting in low efficiency in P value calculations and inability to perform multiple test corrections. In addition, most methods fail to effectively control the false discovery rate (FDR), resulting in a decrease in statistical power when the number of tests increases, and may even fail to find any significant results. To this end, PyNetCor also provides a method to maintain statistical power and prevent FDR for large-scale data analysis, as detailed below: ①Based on linear interpolation P Value Estimation P The calculation of the value depends on the relationship between the correlation statistic and the cumulative distribution function (CDF). r For example, assuming the sample size is n ,butr Corresponding P The value can be passed t The CDF of the distribution ( pt ) is determined as shown in the following formula:

[0085] in, t Calculated by the following formula: .

[0086] To speed up the calculation of complex CDFs, the inventors constructed a pre-calculated P Value and t Statistics (or related systems r ). This table contains L (>10,000) P value( , ,…,1) and its corresponding square correlation coefficient ( ,…, ). Given a new squared correlation coefficient, its approximation can be estimated by linear interpolation from this table P value( ).

[0087]

[0088] in,

[0089] ② P Multiple hypothesis testing correction This embodiment further uses a similar method as above to optimize P Multiple hypothesis testing correction of the value. To this end, five correction methods are provided, namely: Bonferroni, Holm, Hochberg, Benjamini-Hochberg and Benjamini-Yekutieli. Specifically, randomly select M (greater than 10,000, 500,000 is selected in this example) tests, calculate P values ​​and sort them in ascending order ( … ). Then, use the selected P The value correction method is for each original Calculate the corresponding approximation P value( )for( … ).

[0090] It is worth noting that this approach uses a smaller sample size MTo approximate the overall size N Most methods require P The Bonferroni, Holm, and Hochberg methods require that the P Multiply the value by , while the Benjamini-Hochberg method requires P Multiply the value by The Benjamini-Yekutieli method requires P Multiply the value by , and N Similar to the previous method, linear interpolation can be used to calculate any new given true P Approximate correction of values P value.

[0091]

[0092] in, .

[0093] To compare the approximation based on interpolation P Value Method and Classic P To improve the running time of the proposed method, the inventors processed the Pearson correlation matrix of the MHGF dataset and created 6 new datasets (containing 10, 30, 50, 70, 90 and 110 million correlation coefficients respectively) by downsampling.

[0094] The inventors used derive_pvalues ​​in CorALS as a classic P-value method for comparison. With the help of an efficient parallelized linear interpolation engine, pyNetCor can complete the calculation of 110 million approximate P values ​​in 3.36 seconds, which is 118 times faster than the CorALS implementation that does not support multithreading ( Figure 5 (a).

[0095] In addition to its superior performance, pyNetCor also shows high accuracy. On the same dataset, pyNetCor computes approximations of the Pearson, Spearman, and Kendall correlation coefficients. P Value and truth P The absolute error between values ​​is always below 1e-8 ( Figure 5 b). In addition, the approximate correction P Value and truth P The absolute error between values ​​is less than 5e-4 ( Figure 5 (c).

[0096] Example 2 Multi-omics correlation analysis in IBD Inflammatory bowel disease (IBD) is a chronic, relapsing intestinal inflammatory disease that includes ulcerative colitis (UC) and Crohn's disease (CD). Although its genetic pathology remains elusive, many studies have used multi-omics technologies to investigate the pathogenesis of IBD.

[0097] In order to verify the application of the method of the present application in the correlation analysis of multi-omics data, the inventors analyzed the microbiome and metabolomics data from the IBDMDB database. The data set contains 181 CD samples, 102 UC samples, and 105 non-IBD samples. Using non-IBD samples as controls, the inventors used the method of the present application to obtain the Spearman correlation between the microorganisms / metabolites in the diseased samples and the control group.

[0098] Furthermore, the inventors extracted the top 1% of microorganisms and metabolites with the largest correlation changes for visualization ( Figure 6 Figure (a) shows that 10 microbial species were significantly altered in IBD samples. Roseburia inulinivorans and Faecalibacterium prausnitzii It produces butyrate, a short-chain fatty acid with anti-inflammatory properties. Butyrate plays an important role in maintaining intestinal barrier function and inhibiting the development of IBD.

[0099] To identify the metabolites most closely associated with the 10 differential microorganisms, the inventors calculated the Spearman correlation between these microorganisms and all metabolites ( Figure 6 The results showed that there were 103 pairs of microorganisms and metabolites with significant correlation (FDR < 0.05). Figure 5 As shown in b, more than 15 acylcarnitines Roseburia hominis There is a significant correlation. In particular, it is worth noting that C18:1-OH carnitine, C16-OH carnitine, C14:1 carnitine and Roseburia hominis The negative correlation between them is the most prominent.

[0100] It is worth noting that previous studies have found that Roseburia hominis This means that acylcarnitine may be a key metabolite of inflammatory bowel disease and has the potential to become a new therapeutic target for the disease, providing important clues for further research.

[0101] All documents mentioned in this application are cited as references in this application, just as each document is cited as reference individually. In addition, it should be understood that after reading the above teaching content of this application, those skilled in the art can make various changes or modifications to this application, and these equivalent forms also fall within the scope defined by the claims attached to this application.

Claims

1. A method for analyzing the correlation of biological features in bio-omics data, characterized in that: The following steps are involved: S1, obtaining one or more bio-omics data of multiple samples, and obtaining one or more matrices including biological feature quantitative data; S2, for any two matrices, if the product of the number of biometric features of the two matrices exceeds a first preset threshold, at least one of the two matrices is divided into a plurality of sub-matrices along the biometric feature dimension, and one of the matrices or each of its sub-matrices is subjected to matrix calculation respectively with the other matrix or each of its sub-matrices; S3, repeat step S2 until all matrices are calculated to obtain a correlation matrix, wherein the correlation matrix includes correlation coefficient vectors between each pair of the biometric features.

2. The method according to claim 1, characterized in that: In step S3, the correlation coefficient is selected from one of the Pearson correlation coefficient, the Spearman correlation coefficient and the Kendall correlation coefficient.

3. The method according to claim 1 or 2, characterized in that: The correlation matrix also includes P Value vector and / or correction P The method further comprises the following steps: S4, according to P Value or correction P The correlation matrix is ​​screened by the value to obtain the P Value or correction P The number of biometric pairs whose values ​​are less than the second preset threshold is denoted as d; S5, obtaining the first k biometric pairs with the largest correlation coefficient or obtaining the biometric pairs with a correlation coefficient greater than a third preset threshold, wherein 0 <k<d。 4. The method according to claim 3, characterized in that Said P The value is approximated based on linear interpolation. The specific steps include: Pre-acquisition L indivual P The values ​​are sorted from small to large. , and obtain its corresponding square correlation coefficient: ; For any correlation coefficient r , which is approximately P value Calculated using the following formula: in, .

5. The method according to claim 4, characterized in that The correction P The value is approximated based on linear interpolation. The specific steps include: M The correlation coefficient of the biometric pairs P The values ​​are sorted in ascending order as follows: ,…, ; use P The value correction method obtains M indivual P Correction of value P value : ,…, ; For one P Value, its approximate correction P value Calculated using the following formula: Among them, M>10000, .

6. The method according to claim 5, characterized in that Said P The value correction method is selected from at least one of Bonferroni, Holm, Hochberg, Benjamini-Hochberg and Benjamini-Yekutieli.

7. A system for analyzing the correlation of biological features in bio-omics data, characterized in that: Includes the following modules: An input module, for obtaining one or more bio-omics data of a plurality of samples, and for obtaining one or more matrices including quantitative data of biological characteristics; The block matrix calculation module is connected to the input module and is used to implement the following steps: For any two matrices, if the product of the number of biometric features of the two matrices exceeds a first preset threshold, at least one of the two matrices is divided into a plurality of sub-matrices along the biometric feature dimension, and one of the matrices or each of its sub-matrices is subjected to matrix calculation respectively with the other matrix or each of its sub-matrices; Repeat the previous step until all matrices are calculated to obtain a correlation matrix, which includes correlation coefficient vectors between the biological features.

8. The system according to claim 7, characterized in that The correlation matrix also includes P Value vector and / or correction P Value vector, the system further comprises the following modules: A screening module is connected to the block matrix calculation module and is used to P Value or correction P The correlation matrix is ​​screened by the value to obtain the P Value or correction P The number of biometric pairs whose values ​​are less than the second preset threshold is denoted as d; A top-k selection module is connected to the screening module and is used to obtain the first k biometric pairs with the largest correlation coefficient or to obtain the biometric pairs with a correlation coefficient greater than a third preset threshold, wherein 0 <k<d。 9. A computer device, characterized in that: include: Memory for storing computer programs; A processor, configured to implement the steps of the method according to any one of claims 1 to 6 when executing the computer program.

10. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the steps of the method according to any one of claims 1 to 6 are implemented.

Images