Trend analysis method based on identification result of microorganisms in production environment

Abstract

The invention provides a trend analysis method (MBE) based on identification results of microorganisms in a production environment. The MBE platform can manage meta-information (including room, clean room grade, sampling mode and plan category) and identification results of microorganism samples detected in medicine, food and cosmetic production environments, and performs data aggregation according to specified unit time (year, month, day and batch). The genus-species proportion of the environmental microbial flora in a given time range and the distribution of the flora in each classification meta-information are analyzed, the time-varying rule of the flora is analyzed from four dimensions of phenotypic classification, bacterial genus-species level, concerned bacteria and community diversity, and the time and space distribution rule of the flora is analyzed through frequency analysis; and analyzing the correlation between different classifications or between different bacterial genus and species through multivariate statistical analysis. The method provided by the invention can help to identify microorganisms possibly having significant influence on product quality and safety, assess and control microbial pollution risks, and ensure quality and safety.

Description

Technical Field

[0001] The present invention relates to the field of bioinformatics analysis, and mainly to a phenotype database of microorganisms and a method for analyzing the trend of microbial flora in a production environment. Background Art

[0002] In production environments such as medicine and food, process control of microorganisms is an important part of current cGMP and quality risk management. Monitoring of environmental microorganisms requires not only control of the microbial load, but also identification of the types of microorganisms. Because the trend of the isolated microbial species provides additional valuable information for the management and control of the production environment, for example, it can help confirm the common flora in the environment for risk management, discover microorganisms that appear repeatedly in specific areas, or reveal changes in the number of microorganisms of concern, thereby indirectly reflecting the stability of the production environment and the effectiveness of control measures. Therefore, the pharmacopoeia and related regulations also support the inclusion of relevant identification requirements in environmental microbial monitoring to determine the type information of microorganisms and provide more effective data support for environmental control.

[0003] However, the current pharmacopoeia and related regulations do not provide trend analysis methods for microbial species. At present, the trend analysis of manufacturers is often based on load data, corresponding to monitoring technologies such as HPC and even ATP. Some commercial software such as SPASS, MiniTab, Excel, etc. provide statistical models for trend analysis of microbial load data. However, these software also have some shortcomings. First, these software are not specially designed for microbial monitoring data, and the applicability of the model has not been verified. Second, these commercial software only analyze the quantity and cannot combine the microbial species data for trend analysis. Summary of the invention

[0004] In order to overcome the deficiencies of the prior art, the present invention provides a trend analysis method based on the identification results of production environment microorganisms.

[0005] A trend analysis method based on the identification results of production environment microorganisms, the steps are as follows: S1 summarizes the sampling meta-information, identification results, and phenotypic characteristics of the strain to form one data record for one identification result; S2 performs matrix transformation on the summary information obtained by S1, aggregates the data according to the specified time unit, and forms an aggregate matrix; S3 performs genus / species-level proportion statistics on the aggregation matrix obtained in S2, and displays the top 10 species / genera in abundance; S4 The aggregation matrix obtained in S2 is analyzed by source and visualized using a Sange diagram to show the proportional distribution relationship between different classifications and genera / species; S5 uses the aggregate matrix obtained in S2 to perform trend analysis on the abundance of bacterial communities over time and displays it in a bubble chart; S6 Draw a trend analysis line graph for each genus / species of concern over time, and calculate the warning line and correction line based on the model fitting; S7 performs trend analysis on the aggregation matrix obtained in S2 according to phenotype classification, calculates the average value, and displays it in the form of a line graph; S8 uses the aggregation matrix obtained in S2 to calculate α-diversity indices such as Shannon index, Simpson index and coverage, and displays the trend of these indices over time through line graphs; S9 further converts the aggregation matrix obtained in S2 into detection frequency, analyzes the co-detection phenomenon of different bacterial genera between different classifications and the properties of their residence in the environment, and visualizes them through a collection diagram; S10 further converts the aggregation matrix obtained in S2 into detection frequency, analyzes the continuous detection of different bacterial genera within the analysis time range and under the specified source, displays it through a heat map, and identifies the three types of bacterial communities: core, persistent, and transient; S11 performs dimensionality reduction analysis on the aggregation matrix obtained in S2 through multivariate statistics to obtain the correlation and distance information between bacterial genera or species and between different classifications, and displays them visually.

