Clinical detection of human probiotics and human pathogenic bacteria by using a novel high-throughput platform based on next generation sequencing

Abstract

Background

The human body plays host to a vast array of bacteria, found in oral cavities, skin, gastrointestinal tract and the vagina. Some bacteria are harmful while others are beneficial to the host. Despite the availability of many methods to identify bacteria, most of them are only applicable to specific and cultivable bacteria and are also tedious. Based on high throughput sequencing technology, this work derives 16S rRNA sequences of bacteria and analyzes probiotics and pathogens species.

Results

We constructed a database that recorded the species of probiotics and pathogens from literature, along with a modified Smith-Waterman algorithm for assigning the taxonomy of the sequenced 16S rRNA sequences. We also constructed a bacteria disease risk model for seven diseases based on 98 samples. Applicability of the proposed platform is demonstrated by collecting the microbiome in human gut of 13 samples.

Conclusions

The proposed platform provides a relatively easy means of identifying a certain amount of bacteria and their species (including uncultivable pathogens) for clinical microbiology applications. That is, detecting how probiotics and pathogens inhabit humans and how affect their health can significantly contribute to develop a diagnosis and treatment method.

Background

High throughput sequencing can analyze a large amount of sequences, enabling sequencing of 16S rRNA to identify complex bacteria species of pathogens and probiotic bacteria. Many naturally occurring bacteria form a complex population in the environment. The human body plays host to a vast array of bacteria, found in oral cavities, skin, gastrointestinal tract and the vagina. Some bacteria are harmful while others are beneficial to the host.

A pathogen is a microorganism that causes disease in its host. For example, bacterial pathogen include Corynebacterium diphtheria (causes diphtheria), Listeria monocytogenes (causes food poisons), and Legionella pneumophila (causes Legionnaires’ disease). Probiotics, another microorganism, benefit the host and has received considerable attention in recent years. A FAO report in 2001 [1] cited the advantages of probiotics as increasing immunity [2, 3], reducing gastrointestinal discomfort [4, 5], and protecting the flora within urogenital tract [6]. As is well known, probiotics can ameliorate symptoms of diseases [7] and reduce the risk of suffering from diseases [8, 9].

Despite the availability of many approaches to identify probiotics and pathogens, most of them are only applicable to specific and cultivable bacteria but time consuming. For instance, conventional methods detect growth of cultured bacteria in approximately two days, or an additional five days to obtain no-growth culture results [10], which is laborious. Besides, some bacteria cannot be cultured [11], subsequently increasing the difficulty of specifying pathogenic bacteria. Moreover, it is hard to determine whether an infection is caused by one or more bacteria types.

16S rRNA sequences, capable of identifying bacteria on a molecular level, can detect uncultivable bacteria [12]. Use of 16S rRNA sequencing can overcome some problems of conventional culture method [13]. Although 16S rRNA sequencing is a more effective means of identifying bacteria than conventional culture method, 16S rRNA sequencing takes a considerable amount of time in amplifying DNA sequences [14]. Sanger sequencing known as “first-generation” or “conventional” sequencing has been used for DNA sequencing for almost two decades. Next generation sequencing (NGS) can analyze large-scale sequences quicker, enable massively parallel analysis, reduce reagent costs and the size of sample components, and perform high throughput [15]. Thus NGS is more efficient than the Sanger method, which generates one read per sample. In addition, NGS of 16S rRNA more easily identify cultivable or uncultivable bacteria [12].

Because of the improvement of sequencing technology and Bioinformatics approaches, the accuracy in distinguishing bacteria with those methods has been increased. Based on high throughput sequencing technology, this work identifies 16S rRNA sequences of bacteria and analyzes bacteria species. High-throughput sequencing can sequence a large number of 16S rRNA sequence more efficiently; with high-throughput sequencing, researchers can acquire information to identify pathogens and probiotic bacteria [16–18].

Results

Platform application: gut probiotics and pathogens detection

The read statistics of quality filtering and taxonomy assignment are demonstrated in Table 1. Figure 1A illustrated the percentage of probiotics detected by the proposed platform. Table 2 listed the quantities (matched sequenced reads) of probiotics identified in the samples in the case study. The top three identified probiotics in 12 samples are Lactococcus salivarius, Streptococcus thermophilus, and Bifidobacterium longum. Figure 1B and Table 3 listed the proportion and quantities of pathogens, of which top three pathogens are Escherichia coli, Salmonella enteric, and Haemophilus influenza.

