A methylome and metagenomics analysis method

Abstract

The present invention relates to the field of bioinformatics, and in particular to a methylome and metagenomics analysis method. The method comprises: obtaining data to be analyzed after preprocessing the whole genome sulfite sequencing data of the sample to be tested; and performing methylation analysis, bacterial and fungal microecological analysis, and phage analysis on the data to be analyzed. The present invention combines methylome and microbial metagenomics to analyze the whole genome sulfite sequencing data, and can effectively screen biological factors related to sepsis identification and related analysis. These biological factors can be used in sepsis identification, sepsis infection site prediction, and mortality analysis, and are of great value in the field of sepsis detection.

Description

Technical Field

[0001] The present invention relates to the field of bioinformatics, and in particular to a methylome and metagenomics analysis method. Background Art

[0002] Sepsis is a systemic inflammatory response syndrome caused by infection, which is common in patients with severe trauma or infectious diseases. The causes of the disease include infections caused by bacteria, fungi, viruses and parasites, which lead to the body's inflammatory response and imbalance of immune regulation. Currently, there are limited treatment strategies for sepsis, which are mainly supportive rather than curative treatments, focusing on anti-infection, maintaining organ perfusion and hemodynamic stability, and supporting organ function.

[0003] Currently, due to the heterogeneity of sepsis, it is still extremely difficult to accurately diagnose sepsis and judge its prognosis, and delays in diagnosis and intervention can lead to high mortality. Currently, sepsis is identified by culture, quantitative PCR (qPCR), and various inflammatory markers such as procalcitonin and C-reactive protein. However, these methods are limited by a limited pathogen spectrum and low specificity. The most common assessment system currently used to evaluate the severity of organ dysfunction is the Sequential Organ Failure Assessment (SOFA) score.

[0004] Detecting sepsis by plasma cell-free DNA (cfDNA) is a sensitive and specific method. However, simultaneously capturing host methylation and plasma microecological information for sepsis diagnosis and patient mortality prediction remains a technical challenge. Summary of the invention

[0005] In order to solve the problems existing in the prior art, the present invention provides a methylomics and metagenomics analysis method (Multimodal Analysis of Methylomics and Metagenomics in One Tube, MAMMOT).

[0006] The present invention is based on single-sample whole genome bisulfite sequencing (WGBS), integrates methylomes and microbial metagenomics, and provides an analysis method that is particularly suitable for screening sepsis-related biomarkers. Based on this method, the present invention can simultaneously capture host methylation and plasma microecological information, and based on this screening, a large number of sepsis-related cfDNA host methylation sites and microbial species are obtained, and its accuracy has been verified experimentally.

[0007] In a first aspect, the present invention provides a methylome and metagenomics analysis method, comprising:

[0008] The data to be analyzed are obtained after preprocessing the whole genome bisulfite sequencing data of the sample to be tested;

[0009] Methylation analysis, bacterial and fungal microecological analysis, and bacteriophage analysis are performed on the data to be analyzed; the methylation analysis includes:

[0010] Using 400-600 bp sliding window and 150-250 overlap analysis, sites that meet the following conditions are identified as methylation sites: i) located in the promoter or exon region; ii) after deduplication, more than 3 times the number of reads confirm the existence of the site; iii) appearing in 40% of the samples in the same group; iv) the difference in methylation level is greater than 20%.

[0011] Furthermore, in the methylation analysis, for the measurement data that conforms to the normal distribution, the independent sample t-test is used; for the measurement data that does not meet the requirements of the parametric statistical method, the Wilcoxon rank sum test is used; for the categorical data, the comparison of two rates is used; for the correlation analysis of the statistical data that conforms to the normal distribution, the Pearson test is used, and for the correlation analysis of the statistical data that does not conform to the normal distribution, the Spearman test is used.

[0012] Furthermore, the preprocessing includes: removing low-quality reads, GC correction, and obtaining methylation data and microecological data;

[0013] The methylation data are used for subsequent methylation analysis, and the microecological data are used for subsequent bacterial and fungal microecological analysis and phage analysis after annotation and classification.

[0014] Furthermore, the bacterial and fungal microecological analysis includes:

[0015] The data to be analyzed are classified into species, and the composition and abundance of each species are estimated; then alpha diversity analysis and beta diversity analysis are performed.

[0016] Furthermore, the composition and abundance of each species were estimated using a Bayesian algorithm;

[0017] The α diversity analysis estimates the abundance and diversity of the microbial community, and the β diversity analysis estimates the differences in the composition of the microbial community between different samples and different component samples.

[0018] Furthermore, the phage analysis includes:

[0019] Constructing a phage database; aligning the data to be analyzed to the phage database, and retaining the aligned reads;

[0020] After generating Hifi reads, they were aligned to the methylated phage database, methylation detection was performed, duplicates were removed, and the sequences were restricted to appear in at least 100 phage genomes to obtain non-redundant motifs;