[0006] The data record input in step S1 is based on the identification result of a single colony or pure culture, which is expressed in a Latin name. The sampling meta-information includes the sampling point, clean room level, sampling method, plan category and sampling time information of the sample. The identification result information and the sampling meta-information are associated through the sample number, and the phenotypic information is obtained by querying the phenotypic database through the Latin name of the identification result.

[0007] In step S2, data is aggregated for different meta-information or phenotypic information classification fields using the dict and pandas functions in Python according to the specified time unit to generate an aggregation matrix of bacterial species abundance under different classifications.

[0008] In step S3, statistics are performed on the aggregation matrix according to all samples and according to the classification methods of room, clean room level, sampling method, and plan category, and the top 10 bacterial genera / species with abundance and proportion under different classifications are displayed using a Nightingale rose diagram.

[0009] In step S5, the trend analysis of bacterial abundance over time is performed, with time as the X-axis and bacterial genus / species as the Y-axis. The abundance and proportion of bacterial genus / species are mapped to the size and color depth of the bubble for visualization.

[0010] In step S6, a trend analysis line graph of changes over time is drawn for a single genus / species of concern. Warning limits and correction limits need to be calculated. The warning limits and correction limits are obtained by using the fitdistr function to perform maximum likelihood fitting of the data distribution and using the chisq.test method to perform a chi-square goodness of fit test. The best model distribution of the data is determined based on the p-value of the test and calculated using the QC_Lines or nptol.int method.

[0011] In step S7, trend analysis is performed according to phenotypic classification. It is necessary to classify and summarize the bacterial genera / species according to the phenotypic fields, calculate the mean, and draw a line graph of the change of phenotypic group abundance over time with time as the X-axis.

[0012] In step S8, the diversity index is calculated by using the R package vegan to calculate the community α-diversity index under the specified time unit and specified classification, including the Shannon index, Simpson index and good's coverage, and their means, and visualized through a line graph.

[0013] In step S9, the genus / species abundance in the aggregation matrix is ​​normalized to "detected" or "not detected", thereby converting the aggregation matrix into a detection frequency matrix. The pandas function is used to perform a set analysis on different classifications, and the R package VennDetail is used to obtain the common bacterial genera / species between different classifications and display them in a combination of graphs and tables.

[0014] In step S10, the aggregation matrix is ​​converted into a frequency matrix as in step S9, with time as the X-axis and genus / species classification as the Y-axis. The abundance is mapped to color, and a heat map is used to display the temporal continuity of genus / species.

[0015] In step S11, the cor function in the R package corrplot is used to perform correlation analysis on the bacterial genera / species of the aggregation matrix, analyze the correlation between bacterial genera / species, and the similarity between bacterial communities under different classifications, use the principal component analysis method to perform dimensionality reduction analysis on bacterial genera of different classifications, and use the R package FactoMineR to draw scatter plots and confidence ellipses for visualization to identify outliers.

[0016] Beneficial effects of the present invention: At present, the pharmacopoeia and relevant regulations do not provide trend analysis methods for microbial species. At present, the trend analysis of production enterprises is often based on load data. The system that only performs load analysis only provides the results and change trends of the total number of microorganisms, and cannot provide detailed information on the types of microorganisms. There is a lack of specific assessment of potential risks and it is impossible to provide detailed data for targeted control measures. The present invention provides an analysis framework for trend analysis of microbial monitoring in pharmaceutical production environments that combines species information, which is multi-dimensional and comprehensive; the present invention provides a bacterial phenotypic classification database, which provides important data for judging the risks and sources of microorganisms in pharmaceutical production environments; the present invention integrates detailed information on microbial species, quantity and trend changes, which can help identify microorganisms that may have a significant impact on product quality and safety, and can help pharmaceutical companies better assess and control microbial contamination risks and ensure product quality and safety. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 A flow chart of the present invention.

