Rhizosphere Microbial Community Analysis Based on Metabolomics Metagenomics

Abstract

This invention discloses a rhizosphere microbial community analysis method based on metabolomics metagenomics, belonging to the field of microbial ecological data analysis technology. The method includes the following steps: obtaining a metagenomic microbial abundance matrix, a functional pathway abundance matrix, and a metabolite content matrix, and generating three types of standardized matrices and rhizosphere state characteristics; identifying boundary samples to form an initial corrected rhizosphere stage; constructing an initial stage association network based on the initial corrected rhizosphere stage and the three types of standardized matrices; correcting the attribution of boundary samples to obtain a final corrected rhizosphere stage; and generating the final rhizosphere microbial community analysis result based on the final corrected rhizosphere stage. This invention improves the accuracy of grouping and association identification by correcting stages based on rhizosphere state characteristics and identifying synchronous and lagging associations, solving the problems of inaccurate stage boundary division, insufficient identification of lagging associations, and averaging of stage-specific relationships in existing technologies.

Description

Technical Field

[0001] This invention relates to the field of microbial ecological data analysis technology, and in particular to a method for analyzing rhizosphere microbial communities based on metabolomics metagenomics. Background Technology

[0002] With the development of omics detection and data analysis technologies, there is a growing number of studies on plant rhizosphere microbial communities with continuous sampling, multi-omics paired data, and stage records. The combined use of metabolomics and metagenomics provides a data foundation for elucidating the composition, functional potential, and metabolite changes of rhizosphere microorganisms.

[0003] For example, in the scenario of rhizosphere microbial community analysis in reed wetlands, existing technologies typically begin by assigning a unique number to the rhizosphere samples, performing quality control, clustering, and annotation on metagenomic data to obtain microbial abundance tables and functional pathway abundance tables; and identifying and quantifying metabolomics data to obtain metabolite content tables. Subsequently, multi-omics matrices are merged according to sample number, low-quality features are filtered out, differences are compared among preset sample groups, associations between microorganisms and metabolites are calculated, and the rhizosphere microbial community results are interpreted in conjunction with functional genes, pathways, and planting parameters.

[0004] For example, the Chinese invention patent with announcement number CN121237231B discloses a method and system for regulating soil microbial communities, which includes: continuously collecting real-time environmental data based on a target soil region, combining it with future climate prediction data, and generating a predicted scenario sequence using Monte Carlo simulation; constructing a microbial community composition table, screening core bacterial species, and performing metagenomic sequencing; using the core bacterial species as game participants, constructing a participant strategy set, calculating a payoff matrix library, setting environmental correction coefficients for each scenario, performing scenario-specific game solving, and generating an equilibrium summary table; constructing a regulation response prediction model, setting multiple optimization objectives, introducing robustness constraints and expert rules, generating candidate regulation schemes, and solving for the optimal scheme.

[0005] The above-mentioned technology has at least the following technical problems: In existing technologies, when conducting combined metabolomics and metagenomics analysis of rhizosphere microbial communities in reed wetlands, samples from the same stage are typically grouped by sampling date, season, or manual recording stage. However, the rhizosphere is influenced by water level, root zone temperature, root oxygen release, root exudates, and the redox state of sediments, which can easily lead to samples with different actual conditions being grouped together. Furthermore, the release, accumulation, or consumption of metabolites may lag behind changes in microbial abundance and functional pathways, making it easy to miss associations when only synchronous pairing is used. Similarly, in the analysis of plant rhizosphere microbial communities with continuous sampling, multi-omics pairing data, and stage records, there are also problems such as inaccurate stage boundary delineation, insufficient identification of lagging associations, and averaging of stage-specific relationships, which in turn lead to bias in the selection of key objects and reduced stability of results. Summary of the Invention

[0006] To address the technical problems of inaccurate stage boundary delineation, insufficient identification of lag associations, and the averaging of stage-specific relationships in existing technologies, this invention provides a rhizosphere microbial community analysis method based on metabolomics metagenomics. The technical solution is as follows: The metagenomic microbial abundance matrix, functional pathway abundance matrix, and metabolite content matrix were obtained by aligning the sample numbers. The sampling dates were converted into sampling sequence numbers, generating three types of standardized matrices and rhizosphere state characteristics. Boundary samples were identified based on the rhizosphere state characteristics, sampling sequence numbers, and artificially labeled rhizosphere sampling stages to form an initial corrected rhizosphere stage. An initial stage association network was constructed based on the initial corrected rhizosphere stage and the three types of standardized matrices. The boundary samples were assigned to corrected based on the initial stage association network, boundary samples, and the initial corrected rhizosphere stage to obtain the final corrected rhizosphere stage. The final rhizosphere microbial community analysis results were generated based on the final corrected rhizosphere stage and the three types of standardized matrices.

[0007] The beneficial effects of the technical solutions provided in the embodiments of the present invention include at least the following: 1. The rhizosphere microbial community analysis method based on metabolomics metagenomics provided by this invention identifies boundary samples and forms an initial corrected rhizosphere stage based on rhizosphere state characteristics in the rhizosphere microbial community analysis scenario of reed wetlands. Then, it reversely corrects the boundary sample assignment based on the drift of synchronous association network and lagging association network, so that the final corrected rhizosphere stage corresponds to the actual rhizosphere state of reed rhizosphere samples. This reduces the misclassification of stage boundary samples and reduces the disturbance of misclassified samples to the association network. It effectively solves the problem in the prior art that samples with different true states are classified into the same group due to merging samples according to sampling date, season or manual recording stage.

