*2.5. Pyrosequencing Data Analysis*

Obtained raw data were categorized into samples according to the barcode information. The barcode, linker and primers were trimmed off from the reads and reads showing ambiguous nucleotides (2≤), a low-quality score (<25) or a short length (<300 bp) were removed. Potential chimera sequences detected by the Bellerophon method were filtered. The final reads were compared with sequence and taxonomic information acquired from the EzTaxon-e database. This database has curated 16S rRNA gene sequences of species and phylotypes of either cultured or uncultured entries in the GenBank database with taxonomy information. The read sequences were queried with BLASTN, and a maximum of five best-hit sequences were compared by pairwise alignment to determine species or phylotypes. Operational taxonomic units (OTUs) were obtained at a 97% similarity level by the CL community program (ChunLab, Inc., Seoul, Korea) after comparing with the EZ taxon database. Alpha diversity was determined using the "phyloseq" R package. A heatmap analysis was performed based on microbial abundance in the "heatplus" R package. Principal component analysis (PCA) was performed using "devtools" and "ggbiplot" R packages.

## *2.6. Statistical Analysis*

In clinical findings, dmft indices between the BS and control groups were statistically analyzed by Mann–Whitney test (*p* = 0.05), and the association between BS formation and caries development were statistically analyzed by Chi-squared test (*p* = 0.05). After the pyrosequencing, the diversity of the oral microbiome of the subjects was statistically analyzed by the *t*-test (*p* = 0.05). In the heatmap, the Bray–Curtis method was used for calculating cluster dissimilarity in the clustering analysis (*p* = 0.05). The abundance of the predominant bacterial groups was analyzed by *t*-test (*p* = 0.05). Jonckheere's trend test was used for the trend analysis (*p* = 0.05).