Table 1 Results of quality filtering and taxonomy assignment

Figure 1

Relative abundance of probiotics and pathogenic bacteria from human gut of all samples. (A) The percentage of probiotics was identified in the samples. (B) The proportion of pathogenic bacteria was identified in the samples in the case study.

Table 2 The quantities (matched sequenced reads) of probiotics identified in the samples in the case study

Table 3 The quantities (matched sequenced reads) of pathogens identified in the samples in the case study

Table 4 listed the results of disease risk evaluations. It showed that three diseases of two samples (B031 and B034) had similar distributions in the control group. The three diseases are obesity, colorectal cancer, and constipation. Sample B031 had reached the significance level with P-value 0.0333 and 0.0121 < 0.05 of distribution in constipation and colorectal cancer respectively compared to 98 samples control group using binomial test. Sample B034 had reached the significance level with P-value 0.00257 and 0.0121 < 0.05 of distribution in obesity and colorectal cancer. Evaluated by the association of bacterial risk markers and disease, the results suggested that these two samples had higher risk than 98 samples control group in constipation, colorectal cancer, and obesity. Their enterotypes of gut probiotics and pathogens may be one of risk factors which would cause disease.

Table 4 The result of disease risk evaluations of 12 samples

Reproducibility and accuracy evaluation of proposed platform

Two replicated experiments of four samples were performed to estimate the reproducibility of the proposed platform. The results of repeated experiments were consistent. The similarity between two repeated experiments were calculated by using UniFrac [19]. Results of each sample pair (replicate 1 and 2) closely resemble each other. The similarity of UniFrac distance of each sample pair is higher than 0.96 (0.9617 for B014, 0.9872 for B018, 0.9914 for B020, 0.9722 for B033). This implies that the analysis results are reproducible.

Next, accuracy of the platform is evaluated by adding Lactobacillus reuteri to a stool sample (B050). Sample B050 contains 24,408 assigned taxons, and Lactobacillus reuteri has no detected count. Whether the counts of this species in positive control sample (B050S_L) are elevated must be determined. Analysis results indicate that 27,113 taxons are detected in sample B050S_L. In fact, the detected counts of Lactobacillus reuteri in sample B050S_L are 1,430, and the percentage of Lactobacillus reuteri markedly increases from 0% to 5%.

In short, our platform is accurate and reproducible in terms of detecting the quantities of bacterial species of the proposed platform. The results evaluate the accuracy and feasibility of proposed platform in order to identify probiotics and pathogens. While requiring only about one day for detection, not limited in identifying certain bacteria, the proposed platform can detect and quantify multiple bacteria simultaneously.

Discussion

Because of the constraint of costs and technical limitations, 16S rRNA sequences obtained in most databases are partial sequences. Many studies thus assign taxonomy by using partial 16S rRNA sequences. In our probiotics and pathogens 16S rRNA sequence database, 17,964 sequences are collected from NCBI nucleotide database, NCBI 16S microbial rRNA database, Greengenes database, and SILVA. Our probiotics and pathogens 16S rRNA database contain less than 39% of 16S rRNA sequences which are longer than 1400 bps. Only 9% of the sequences are close to full length.

This work extracts the V4 region from full length 16S rRNA of microbiome in the human gut as a platform application. Some 16S rRNA variable regions are more dependable than other regions in assigning taxonomy like V3 and V4 [20, 21]; in addition, some 16S rRNA variable regions are much conserved. The proportion and diversity of probiotics and pathogens may be made diverse by using different 16S rRNA variable regions. The proposed platform is also applicable to other 16S rRNA variable regions for taxonomy assignment. Importantly, a more appropriate region than others must be selected to produce an outcome that is close to full length 16S rRNA sequence.

This work further attempt is to collect common probiotics and pathogens from the literature. Although it may be incomplete, recent advances in sequencing technology make it possible to identify and define an increasing number of bacteria, implying an obvious increase in the number of identified probiotics and pathogens in the future. Efforts are underway in our laboratory to update the list of used probiotics and pathogens.