[0008] 2. This invention generates synchronous and lagging associations separately within the final modified rhizosphere stage, and labels the association edges as stable associations, stage-specific associations, and stage-migrating associations. In the context of plant rhizosphere microbial community analysis with continuous sampling, multi-omics paired data, and stage records, this invention effectively distinguishes the synchronous and sequential changes among metabolites, microbial abundance, and functional pathways. This improves the completeness of lagging association identification and the distinguishability of stage-specific relationships, and solves the problems of missed associations due to synchronous pairing alone, stage-specific relationships being averaged across stages by a unified network, leading to bias in the selection of key objects and reduced result stability. Attached Figure Description

[0009] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0010] Figure 1 A flowchart illustrating the rhizosphere microbial community analysis method based on metabolomics metagenomics provided in this application embodiment; Figure 2 A flowchart of boundary sample network drift reverse correction provided in this application embodiment; Figure 3 A network diagram showing the final stage of association provided in the embodiments of this application; Figure 4 A comparison diagram of network drift between adjacent stages provided in an embodiment of this application. Detailed Implementation

[0011] Embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While some embodiments of the present disclosure are shown in the drawings, it should be understood that embodiments of the present disclosure may be implemented in various forms and should not be construed as limited to the embodiments set forth herein. Rather, these embodiments are provided to provide a more thorough and complete understanding of the present disclosure.

[0012] It should be understood that the accompanying drawings and embodiments of this disclosure are for illustrative purposes only and are not intended to limit the scope of protection of this disclosure. In the description of the embodiments of this disclosure, the term "comprising" and similar terms should be understood as open-ended inclusion, i.e., "including but not limited to". The term "based on" should be understood as "at least partially based on". The term "one embodiment" or "this embodiment" should be understood as "at least one embodiment". The terms "first", "second", etc., may refer to different or the same objects.

[0013] It should be noted that this embodiment will use the combined metabolomics and metagenomics analysis of the rhizosphere microbial community in reed wetlands as an example. Reed rhizosphere samples were collected sequentially over time from the surface of reed roots and near-root soil within the same wetland plot. During sampling, a mixed rhizosphere sample representing the rhizosphere environment at each sampling point was obtained, and a unique sample number was assigned to each sample. The sampling date and artificially marked rhizosphere sampling stage were recorded. In the reed wetland scenario, artificially marked rhizosphere sampling stages could include the greening stage, rapid growth stage, heading and flowering stage, and maturity and yellowing stage. The same rhizosphere sample was divided into two parts: one for metagenomic sequencing and the other for metabolomics analysis. For the sample used for metagenomic sequencing, total DNA was first obtained using soil or rhizosphere microbial total DNA extraction methods, followed by sequencing library construction and high-throughput sequencing. The resulting metagenomic sequencing data, after species annotation, was statistically analyzed by sample number to form a metagenomic microbial abundance matrix. Then, after functional gene or metabolic pathway annotation, functional pathway abundance was statistically analyzed by sample number to form a functional pathway abundance matrix. For samples used in metabolomics analysis, rhizosphere metabolites were first extracted using methanol, water, or other metabolite extraction solvents. Metabolite detection data were then obtained using liquid chromatography-mass spectrometry (LC-MS), gas chromatography-mass spectrometry (GC-MS), or other mass spectrometry methods. The obtained metabolomics data were then processed by peak identification, peak area analysis, or relative content analysis to form a metabolite content matrix based on sample number. The functional pathway abundance matrix represents the pathway abundance results obtained after annotating functional genes or metabolic pathways in metagenomic sequences. It reflects the functional responses of the microbial community in the sample in carbon metabolism, nitrogen metabolism, sulfur metabolism, methane metabolism, or organic matter degradation. The microbial abundance matrix, functional pathway abundance matrix, and metabolite content matrix all use sample number as row identifiers and microorganisms, functional pathways, and metabolites as column features, respectively.

[0014] To make the technical problems, technical solutions and advantages of the present invention clearer, a detailed description will be given below in conjunction with the accompanying drawings and specific embodiments.

[0015] like Figure 1 The diagram shows a flowchart of a rhizosphere microbial community analysis method based on metabolomics metagenomics provided in this application. The method includes the following steps: S1: Read the sample IDs from the metagenomic microbial abundance matrix, functional pathway abundance matrix, and metabolite content matrix, and use the sample IDs that exist in all three types of matrices as the common sample ID set. Build an aligned sample table based on the common sample ID set. In the aligned sample table, each row corresponds to a common sample ID, and the sampling date, manually labeled rhizosphere sampling stage, and row position in the three types of matrices are written for that sample ID. All three types of matrices are reordered according to the common sample ID set, so that the same row corresponds to the same sample ID. If any sample ID lacks a corresponding data row in the metagenomic microbial abundance matrix, functional pathway abundance matrix, or metabolite content matrix, the data row corresponding to that sample ID is deleted, and the sample is not included in subsequent calculations. Arrange the aligned sample table in ascending order of sampling date; the earliest sampling date corresponds to sampling sequence number 1, and the sampling sequence number increments by 1 for each subsequent change in sampling date; when multiple duplicate samples exist for the same sampling date, the multiple duplicate samples use the same sampling sequence number and are retained as different data rows according to the sample ID.

[0016] Logarithmic transformation and standardization were performed on the metagenomic microbial abundance matrix, functional pathway abundance matrix, and metabolite content matrix, respectively. Before logarithmic transformation, zero and missing values ​​were checked column by column. For a given feature column, if a positive value existed in the column, the zero and missing values ​​were replaced with half of the smallest positive value in the column; if no positive value existed in the column, it meant that the feature had no effective observation in the current dataset, and the feature column was deleted. After the replacement, this embodiment used natural logarithm for logarithmic transformation, that is, the natural logarithm was taken for each value after replacement. Natural logarithm is a specific method of logarithmic transformation, used to compress the order-of-magnitude differences in metagenomic abundance, functional pathway abundance, and metabolite content. This process is used to avoid zero values ​​not being able to be logarithmized and to convert abundance or content changes across orders of magnitude into values ​​that are easy to compare. Standardization was performed column by column. For any feature column, the values ​​of the feature column in all retained samples were added together and then divided by the number of retained samples to obtain the mean of the feature column; then the mean was subtracted from the value of each sample in the feature column. Subsequently, the square of the difference between each sample value and the mean is calculated. All squares are summed and divided by the number of samples to be retained. The square root of the result is then taken to obtain the standard deviation of that feature column. Finally, the result of subtracting the mean from each sample is divided by the standard deviation. If the standard deviation of a feature column is 0, it indicates that the feature remains unchanged across samples. In this case, the feature cannot distinguish between different rhizosphere sampling stages and would also result in a denominator of 0 in the standardization step; therefore, this feature column is deleted. After processing, the standardized microbial abundance matrix, standardized functional pathway abundance matrix, and standardized metabolite content matrix are obtained, collectively referred to as the three types of standardized matrices.