[0021] The non-redundant motifs were aligned to a bioinformatics database, the aligned reads were retained, and the phage families and abundances were analyzed.

[0022] Furthermore, the phage database is prepared by the following method:

[0023] The viral sequences were searched using phage and bacteriophage as search terms from large databases to obtain sequences of open reading frames, genes, proteins, regions and modified genomes, which were combined as a bacteriophage database.

[0024] Furthermore, the comparison with the phage database includes a first Blast and a second Blast;

[0025] The first Blast includes:

[0026] Using the phage database as a Blast database, the non-human reads in the data to be analyzed are compared, and the Blast conditions include: word size is 28, e value is 0.0005, elimination limit is 1, and matching reads are retained; then, reads that can match human sequences in these reads are removed to obtain the first comparison read combination;

[0027] The second Blast includes: using the phage database as the Blast database, aligning the first comparison reads combination, and the Blast conditions include: word size is 28, using 8 threads, e-value cutoff value is less than 0.0005, elimination limit is 10, and obtaining reads covering highly conserved sequences between phages.

[0028] The phage database mentioned above can use a self-constructed database, or a phage database in the prior art (such as the publicly available CPD at https: / / doi.org / 10.5281 / zenodo.7154236).

[0029] For the first time, the present invention generated a local phage BLAST database by self-building the database, then removed any sequences with homology to the human genome, and used secondary BLAST to capture reads covering highly conserved sequences between bacterial phages, and aligned the consistent reads (HIFI reads) generated by the PacBio CCS mode to 8,848 high-quality phage genomes. This largely eliminated the impact of methylation on the data.

[0030] In a second aspect, the present invention provides use of the method in screening biological factors for identifying sepsis.

[0031] Furthermore, the biological factors include: methylation sites, bacteria or bacteriophages.

[0032] The present invention has the following beneficial effects:

[0033] The present invention is based on single-sample whole-genome bisulfite sequencing data, integrates methylome and microbial metagenomics methods, and provides a methylome and metagenomics analysis method through a specific comparison method. The method has high practicality and accuracy. Based on this method, the present application has discovered a considerable number of sepsis-related cfDNA host methylation sites, as well as sepsis-related bacteria and phage species, which is of great significance in the field of sepsis detection. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] In order to more clearly illustrate the technical solutions in the present invention or the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0035] Figure 1 It is a schematic diagram of the process flow of the methylome and metagenomics analysis method provided in Example 1 of the present invention.

[0036] Figure 2 It is the sepsis-related plasma cfDNA methylation feature provided in Example 1 of the present invention; the figure is a heat map of the top 200 differentially expressed methylated regions from cfDNA, involving Sepsis vs Non-Sepsis vs Health; Death vs. Live; abdominal cavity, lungs, bloodstream, urinary tract and other infection sites; the color range of the heat map in the figure represents the normalized differential expression methylation degree Z value, ranging from +4 (red for high methylation) to -4 (blue for low methylation).

[0037] Figure 3 This is a volcano plot of differentially methylated regions provided in Example 1 of the present invention; in the figure, A: Sepsis vs. Health; B: Sepsis vs. Non_Sepsis; C: Abdominal vs. Pulmonary; D: Death vs. Live; red in the figure indicates high-methylated DMRs between groups, and blue indicates low-methylation.

[0038] Figure 4It is a graph of WGBS and mNGS results of microbial cfDNA for clinical pathogen identification in 48 validation cohorts provided in Example 1 of the present invention; wherein A: genomic abundance of microorganisms in WGBS identified by blood culture; B: genomic abundance of microorganisms in mNGS identified by blood culture; C: WGBS and BC positive overlap of different pathogens; D: mNGS and BC positive overlap of different pathogens; E: WGBS and mNGS positive overlap of different pathogens, the numbers in the histogram represent the corresponding cases of the pathogen, and Double means that the pathogen is identified by both methods at the same time.

[0039] Figure 5 : This is a graph of the alpha diversity of bacteria and fungi in plasma provided in Example 1 of the present invention; the graph compares the differences in the parameters of bacterial and fungal richness and uniformity levels at the six evaluation genus levels among the three groups, the colors represent the corresponding groups, and the ordinate values ​​represent the degree of the parameters; in the graph, A: Sepsis vs. Non_Sepsis vs. Health; B: Death vs. Live; C: Abdominal vs. Pulmonary vs. Blood vs. Urine.

[0040] Figure 6 This is an NMDS analysis diagram of β diversity provided in Example 1 of the present invention; the colors in the figure represent the corresponding groups, R2 represents the difference explanation, and the close spacing of dimensions indicates that the discrete degree of group samples is small; in the figure, A: Sepsis vs. Non_Sepsis vs. Health; B: Death vs. Live; C: Abdominal vs. Pulmonary vs. Blood vs. Urine.