Previous studies [22–24] identified pathogen or probiotic bacteria by using antibody, 16S rRNA gene microarrays, fluorescence in situ hybridization (FISH), and proteomic methods. In this work, the proposed platform can detect various pathogens and probiotics based on 16S rRNA (rDNA) sequences of bacteria using NGS and Bioinformatics method. An average of 126,451 reads was acquired per sample in this work. It is doubt that the sequencing depth is enough to detect a small amount of probiotics and pathogens. Although increasing the coverage of sequencing can advance the sensitivity of detecting probiotics and pathogens, the sequencing cost will increase. It is important to work out an appropriate coverage of sequencing for detecting probiotics and pathogens.

The results of disease risk evaluations revealed that most of 12 samples did not have resembled distributions of bacteria markers with control group. Only two samples had reached the significance level of distributions. The reason for the phenomenon may be the overlapped bacteria markers between diseases. 28 markers are used in colorectal cancer, and 17 markers are used in irritable bowel syndrome. Six markers are overlapped. For sample B031, the significant distributions in colorectal cancer were partly contributed to the significance in irritable bowel syndrome owing to the overlapped markers. Similarly, two overlapped markers for sample B034 were in colorectal cancer and obesity. In this kind of speculation, the influence of colorectal cancer to irritable bowel syndrome would be six (overlapped markers of CC and IBS) over seventeen (markers of IBS), and the influence of colorectal cancer to obesity would be two (overlapped markers of CC and obesity) over nine (markers of obesity). In addition, the influence of colorectal cancer to constipation and ulcerative colitis would be one over six and two over ten, respectively.

In addition to that some bacteria markers in species level are belong to the marker of genus level and species level, both genus marker and species markers may have associated with affecting the distributions mutually. Continually, collecting more markers and evaluating the distributions with markers in the same level are required for constructing a global prediction model in Taiwanese.

Conclusions

This work constructed a bacterial disease risk evaluation model for seven diseases and developed a novel platform by using NGS and Bioinformatics approach. Compared with the traditional bacteria culture method, our proposed platform can reduce experiment time. Besides, the proportion of probiotics and pathogens (including uncultivable pathogens) in the human body can be detected rapidly with 16 s RNA database of probiotics and pathogens. Furthermore, the proposed platform provides further insight into the cause of disease based on the relation of probiotics, pathogens, and disease. For instance, the type of antibiotics can be adjusted if the pathogens of disease are identified from infected patients. In addition, the proposed platform allows researchers to determine whether the intake of probiotics impacts the human body [25–29]. In the future, this preliminary study will be continuously extended for more bacterial disease markers. For more comprehensive applications, this work will also collect bacteria from other parts of human body as control group data. In fact, a detective method of how the probiotics and pathogens inhabit human can provide new insight for human health. It could improve diagnosis and treatment method.

Methods

Figure 2 illustrates the bioinformatics system flow of the proposed platform, which includes analysis pipeline of NGS. The Figure 2 contains four parts: sequence quality filtering, construction of bacteria sequence database, taxonomy assignment, and disease risk model evaluation. The detailed components in the proposed platform are described below.

Figure 2

System flow of bioinformatics analysis in the proposed platform. The proposed platform comprises the analysis pipeline of NGS, construction of probiotics and pathogens database, bacterial disease risk model evaluation and the application of individualized bacteria sequencing profile.

Sample collection

In this study, stool samples of 98 Taiwan volunteers were gathered. The samples were collected by Sigma-transwab (Medical Wire) into a tube with Liquid Amies Transport Medium, and stored at 4°C until processing.

DNA extraction

In the case study, fresh faeces were obtained from participants. DNA was extracted directly on stool samples by using a QIAamp DNA Stool Mini Kit (Qiagen). A swab was vortexed vigorously and incubated at room temperature for 1 min. The sample was then transferred to microcentrifuge tubes containing 560 μl Buffer ASL, vortexed, and incubated at 37°C for 30 min. In addition, the suspension was incubated at 95°C for 15 min, vortexed, and centrifuged at 14,000 rpm for 1 min into pellet stool particles. Extraction was performed following the protocol of the QIAamp DNA Stool Mini Kit. The DNA was eluted with 50 μl Buffer AE, and centrifuged at 14,000 rpm for 1 min. Moreover, the DNA extract was stored at-20°C until further analysis. Finally, DNA extraction was performed, depending on the sample collected.