[0017] For each sample number, the values ​​in the standardized microbial abundance matrix are first read according to a fixed order of microbial characteristic columns, and then the values ​​in the standardized metabolite content matrix are read according to a fixed order of metabolite characteristic columns. These two values ​​are then concatenated to form a sample vector, which serves as the rhizosphere state feature, characterizing the microbial community structure and overall metabolite distribution in the reed rhizosphere at a given sampling time. It should be noted that the standardized functional pathway abundance matrix is ​​not involved in the generation of the rhizosphere state feature; it is only retained for subsequent interpretation of the functional response under this rhizosphere state, thus avoiding the functional pathway data simultaneously serving the dual functions of stage division and functional interpretation.

[0018] After the above operations, S1 outputs the aligned sample table, three types of standardized matrices, rhizosphere state characteristics, sampling sequence number, and manually labeled rhizosphere sampling stage. The aligned sample table includes at least the sample number, sampling date, sampling sequence number, manually labeled rhizosphere sampling stage, initial corrected rhizosphere stage field, and boundary sample label field. Furthermore, to ensure traceability of subsequent steps, the matrix column names maintain a one-to-one correspondence before and after preprocessing, and deleted feature columns record the reason for deletion, including the absence of positive values ​​or a standard deviation of 0. If the research object lacks a functional pathway abundance matrix, the association between functional pathways and metabolites cannot be output, but the association between microorganisms and metabolites can still be output following the same process. However, this embodiment assumes that all three types of matrices exist, and that each sample entering the common sample number set has metagenomic and metabolomics paired data, and does not consider the above situation.

[0019] S2: Based on rhizosphere state characteristics, sampling sequence numbers, and artificially labeled rhizosphere sampling stages, boundary samples are identified to form initial corrected rhizosphere stages. First, sample numbers are grouped according to artificially labeled rhizosphere sampling stages. The arithmetic mean of the rhizosphere state characteristics corresponding to sample numbers belonging to the same artificially labeled rhizosphere sampling stage is calculated column by column. The resulting average vector is the rhizosphere state characteristic center of that artificially labeled stage, representing the average state of the microbial community structure and overall metabolite distribution within that stage. Specifically, for the j-th feature dimension within a certain artificially labeled stage, the values ​​of all samples in that stage on the j-th feature dimension are summed, and then divided by the number of samples in that stage to obtain the center value of the characteristic center of that stage on the j-th feature dimension. After performing this operation on all feature dimensions, the rhizosphere state characteristic center of that artificially labeled stage is obtained. For each sample number, calculate the distance from its root state feature to the feature center of the current artificial labeling stage, and the distance to the feature center of the adjacent artificial labeling stage. Specifically, the distance is calculated using the Euclidean distance between standardized feature vectors. Assume that the root state feature of a sample includes the 1st to mth feature values, and the feature center of the artificial labeling stage to which the sample belongs also includes the 1st to mth center values. First, subtract the 1st center value from the 1st feature value and square it; then subtract the 2nd center value from the 2nd feature value and square it, and so on up to the mth feature dimension. Then, sum the m squared values ​​and take the square root of the sum to obtain the distance from the sample to the feature center of the current artificial labeling stage. The distance from the sample to the feature center of the adjacent artificial labeling stage is calculated in the same way, except that the center value is replaced with the corresponding center value from the adjacent artificial labeling stage feature center. The first artificial labeling stage is compared only with the adjacent next artificial labeling stage, the last artificial labeling stage is compared only with the adjacent previous artificial labeling stage, and the intermediate artificial labeling stages are compared with the adjacent preceding and following artificial labeling stages respectively. That is, for the same sample number, the distance of the sample to the feature center of the current artificial labeling stage is compared with the distance of the sample to the feature center of the adjacent artificial labeling stage. The smaller the distance of the sample to the feature center of a certain artificial labeling stage, the closer the overall root state of the sample is to the average state of that stage.

[0020] For example, if the abundance of organic acids, sugars and related rhizosphere microorganisms in a rhizosphere sample manually recorded as being in the rapid growth stage is close to the overall state of the sample in the heading and flowering stage, then the distance of this sample to the characteristic center of the heading and flowering stage may be less than the distance to the characteristic center of the rapid growth stage. In this case, this sample is more suitable to participate in the subsequent association network construction according to the rhizosphere state of adjacent stages.