[0041] Figure 7 It is a relative abundance characteristic diagram of the phage family provided in Example 1 of the present invention; wherein the ordinate is the relative abundance, and the abscissa is the different groups; in the figure A: the top 15 distribution of different groups; B: the top 15 distribution of different infection sites.

[0042] Figure 8 Schematic diagram of the changes in Escherichia phage characteristics in plasma to distinguish Sepsis from Non-Sepsis provided in Example 1 of the present invention; wherein A: average distribution of unique phages in Sepsis and Non-Sepsis groups by bacterial host genus; B: distribution of Escherichia phage families; C: top four Escherichia phage families with significant differences.

[0043] Fig. 9This is a schematic diagram of the changes in the characteristics of Enterobacter phages in plasma to distinguish sepsis from non-sepsis provided in Example 1 of the present invention; wherein A: the distribution of Enterobacter phage families; B: the top four phage families with significant differences are Escherichia coli.

[0044] Fig.10 Schematic diagram of the changes in phage characteristics in plasma for distinguishing lung and abdominal infection sites provided in Example 1 of the present invention; wherein A: average distribution of unique phages in lung infection and abdominal infection groups by bacterial host genus; B: distribution of phage families of the genus Klebsiella; C: top four phage families with significant differences in Klebsiella phage family.

[0045] Fig.11 Schematic diagram of the changes in phage characteristics in plasma provided in Example 1 of the present invention for distinguishing between survival and death groups; wherein, A: average distribution of unique phages in survival and death groups by bacterial host genus; B: distribution of phage families of Klebsiella genus; C: top four phage families with significant differences in Staphylococcus genus.

[0046] Fig.12 It is a scatter plot of the percentage of GC content of bacteriophages between groups provided in Example 1 of the present invention; in the figure, methylation represents the methylation treatment group, and nonmethylation represents the non-methylation control group.

[0047] Fig.13 is a scatter plot of the correlation between two groups of bacteriophages provided in Example 1 of the present invention; in the figure, AI: are Escherichia-virus-HK630 , Escherichia-virus-Lambda , Klebsiella-phage-KP179 , Klebsiella-phage-Mineola , Komagataeibacter-phage-phiKM1 , Mycobacterium-phage- Xula , Streptococcus-phage-142 , Tenacibaculum-virus-PTm1 , Uncultured-Phage ; The relative abundance of the two groups of samples is zero and no display is shown. The size of the correlation coefficient R reflects the consistency of the relative abundance between the phage groups.

[0048] Fig.14 is a scatter plot of the correlation between two groups of bacteriophage hosts provided in Example 1 of the present invention; in the figure, AI: is Escherichia , Enterobacter , Aeromonas , Klebsiella , Streptococcus , Staphylococcus , Lactobacillus , Campylobacter , Prochlorococcus ; The relative abundance of the two groups of samples is zero and no display is shown. The size of the correlation coefficient R reflects the consistency of the relative abundance between the phage groups. DETAILED DESCRIPTION

[0049] In order to make the purpose, technical solution and advantages of the present invention clearer, the technical solution of the present invention will be clearly and completely described below in conjunction with the drawings of the present invention. Obviously, the described embodiments are 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.

[0050] The experimental methods involved in the following examples, unless otherwise mentioned, are all conventional methods in the art, for example, reference may be made to experimental manuals in the art, or the conditions recommended in the manufacturer's instructions.

[0051] Unless otherwise specified, the experimental materials and reagents involved in the following examples can be obtained from commercial sources.

[0052] Example 1

[0053] The present invention provides a highly sensitive, highly specific, and accurate single-tube multimodal methylome and metagenomics analysis method for sepsis (process as shown in the following example): Figure 1 shown), including:

[0054] 1. The whole genome bisulfite sequencing (WGBS) data of the existing cohort (including 12 healthy subjects, 40 ICU non-sepsis patients and 134 sepsis patients) were analyzed as follows:

[0055] 1.1 Quality analysis (data cleaning)

[0056] The tool used for quality analysis is Trim Galore 0.4.1 software. This software can remove reads with too low quality, unrecognized, and reads with short insert fragments that lead to the detection of adapter sequences. At the same time, the base distribution and quality of the data before and after screening, as well as the sequence length, are evaluated to delete low-quality reads (the number of bases with a quality value Q≤20 accounts for more than 50% of the total reads).

[0057] 1.2 GC correction and sequence alignment

[0058] (1) After aligning the sequencing sequence of each sample with hg38, 100-120M reads were obtained, and low-quality, multiple-matching, and incompletely matching sequences were removed again.

[0059] (2) The chromosome was divided into windows of 500 bp each, and the start and end positions of each window were recorded. The GC values ​​of all sequences in each window were counted, and the number of sequences in each window was calculated in 10 samples to obtain the average number of sequences in each window. Windows with an average number of sequences exceeding three standard deviations were removed.

