Calculation acceleration method and system for microbial metagenome function prediction and application
Through the method of sparse matrix storage and calculation, combined with the optimization processing of the armadillo library and Perl language, the computational efficiency problem of predicting the sequence function of microbial metagenomic marker genes is solved, and more efficient data processing and less memory usage are achieved.
Patent Information
- Application Number
- CN202410060875.4
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2024-01-16
- Publication Date
- 2025-07-18
AI Technical Summary
The existing microbial metagenomic marker gene sequence function prediction methods are ineffective in the face of the surge in data volume caused by the rapid development of second-generation sequencing technology, which easily leads to memory overflow and cannot complete functional prediction within a reasonable time range.
The method of sparse matrix storage and calculation is adopted, and the feature-function unit mapping matrix is read and calculated through the armadillo library, the path abundance is inferred in combination with the minpath algorithm, and the description information is added line by line using Perl language to optimize the generation process of samples, features, and functional unit mapping tables.
It significantly reduces computing time and memory usage, improves computing efficiency, can process larger data volumes of microbial metagenomic data, reduces the risk of program crashes, and simplifies code maintenance.
Smart Images
Figure HDA0004666442470000011 
Figure HDA0004666442470000021
Abstract
Description
Technical Field
[0001] The present invention belongs to the technical field of microbial metagenomic data processing, and relates to a method, system and application for accelerating the prediction of microbial metagenomic functions based on marker genes. Background Art
[0002] With the development of next-generation sequencing (NGS) technology and the continuous innovation of technology, the cost of gene sequencing has further declined. In 2021, the average cost of gene sequencing per megabyte of data was only $0.006. At the same time, the data of gene sequencing is showing explosive growth. Taking Illumina NovaSeq 6000 as an example, a single run can produce outputs of up to 6TB and 20 billion reads. In recent years, the growth rate of genomic sequencing data worldwide has exceeded Moore's law. Between 2008 and 2016, the amount of gene sequencing data worldwide doubled every seven months. How to quickly process the ultra-high-speed growing data faces severe challenges.
[0003] Amplicon sequencing is a highly targeted method for analyzing gene variations in specific genomic regions. Ultra-deep sequencing of PCR products (amplicons) can effectively identify variations and characterize them. The general idea is to target and capture the target region, then perform next-generation sequencing (NGS), analyze the sequencing results, and obtain corresponding information. It mainly includes 16S rDNA sequencing, 18S rDNA sequencing, ITS sequencing, and target region amplicon sequencing, etc. Sequences of a certain hypervariable region of 16S / 18S / ITS measured by the second-generation high-throughput sequencing platform are used to reflect the differences between species in the classification of bacteria, fungi, and archaea in environmental samples, which has important guiding significance for studying the microbial composition in environments such as the ocean, soil, and intestinal feces. At the same time, it is also a widely used method in phylogenetic and taxonomic research, especially applied to different microbial metagenomics samples.
[0004] The results obtained by using amplicon sequencing method, that is, when targeting and detecting a certain marker gene sequence, are more likely to be affected by SNVs caused by sequencing errors, resulting in incorrect sequence classification, and ultimately detecting similar but incorrect microorganisms, or mistakenly thinking that new microorganisms have been discovered. To address this problem of amplicon sequencing, there are currently two commonly used analysis strategies - OTU and ASV - to reduce the impact caused by sequencing errors. OTU is essentially a clustering method based on its own data for clustering. And ASV simply means that after removing the incorrect sequences, the Identity standard is set to 100% for clustering. In recent years, more and more articles have started to adopt ASV and abandon OTU.
[0005] One limitation of microbial community marker gene sequencing is that it does not provide information on the functional composition of the sampled community. Currently, a variety of tools have been developed to predict the functional potential of bacterial communities based on marker gene sequencing profiles. Commonly used tools for functional prediction of marker gene profiles include PICRUSt, PICRUSt2, Tax4Fun, etc. Among them, the PICRUSt2 software is the most widely used in marker gene function prediction. These mature functional prediction tools for marker gene sequences share a common feature. Their basic idea is to obtain the mapping relationship between feature sequences (OTU / ASV) and functions through the mapping relationship between feature sequences (OTU / ASV) and species and the mapping relationship between species and functions. With the rapid development of next-generation sequencing technology, the amount of sequences obtained in a single experiment is larger. At the same time, amplicon denoising methods generate more feature sequences compared to traditional taxonomic unit clustering methods. Facing the increasing computational requirements for function prediction, users hope that the function prediction tool has higher running efficiency so that the function prediction of marker gene sequences can be completed within a reasonable time range. The current solution of PICRUSt2 to this problem is to use the Python language to call the pandas cache matrix and perform multi-process parallel computing by slicing the matrix. Although the efficiency of block parallelism is higher than that of single-threading, it still needs to read all the data into the cache at once and then slice it. When facing huge amounts of data, this process is slow and requires a large amount of memory, and the computing server is prone to memory overflow, resulting in the failure of the entire task. Summary of the Invention
[0006] In order to solve the deficiencies of the existing technology, the purpose of the present invention is to provide a computational acceleration method for predicting the functions of microbial metagenomic marker gene sequences, predicting marker gene functions, gene families, reading sequence abundance tables, and finally outputting predicted function abundance files for each sample. This method is based on sparse matrices to improve the problem of the long time-consuming of current gene sequence function prediction and lay a solid foundation for subsequent data analysis. The present invention adopts a simpler way to store and calculate the sample-feature matrix and the feature-function unit mapping matrix, which can not only make good use of the characteristics of sparse matrices, but also has simpler code maintenance and better performance compared with the implementation method of pandas. In addition, the present invention optimizes the steps of generating sample, feature, and function unit mapping tables, greatly reducing the computational time compared with the original PICRUSt2. The present invention also uses some tools used by the original PICRUSt2 itself for sequence alignment, phylogenetic tree construction, prediction, etc. These tools can greatly simplify the operation, lower the usage threshold, and make function prediction easier to implement.
[0007] To achieve the above object, the technical solution provided by the present invention is: A computational acceleration method for predicting the functions of microbial metagenomic marker gene sequences, as Figure 2As shown, it includes the following steps:
[0008] 1) Preprocess the characteristic abundance data and characteristic sequences of microbial metagenomic marker genes to obtain the mapping relationship between characteristic sequences and functional units, that is, the characteristic-functional unit mapping matrix; if the characteristic sequence is the 16S rRNA gene, correct the characteristic abundance data to obtain the copy number of the 16S rRNA gene of the species corresponding to the characteristic sequence and the characteristic abundance matrix corrected by the 16S rRNA copy number; in the present invention, since the 16S rRNA characteristic sequences obtained by sequencing come from different species, and the copy numbers of the 16S rRNA genes of different species are not consistent. Therefore, the number of 16S rRNA sequences obtained by sequencing cannot directly reflect the proportion of species composition in the community (the number of sequences obtained by sequencing is equivalent to the result of the actual species individuals in the community weighted by the copy number), and subsequent analysis needs to eliminate the influence of the 16S rRNA copy numbers of different species through 16S rRNA copy number correction;
[0009] 2) Read and calculate the characteristic abundance matrix corrected by the 16S rRNA copy number obtained in step 1) or the externally input characteristic abundance matrix and the characteristic-functional unit mapping matrix obtained in step 1) through the armadillo library to obtain the functional unit abundance matrix and the characteristic-functional unit abundance matrix;
[0010] 3) For the functional unit abundance matrix obtained in step 2), infer the pathway abundance through minpath and calculate the pathway abundance matrix;
[0011] 4) Add description information to the functional unit abundance matrix, the characteristic-functional unit abundance matrix obtained in step 2), and the pathway abundance matrix obtained in step 3);
[0012] 5) Output the predicted functional abundance file for each sample.
[0013] In step 1), the preprocessing of the characteristic abundance data and characteristic sequences of microbial metagenomic marker genes includes reading the characteristic sequences and characteristic abundance data, saving the characteristic abundance data in text form, and establishing the mapping relationship between characteristic sequences and functional units;
[0014] The establishment of the mapping relationship between characteristic sequences and functional units includes the following steps (consistent with the Picrust2 method):
[0015] 1.1) For the characteristic sequences of the functions to be predicted, align them to the reference marker gene sequences through hmmalign; hmmalign is a multiple sequence alignment tool using the hidden Markov model (HMM), which is used to perform HMM alignment of all characteristic sequences with a pre-set configuration file (the configuration file contains reference sequences) respectively, and output the alignment results;
[0016] 1.2) For the aligned feature sequence to be predicted and the reference sequence in the preset configuration file, obtain the evolutionary distance between the feature sequence to be predicted and the reference sequence (reference marker gene sequence) through epa-ng, construct an evolutionary tree, and save it in the jplace format; epa-ng is a completely rewritten program of the evolutionary placement algorithm (EPA) implemented in RAxML, used to perform maximum likelihood-based phylogenetic placement of genetic sequences on the pre-constructed reference tree provided and the alignment results obtained in step 1; the evolutionary distance refers to the differences between two sequences and can be used to measure their similarity; the reference tree is constructed based on the 16S rRNA of a large number of microorganisms with known genomes. Since the genomes of these microorganisms are known, the gene distribution on these genomes and the potential metabolic functions of these genes are also known. The purpose of Picrust2 is to predict the metabolic functions of microorganisms / species of genomes of 16S rRNA from unknown sources. The basic principle is that species with closer genetic relationships have more similar genetic information in evolution, so they usually share more metabolic pathways and genes. This results in greater similarity in their metabolic functions. By calculating the similarity between the 16S rRNA sequence from an unknown source and the 16S rRNA sequences in the reference tree, the 16S rRNA sequence from an unknown source is "inserted" into the reference tree, and its genetic relationship with known genomes is inferred based on its position in the tree. Finally, the functional potential of the 16S rRNA from an unknown source is predicted using the gene distribution of known genomes.
[0017] 1.3) For the evolutionary tree file in the jplace format, convert it to the newick format through gappa; gappa is a set of toolkits for processing phylogenetic data, used to analyze and compare different jplace files, edit, operate on, and convert phylogenetic data files in different formats, etc.;
[0018] 1.4) For the newick evolutionary tree containing the evolutionary relationship between the feature sequence to be predicted and the reference sequence and the functional unit composition matrix obtained from the genomes corresponding to the reference sequence, perform hidden state prediction through castor, calculate the functional unit composition of each feature sequence to be predicted, obtain the feature-functional unit mapping matrix, and save it in text form; castor is an R language package for efficient algorithms for large-scale phylogenetic analysis, based on the classical ancestral state reconstruction algorithm, used for feature reconstruction of phylogenetic branches with unknown metabolic diversity and phenotypes;
[0019] 1.5) Calculate the distance NSTI between each feature sequence to be predicted and the reference sequence with the closest evolutionary distance through castor. NSTI (Nearest Sequenced Taxon Index) is an index in Picrust2 used to evaluate the matching degree or similarity of 16S rRNA sequence data of the microbiome on the phylogenetic tree. Specifically, NSTI represents the similarity between the 16S rRNA sequences in a given microbial sample and the sequenced species on the phylogenetic tree of known bacteria and archaea. The lower the NSTI value, the closer the position of the 16S rRNA sequences in the sample is to that of the known species on the phylogenetic tree, so the functional prediction is usually more reliable. The higher the NSTI value, the farther the position of the 16S rRNA sequences in the sample is from that of the known species on the phylogenetic tree, so the reliability of the functional prediction may be reduced.
[0020] In step 2), read the feature abundance matrix corrected by 16S rRNA copy number obtained in step 1) or the externally input feature abundance matrix, and the feature-functional unit mapping matrix saved in text form obtained in step 1.4); filter the feature sequences to be predicted whose distance (NSTI) between each feature sequence to be predicted and the reference sequence with the closest evolutionary distance is greater than the set threshold (refer to the default value of Picrust2, set to 2); perform multiplication calculation on the feature abundance matrix and the feature-functional unit mapping matrix to obtain the functional unit abundance matrix, and obtain the feature-functional unit abundance matrix through an iterator. It includes the following steps:
[0021] 2.1) Matrix reading: According to the characteristics of the feature abundance matrix and the feature-functional unit mapping matrix, read the matrix data by using the integer sparse matrix class through the armdillo library to obtain the feature abundance sparse matrix and the feature-functional unit sparse matrix respectively. A sparse matrix is a matrix with most elements being zero, used to store very large matrices; the sparse matrix class uses a hybrid storage framework, automatically and seamlessly switching between three data storage formats according to the format most suitable for a specific operation: (i) compressed sparse column, used for efficient basic arithmetic operations such as matrix multiplication and addition, and efficient reading of individual elements; (ii) coordinate list, used to facilitate operations involving batch coordinate conversion; (iii) red-black tree, used for robust and efficient incremental construction of sparse matrices (i.e., constructed by setting individual elements one by one). To further improve the execution efficiency, this class uses C++ features such as template metaprogramming to provide a compile-time expression evaluator, which can automatically detect and optimize common mathematical expression patterns;
[0022] 2.2) Iterate through the non-zero values in the feature abundance matrix using the iterator of the sparse matrix class, and quickly filter out the to-be-predicted feature sequences with an overly large evolutionary distance from the nearest reference sequence, which are considered to have inaccurate prediction results, to obtain a corrected feature-functional unit sparse matrix; the filtering operation can basically eliminate the "garbage" sequences that cannot be located in the reference tree, and these sequences are usually off-target sequences or cannot be classified at the phylum level. Use the remaining sequences for the subsequent steps;
[0023] 2.3) Through the matrix multiplication operation function of the matrix class, perform matrix multiplication on the feature abundance sparse matrix in step 2.1) and the corrected feature-functional unit sparse matrix in step 2.2) to obtain a functional unit abundance matrix; iterate through the non-zero rows in the functional unit abundance matrix using the iterator of the sparse matrix class to obtain a feature-functional unit abundance matrix; through this method, the functional composition of different samples can be predicted from the sequencing results and used in downstream bioinformatics analysis.
[0024] In step 3), for the functional unit abundance matrix obtained in step 2), infer the existing functional pathways through minpath and calculate the pathway abundance matrix. (Consistent with the Picrust2 method)
[0025] First, obtain the mapping relationship between gene families (functional units) and pathways from known databases such as KEGG. This mapping is usually a many-to-many relationship, that is, a gene family may belong to multiple pathways, and a pathway may contain multiple gene families. Then, conservatively reconstruct the pathways through minpath, select the smallest set of pathways that can cover all gene families in the functional unit abundance matrix (this set of pathways is sufficient to explain the existence of all gene families in the dataset), and calculate the pathway abundance matrix; this conservative estimate provides a more reliable estimate of the functional diversity of the samples.
[0026] In step 4), for the functional unit abundance matrix, feature-functional unit abundance matrix obtained in step 2), and the pathway abundance matrix obtained in step 3), read the content of the file line by line using the perl language, store each line of data in a variable, add descriptive text to the functional units and pathways of each line of data, and output it to the file line by line. Read the matrix data and the description files of the functional units and pathways, and add the descriptions of the functional units and pathways to each abundance matrix.
[0027] The present invention also provides a system for implementing the above method for accelerating the functional prediction calculation of microbial metagenomes, and the system includes: a feature sequence processing module, a feature abundance data processing module, a sparse matrix calculation module, a functional unit abundance matrix inference module, and a description information adding module; wherein,
[0028] Feature sequence processing module: It is used to process the feature sequences of the to-be-predicted functions of microbial metagenomic marker genes, including steps such as aligning to reference marker gene sequences, calculating evolutionary distances, and constructing phylogenetic trees.
[0029] Feature abundance data processing module: It is used to read the feature abundance data of microbial metagenomic marker genes, save it in text form, and establish the mapping relationship between feature sequences and functional units.
[0030] Sparse matrix calculation module: It uses the armadillo library to read and calculate the feature abundance matrix corrected by 16S rRNA copy number or externally input and the feature-functional unit mapping matrix to obtain the functional unit abundance matrix and the feature-functional unit abundance matrix.
[0031] Functional unit abundance matrix inference module: For the obtained functional unit abundance matrix, it infers the pathway abundance through the minpath algorithm and calculates to obtain the pathway abundance matrix.
[0032] Description information adding module: It adds description information to the obtained functional unit abundance matrix, feature-functional unit abundance matrix, and pathway abundance matrix, and outputs the predicted functional abundance file for each sample.
[0033] The joint operation of these modules realizes the calculation acceleration method of the present invention.
[0034] The present invention also provides an application of the above method in the research field of microbial community composition, such as environmental microbiology, medical microbiology, food microbiology, etc. Through this method, the functions and gene families of marker genes in microbial samples can be predicted, and the predicted functional abundance file for each sample can be output, providing a basis for subsequent data analysis. This method stores and calculates the sample-feature matrix and the feature-functional unit mapping matrix in a simpler way, with simple code maintenance and better performance. In addition, this method optimizes the steps of generating the sample, feature, and functional unit mapping table, greatly reducing the calculation time-consuming.
[0035] Compared with the prior art, the present invention has the following advantages and beneficial effects:
[0036] The method of the present invention uses the iterator and operation functions of the sparse matrix class in Armadillo to calculate and generate the predicted functional abundance matrix of the feature marker genes. At the same time, through the Perl language, it realizes the method of reading the matrix line by line and adding function and pathway descriptions. Compared with the existing Picrust2 program, the present invention can process larger amounts of data when the memory of the computing server is limited; compared with the implementation methods such as matrix slicing and multi-process parallelism, the scheduling and fault tolerance mechanism of the present invention are easier to develop and maintain; compared with the existing method of using the Python language to call pandas to read and calculate matrices, the implementation of the present invention is simpler and the code is easier to maintain; compared with the implementation of the existing Picrust2 program, the performance of the present invention is higher, and the use of fewer processes, computing time, and memory usage are all significantly reduced. Native Picrust2 uses the pandas extension library of the Python language to operate on matrix data (modifying, adding, deleting rows and columns, and numerical operations between matrices). With the development of sequencing technology and analysis technology (the increase in sequencing volume, and more features are obtained by the ASV method compared to the OTU method), the number of feature sequences for a single functional prediction analysis is increasing. At the same time, the number of functional units in each publicly available functional database is also increasing. For large-scale data, Picrust2 may generate vectors or matrices with lengths in the hundreds of millions or billions when calculating intermediate data using Pandas, and read and write on the hard disk. The present invention rewrites this function using the Armadillo library in C++. Compared with pandas, the advantage of Armadillo is that it is very fast when performing large-scale sparse matrix multiplication calculations, using efficient algorithms and data structures and less memory.
[0037] At the same time, the method of the present invention also optimizes the function of adding descriptions to large matrices. Compared with the existing method of reading the matrix into the cache at one time, the method of reading and processing line by line is adopted, which greatly reduces the time consumption and memory usage of the description adding step. When modifying and adding information to the prediction results, pandas needs to read the huge prediction output results into the memory at one time. Reading into the memory at one time will cause the problem of insufficient memory, while processing the text line by line only occupies a small amount of memory when reading each line. The present invention uses the Perl language to read and process the text line by line, which can reduce the memory occupancy; at the same time, it can respond to the user's request faster, without waiting for the entire file to be read before starting to process, reducing the user's waiting time; it can also reduce the risk of program crashes, because it can avoid memory leaks and other errors that may occur when reading the entire file into the memory. If the program crashes, only the currently processed line will be lost, without affecting the processing of the entire file. Brief Description of the Drawings
[0038] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the accompanying drawings required for the description of the embodiments or the prior art. Obviously, the accompanying drawings in the following description are only some embodiments of the present invention. For those skilled in the art, without creative efforts, other drawings can be obtained based on these drawings.
[0039] Figure 1 is the step flow chart in the embodiment of the present invention.
[0040] Figure 2 is the flow chart of the present invention. Detailed implementation manners
[0041] Combined with the following specific embodiments and the accompanying drawings, the present invention will be further described in detail. The processes, conditions, experimental methods, etc. for implementing the present invention, except for the specifically mentioned content below, are all common knowledge and well-known common sense in the art, and the present invention has no special restricted content.
[0042] The following further illustrates the present invention with specific embodiments.
[0043] The computational acceleration method for predicting the functions of microbial metagenomic marker gene sequences provided in this embodiment, as Figure 1 shown, includes the following steps:
[0044] 1. Preprocess the feature abundance biom file and the feature sequence fasta file obtained from the marker gene sequencing off-machine data through methods such as clustering and denoising; in the present invention, through clustering, denoising and other processing, sequences with biological significance can be obtained, sequences that can represent actual species taxa can be obtained, and redundant sequences can be excluded.
[0045] Armadillo is a high-quality linear algebra library for the C++ language, providing syntax and functions similar to Matlab; used for directly developing algorithms in C++, providing efficient classes for vectors and matrices; can be used in machine learning, pattern recognition, computer vision, signal processing, bioinformatics, statistics, finance, etc.
[0046] Pandas is an extension library of the Python language and is a fast, flexible and easy-to-use data analysis and operation tool for the Python language.
[0047] DataFrame is a two-dimensional labeled data structure and is the main data structure and the most commonly used object of pandas.
[0048] MinPath is a parsimonious method for reconstructing biological pathways using gene family prediction, which realizes a more conservative and reliable estimation of the biological pathways of the query dataset.
[0049] Step 1 includes the following sub - steps:
[0050] 1.1. Read data: Read two files, namely the clustered / denoised feature.biom and rep.fasta, from the local file system;
[0051] 1.2. Use biom convert to convert the biom - formatted abundance data into a tab - delimited feature.txt text file and store it in the local file system;
[0052] 1.3. Generate a feature - function unit mapping matrix fun_predict.txt for the genome corresponding to rep.fasta. The process is as follows:
[0053] 1.3.1. For the rep.fasta file, use hmmalign for multiple sequence alignment to the reference sequence ref.fasta so that the bases without mutations, insertions, or deletions between pairwise sequences are in the same position, obtaining align.fasta;
[0054] 1.3.2. For align.fasta, calculate the evolutionary distance representing the differences in aligned bases between pairwise sequences through epa - ng to obtain an evolutionary tree epa.jplace, and then convert it to epa.newick through gappa;
[0055] 1.3.3. For epa.newick and the functional mapping matrix fun_ref.txt of the reference sequence, predict the functional mapping of all unknown nodes in epa.newick through castor to generate a feature - function unit mapping matrix fun_predict.txt corresponding to rep.fasta.
[0056] 2. For the feature abundance matrix feature.txt and the feature - function unit mapping matrix fun_predict.txt, use the armadillo sparse matrix class to read, fuse data, and perform matrix multiplication. Use a non - zero value iterator to operate on the sparse matrix to reduce calculation time, obtaining a functional unit abundance matrix fun_unstrat.txt, which includes a feature - function unit abundance matrix fun_strat.txt. Figure 1 In this context, sp_fmat represents using the floating - point sparse matrix class in the armadillo library for matrix caching and calculation.
[0057] 3. For the functional unit abundance matrix fun_unstrat.txt, the total number of functions is a finite number. The performance of pandas is close to that of armadillo. Use the DataFrame method of pandas to read the data and extract all the predicted functional unit sets. Based on the known function-pathway mapping relationship, calculate the minimum set of pathways mapped by the function set through minpath as the set of predicted pathways, and sum the functional abundances based on the mapping relationship to obtain the pathway abundance matrix pathway.txt.
[0058] DataFrame uses numpy as its underlying data container, which is a collection of one-dimensional numpy arrays with different data types, along with two indexes (one for rows and one for columns), providing many functions, but all these functions come at the cost of performance. And in the native Picrust2, the apply function provided by pandas is used for matrix multiplication calculation, which essentially iterates over rows or columns and applies the multiplication function without skipping elements with a value of 0.
[0059] 4. Further, for the various abundance matrices fun_unstrat.txt, fun_strat.txt, and pathway.txt obtained in the above steps, use the perl language to establish a hash table for the known function-description mapping, read each function and pathway of the abundance matrix line by line, add the description information, and then output line by line. This is to achieve the purpose of reducing calculation time and memory usage.
[0060] Under the condition of using the same input file, the output result of the present invention is exactly the same as that of Picrust2, but the calculation time is shorter and the memory usage is less. When the data volume is huge, the effect is more obvious, more computing resources are saved, more analysis tasks can be completed per unit time, and the analysis cycle is compressed.
[0061] Example
[0062] The test computing server has a CPU model of Intel(R) Xeon(R) Gold 5218 CPU @ 2.30 GHz and a maximum available memory of 251 GB.
[0063] Use the sequence abundance table (feature.biom) of 42 samples and 8211 feature sequences (a matrix of size 42×8211) as the input data:
[0064] Output the feature sequence-functional unit stratification table (fun_strat.txt):
[0065] ● The CPU running time of the method of the present invention is 255.793 s, and the maximum virtual memory used during operation is 2.493 GB.
[0066] Add a description to the feature sequence - functional unit stratification table (fun_strat_desc.txt):
[0067] · The CPU running time of the method of the present invention is 170.349 s, and the maximum virtual memory used during operation is 10.594 MB;
[0068] Use the sequence abundance table (feature.biom) with 152 samples and 52,228 feature sequences (a matrix of size 152×52,228) as input data:
[0069] Output the feature sequence - functional unit stratification table (fun_strat.txt):
[0070] · The CPU running time of the method of the present invention is 1736.505 s, and the maximum virtual memory used during operation is 12.160 GB;
[0071] Add a description to the feature sequence - functional unit stratification table (fun_strat_desc.txt):
[0072] · The CPU running time of the method of the present invention is 170.349 s, and the maximum virtual memory used during operation is 10.594 MB.
[0073] Comparative example
[0074] The test computing server has a CPU model of Intel(R) Xeon(R) Gold 5218 CPU @ 2.30 GHz and a maximum available memory of 251 GB.
[0075] Use the sequence abundance table (feature.biom) with 42 samples and 8,211 feature sequences (a matrix of size 42×8,211) as input data:
[0076] Output the feature sequence - functional unit stratification table (fun_strat.txt):
[0077] · The CPU running time of native Picrust2 is 523.167 s, and the maximum virtual memory used during operation is 52.505 G;
[0078] Add a description to the feature sequence - functional unit stratification table (fun_strat_desc.txt):
[0079] · The CPU running time of native Picrust2 is 538.646 s, and the maximum virtual memory used during operation is 10.368 GB;
[0080] Use the sequence abundance table (feature.biom) with 152 samples and 52,228 feature sequences (a matrix of size 152×52,228) as input data:
[0081] Output the feature sequence - functional unit stratification table (fun_strat.txt):
[0082] · The native Picrust2 running memory exceeded the maximum available memory of the test server (251GB) and could not be calculated;
[0083] Add a description to the feature sequence - functional unit stratification table (fun_strat_desc.txt):
[0084] · The CPU running time of the native Picrust2 was 538.646s, and the maximum virtual memory used during operation was 10.368GB.
[0085] The protected content of the present invention is not limited to the above embodiments. Without departing from the spirit and scope of the inventive concept, changes and advantages that can be conceived by those skilled in the art are included in the present invention, and the scope of protection is defined by the appended claims.
Claims
1. A computational acceleration method for microbial metagenomic function prediction, characterized in that The method includes the following steps: Step 1: Preprocess the characteristic abundance data and characteristic sequences of microbial metagenomic marker genes to obtain the mapping relationship between characteristic sequences and functional units, i.e., the characteristic-functional unit mapping matrix; Step 2: Read and calculate the characteristic abundance matrix and the characteristic-functional unit mapping matrix obtained in Step 1 through the armadillo library to obtain the functional unit abundance matrix and the characteristic-functional unit abundance matrix; Step 3: For the functional unit abundance matrix obtained in Step 2, infer the pathway abundance through minpath and calculate to obtain the pathway abundance matrix; Step 4: Add description information to the functional unit abundance matrix, the characteristic-functional unit abundance matrix obtained in Step 2, and the pathway abundance matrix obtained in Step 3; Step 5: Output the predicted functional abundance file for each sample.
2. The method according to claim 1, wherein In Step 1, if the characteristic sequence is the 16S rRNA gene, correct the characteristic abundance data to obtain the copy number of the 16S rRNA gene of the species corresponding to the characteristic sequence, and eliminate the influence of the 16S rRNA copy numbers of different species.
3. The method according to claim 1, characterized in that, In Step 1, the preprocessing includes reading the characteristic sequences and characteristic abundance data, saving the characteristic abundance data in text form, and establishing the mapping relationship between the characteristic sequences and functional units.
4. The method according to claim 1, wherein In Step 1, the establishment of the mapping relationship between the characteristic sequences and functional units further includes: Step 1.1: For the characteristic sequences of the functions to be predicted, perform hidden Markov model alignment of all characteristic sequences with the reference sequences in the preset configuration file respectively, align the characteristic sequences of the functions to be predicted to the reference sequences, and output the alignment results; Step 1.2: For the aligned characteristic sequences to be predicted and the reference sequences in the preset configuration file, obtain the evolutionary distance between the characteristic sequences to be predicted and the reference sequences through epa-ng, construct an evolutionary tree, and save it in the jplace format; Step 1.3: Convert the evolutionary tree file in the jplace format to an evolutionary tree file in the newick format through gappa; Step 1.4: For the newick evolutionary tree containing the evolutionary relationship between the characteristic sequences to be predicted and the reference sequences and the functional unit composition matrix obtained from the genomes corresponding to the reference sequences, perform hidden state prediction through castor, calculate the functional unit composition of each characteristic sequence to be predicted, obtain the characteristic-functional unit mapping matrix, and save it in text form; Step 1.5: Calculate the distance between each characteristic sequence to be predicted and the reference sequence with the closest evolutionary distance through castor.
5. The method according to claim 1, characterized in that, In Step 2, the characteristic abundance matrix is the characteristic abundance matrix corrected by the 16S rRNA copy number or externally input; Read the characteristic abundance matrix and the characteristic-functional unit mapping matrix obtained in Step 1, filter out the characteristic sequences to be predicted whose distance between each characteristic sequence to be predicted and the reference sequence with the closest evolutionary distance is greater than the set threshold; perform multiplication calculation on the characteristic abundance matrix and the characteristic-functional unit mapping matrix to obtain the functional unit abundance matrix, and obtain the characteristic-functional unit abundance matrix through iteration of the iterator.
6. The method according to claim 5, wherein The said Step 2 further includes: Step 2.1: According to the characteristics of the feature abundance matrix and the feature-functional unit mapping matrix, read the matrix data by using the integer sparse matrix class through the armdillo library; Step 2.2: Iterate through the non-zero values in the feature abundance matrix by the iterator of the sparse matrix class, filter the to-be-predicted feature sequences whose evolutionary distance from the reference sequence is greater than the set threshold and are considered to have inaccurate prediction results, and obtain the corrected feature-functional unit sparse matrix; Step 2.3: Through the matrix multiplication operation function of the matrix class, perform matrix multiplication on the feature abundance matrix in Step 2.1 and the corrected feature-functional unit mapping matrix in Step 2.2 to obtain the functional unit abundance matrix; iterate through the non-zero rows in the functional unit abundance matrix by the iterator of the sparse matrix class to obtain the feature-functional unit abundance matrix.
7. The method according to claim 1, wherein In Step 3, for the functional unit abundance matrix obtained in Step 2, infer the existing functional pathways through minpath and calculate the pathway abundance matrix; the steps are as follows: Obtain the mapping relationship between gene families and pathways from the known database, then conservatively reconstruct the pathways through minpath, select the smallest set of pathways that can cover all gene families in the functional unit abundance matrix, and calculate to obtain the pathway abundance matrix.
8. The method according to claim 1, wherein In Step 4, for the functional unit abundance matrix, the feature-functional unit abundance matrix obtained in Step 2, and the pathway abundance matrix obtained in Step 3, read the matrix data and the description files of functional units and pathways through the perl language, and add the descriptions of functional units and pathways to each abundance matrix.
9. A system for implementing the method according to any one of claims 1-8, characterized in that, The system includes a feature sequence processing module, a feature abundance data processing module, a sparse matrix calculation module, a functional unit abundance matrix inference module, and a description information adding module; among them, The feature sequence processing module is used to process the feature sequences of the to-be-predicted functions of microbial metagenomic marker genes, including aligning to the reference marker gene sequences, calculating the evolutionary distance, and constructing an evolutionary tree; The feature abundance data processing module is used to read the feature abundance data of microbial metagenomic marker genes, save it in text form, and establish the mapping relationship between feature sequences and functional units; The sparse matrix calculation module uses the armdillo library to read and calculate the feature abundance matrix corrected by the 16S rRNA copy number or the externally input feature abundance matrix and the feature-functional unit mapping matrix to obtain the functional unit abundance matrix and the feature-functional unit abundance matrix; The functional unit abundance matrix inference module infers the pathway abundance for the obtained functional unit abundance matrix through the minpath algorithm and calculates to obtain the pathway abundance matrix; The description information adding module adds description information to the obtained functional unit abundance matrix, feature-functional unit abundance matrix, and pathway abundance matrix, and outputs the predicted functional abundance file for each sample.
10. The application of the method according to any one of claims 1-8, or the system according to claim 9 in the study of microbial community composition.