Library construction and sequencing for V4 region of 16S ribosomal DNA

The PCR primers, F515 (5′-GTGCCAGCMGCCGCGGTAA-3′) and R806 (5′-GGACTACHVGGGTWTCTAAT-3′), were designed to amplify the V4 domain of bacterial 16S ribosomal DNA as described previously [30]. PCR amplification was performed in a 50 μl reaction volume containing 25 μl 2X Taq Master Mix (Thermo Scientific), 0.2 μM of each forward and reverse primer, and 20 ng DNA template. The reaction conditions consisted of an initial 95°C for 5 min, followed by 30 cycles of 95°C for 30 sec, 54°C for 1 min, and 72°C for 1 min, as well as a final extension of 72°C for 5 min. Next, amplified products were checked by 2% agarose gel electrophoresis and ethidium bromide staining. Amplicons were purified using the AMPure XP PCR Purification Kit (Agencourt), and quantified using Qubit dsDNA HS Assay Kit (Qubit) on Qubit 2.0 Fluorometer (Qubit)–all according to respective manufacturer instructions. For V4 library preparation, Illumina adapters were attached to the amplicons using the Illumina TruSeq DNA Sample Preparation v2 Kit. Purified libraries were applied for cluster generation and sequencing on the MiSeq system. The raw sequence files are available for download at http://clinic.mbc.nctu.edu.tw/.

16S rRNA (rDNA) sequence data quality filtering

The raw fastq files obtained by Illumina sequencing machine were quality-filtered using the FASTX-Toolkit^a. The paired-end 150 bp reads were performed using the minimum acceptable phred quality score of 20, as well as the 70% of bases that must exceed 20 phred quality score. Sequence shorter than 100 nucleotides would be omitted after quality trimming from reads tail. Notably, reads containing ambiguous characters were discarded.

Construction of probiotics and pathogens database

The list of probiotics and pathogens were obtained from literatures or the claims of official departments. Additional file 1: Table S1 lists species of probiotics which were adapted from both literatures [7, 9] and the claims of official departments, such as Taiwan Food and Drug Administration [31] and Health Canada [32]. 99 bacterial pathogens were collected from literature [25, 26, 33–42] and Taiwan Food and Drug Administration [31] (Additional file 1: Table S2).

The 16S rRNA sequences of probiotics and pathogens used for taxonomy mapping were retrieved from the NCBI nucleotide database, NCBI 16S microbial rRNA database, Greengenes database [43] and SILVA [44]. Following sequence data collection, we assemble partial sequences which used the same species classification and removed redundant sequences. Additionally, we also removed the unique sequence from only one research support with 3% similarity which shared the same species classification with other sequence.

Taxonomy mapping

To generate taxonomy assignments, the proposed platform invoked a modified Smith-Waterman algorithm from miRExpress [45], which can compare pairs of sequences in parallel, for mapping reads to taxons. miRExpress was designed for identifying the best similarity between sequencing reads and miRNA precursor sequences. In our model, it was modified for identifying multiple hits of 16S rRNA sequence mapping results with similarity threshold 0.97. In order to reduce the storage space of output, the SAM format [46] was used to replace the original miRExpress output format for storing alignment results. Furthermore, two kinds of output format were designed. One format records whole mapped sequencing reads based on taxons. The other one records which taxons could be assigned based on sequencing reads. These two kinds of output could support the important information for assigning sequencing reads to suitable taxon. miRExpress was originally designed for dealing with single-end sequencing data. Therefore, the additional program was added for processing paired-end sequencing data. In this part, both end sequencing reads need to be assigned to the same taxon. If paired-end sequencing reads were mapped to different taxons, this paired sequence would be dropped. The probiotics and pathogens 16S rRNA sequence from our database were built in FASTA format. Following quality filtering, all paired-end sequences were aligned to the probiotics and pathogens database with whole read aligned from one end to the other end. Reads were then truncated with an identity lower than 97%, according to previous research in order to achieve a better compromise between sequences from PCR sequencing errors and taxonomic relatedness [27].

The construction of Bacterial disease risk evaluation model (BDREM)