[0060] (3) Finally, perform Loess fitting on the sequence number and GC value of all windows. The calculation formula is as follows: f(GC, n) = Loess(n ik ~GC ik ), where i=1,2,3......22,X,Y; n ik is the original sequence number of the kth window of the i-th chromosome, GC ik represents the GC value of the kth window of the ith chromosome, and f(GC, n) represents the fitting relationship between the window sequence number and the GC value.

[0061] According to the above fitting relationship, the GC value of each window is used to predict the number of new sequences, which is recorded as the predicted value. The calculation formula is as follows: ik = f(GC ik , n ik ). Where i=1,2,3.....22,X,Y; pre ik represents the predicted value, f(GC ik , n ik ) represents the fit between the sequence number and GC value of the kth window on the i-th chromosome.

[0062] According to the predicted values ​​obtained above, the number of sequences after correction in each interval is calculated. The calculation formula is as follows: adjust_count=raw_count pre ik / n ik . Among them, i=1,2,3.....22, X, Y; adjust_count is the number of sequences after correction in the interval to be tested, raw_count is the number of original sequences in the interval to be tested, pre ik represents the predicted value, n ik Indicates the number of original sequences in the window where the interval to be tested is located.

[0063] (4) Use MethylDackel software to obtain methylation data. Use Bowtie2 to align the sequence to the host genome, and then remove the sequence aligned to the host to obtain microecological data.

[0064] 1.3 Data Processing

[0065] CpG sites are randomly distributed and discrete in the genome. A 500 bp window was used to slide and the methylation percentage was calculated to determine the methylation status of the site.

[0066] In order to study the species composition and diversity information of the sample microecology, Kraken2 software was used to annotate and classify the valid sequences of all samples.

[0067] 1.4 Statistics and methylation analysis

[0068] SPSS22.0 and GraphPad Prism8.0 software were used, and independent samples t The test is used for the statistics of normally distributed quantitative data. If it meets the non-parametric conditions, the Wilcoxon rank sum test is used. In the statistical inference of categorical data, such as the comparison of two rates, the chi-square test is used. The Person test is used for the correlation analysis of statistical data that meets the normal distribution, and the Spearman test is used for non-normal distribution. P <0.05 indicates a significant difference. Differential methylation site detection logic (sites that meet the following conditions are considered differential methylation sites):

[0069] (1) After deduplication, more than 3 times the number of reads confirms the existence of the site;

[0070] (2) 40% of the sites appeared in the same group and were confirmed to be included;

[0071] (3) 500 bp window sliding and 200 bp overlap analysis; methylation degree difference > 20%;

[0072] (4) Only methylation occurred in promoter and exon regions.

[0073] Heatmaps were created using Pheatmap software, and volcano maps were visualized using the R package ggplot2.

[0074] 1.5 Bacterial and fungal microecological analysis

[0075] (1) Species classification: Bracken software classified the classification results obtained by Kraken2 and used the Bayesian algorithm to re-estimate the species composition abundance of the metagenomic samples. The top 15 species in terms of abundance were visualized using the R package ggplot2, and the statistical test was the Wilcoxon rank sum test.

[0076] (2) α diversity analysis: α diversity reflects the abundance and diversity of microbial communities, including a series of statistical analysis indices to estimate the species abundance and diversity of each sample. Richness: refers to the number of species in a community. Chao: is an index that uses the Chao1 algorithm to estimate the number of species contained in a sample. The Chao index is used to evaluate the number of species in a sample. The larger the Chao index, the more species there are in the sample. Ace: is one of the commonly used indices for estimating the total number of species in ecology, which is different from the Chao1 algorithm. Shannon: is used to estimate one of the microbial diversity indices in a sample. The larger the Shannon value, the higher the community diversity. Simpson: is used to estimate one of the microbial diversity indices in a sample, and is often used in ecology to quantitatively describe the biodiversity of a region. The present invention uses the R package Vegan package to calculate the α diversity index of each sample.

[0077] (3) β diversity analysis: β diversity is a comparative analysis of the microbial community composition of different samples and samples from different groups. The dimensionality reduction analysis used in this study is non-metric multidimensional scaling (NMDS) analysis. The Bray Curtis distance between samples was calculated using the R package vegan and the distance between samples was calculated using the phyloseq package. The principal coordinate combination with the largest contribution rate was selected for plotting. The closer the distance between the samples in the figure, the more similar the species composition structure of the samples. The calculation was performed using adonis P The statistical method was Wilcoxon rank sum test.

[0078] 1.6 Phage analysis

[0079] (1) Construction of bacteriophage database: The comprehensive phage database (CPD) was constructed based on phage genomes retrieved from the nucleotide database of the National Center of Biotechnology Information (NCBI) (the search was described below); the search was limited to viral sequences, and "phage" or "bacteriophage" was used as the search term. Repeated sequences or sequences containing only open reading frames, genes, proteins, regions, or modified genomes were excluded. A total of 26,159 phage sequences were generated, which were merged into a single FASTA file and used to generate a local phage BLAST database.