[0021] For each sample number, if its distance to the feature center of the nearest adjacent artificially labeled stage is less than its distance to the feature center of the current artificially labeled stage, then the sample number is assigned to the nearest adjacent stage. If the absolute difference between its distance to the feature center of the current artificially labeled stage and its distance to the feature center of the nearest adjacent artificially labeled stage does not exceed the stage boundary tolerance, then the sample number is marked as a boundary sample and temporarily assigned to the side with the closer distance. If its distance to the feature center of the current artificially labeled stage is equal to its distance to the feature center of the nearest adjacent artificially labeled stage, then the sample number is marked as a boundary sample and temporarily retained in the root sampling stage of the artificial label to which the sample number belongs. The remaining samples are retained in the current artificially labeled stage. The stage boundary tolerance is obtained from the data itself. Specifically, the distances from samples within each artificially labeled stage to the feature center of the current stage are first summarized to obtain the internal distance distribution of the artificial stage. Then, all distance values ​​in this internal distance distribution are arranged in ascending order, and the distance value located in the upper quartile is taken as the stage boundary tolerance. This embodiment uses the distance value corresponding to the upper quartile position because this value comes from most of the fluctuation range within the artificial stages in the current dataset, which can preserve the common root zone state differences within the artificial stages, while identifying samples close to the center of adjacent stages as boundary samples. The stage boundary tolerance represents the upper limit of acceptable root zone state fluctuation within the artificial stage.

[0022] After completing the above division, the number of samples and the number of sampling numbers covered in each provisional stage are counted. If a provisional stage does not reach the minimum number of samples, it is merged into a provisional stage that is temporally adjacent and closer to the center of the rhizosphere state features. The minimum number of samples is determined based on the minimum data requirements for subsequent stage mean comparison and leave-one-out cross-validation. Specifically, in S3, the sample mean of each feature within the stage needs to be calculated, and the mean changes between adjacent stages need to be compared. When there are fewer than 3 samples in a stage, the fluctuation of a single sample will significantly affect the stage mean, making it difficult to reflect the overall rhizosphere state of that stage. At the same time, leave-one-out cross-validation also needs to be performed on the samples within the stage in S3. Leave-one-out cross-validation requires that 1 sample be left out as validation data each time, and the remaining samples as training data. Therefore, after leaving out 1 sample, the training data still needs to retain no less than 2 samples to complete the regression model establishment and parameter selection. Based on the above requirements, in this embodiment, the minimum number of samples is no less than 3 samples, and it covers no less than 2 sampling numbers. If a provisional stage does not meet the minimum sample size requirement, a separate association network is not constructed for that provisional stage. Instead, it is first merged into an adjacent provisional stage according to the aforementioned rules, and then the merged stage is used as the initial corrected root stage to participate in the S3 network construction. After merging, the initial corrected root stage is formed, and the initial corrected root stage corresponding to each sample number and the boundary sample label are output. Non-boundary samples are not re-partitioned in the subsequent S4, and boundary samples are only allowed to undergo one assignment correction between temporally adjacent stages.

[0023] S3: Within each initial modified rhizosphere stage, sample data belonging to that stage are extracted from the three types of standardized matrices. For each microbial feature within that stage, using the microbial abundance as the object to be explained and the total metabolite content within that stage as the input variable, candidate relationships between microorganisms and metabolites are established. For each functional pathway feature within that stage, using the functional pathway abundance as the object to be explained and the total metabolite content within that stage as the input variable, candidate relationships between functional pathways and metabolites are established. These candidate relationships serve as the modeling objects for the synchronous association network.

[0024] Specifically, the synchronous association network is established using a regression model with sparse constraints. This embodiment uses L1-constrained linear regression, also known as Lasso regression. When building a model for a given object, the input is the standardized metabolite content of each sample within the same initial corrected rhizosphere stage, and the output is the corresponding microbial abundance or functional pathway abundance within that stage. L1 constraints reduce the regression coefficients of some metabolites to 0, and metabolites with non-zero coefficients are considered candidate association objects. In L1-constrained linear regression, the Lasso sparse constraint parameter controls the degree to which regression coefficients are compressed to 0; the larger this parameter is, the fewer non-zero coefficients are usually retained, therefore this parameter needs to be determined from the sample data within the same stage. The Lasso sparse constraint parameter is obtained through cross-validation within the same stage, where the number of folds in the cross-validation does not exceed the number of samples within that stage. If the number of samples in a phase is no less than 5, the samples in that phase are divided into 5 parts, with one part used as validation data and the remaining 4 parts used as training data. If the number of samples in a phase is 3 or 4, one sample is reserved as validation data each time, and the remaining samples are used as training data, until each sample has been reserved for validation once. The parameter value corresponding to the minimum mean squared error of cross-validation is selected. If the mean squared errors of multiple parameter values ​​are the same, the parameter value that results in fewer non-zero coefficients is selected.

[0025] To reduce the impact of sample perturbation on single modeling, repeated sampling modeling is performed within the same initial corrected root zone stage. Specifically, sampling with replacement is performed according to the sample size of that stage. If there are n sample numbers in that stage, n random sample numbers are drawn from these n sample numbers. After each draw, the sample number is returned to the original sample so that it can still be drawn in the next draw. After n draws, a resampled sample set containing n sampling locations is formed. The same sample number may appear multiple times in a resampled sample set, or it may not appear at all. Subsequently, a Lasso regression model is rebuilt on this resampled sample set according to the same modeling rules. In this embodiment, the number of repeated samplings is set to 100. For an edge between a microorganism and a metabolite, or an edge between a functional pathway and a metabolite, its retention rate is equal to the number of times the edge has a non-zero coefficient in 100 repeated modelings divided by 100. Each time a non-zero coefficient appears on an edge, it is recorded whether the coefficient is positive or negative. If all non-zero coefficients on the edge are positive or negative in all repeated modeling results, the regression coefficients of that edge are considered to have consistent directions. If the same edge has both positive and negative values, the regression coefficients are considered to have inconsistent directions. The determination of consistent regression coefficient directions is used to decide whether to include the edge in the synchronous association network. It should be noted that the above number of repeated samplings is only an example, and this embodiment does not impose any constraints on it. Technicians can adjust it freely according to the actual situation.