[0018] Figure 2 This is a sample data for the input data summary table.

[0019] Figure 3 This is an example of aggregated data after the input data is processed.

[0020] Figure 4 This is the Nightingale rose diagram of the top 10 bacterial genera / species.

[0021] Figure 5 This is a Sankey diagram of the bacterial community composition distribution.

[0022] Figure 6 This is a bubble chart showing the top 10 bacterial genera / species changing over time.

[0023] Figure 7 Line chart for trend analysis of bacteria of concern.

[0024] Figure 8 Line graph for phenotypic trend analysis.

[0025] Fig. 9 Line graph for diversity index zone analysis.

[0026] Fig.10 A collection diagram for analyzing common and unique bacterial genera / species from different sources.

[0027] Fig.11 Heat map for analysis of core, persister and transient microbiota.

[0028] Fig.12 PCA diagram for analyzing the distance between bacterial communities at different sites.

[0029] Fig.13 This is a genus / species correlation analysis. DETAILED DESCRIPTION

[0030] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention. It should be noted that the embodiments of the present invention and the features in the embodiments can be combined with each other without conflict.

[0031] Example 1 This embodiment provides a method for analyzing bacterial flora trends based on the identification results of microorganisms collected in the production environment. Figure 1 , the steps include: S1 Organize the summary table of identification results according to the sampling meta-information (including room, clean room level, sampling method, plan category) and the Latin name of the identification result, and query the phenotype database with the Latin name of the identification result to establish a summary table containing phenotype classification. Some data in the summary table are as follows Figure 2 shown.

[0032] S2 performs matrix transformation on the summary information obtained in S1, aggregates the data according to the specified time unit, and aggregates the data to form an aggregate matrix containing species / genus abundance. This example aggregates the abundance at the genus level with the month as the time unit. The matrix style is as follows: the first column is the Latin name of the species genus, the first row is the classification information, including sources such as sedimentation bacteria, floating bacteria, PW, etc., and each cell value is the number of detected bacterial genera under the specific classification of the month. Figure 3 For example Figure 3 The 44 in the fourth column of the second row indicates that species of the genus Afpia under the PW classification were detected 44 times in that month.

[0033] S3 aggregates the data in the summary information table obtained in S1 at the genus and species level, calculates the proportion of the top 10 bacterial genera / species, and draws a Nightingale rose diagram for visualization ( Figure 4 ). The example shows the Latin names and proportions of the top 10 genera. For example, the first one is Phyllobacterium, which was detected 46 times, accounting for 9.39% of the detection ratio of this category. They are arranged from high to low (the size of the color block corresponds to the proportion). Species outside the top 10 are summed up in Other.

[0034] S4 aggregates the data of the summary information table obtained in S1 according to the genus / species of the identification results, and combines the classification of meta-information to perform distribution analysis on the multi-level meta-information and the summarized identification results, and uses the Sankey function in python to draw a Sankey diagram for visualization. The values ​​of meta-information and species / genus are displayed in rectangular boxes, the size of the rectangular boxes represents the quantity distribution, and the lines between the rectangular boxes represent the flow relationship between the two variables. Figure 5 It shows that Gram-negative bacteria mainly come from PW and WFI, and from PW Station 2, Workshop 209, Workshop 701, and Workshop NXM. In addition, it also shows which genera the negative bacteria, positive bacteria, and fungi in the phenotype are composed of. For example, the negative bacteria in the example are mainly composed of Phyllobacterium, Afipia, Sphingomonas, etc.