[0080] The cleaned and processed reads were annotated by BLAST alignment search against non-human reads, using CPD as BLAST database (used as '-db'), word size of 28, e-value of 0.0005 and rejection limit of 1. The output format was set to '6 qseqid sseqid pident length evalue stitle'. Reads that had matches to phage sequences were then extracted from the FASTA file. This was done using the command 'awk '{print $1}' blast_output.txt>tokeep.txt' which extracted the sample read ID ('qseqid') from the BLAST output table. These reads were then retained using Seqkit, resulting in a FASTA file containing only potential phage hits, using the command 'seqkitgrep --pattern-file tokeep.txt sample.fasta>subset_sample.fasta'. A BLAST alignment search was then performed for all reads with bacteriophage hits below the significance threshold, restricting the subset of the NCBI nucleotide database to those with human taxonomic IDs to identify bacteriophage reads that matched human homology. This was done using a BLAST database created from sequences retrieved on November 28, 2021 with taxid 9606. The BLAST search was run with an e-value cutoff of 0.0005, a rejection limit of 1, a number of threads of 8, and a word size of 28. The output format was set using '6qseqid sseqid pident length evalue stitle'. Reads with both bacteriophage hits and matches to human sequences with an e-value less than or equal to 0.0005 were then removed from the FASTA file using the command 'awk '{print $1}' blast_output.txt>toremove.txt', which extracted the sample read ID ('qseqid') from the BLAST output table. These reads were then removed using Seqkit, resulting in a FASTA file containing only sequences that aligned to phage sequences, removing any sequences with human genome homology using the command 'seqkit grep -v --pattern-file toremove.txt subset_sample.fasta>depleted_subset_sample.fasta'.

[0081] (2) Second BLAST: The remaining sequences were subjected to a second BLAST search against the CPD with a word size of 28, using 8 threads, an e-value cutoff of less than 0.0005, and a rejection limit of 10 to capture potential identities of reads covering highly conserved sequences among bacterial phages. The output includes qseqid, sseqid, stitle, pident, length, and e-value. The CPD includes phages and their bacterial hosts and was created by annotation from the Virus-Host DB40, phage NCBI classification entries, phage names, and manual curation of manuscripts associated with phage sequence submissions. Phages with only genus-level host details were annotated. Host bacterial identity was characterized by literature search. Accession numbers of enteric phages from the Cenote Human ViromeDatabase32 were used to identify the proportion of annotated, uncharacterized phages that have been found to be enteric phages.

[0082] (3) Methylation motifs: The construction of a methylated bacterial phage database refers to the methylation sites of 8,848 human phage genomes sequenced by single molecule real-time sequencing (SMRT).