[0026] The association retention ratio threshold is obtained through scrambled data. Specifically, while maintaining the numerical distributions of microbial abundance, functional pathway abundance, and metabolite content, the sample numbers of the standardized metabolite content matrix are randomly shuffled to disrupt the true sample correspondence between metabolites and microorganisms or functional pathways. On this scrambled data, the retention ratio of each candidate edge is calculated according to the same repeated sampling modeling rules as the real data. In this embodiment, the scrambling is repeated 100 times, and the retention ratios of all scrambled candidate edges are summarized. The retention ratio corresponding to the highest quantile position is taken as the association retention ratio threshold for this initial corrected rhizosphere stage to exclude most candidate associations generated by random pairing. This embodiment uses the retention ratio corresponding to the 95th percentile as the specific implementation method. When the retention ratio of an edge in the real data reaches this threshold and the coefficient direction is consistent, the edge is written into the synchronous association network. The synchronous association network consists of microorganisms, functional pathways, and metabolites retained in the same stage as nodes, and non-zero associations with consistent directions as edges.

[0027] For candidate association edges between microorganisms and metabolites that have not entered the synchronous association network, and for candidate association edges between functional pathways and metabolites, the changes in their stage mean values ​​between adjacent initial modified rhizosphere stages are further compared. Specifically, for a certain feature, the sample mean value of that feature in two adjacent initial modified rhizosphere stages is calculated. Then, the sample mean value of the feature in the later initial modified rhizosphere stage is subtracted from the sample mean value of the earlier initial modified rhizosphere stage. The difference is the stage mean value change result. When the difference is greater than 0, the stage mean value change result is an increase; when the difference is less than 0, the stage mean value change result is a decrease; when the difference is equal to 0, the stage mean value change result is no change. A stage mean value change result not equal to 0 indicates that the feature has changed. If one endpoint of a candidate edge undergoes a non-zero change in the previous stage, and the other endpoint of the candidate edge undergoes a non-zero change in the later stage, it is included in the lag candidate set.

[0028] When a metabolite changes in the preceding initial modified rhizosphere stage, and a microorganism or functional pathway shows a non-zero change in the following initial modified rhizosphere stage, and the interval does not exceed the maximum lag stage number, a metabolite-first lag candidate edge is established; when a microorganism or functional pathway changes first, followed by a metabolite change, a microorganism or functional pathway-first lag candidate edge is established. The maximum lag stage number represents the maximum number of modified stages between the object that is allowed to change first and the other endpoint of the same candidate edge. In the reed wetland rhizosphere example, the maximum lag stage number is fixed at 1, meaning only adjacent preceding and following stages are checked to avoid incorrectly pairing seasonal changes that are too far apart. For example, if a certain type of organic acid shows a stage mean change first between the rapid growth period and the heading and flowering period, while a certain microorganism or functional pathway involved in the nitrogen cycle shows a stage mean change only in the following adjacent stage, then this combination enters the metabolite-first lag candidate edge; conversely, if a rhizosphere microorganism or functional pathway changes first, and the related amino acid or sugar metabolite changes later, then this combination enters the microorganism or functional pathway-first lag candidate edge.

[0029] Lagged candidate edges are repeatedly modeled in the corresponding preceding and following stage samples according to the sampling with replacement rule, and the retention ratio and coefficient direction of the lagged candidate edges are calculated. Lagged candidate edges that reach the association retention ratio threshold and have the same preceding direction are written into the initial stage association network, and the lagged stage number and preceding object are recorded. When both the metabolite-preceding and microbial or functional pathway-preceding directions reach the association retention ratio threshold, they are recorded as bidirectional lagged edges; bidirectional lagged edges only indicate that a unique order cannot be determined at the current sampling interval, and are not interpreted as the two objects being mutually causal. Thus, S3 outputs the initial stage association network, where each edge records at least two endpoints, the stage it belongs to, the direction of the regression coefficient, the synchronous or lagged type, the lagged stage number, the preceding object, and the retention ratio.

[0030] To ensure comparability of network results across different stages, the same feature column order is used for the same type of modeling object across all stages. If sample composition changes due to sample merging or boundary sample correction at a particular stage, only the rows of samples involved in the modeling are changed; the already determined order of microbial, functional pathway, and metabolite feature columns remains unchanged. If all samples of a particular object have the same value in the current stage, making it impossible to establish an effective regression model, then no related edges are generated for that object in that stage, and the reason is noted in the stage's network construction record.

[0031] S4: Based on the initial stage association network, boundary samples, and the initial corrected root stage, the boundary samples are assigned a specific category to obtain the final corrected root stage. It should be noted that this step only corrects the category of boundary samples; non-boundary samples retain the initial corrected root stage output from S2 and are used as fixed samples in the modeling process when reconstructing adjacent stage networks under both candidate categories. Each boundary sample undergoes only one candidate category comparison, without iterative iteration; candidate categories are limited to the preceding and following corrected stages in time and are not allowed to move across stages.

[0032] like Figure 2 The diagram shown is a flowchart of the boundary sample network drift reverse correction process provided in this application embodiment. For each boundary sample, two candidate assignments are constructed: one for the previous correction stage and one for the next. Under each candidate assignment, only the association networks of the adjacent two stages of that boundary sample are reconstructed. The specific reconstruction rules are the same as in S3, including synchronous association network establishment, lag candidate edge establishment, repeated sampling modeling, association retention ratio threshold judgment, and coefficient direction consistency judgment.