[0035] S5 The summary information table obtained in S1 is aggregated by genus, species, and time unit (e.g., by "month"), with time as the x-axis and the top 10 genera / species as the y-axis. The abundance of each month is mapped to the size of the bubble for visualization ( Figure 6 ).For example Figure 6 Afipia was detected from March to July, with the highest detection amount in April, and the bubbles were also the largest and darkest; Staphylococcus had the highest detection amount in July; Phyllobacterim showed a gradual increasing trend from January to June, and the bubble chart can better display the changes in detection trends.

[0036] S6 Aggregate the bacteria of concern by genus, species and time unit (e.g., by "month"), plot the time trend of individual bacteria with time as the x-axis, and draw warning lines and correction limits, marking the points that exceed the limit ( Figure 7 ). Figure 7 To extract the Phyllobacterium alone, the test results from January to December 2022 are summarized and plotted on a monthly basis, with the abscissa representing the month and the ordinate representing the abundance of detection. The yellow warning line is the upper limit of the 95% interval calculated based on the negative binomial distribution model. The month that exceeds the warning line is identified through this graph as June 2022. Months that exceed the warning line represent the loss of control of the environmental flora, that is, the destruction of the stable state of the environment, and are time nodes that require attention and measures to be taken.

[0037] S7 summarizes the data obtained in S1 using phenotypic classification information and time units, visualizes the proportion of each component in the form of a line graph, and calculates the average proportion of each classification, with arrows indicating the direction of change over time ( Figure 8 ). Figure 8Negative means Gram-negative bacteria, positive means Gram-positive bacteria, fungi means eukaryotic microorganisms such as molds and yeasts, and the rest that cannot be classified are unknown. The aggregation matrix containing the Latin names of the identification results is further classified and summarized according to the phenotype after querying the phenotype to obtain the number of microorganisms detected in the corresponding phenotype each month, and a line graph is drawn to analyze the trend of different phenotypes over time. For example Figure 8 The proportion of Gram-negative bacteria (blue line) continued to be higher than that of Gram-positive bacteria (red line) from March 2022 to June 2022. From September to December, the proportion of Gram-positive bacteria exceeded that of Gram-negative bacteria. In the clean room, the ratio of Gram-negative bacteria to Gram-positive bacteria reflects the quality of clean room control. Negative bacteria mean a higher risk of pathogenicity and a humid clean room environment.

[0038] S8 uses the matrix of genus, species and time unit (e.g., "monthly") to aggregate the summary information table obtained in S1, and calculates the Shannon index of the bacterial community for each month. The calculation formula is as follows: H = -Σ(p(i) * log(p(i))), where p(i) is the relative abundance of genus / species i in the microbial community. The calculation formula of Simpson index is: D = 1 -Σ(p(i) * p(i)), and the calculation method of goods_coverage is 1-Σ(x==1) / -Σ(x), where x is the abundance of genus / species in the microbial community. Visualize the trend of changes over the month in the form of a line graph ( Fig. 9 ). Fig. 9 The example shows that the Shannon index (red line) is between 2 and 5 from January to December, which is relatively stable. It drops significantly in August 2022, when goods_coverage also drops, indicating a decrease in diversity and coverage index due to a decrease in the number of samples.

[0039] S9 classifies the summary information table obtained in S1 according to different meta-information (here taking the room field as an example), calculates the species / genus intersection between different classifications, and helps us understand the connection between different groups. For example, the common and unique detected microorganisms between different cleanliness level areas, different rooms, and different sources are visualized in the form of a collection diagram ( Fig.10 ). In the figure, the connected points below the X-axis indicate that these categories have commonly detected microorganisms. For example, in the example, there are 6 unique microorganisms in CXC, 3 microorganisms that are only shared by CXC and DN, 1 microorganism shared by CXC and LABQ, and 1 microorganism that is detected in 6 areas including Workshop 209, QC2, CXC, NXM Workshop, DN, and R300PW Station 2. Set analysis can reflect the wide distribution of genera and the similarity of bacterial communities between different areas, thereby reflecting the contamination path.