[0083] Consensus reads (HIFI reads) generated by PacBio CCS mode were aligned to 8848 high-quality phage genomes (Chinese Human Gut Virome-High-Quality, CHGV-HQ), and then methylated bases were detected using the SMRTlink tool, which compares the reference genome with the unmethylated bases and then the average pulse interval ratio (IPDr) of all subreads belonging to the aligned HIFI reads at a certain position. In addition, methylation motifs were identified using MEME-ChIP[25b (v5.4.1; default parameters). The motifs were then deduplicated using a custom R pipeline (https: / / github.com / whchenlab / gutphagemethylome). To further remove false positives, these motifs were also required to be found in at least 100 phage genomes. A total of 157 non-redundant motifs were identified. These motifs were compared with those in REBASE using a custom R pipeline. Motifs that were considered to be one of the possibilities of motifs in public databases were considered to overlap. A total of 155 were identical to motifs in REBASE and were finally identified as the human phage methylation reference genome.

[0084] (4) Phage and host annotation: Phage sequences were analyzed using R. Phage access IDs were first extracted and the corresponding taxonomic lineage information was assigned using the taxize v0.9.99 package. The R package vegan v2.5-7 calculated phage diversity using the Simpson's diversity index (SDI). The R code output from CPD and blast was then used to interpret the representativeness of the phage group for phage families and host genera / characteristics (the code is available at https: / / doi.org / 10.5281 / zenodo.7734114). The R package factoextra v1.0.6 was used to generate the Pearson dissimilarity matrix. Genetically defined taxonomic lines were used from publicly available databases for alternative analysis. For example, a subset of E. coli phages was analyzed using gene-defined lines from the INPHARED database v1.7 instead of the phage family identities present in CPD. The sample composition abundance was re-estimated using a Bayesian algorithm. The top 15 species in terms of abundance were visualized using the R package ggplot2, and the statistical test was the Wilcoxon rank sum test.

[0085] Table 1 Clinical characteristics of 12 healthy subjects, 40 ICU non-sepsis patients, and 134 sepsis patients

[0086]

[0087] Note: COPD: chronic obstructive pulmonary disease, SOFA: Sequential Organ Failure Assessment, : Indicates Sepsis vs. Non_Sepsis, # : stands for Sepsis vs. Health.

[0088] Table 2 Clinical characteristics of 48 patients validated by MAMMOT

[0089]

[0090] 2. Experimental results

[0091] (1) Differentially methylated regions of plasma cfDNA

[0092] The heat map of the top 200 DMRs showed a clustering relationship between the grouped samples. Compared with Health and Non_Sepsis, the methylation pattern of plasma cfDNA in the sepsis group changed. At the same time, the methylation patterns of plasma cfDNA in each infection site of sepsis and the clinical outcomes of patients within 28 days were also different ( Figure 2 ). Compared with healthy subjects, a total of 9901 hypermethylated DMRs and 23803 hypomethylated DMRs were identified in the sepsis group; compared with ICU non-sepsis patients, a total of 3112 hypermethylated DMRs and 5386 hypomethylated DMRs were identified in the sepsis group; compared with the lung infection site group, a total of 22479 hypermethylated DMRs and 8549 hypomethylated DMRs were identified in the abdominal infection site group; compared with the survival group, a total of 10238 hypermethylated DMRs and 14144 hypomethylated DMRs were identified in the death group ( Figure 3 ).

[0093] (2) Sepsis diagnosis and mortality prediction based on bacterial and fungal characteristics in plasma

[0094] The present invention uses 12 healthy cohorts as controls and seeks the optimal auc threshold for all 48 experimental sample blood culture positive species. This threshold is used as the standard, and microbial detection is positive when it is greater than this threshold. The present invention compares the genomic abundance of bacterial cfDNA in patients in the methylation group (48 samples matched with a positive clean-catch blood culture, WGBS group) ( Figure 4 Similarly, the present invention compared the genomic abundance of bacterial cfDNA in patients in the non-methylated group (48 samples matched with a positive clean-catch blood culture, mNGS group) ( Figure 4 B in the figure). When blood culture was used as the gold standard, both the WGBS and mNGS methods of the present invention identified a large number of Escherichia coli. However, when it comes to consistency with blood culture, WGBS is superior to mNGS (87.50%>62.50%, chi-square test, two-sided P=0.0047). The present invention further evaluated the consistency of the number of samples corresponding to blood culture pathogens by different methods. Bacillus megaterium, Bacteroides vulgatus and Burkholderia contaminans were not identified by the WGBS and mNGS methods. Escherichia coli and Klebsiella pneumoniae were identified the most in both methods ( Figure 4 C and D in ). This may be due to the limitations of sensitivity and specificity of conventional blood culture. The method of the present invention is better than mNGS, but the consistency is poor. Although WGBS and mNGS are both from the same batch of samples and constructed using single-stranded libraries, the present invention changes the sequencing depth based on the needs of data mining.

[0095] (3) Sepsis diagnosis and mortality prediction based on bacterial and fungal characteristics in plasma

[0096] i) Alpha diversity analysis

[0097] The study found that richness, ACE and Chao1 indicators can distinguish sepsis from ICU non-sepsis or sepsis from healthy physical examination groups, and the differences are statistically significant ( P <0.0001)( Figure 5 A in the figure), and these three indicators can also distinguish the survival group from the death group ( P <0.05)( Figure 5 Simpson was able to distinguish between sepsis, ICU non-sepsis, and healthy subjects ( P <0.05, P <0.0001). However, for the classification of abdominal and lung infection sites, all six indicators showed no significance ( Figure 5 C in.

[0098] ii) β-diversity analysis

[0099] NMDS analysis explained the difference R in sepsis, ICU non-sepsis and healthy subjects 2 The difference was statistically significant (7%) P =0.001)( Figure 6 However, in the death and survival groups, R 2 The difference was 0.7%, which was not statistically significant ( P =0.36). In the abdominal, lung, bloodstream, and urinary infection site groups, R 2 The difference was statistically significant ( P =0.002).

[0100] (4) Sepsis diagnosis and mortality prediction based on phage characteristics in plasma

[0101] i) Distribution of phages in plasma among different groups

[0102] The present invention explores the distribution characteristics of phages in different populations and infection sites in plasma in order to subsequently find appropriate biomarkers. Figure 7 The top 15 representative phage families in terms of relative abundance in each population group are shown. In each population, abundant phage communities were identified, and the phage family with the highest relative abundance was Siphoviridae, which is consistent with previous studies.

[0103] ii) Alterations in phage signatures in plasma distinguish Sepsis from Non-Sepsis