[0033] Using the initial stage association network of the adjacent two sides of the boundary sample in S3 as a control, three types of conflict edges are counted. The first type is the number of synchronous association edges that disappear, i.e., when a synchronous edge exists in the control network but does not exist in the corresponding stage after the candidate attribution reconstruction, the number of synchronous association edges that disappear is incremented by 1. The second type is the number of lagging association edges that disappear, i.e., when a lagging edge exists in the control network but does not exist in the corresponding stage after the candidate attribution reconstruction, the number of lagging association edges that disappear is incremented by 1. The third type is the number of lagging edge predecessor object reversals, i.e. when a lagging edge exists both before and after reconstruction, but the predecessor object changes from metabolite to microorganism or functional pathway, or from microorganism or functional pathway to metabolite, the number of lagging edge predecessor object reversals is incremented by 1. In this embodiment, only edges in S3 that have reached the association retention ratio threshold and whose coefficient directions are consistent are counted. Candidate edges that have not reached the stable retention condition are not included in the conflict count, and the same edge is counted only once under the same candidate attribution. If the same edge meets both the disappearance and reversal conditions, it is counted as disappearance and is not counted again for reversal. The number of conflicting edges to which a candidate belongs is obtained by adding the three types of numbers together.

[0034] Compare the number of conflicting edges corresponding to the two candidate assignments. If the number of conflicting edges for the candidate assignment in the previous correction stage is less, then the candidate assignment from the previous correction stage is adopted; if the number of conflicting edges for the candidate assignment in the subsequent correction stage is less, then the candidate assignment from the subsequent correction stage is adopted. If the total number of conflicting edges corresponding to the two candidate assignments is the same, then the assignment determined by the root-state feature distance in S2 is retained, that is, the assignment from the boundary sample to the feature center of the previous stage and the feature center of the subsequent stage with the smaller distance is retained. After all boundary samples have undergone one assignment correction, the final corrected root-state stage is obtained. If there are no boundary samples, then the final corrected root-state stage is equal to the initial corrected root-state stage output by S2.

[0035] S5: Based on the final revised rhizosphere stage and the three types of standardized matrices, generate the final rhizosphere microbial community analysis results. Specifically, within each final revised rhizosphere stage, synchronous and lag associations are regenerated according to the processing method in S3. When rebuilding the network, microbial abundance and functional pathway abundance are still treated as objects to be explained, and metabolite content is still used as input variables. Among them, the repeated sampling modeling, the determination of the association retention ratio threshold, the judgment of coefficient direction consistency, the establishment of lag candidate edges, the limitation of the maximum number of lag stages, and the rules for recording bidirectional lag edges are all consistent with S3.

[0036] like Figure 3 The diagram shown is a network diagram of the final stage of the rhizosphere provided in this application embodiment. The network of the final modified rhizosphere stage includes three types of nodes: microbial nodes, metabolite nodes, and functional pathway nodes. Microbial nodes represent microbial features in the standardized microbial abundance matrix, metabolite nodes represent metabolite features in the standardized metabolite content matrix, and functional pathway nodes represent functional pathway features in the standardized functional pathway abundance matrix. In the diagram, rhizobia, methanogenic bacteria, Pseudomonas, actinomycetes, Bacillus, and sulfur-reducing bacteria are used as example microbial nodes in the rhizosphere samples of reed wetlands; organic acids, sugars, amino acids, and phenolic acids are used as example metabolite nodes; and methane metabolism, carbon metabolism, nitrogen metabolism, sulfur metabolism, and organic matter degradation are used as example functional pathway nodes. Solid lines represent synchronous associations written into the synchronous network within the same final modified rhizosphere stage; dashed lines with arrows represent lag associations that have reached the association retention ratio threshold and have a defined preceding object, with the arrow pointing from the preceding object to the object that changes later; bidirectional dashed lines represent bidirectional lag edges that have reached the association retention ratio threshold in both directions at the current sampling interval. The thickness of the edge is used to indicate the retention rate. A thick line indicates that the edge has a higher retention rate in repeated sampling modeling, while a thin line indicates that the edge has a lower retention rate in repeated sampling modeling.

[0037] After the final network is completed, the associated edges in adjacent final modified rhizosphere stages are compared. If two edges connect endpoints with the same type and name, they are considered to be the same comparable edge. Endpoint types include microorganisms, functional pathways, and metabolites, and endpoint names are the corresponding column names in the three types of standardized matrices. That is, when two edges connect the same microorganism column name and the same metabolite column name, or both connect the same functional pathway column name and the same metabolite column name, they are considered to be the same comparable edge. For the same comparable edge, the regression coefficient direction, lag number, and preceding objects are further compared.

[0038] Specifically, if the comparable edge is a synchronous correlation edge, then its existence and regression coefficient direction are compared in adjacent final modified root stages; when the same synchronous correlation edge exists in both adjacent final modified root stages and the regression coefficient direction is the same, the synchronous correlation edge is marked as a stable correlation; when the synchronous correlation edge exists only in one final modified root stage, the synchronous correlation edge is marked as a stage-specific correlation; when the same synchronous correlation edge exists in both adjacent final modified root stages but the regression coefficient direction is different, the synchronous correlation edge is marked as a stage-transitional correlation. If the comparable edge is a lagged correlation edge, then compare its existence, regression coefficient direction, number of lag stages, and preceding objects in adjacent final modified root stages. When the same lagged correlation edge exists in both adjacent final modified root stages and has the same regression coefficient direction, number of lag stages, and preceding objects, the lagged correlation edge is marked as a stable correlation. When the lagged correlation edge exists only in one final modified root stage, the lagged correlation edge is marked as a stage-specific correlation. When the same lagged correlation edge exists in both adjacent final modified root stages but has different regression coefficient directions, number of lag stages, or preceding objects, the lagged correlation edge is marked as a stage-transitional correlation.