To study the associations between bacteria and diseases, we collected related information from literatures. We concerned bacteria that are associated with seven diseases: constipation [28, 47, 48], obesity [29, 49–52], irritable bowel syndrome (IBS) [28, 53–58], ulcerative colitis (UC) [53, 59–61], colon cancer (CC) [62–64], Atopic Dermatitis (AD) and Allergic rhinitis (AR), were collected positive correlation and negative correlation data, and the individual risk of disease was evaluated.

The association data were majorly collected from case–control studies which the quantities of bacteria were obtained from NGS data, and few well-known bacteria validated by multiple studies through cultural experiments were also included. We further eliminated some conflicted data with both positive and negative correlation between bacteria and disease in different studies.

Health Asians stool samples of 98 Taiwan volunteers were gathered. Following deep sequencing and sequencing data processing, the proportion of 78 bacteria from control group was applied as risk markers (constipation: 6, obesity: 9, IBS: 17, UC: 10, CC: 28, AD: 4, AR: 4) to predict disease risk to seven diseases in this study (Table 5).

Table 5 Disease-related biomarkers of seven diseases

The mathematical formula of BDREM in this study was developed as the following steps. Let λ be a N × S matrix, where N is the number of markers selected in the prediction model of constipation and S is the number of health subjects in 7 prediction models. T _i was defined as one of the two notches of median for each row of λ[65]. T _i is a threshold to distinguish λ _ij from normal proportion level to abnormal (fail to success in one trail of binomial distribution). Smaller notch was selected to T _i when each marker was recorded as a negative association to the disease, and a success trail was identified when λ _ij is smaller than T _i . On the opposite, larger notch was selected when association was positive, and a success trail was identified when λ _ij is larger than T _i .

$$\mathit{Ti}=\left\{\begin{array}{ll}\mathit{\text{Medain}}\mathit{\text{of}}& \hfill \left\{{\lambda }_{i1},{\lambda }_{i2},\dots ,{\lambda }_{\mathit{iS}}\right\}+\frac{1.58\times \mathit{\text{IQR}}}{\sqrt{S}},\hfill \\ \mathit{\text{positive}}\mathit{\text{association}}\\ \mathit{\text{Medain}}\mathit{\text{of}}& \hfill \left\{{\lambda }_{i1},{\lambda }_{i2},\dots ,{\lambda }_{\mathit{iS}}\right\}-\frac{1.58\times \mathit{\text{IQR}}}{\sqrt{S}},\hfill \\ \mathit{\text{negative}}\mathit{\text{association}}\end{array}\right.$$

Let P _j be the probability of successful trails in the j^th column of λ. The meaning of P _j is the personal probability that abnormal proportion level happened.

$$P_j=\frac{\#\mathit{\text{success}}\mathit{\text{trails}}\mathit{\text{in}}\text{the}\mathit{\text{jth}}\text{column}\text{of}\lambda .}{N}$$

Let P _h be the mean of P _j . It represents how frequent the abnormal proportion level happened to all P _j in average, regarded as the hypothesized probability of success in each P _j .

$$P_h=\frac{{\displaystyle {\sum }_{j=1}^S}{P}_j}{N}$$

Assume P _j obey a binomial distribution, and let P _h be the hypothesized probability (0.05051 for constipation, 0.07239 for obesity, 0.06952 for IBS, 0.05227 for UC, 0.09280 for CC, 0.04924 for AD, 0.05114 for AR). A binomial test was used to P _j and P _h . Alpha = 0.05 was choose to judge if a subject is significantly differently from the others in λ.

Figure 3 illustrated an example for evaluating the risk of obesity of B034 and B031. The model used lower and upper proportion bound of 9 markers from 98 control samples to define risk markers of these two samples following by using binomial test. Four markers of B034 exceed the lower bound and upper bound of obesity. The binomial test P-Value of B034 is 2.572e-03 < 0.05, Since P-Value < = hypothesized probability 0.07239, this case is specifically associated (significantly) with disease than random chance. There are two markers of B031 exceed lower bound of obesity. The P-Value of B031 is 0.1344 > 0.05, the case is no more associated with disease than random chance. As the results, we can assume that B034 had higher probability to cause Obesity.

Figure 3

An example for evaluating the risk of obesity by using bacterial disease risk evaluation model. The model used lower and upper proportion bound of 9 markers from 98 control samples to define risk markers of these two samples (B034 and B031) following by using binomial test.

Endnote

^ahttp://hannonlab.cshl.edu/fastx_toolkit/index.html.