[0104] In the previous study, this method has analyzed a large amount of fungal, host bacteria and phage information. Therefore, the unique phages in the Sepsis and Non_Sepsis groups were first distributed evenly according to the bacterial host genus ( Figure 8 Brown-Forsythe and Welch two-sided Dunnett's T3 multiple comparison tests were used to test the relative abundance of phage families between the two groups. It was found that the most significant difference between the two groups was in the genus Enterobacter ( Enterobacter ) and Escherichia coli ( Escherichia ) and were enriched in the Non_Sepsis group. Then, further analysis was performed on the genera Enterobacter and Escherichia. Siphoviridae and Myoviridae The enrichment number was the highest in both groups, and Siphoviridae, Myoviridae, Podoviridae and Autographiviridae There was a statistically significant difference in relative abundance between the two groups ( Figure 8 B, C in ). Similarly, this has been found in the Enterobacter phage family ( Fig. 9 ), P <0.001; P <0.0001; P <0.00001.

[0105] iii) Alterations in phage signatures in plasma distinguish between pulmonary and peritoneal infection sites

[0106] In the above, the present invention comprehensively evaluates the characteristics of phage families in severe infection populations. In the early stage of infection, determining the correct infection site is very important for the treatment of severe sepsis and septic shock. Therefore, the unique phages in the lung infection group and the abdominal infection group were evenly distributed according to the bacterial host genus ( Fig.10 A in Figure 1). Multiple comparisons were performed to test the abundance of phage families between the two groups. The genus with statistically significant differences was Klebsiella ( Klebsiella ) Siphoviridae , P＜ 0.05 。

[0107] iv) Distribution of phages in plasma in the survival and death groups at 28 days

[0108] Finally, the patients were divided into two groups according to their clinical outcomes, and the unique phage families were evenly distributed according to the bacterial host genus ( Fig.11 A in Figure 1). Multiple comparisons were performed to test the relative abundance of phage families between the two groups. The difference was statistically significant for Staphylococcus Herelleviridae, P < 0.05.

[0109] (5) Verification of the composition of bacteriophages and their host bacteria in microbial metagenomes

[0110] i) GC content of sample sequence

[0111] In WGBS, sulfite treatment converts unmethylated cytosine (C) into uracil (U), thereby converting the cytosine band into the uracil band, and then detecting the DNA methylation level, while the methylated cytosine (C) remains unchanged. Comparing the overall phage GC content of the two groups of samples, the unmethylated group was significantly higher than the methylated group, and the difference was significant (P < 5.9e-5), see Fig.12 The few individual differences may be due to the bias of the nucleotide database reflected by CPD in MAMMOT in this study towards the enrichment of human-related pathogens and their phages.

[0112] ii) Correlation of sample abundance between the two groups

[0113] a) Correlation of relative abundance of phages

[0114] like Fig.13 As shown, the Random Forest Importance with good classification performance in the second part of phages is: Escherichia-virus-HK630 ( P =4.2e -5 )、 Escherichia-virus-Lambda ( P =1.6e -14 )、 Klebsiella-phage-KP179 ( P <2.2e -16 )、 Klebsiella-phage-Mineola ( P =0.00042)、 Komagataeibacter-phage-phiKM1 ( P =2.2e -8 )、 Mycobacterium-phage-Xula ( P =1.7e -8 )、 Streptococcus-phage-142 ( P =2.2e -16 )、 Tenacibaculum-virus-PTm1 ( P =3.4e -15 )、 Uncultured-Phage ( P =3e -19 The relative abundance of phages that have been reported but not yet cultured in pure culture was tested for linear correlation (Pearson test), and there was a linear relationship between the groups (the minimum R was 0.49 and the maximum was 0.87).

[0115] b) Correlation between relative abundance of host bacteria

[0116] like Fig.14 As shown, at the same time, the annotated phage host pathogens are: Escherichia , Enterobacter , Aeromonas , Klebsiella , Streptococcus, Staphylococcus , Lactobacillus , Campylobacter , Prochlorococcus There was a linear relationship between the two groups, and the consistency was good.

[0117] Example 2

[0118] Based on the method provided in Example 1, the present invention further obtains multiple biological factors that can be applied to sepsis diagnosis and patient mortality through machine learning such as forward feature selection, RF and LR, specifically:

[0119] 1. Based on MAMMOT, biological factors related to the diagnosis of sepsis and non-sepsis were obtained, including:

[0120] Escherichia-virus-Lambda、Escherichia-virus-HK630、Citrobacter- werkmanii, chr7:67292401-6729290 0 ( PMS2P4 ) 、 chrX:149525801-14952630 0 ( IDSP1 ) 、 Streptococcus-oralis, chr20:63537201-63537700 ( PTK6 ) 、Enterobacteriaceae- endosymbiont-of-Donacia-fulgens, chr10:42333001-42333500 ( LOC441666 ) 、 Paenibacillus-psychroresistens、 chr3:48566001-48566500 ( COL7A1 ), chr12: 6772401-6772900 ( LAG3 )、 Campylobacter-phage-A11a、 chr9:13223601-13224100 ( MPDZ ).

[0121] The above DMSs and microbial relative abundances were used as the classifier training set: validation set = 7 / 3. The sensitivity was 86.77%, the specificity was 88.21%, the area under the curve (AUC) of the training set was 0.954 (95% CI: 0.944-0.965), and the area under the curve of the validation set was 0.948.