[0039] like Figure 4The image shown is a comparison diagram of adjacent stage network drift provided in an embodiment of this application, used to illustrate the continuation, disappearance, and migration of association edges between different final revised rhizosphere stages. The left side of the image shows the stage association network of final revised stage 1, and the right side shows the stage association network of final revised stage 2. Both stages include microbial nodes, metabolite nodes, and functional pathway nodes, and the association edges are compared according to the same endpoint type and endpoint name. If an association edge exists in both final revised stage 1 and final revised stage 2, and the regression coefficient direction, lag stage number, and preceding object are the same, then the association edge is marked as a stable association and indicated by a solid arrow. If an association edge appears only in one of the final revised rhizosphere stages, then the association edge is marked as a stage-specific association and indicated by a dashed arrow. If an association edge exists in two adjacent final revised rhizosphere stages, but its regression coefficient direction, lag stage number, or preceding object changes, or the association object shifts from the original metabolite-functional pathway relationship to another metabolite-functional pathway relationship, then the association edge is marked as a stage-migrating association and indicated by a hollow dashed arrow. The associations between rhizobia, organic acids, and carbon metabolism remain consistent across both stages, serving as examples of stable associations. The association between actinomycetes and carbohydrates exhibits stage-specific differences between adjacent stages, serving as examples of stage-specific associations. Furthermore, the associations related to actinomycetes shift from carbohydrate or nitrogen metabolism towards organic matter degradation, serving as examples of stage-migrating associations. This network drift comparison visually demonstrates the stage-specific changes in the rhizosphere microbial community of reed wetlands across different final revised rhizosphere stages, avoiding the averaging of cross-stage differences into a single unified network.

[0040] The final rhizosphere microbial community analysis results include associations between microorganisms and metabolites, associations between functional pathways and metabolites, synchronous or lag types, number of lag stages, preceding objects, regression coefficient direction, retention ratio, and markers for stable associations, stage-specific associations, and stage-migrating associations within each final revised rhizosphere stage. Artificially labeled rhizosphere sampling stages are only added as annotations to the final rhizosphere microbial community analysis results and do not alter the final revised rhizosphere stage. This approach allows for the analysis of stage-specific synchronous or lag relationships between rhizosphere metabolites such as organic acids, phenolic acids, sugars, or amino acids in reed wetland rhizosphere samples and rhizosphere microorganisms and functional pathways involved in carbon metabolism, nitrogen metabolism, sulfur metabolism, methane metabolism, or organic matter degradation. For example, in continuous rhizosphere samples from the rapid growth stage to the heading and flowering stage of reeds, it is possible to identify whether there is an adjacent stage lag relationship between changes in organic acid abundance and changes in nitrogen cycle-related microbial abundance, and whether changes in sugars or amino acids are stable, stage-specific, or stage-migrating associations with changes in carbon metabolism functional pathways. Meanwhile, by relying on boundary sample identification and network drift-based attribution correction, the impact of inaccurate manual stage boundary division on the associated network is reduced, the lag association omission caused by pairing based solely on synchronous sampling is decreased, and the stability and reproducibility of stage-specific relationships in continuous rhizosphere sampling data are improved.

[0041] Through the above description of the implementation methods, those skilled in the art can clearly understand that, for the sake of convenience and brevity, only the division of the above functional modules is used as an example. In actual applications, the above functions can be assigned to different functional modules as needed, that is, the internal structure of the above functions can be divided into different functional modules to complete all or part of the functions described above.

[0042] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions within the technical scope disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for analyzing rhizosphere microbial communities based on metabolomics and metagenomics, characterized in that, Includes the following steps: S1: Obtain metagenomic microbial abundance matrix, functional pathway abundance matrix, and metabolite content matrix aligned by sample number, convert sampling date to sampling sequence number, and generate three types of standardized matrices and rhizosphere state characteristics; S2: Based on the rhizosphere state characteristics, the sampling sequence number, and the manually marked rhizosphere sampling stage, identify boundary samples to form an initial corrected rhizosphere stage; S3: Construct the initial stage association network based on the initial modified root stage and the three types of normalized matrices; S4: Based on the initial stage association network, the boundary samples, and the initial corrected root stage, the boundary samples are assigned to correct their affiliation to obtain the final corrected root stage; S5: Based on the final corrected rhizosphere stage and the three types of standardized matrices, generate the final rhizosphere microbial community analysis results.

2. The rhizosphere microbial community analysis method based on metabolomics metagenomics as described in claim 1, characterized in that, The metagenomic microbial abundance matrix, functional pathway abundance matrix, and metabolite content matrix obtained by alignment include: Read the sample numbers from the metagenomic microbial abundance matrix, the functional pathway abundance matrix, and the metabolite content matrix, and write the data rows with the same sample number into the aligned sample table accordingly; If any sample number is missing a data row corresponding to that sample number in any of the metagenomic microbial abundance matrix, the functional pathway abundance matrix, and the metabolite content matrix, then delete the data row corresponding to that sample number. Sample numbers are generated in chronological order based on the sampling date.

3. The rhizosphere microbial community analysis method based on metabolomics metagenomics as described in claim 2, characterized in that, include: Logarithmic transformation and standardization were performed on the metagenomic microbial abundance matrix, the functional pathway abundance matrix, and the metabolite content matrix to obtain three types of standardized matrices. The three types of standardized matrices include a standardized microbial abundance matrix, a standardized functional pathway abundance matrix, and a standardized metabolite content matrix; The standardized microbial abundance and standardized metabolite content corresponding to each sample number are connected in a fixed order to generate the rhizosphere state characteristics.

4. The rhizosphere microbial community analysis method based on metabolomics metagenomics as described in claim 1, characterized in that, The formation of the initial modified root zone stage includes: According to the artificially labeled rhizosphere sampling stage, the rhizosphere state features corresponding to the sample numbers belonging to the same artificially labeled rhizosphere sampling stage are summarized to generate the rhizosphere state feature center for each artificially labeled rhizosphere sampling stage. For each sample number, calculate the distance from its rhizosphere state feature to the feature center of the current artificial labeling stage, and the distance to the feature center of the adjacent artificial labeling stage; Wherein, the distance to the feature center of this artificial labeling stage is the sum of the differences between the rhizosphere state feature corresponding to the sample number and the feature center of this artificial labeling stage on the same feature dimension, and the distance to the feature center of the adjacent artificial labeling stage is the sum of the differences between the rhizosphere state feature corresponding to the sample number and the feature center of the adjacent artificial labeling stage on the same feature dimension.