[0040] S10 summarizes the genera and species by month, and classifies the detected bacterial communities according to the distribution and continuity of the detection time. The core bacterial community refers to the genera that appear at more than 90% of the time points and are more than 80% continuous; the persistent bacterial community is defined as the genera that appear at more than 60% of the time points and are more than 80% continuous; the transient bacterial community is defined as the genera detected except the core bacterial community and the persistent bacterial community. The detection situation and the type of bacterial genus are visualized in the form of a heat map ( Fig.11 ). Fig.11 It shows that the Ralstonia genus was detected continuously from March to August 2022, accounting for 6 / 10*100%=60% of the months when samples were sent, and the detection time was 100% continuous, which meets the definition of persistent flora and is marked as "Persistent" on the far left of the heat map. The continuous detection of flora indicates that colonization has occurred in the environment, or that the production environment has failed to resolve the source of pollution. Transient flora may be caused by accidental factors entering the system, not by system control problems. This analysis is used to distinguish between system problems and accidental factors.

[0041] S11 aggregates the summary information table obtained in S1 by genus, species, and time unit (e.g., by "month"), and further filters it according to the meta-information. The data after classification and screening are reduced in dimension through principal component analysis and visualized in the form of scatter plots and confidence ellipses ( Fig.12 ). Fig.12 The Euclidean distance of the bacterial communities between different sources (ID, PW, WFI, sedimentation bacteria, floating bacteria, bioindicators, and NA classification) was analyzed, reflecting the differences in bacterial communities from different sources, and the differences in bacterial communities from different months from the same source. For example, the purple diamond letter ellipse representing sedimentation bacteria is small, indicating that there is little difference between months, but there are two points on the vertical axis that are far away from the central point, indicating that there are two months with relatively large differences in bacterial community structure, suggesting that the investigation of related environmental factors is needed. In addition, the R package corrplot was used to analyze the similarity of the detected genera / species between different months, sources, classifications, and environments, and the correlation between different species / genera was analyzed by calculating the correlation coefficient, and a correlation coefficient heat map was drawn to display ( Fig.13 ). Fig.13 The results showed that the correlation coefficients of Phyllobacterium, Pelomonas, Afipia, and Sphingomonas genera were close to 1, and the significance test was less than 0.001, indicating that the detection of these genera was significantly positively correlated, that is, they were detected together, representing the pattern of these bacterial communities in the environment, that is, they may coexist in biofilms or have common carriers. The information is used for inference of environmental pollution investigations.

[0042] The above description is only a preferred embodiment of the present invention, and does not limit the implementation mode and protection scope of the present invention. For those skilled in the art, it should be aware that all solutions obtained by equivalent substitutions and obvious changes made using the description and illustrations of the present invention should be included in the protection scope of the present invention.

Claims

1. A trend analysis method based on the identification results of production environment microorganisms, characterized in that: Here are the steps: S1 summarizes the sampling meta-information, identification results, and phenotypic characteristics of the strain to form one data record for one identification result; S2 performs matrix transformation on the summary information obtained by S1, aggregates the data according to the specified time unit, and forms an aggregate matrix; S3 performs genus / species-level proportion statistics on the aggregation matrix obtained in S2, and displays the top 10 species / genera in abundance; S4 The aggregation matrix obtained in S2 is analyzed by source and visualized using a Sange diagram to show the proportional distribution relationship between different classifications and genera / species; S5 uses the aggregate matrix obtained in S2 to perform trend analysis on the abundance of bacterial communities over time and displays it in a bubble chart; S6 Draw a trend analysis line graph for each genus / species of concern over time, and calculate the warning line and correction line based on the model fitting; S7 performs trend analysis on the aggregation matrix obtained in S2 according to phenotype classification, calculates the average value, and displays it in the form of a line graph; S8 uses the aggregation matrix obtained in S2 to calculate α-diversity indices such as Shannon index, Simpson index and coverage, and displays the trend of these indices over time through line graphs; S9 further converts the aggregation matrix obtained in S2 into detection frequency, analyzes the co-detection phenomenon of different bacterial genera between different classifications and the properties of their residence in the environment, and visualizes them through a collection diagram; S10 further converts the aggregation matrix obtained in S2 into detection frequency, analyzes the continuous detection of different bacterial genera within the analysis time range and under the specified source, displays it through a heat map, and identifies the three types of bacterial communities: core, persistent, and transient; S11 performs dimensionality reduction analysis on the aggregation matrix obtained in S2 through multivariate statistics to obtain the correlation and distance information between bacterial genera or species and between different classifications, and displays them visually.