[0122] 2. Biological factors related to the abdominal and pulmonary infection sites of sepsis based on MAMMOT, including:

[0123] Bradyrhizobium-sp.-PSBB068、Uncultured-Phage (Phages that have been reported but not yet cultured in pure form) 、 chr11-65492201-65492700( LINC02736 ), chr8: 144052201-144052700 ( OPLAH ), chr13: 25095601-25096100 ( PABPC3 ), chr17: 81862801-81863300 ( LOC124904086 ), chr17: 10683601-10684100 ( SCO1 ), chr17: 10683801-10684300 ( SCO1 ).

[0124] The above DMSs and microbial relative abundance were used as classifiers. Training set: validation set = 7 / 3. The sensitivity was 91.00%, the specificity was 80.39%, the area under the curve (AUC) of the training set was 0.939 (95% CI: 0.924-0.955), and the area under the curve of the validation set was 0.913.

[0125] 3. Biological factors related to patient mortality prediction based on MAMMOT include:

[0126] chr12:10052601-10053100 ( CLEC9A) , chr2: 127566601-127567100 ( MYO7B ), chr1: 25240201-25240700 ( LOC124903881 ), chr7: 149866801-149867300 ( ZNF862 ) 、 Mycobacterium-phage-Hail、Winogradskyella-forsetii、 chr17:44364001-44364500 ( FAM171A2 ) 、Lachnoclostridium-phytofermentans、Streptococcus-phage-142、 chr14:3911360139114100 ( GEMIN2 ), chr22: 29882601-29883100 ( MTMR3 ) 、Mycobacterium-phage- Xula、Tenacibaculum-virus-PTm1、Komagataeibacter-phage-phiKM1 .

[0127] The above DMSs and microbial relative abundance were used as classifiers. Training set: validation set = 7 / 3. The sensitivity was 83.50%, the specificity was 85.34%, the area under the curve (AUC) of the training set was 0.921 (95% CI: 0.904-0.938), and the area under the curve of the validation set was 0.915.

[0128] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. Application of a methylome and metagenomics analysis method in screening biological factors for sepsis identification, characterized in that: include: The data to be analyzed are obtained after preprocessing the whole genome bisulfite sequencing data of the sample to be tested; Performing methylation analysis, bacterial and fungal microecological analysis, and bacteriophage analysis on the data to be analyzed; The methylation analysis includes: Using 400-600 bp window sliding and 150-250 bp overlap analysis, sites that meet the following conditions are identified as methylation sites: i) located in the promoter or exon region; ii) after deduplication, more than 3 times the reads confirm the existence of the site; iii) appearing in 40% of the samples in the same group; iv) the difference in methylation level is greater than 20%; The preprocessing includes: removing low-quality reads, GC correction, and obtaining methylation data and microecological data; the methylation data is used for subsequent methylation analysis, and the microecological data is used for subsequent bacterial and fungal microecological analysis and phage analysis after annotation and classification; In the methylation analysis, for measurement data that conform to normal distribution, an independent sample t test is used; for measurement data that do not meet the requirements of parametric statistical methods, a Wilcoxon rank sum test is used; for categorical data, a comparison of two rates is used; for correlation analysis of statistical data that conform to normal distribution, a Pearson test is used, and for correlation analysis of statistical data that do not conform to normal distribution, a Spearman test is used; The bacterial and fungal microecological analysis includes: classifying the data to be analyzed by species, and estimating the composition and abundance of each species by using a Bayesian algorithm; then performing alpha diversity analysis and beta diversity analysis; the alpha diversity analysis estimates the abundance and diversity of the microbial community, and the beta diversity analysis estimates the differences in the composition of the microbial community of different samples and different component samples; The phage analysis comprises: constructing a phage database; aligning the data to be analyzed to the phage database, retaining the aligned reads; generating Hifi reads and then aligning them to a methylated phage database, performing methylation detection, deduplication and limiting the occurrence in at least 100 phage genomes to obtain non-redundant motifs; aligning the non-redundant motifs to a bioinformatics database, retaining the aligned motifs, and analyzing phage families and abundance; the phage database is prepared by the following method: searching in virus sequences using phage and bacteriophage as search terms from a large database to obtain sequences of open reading frames, genes, proteins, regions and modified genomes, and merging them as a phage database; The comparison of the phage database includes a first Blast and a second Blast; the first Blast includes: using the phage database as a Blast database, comparing the non-human reads in the data to be analyzed, and the Blast conditions include: a word size of 28, an e value of 0.0005, a removal limit of 1, and retaining the matched reads; then removing the reads that can match the human sequence in these reads to obtain a first comparison read combination; the second Blast includes: using the phage database as a Blast database, comparing the first comparison read combination, and the Blast conditions include: a word size of 28, using 8 threads, an e value cutoff value less than 0.0005, a removal limit of 10, and obtaining reads covering highly conserved sequences between phages.

Images