5. The rhizosphere microbial community analysis method based on metabolomics metagenomics as described in claim 4, characterized in that, The formation of the initial modified root zone stage also includes: When the distance from the root state feature corresponding to the sample number to the feature center of the adjacent artificially labeled stage is less than the distance to the feature center of the current artificially labeled stage, the sample number is assigned to the nearest adjacent stage. When the absolute difference between the distance to the feature center of the adjacent artificially labeled stage and the distance to the feature center of this artificially labeled stage does not exceed the stage boundary tolerance, the sample is numbered and marked as a boundary sample, and is included according to the side with the closer distance. When the distance to the feature center of the adjacent artificial labeling stage is equal to the distance to the feature center of the current artificial labeling stage, the sample number is marked as a boundary sample and temporarily retained in the root sampling stage of the artificial label to which the sample number belongs; After the division is completed, if the number of samples in a certain initial modified rhizosphere stage is lower than the minimum number of samples, the initial modified rhizosphere stage is merged into an initial modified rhizosphere stage that is temporally adjacent and closer to the center of the rhizosphere state features.

6. The rhizosphere microbial community analysis method based on metabolomics metagenomics as described in claim 1, characterized in that, The construction of the initial stage association network includes: Within each initial modified root stage, sample data belonging to that initial modified root stage are extracted from the three types of standardized matrices; Using the standardized microbial abundance and standardized functional pathway abundance within the initial modified rhizosphere stage as the objects to be explained, and the standardized metabolite content within the initial modified rhizosphere stage as the input variable, candidate relationships between microorganisms and metabolites, and candidate relationships between functional pathways and metabolites are established, and the candidate relationships are used as the modeling objects of the synchronous association network.

7. The rhizosphere microbial community analysis method based on metabolomics metagenomics as described in claim 6, characterized in that, The construction of the initial stage association network also includes: Regression models with sparse constraints are established for the candidate relationships between the microorganisms and metabolites and the candidate relationships between the functional pathways and metabolites, and repeated sampling modeling is performed. When the retention ratio of candidate association edges between microorganisms and metabolites, or between functional pathways and metabolites, reaches the association retention ratio threshold and the coefficients are in the same direction, the candidate association edges between the microorganisms and metabolites or between the functional pathways and metabolites are written into the synchronous association network.

8. The rhizosphere microbial community analysis method based on metabolomics metagenomics as described in claim 7, characterized in that, The construction of the initial stage association network also includes: For candidate association edges between microorganisms and metabolites that have not entered the synchronous association network, and for candidate association edges between functional pathways and metabolites, compare the changes in their stage mean between adjacent initial modified rhizosphere stages. Wherein, the change in the stage mean is the sample mean of the next initial corrected root stage minus the sample mean of the previous initial corrected root stage. When metabolites change in the previous initial modified rhizosphere stage, and microorganisms or functional pathways change in the subsequent initial modified rhizosphere stage, and the interval does not exceed the maximum number of lag stages, metabolite lag candidate edges are established, wherein the change is the change in the mean of the stage that is not equal to 0. When microorganisms or functional pathways change first and metabolites change later, establish candidate edges for microorganisms or functional pathways that change first and then lag. The lagging candidate edges are repeatedly modeled in the corresponding preceding and following stage samples. The lagging candidate edges that reach the association retention ratio threshold and have the same preceding direction are written into the initial stage association network, and the number of lagging stages and preceding objects are recorded. When both directions reach the associated retention ratio threshold, it is recorded as a bidirectional hysteresis edge.

9. The rhizosphere microbial community analysis method based on metabolomics metagenomics as described in claim 1, characterized in that, The process of obtaining the final modified root stage includes: Only the boundary samples are assigned a new classification; non-boundary samples are not reclassified. For each boundary sample, two candidate assignments are constructed: one for the preceding correction stage with adjacent assignment time and the other for the following correction stage with adjacent assignment time. Under two candidate attributions, only the association network of the adjacent two sides of the boundary sample is reconstructed, and the number of synchronous association edges that have been written in the control network disappearing, the number of lagging association edges that have been written in the control network disappearing, and the number of lagging edge preceding objects being reversed are counted respectively. When a candidate assignment corresponds to fewer conflicting edges, that candidate assignment is adopted; When the number of conflicting edges corresponding to two candidate attributions is the same, the attribution determined by the root state feature distance in S2 is retained; After all boundary samples have undergone one assignment correction, the final corrected root stage is obtained.

10. The rhizosphere microbial community analysis method based on metabolomics metagenomics as described in claim 1, characterized in that, The final rhizosphere microbial community analysis results include: Based on the final revised root stage and the three types of standardized matrices, synchronous associations and lag associations are regenerated in each final revised root stage; Compare the associated edges in adjacent final modified root stages; For the same synchronous correlation edge, if it exists in both adjacent final modified root stages and the regression coefficients are in the same direction, the synchronous correlation edge is marked as a stable correlation. If a synchronous association edge exists only in a final modified root stage, it is marked as a stage-specific association. When a synchronous correlation edge exists in both adjacent final modified root stages but the regression coefficients are in different directions, the synchronous correlation edge is marked as a stage migration correlation. For the same lagged correlation edge, if it exists in both adjacent final modified root stages and the regression coefficient direction, the number of lagged stages, and the preceding objects are all the same, the lagged correlation edge is marked as a stable correlation. When a hysteretic edge exists only in a final modified root stage, it is marked as a stage-specific edge. If it exists in both adjacent final modified root stages, but the regression coefficient direction, lag stage number or preceding object are different, then the lag association edge is marked as a stage migration association. Artificially labeled rhizosphere sampling stages were included as annotations in the final rhizosphere microbial community analysis results.

Images