4.2.4. 16S Microbiome Data Analysis
We analyzed the paired-end 16S V3-V4 raw sequencing read files (i.e., fastq.gz files) using two bioinformatics pipelines: QIIME2 and DADA2 (Version 1.2.0). The executing software for QIIME2 was Python software and for DADA2 was R (version 3.3.2). First, the sequencing reads were sorted using unique barcodes for each PCR product. The barcode, linker, and primer sequences were then removed from the original sequencing reads. The removed reads were merged paired-end reads using FLASH v 1.2.11. The merged reads containing two or more ambiguous nucleotides, those with a low-quality score (average score < 20), or reads shorter than 300 bp, were filtered out. All pre-filtered sequence data quality control and feature table construction were completed using the Divisive Amplicon Denoising Algorithm 2 (DADA2) in the QIIME 2 plugin, which detects and corrects amplicon errors and filters out PhiX chimeric sequences. The pre-processed reads from each sample were used to calculate the number of Amplicon Sequence Variants (ASVs). The number of ASVs was determined by clustering the sequences from each sample using a 99% sequence identity cut-off using DADA2 pipeline-based QIIME2 software. After the denoising step, the microbial diversity analysis (alpha and beta) was performed using the “diversity” category in the QIIME2 plugin. The alpha diversity was calculated by Observed_ASVs, Chao1, Shannon, Simpson, and Pielou_e alpha diversity indices, which were used to estimate microbial richness and evenness scores, and the Kruskal–Wallis nonparametric ANOVA test and the Mann–Whitney U test were applied to determine the statistical significance of each comparison analysis in the present skin microbiome study. Microbial beta-diversity between each comparative group was confirmed through distance dissimilarity calculated using Bray–Curtis and unweighted_UniFrac distance matrices. The permutational multivariate analysis of variance (PERMANOVA) non-parametric statistical test was used to determine the statistical significance of the beta-diversity analysis. The results of the microbial beta-diversity analysis were visualized through a principal coordinate analysis (PCoA) plot using R studio. Finally, taxonomic relative abundance was counted with a naïve Bayes classifier using a confidence threshold of 0.7%, derived from the ASVs data for each sample. The classified microbial composition based on SILVA v138 16S rRNA gene database was normalized using the value calculated from the taxonomy abundance count, divided by the number of pre-processed reads for each sample.