2. The method according to claim 1, characterized in that The data record input in step S1 is based on the identification result of a single colony or pure culture, which is expressed in a Latin name. The sampling meta-information includes the sampling point, clean room level, sampling method, plan category and sampling time information of the sample. The identification result information and the sampling meta-information are associated through the sample number, and the phenotypic information is obtained by querying the phenotypic database through the Latin name of the identification result.

3. The method according to claim 1, characterized in that In step S2, data is aggregated for different meta-information or phenotypic information classification fields using the dict and pandas functions in Python according to the specified time unit to generate an aggregation matrix of bacterial species abundance under different classifications.

4. The method according to claim 1, characterized in that In step S3, statistics are performed on the aggregation matrix according to all samples and according to the classification methods of room, clean room level, sampling method, and plan category, and the top 10 bacterial genera / species with abundance and proportion under different classifications are displayed using a Nightingale rose diagram.

5. The method according to claim 1, characterized in that In step S5, the trend analysis of bacterial abundance over time is performed, with time as the X-axis and bacterial genus / species as the Y-axis. The abundance and proportion of bacterial genus / species are mapped to the size and color depth of the bubble for visualization.

6. The method according to claim 1, characterized in that In step S6, a trend analysis line graph of changes over time is drawn for a single genus / species of concern. Warning limits and correction limits need to be calculated. The warning limits and correction limits are obtained by using the fitdistr function to perform maximum likelihood fitting of the data distribution and using the chisq.test method to perform a chi-square goodness of fit test. The best model distribution of the data is determined based on the p-value of the test and calculated using the QC_Lines or nptol.int method.

7. The method for analyzing bacterial flora trends based on the identification results of a single microorganism collected from a production environment according to claim 1, characterized in that: In step S7, trend analysis is performed according to phenotypic classification. It is necessary to classify and summarize the bacterial genera / species according to the phenotypic fields, calculate the mean, and draw a line graph of the change of phenotypic group abundance over time with time as the X-axis.

8. The method according to claim 1, characterized in that In step S8, the diversity index is calculated by using the R package vegan to calculate the community α-diversity index under the specified time unit and specified classification, including the Shannon index, Simpson index and good's coverage, and their means, and visualized through a line graph.

9. The method according to claim 1, characterized in that: In step S9, the abundance of genus / species in the aggregation matrix was normalized to "detected" or "not detected", thereby converting the aggregation matrix into a detection frequency matrix. The pandas function was used to perform a set analysis on different classifications, and the R package VennDetail was used to obtain the common bacterial genera / species between different classifications and displayed in a combination of graphs and tables.

10. The method according to claim 1, characterized in that In step S10, the aggregation matrix is ​​converted into a frequency matrix as in step S9, with time as the X-axis and the genus / species classification as the Y-axis. The abundance is mapped to color, and a heat map is used to display the continuity of the genus / species over time. In step S11, the aggregation matrix is ​​subjected to correlation analysis on the genus / species using the cor function in the R package corrplot, and the correlation between genera / species and the similarity between bacterial communities under different classifications are analyzed. The principal component analysis method is used to perform dimensionality reduction analysis on genera of different classifications, and the R package FactoMineR is used to draw scatter plots and confidence ellipses for visualization to identify outliers.

Images