Identifying a Correlation among Qualitative Non-Numeric Parameters in Natural Fish Microbe Dataset Using Machine Learning

Abstract

Recent technical innovations and developments in computer-based technology have enabled bioscience researchers to acquire comprehensive datasets and identify unique parameters within experimental datasets. However, field researchers may face the challenge that datasets exhibit few associations among any measurement results (e.g., from analytical instruments, phenotype observations as well as field environmental data), and may contain non-numerical, qualitative parameters, which make statistical analyses difficult. Here, we propose an advanced analysis scheme that combines two machine learning steps to mine association rules between non-numerical parameters. The aim of this analysis is to identify relationships between variables and enable the visualization of association rules from data of samples collected in the field, which have less correlations between genetic, physical, and non-numerical qualitative parameters. The analysis scheme presented here may increase the potential to identify important characteristics of big datasets.

1. Introduction

Technical innovation has enhanced progress in many fields of research and broadened the potential for data collection and analysis. The term “omics” has recently emerged, referring to exhaustive datasets of pools of biological molecules—for example, metabolomics, genomics, and proteomics. Next-generation sequencing is an example of a technical innovation that has contributed to genomics and has been used to elucidate the reasons that many species cannot grow in particular environments, leading to a deeper understanding of microbe communities that grow without cultivation [1,2,3]. Microbes exists in various locations, sometimes in symbiotic relationships with various animals [4]. While living on animals, the microorganisms interact with hosts or items that the host carries—for example, with food residue [5]. Commensal microbes are those which benefit from the host, while the host receives no benefit; they exhibit high cell turnover and produce high levels of secretions, including sweat and mucus [6,7,8]. Commensal microbes can, however, affect the host’s immunity and behavior—for example, via physical attachment and the release of bacterial components and metabolites [9,10,11,12]. Because of the effects that commensal microbes have on the host, exhaustive information on the microbial content of the host’s intestine may provide information relating to the status of health [13,14]. Moreover, to clarify the significance of not only commensal microbes but also the workings of nature for the host or human social and environmental benefits, humans might control nature in pursuit of the sustainable development goals (SDGs), particularly food consumption. However, in nature, a significant variety of parameters exist that cannot be controlled. Although technical innovation provides a relevant amount of information for us, we must develop appropriate methods to overcome the problems in accomplishing the SDGs.

Informatics is a constantly evolving technique with collateral performance advancements in computers, and various methods of calculating the unique characteristics of experimental datasets have been developed—for example, the similarity by clustering methods and order reduction for the overview of a dataset. In recent years, machine learning and deep learning have become essential in order to obtain an overview of big data, such as omics datasets [15,16]. Several complementary methods may be used simultaneously to solve complex data mining issues. Moreover, several methods were improved and their concept and application were extended [17,18].

Numerical data are typically used when the statistical analysis of experimental data is undertaken; when qualitative data are used, appropriate quantification must be performed before statistical analysis is possible. Thus, qualitative data, which cannot be easily converted into numerical data, present a significant challenge to statistical analysis [17].

The present study aimed to mine characteristic variance from crude omics data, and to identify association trends within subjective and qualitative data, which are difficult to convert to numerical data using informatics. We used two machine learning algorithms: K-means clustering and random forest [19,20]. The former is a well-known clustering analysis method that uses centroids to classify data coordinates, but the optimal number of groups for a particular dataset and the validity of separation cannot be known prior to the analysis. Random forest is an ensemble learning method for classification and regression. It is a high-performance method but requires datasets where some class information is already known—for example, control group, study group, and dosage amount—and is thus classified as a supervised learning method. We attempted to analyze non-experimental datasets from the field, which did not have supervised information. Using a combination of the two algorithms, we were able to overcome the limitations of each method and extract the characteristics, which are notable and/or interesting points of the dataset to study, and assign non-numerical classification information for downstream analysis. We propose that our advanced analysis scheme is applicable for data with unsupervised and/or multiple non-numerical variables to extract characteristic information for following the study and association rules between multiple variables, including non-numerical parameters.

2. Related Works

Several researchers have used machine learning to extract important characteristics and classification. In 2020, Li et al. used machine learning to diagnose a disorder [21]. Mesiar and Sheikhi applied machine learning to nonlinear entities [18]. Wei et al. and Tatsumi et al. used machine learning for data from the field, such as those from land and sea [22,23]. Wang et al. attempted to improve calculation performance by combining machine learning [24]. Shiokawa et al. applied market basket analysis for the visualization and selection association of information between human lifestyle factors and experimental measurements [17].

3. Proposed Framework

Our proposed framework consists of three well-known algorithms, K-means, random forest, and Apriori. The algorithms listed below provide a brief but detailed explanation of previous reports [17,25,26,27,28]. In this study, the algorithms were computed by R using the packages “Randomforest” and “arules”.

3.1. K-Means

n-dimensional points are clustered into K (1, … K), which is the arbitrary number of users. Set initial k centers $C=\left\{C_1,C_2,\dots ,C_k\right\}$. Cluster $C_i\left(i\in \left\{1,\dots ,k\right\}\right)$ is set to be a set of points in X that are closer to C_i than C_j in all j (j ≠ i), and we calculate the center of mass of all points in C_i: $C_i=\frac{1}{\left|C_i\right|}\sum _{x\in C_i}x$ Repeat until C is stable.

3.2. Random Forest

Random forest is a machine learning method based on the decision tree algorithm for classification and regression, developed by Breiman [27]. This algorithm displays the training method called bootstraps, and the decision tree uses Gini impurity to generate branches and extract variable importance. Gini impurity for a set of data with J classes ($i\in \left\{1,2\dots ,J\right\}$) labeled p is computed as follows.

$$\mathrm{Gini} \mathrm{impurity}\left(p\right)=1-\sum _{i=1}^J{p}_i^2$$

3.3. Association Rule Mining (Apriori)

We performed association rule mining as previously described [17]. Briefly, we used the Apriori algorithm [28] and applied parameters in R with the packages “arules” with “support”, “confidence”, and “lift”. The meanings are described below using formulas with probability of X and Y.

$$\mathrm{support}\left(X\Rightarrow Y\right)=P\left(X{{\displaystyle \cap }}^{ }Y\right)$$

$$\mathrm{confidence}\left(X\Rightarrow Y\right)=\frac{P\left(X{{\displaystyle \cap }}^{ }Y\right)}{P\left(X\right)}$$

$$\mathrm{lift}\left(X\Rightarrow Y\right)=\frac{P\left(X{{\displaystyle \cap }}^{ }Y\right)}{P\left(X\right)P\left(Y\right)}$$

In this study, the parameters were set as follows: “support” was 0.063, “confidence” was 0.25, and “lift > 1” with “maxlen = 2.” Our our data preparation for Apriori was converted using datasets to zero-one data by ranking and a quarter of the whole variables belonging to “high” or “low” were “one” and three quarters were “zero”. Therefore, when any association rules occurred randomly, support was approximately 0.063 and confidence was 0.25.

The analysis flow details are given in the Materials and Methods section below.

4. Materials and Methods

4.1. Overview

The analytical scheme of our study is presented in Figure 1. We used a two-stage analysis: first, the group was divided using K-means clustering of the applied dataset, and the validity of the result was evaluated using random forest based on the error rate. Important variables were extracted simultaneously with this step. This step was therefore termed “importance-based K-means”, referred to as “i-means” for short. Second, association rules were mined in each dataset and qualitative information was evaluated by Apriori [29].

NMR measurements were normalized by 2,2-dimethyl-2-silapentane-5-sulfonate (DSS). Each dataset was converted as per the ratio to the sum of each sample. K-means clustering was performed on each dataset and the validity checked using the random forest importance-based error rate. This part is referred to as “importance-based K-means”, or “i-means”. The resulting and original datasets were converted into zero-one data by data ranking or class information resulting from i-means. The zero-one data were analyzed by the Apriori algorithm. Finally, we selected meaningful association rules by the importance calculated by random forest from the extracted association rules by the Apriori algorithm.

4.2. Sample Preparation

We used natural fish samples from Japanese hydrosphere cultivation at 118 points in the river or sea. The samples (n = 315 individuals from 21 orders of taxonomy) were collected over the period of 2012–2016. Collected samples were stored at −30 °C or −80 °C until anatomic analysis was performed in the lab, to obtain intestinal content for measurements. Observers assessed and attached qualitative information, feeding behavior, ecology, habitat, body shape, tail type, season of collection, color, and scale type (

Table 1. Added qualitative information and sample numbers.

Qualitative Parameters	Observed Number
Bacterial Class 1 (from i-means)	26
Bacterial Class 2	93
Bacterial Class 3	90
Body1 (Spindle-shaped)	64
Body2 (Compressed)	144
Body3 (Cubic-shaped)	81
Body4 (Flat-shaped)	16
Body5 (Long, slender)	14
Color1 (Red)	85
Color2 (Blue)	16
Color3 (Yellow)	10
Color4 (Brown)	75
Color5 (Black)	131
Ecology1 (Migratory)	36
Ecology2 (Territorial-staying)	188
Ecology3 (Rockfish)	66
Ecology4 (Demersal fish)	26
Food1 (Plankton eater)	42
Food2 (Herbivore)	41
Food3 (Polychaeta eater)	150
Food4 (Crustacea eater)	259
Food5 (Mollusk eater)	82
Food6 (Fish eater)	162
Habitat1 (Freshwater)	24
Habitat2 (Brackish water)	22
Habitat3 (Coast)	136
Habitat4 (Offshore)	92
Habitat5 (Deep sea)	48
Length1 (50–150 mm)	94
Length2 (150–200 mm)	97
Length3 (200–300 mm)	88
Length4 (Over 300 mm)	58
NMR Class 1 (from i-means)	192
NMR Class 2	99
NMR Class 3	7
NMR Class 4	11
Scales1 (Matte)	52
Scales2 (Glossy)	133
Scales3 (Shiny)	131
Season1 (March–May)	90
Season2 (June–August)	68
Season3 (September–November)	115
Season4 (December–February)	46
Tail1 (Two distinct)	122
Tail2 (Others)	191

The sample had a single category; if there were overlaps and missing values, the number did not always reflect the sample number.

). The fish samples were measured by length and divided into four groups by length average and variance. The intestinal content of the fishes was collected in a sample tube, freeze-dried, powdered, and stored at −80 °C for further experiments.

4.3. Nuclear Magnetic Resonance

For nuclear magnetic resonance (NMR) observations, 18 mg of each powdered sample was extracted in 600 μL of KPi buffer containing 90% deuterium oxide and 1 mM sodium 2,2-dimethyl-2-silapentane-5-sulfonate (DSS) at 65 °C for 15 min, and then centrifuged at 17,800 G for 5 min. The entire supernatant was mixed and transferred to a 5-mm NMR tube. Two-dimensional J-resolved (2D J-RES) NMR spectra were acquired at 298 K using a Bruker AVANCE II 700 spectrometer equipped with an ¹H inverse triple-resonance cryogenically cooled probe with Z-axis gradients (Bruker BioSpin GmbH, Rheinstetten, Germany). In brief, 2D J-RES NMR spectra were acquired using the standard Bruker pulse program jresgpprqf, with 16 K (F2) and 16 (F1) points, and were then categorized into 32 transients and 16 dummy scans.

4.4. Data Processing

The 2D J-RES NMR spectra were processed using TopSpin software (Bruker Biospin: https://www.bruker.com/en/products-and-solutions/mr/nmr-software/topspin.html, accessed on 29 April 2022). Tilt correction and symmetrization were performed and projections of the 2D spectra were obtained. We defined 235 regions of interest manually using Revolution R Open software (https://cran.r-project.org/, accessed on 29 April 2022).

4.5. Data Preprocessing for Analysis and Annotation

The NMR measurement outcomes were normalized according to the DSS signal and peaks were annotated by comparison with premeasured standards from our database of metabolites acquired under the same conditions (Figure S1,

Table 2. NMR signal annotation of gut contents of fish extracts using KPi/D₂O.

Peak Number	Δ1H (ppm)	Δ13C (ppm)	Annotation
1	0.93	13.9	Ile
2	0.946	23.7	Leu
3	0.97	19.4	Val
4	1.001	20.6	Ile
5	1.034	20.6	Val
6	1.319	22.2	Lactate, Thr
7	1.471	18.9	Ala
8	1.701	26.7	Leu
9	1.72	29	Lys
10	1.909	25.9	Acetate
11	1.973	38.6	Ile
12	1.999	31.9	Pro
13	2.045	29.7	Glu
14	2.124	29	Gln
15	2.261	31.9	Val
16	2.339	36.2	Glu
17	2.394	36.9	Succinate
18	2.636	31.5	Met
19	2.677	39.4	Asp
20	2.711	37.2	DMA
21	2.807	39.4	Asp
22	2.876	47.4	TMA
23	3.025	39.7	Creatine
24	3.196	56.6	Choline
25	3.233	43.3	Arg
26	3.223	76.9	Glucose
27	3.258	62.3	TMAO, Taurine
28	3.414	38.1	Taurine
29	3.501	70	Choline
30	3.548	44.2	Gly
31	3.645	65.3	Glycerol
32	3.921	56.5	Creatine
33	3.955	63	Ser
34	4.103	71.2	Lactate
35	4.241	68.7	Thr
36	4.642	98.6	Glucose
37	5.224	94.8	Glucose
38	6.091	91.1	inosine
39	6.887	118.6	Tyr
40	7.108	119.8	His
41	7.181	133.4	Tyr
42	7.265	124.8	Trp
43	7.318	131.9	Phe
44	7.365	130.5	Phe
45	7.412	132	Phe
46	7.518	114.6	Trp
47	7.721	121.2	Trp
48	7.999	138.6	His
49	8.178	148.3	Inosine
50	8.223	148.9	Inosine
51	8.338	143	Inosine

). Datasets from NMR and gut microbe analysis (MiSeq) were processed to obtain the composition ratio to the sum of the measurement in fish intestines for each sample. Subjective qualitative information was assigned to samples by observation, as mentioned above (Table 1). These observations were made by multiple persons subjectively. The classification of a sample was difficult if qualitative information belonged to multiple categories.

4.6. Fish Gut Microbe Analysis by MiSeq

Fish intestinal microbial DNA was extracted as described in a previous report [30]. Further, DNA amplification was conducted according to the same report. Briefly, we used universal primers 954 f and 1369 r, targeted to the V6–8 regions of the 16S rRNA cording region.

4.7. Data Analysis

4.7.1. i-Means Analysis

Median points k of the NMR and MiSeq dataset were calculated using K-means clustering in R (https://cran.r-project.org/, accessed on 29 April 2022). R required the number of median points to perform K-means. Median point values were determined by reference to previous reports [22,31,32]. Briefly, the k of the bacterial data comprised three clusters, and the NMR data comprised four. Next, the validity of the calculated median values was assessed by evaluating the random forest error rate using the “Randomforest” package in R. In other words, K-means clustering transformed tentative class information into datasets. Then, the random forest algorithm confirmed that we could correctly classify using tentative class information as supervised data. At the same time, we computed the silhouette index of each K-means result for an internal validation index to evaluate the clustering quality. These steps were repeated at least 10,000 times to improve validity and to identify classification and characteristic variances, which were a result of plateaus on the error rate curve. The i-means analysis was performed at least 10,000 times and importance parameters that appeared at high frequency (over four times in the five times i-means was performed) were adopted as meaningful characteristic variances and used for subsequent analyses. In this step, we measured the elapsed times for reference with differential performance with two computers using MiSeq data ($53\times 209$ matrix, 10,000 repetitions). The elapsed time clearly depended on the data size, repetition time, and computer environment.

4.7.2. Association Rule Mining (Apriori)

Numerical data from NMR and MiSeq were converted to zero-one data and categorized as high, low, or null according to the interquartile range, as per the previous study [17]. “High” was defined as higher than the interquartile range; “low,” lower than the interquartile range; and other data points were classified as “null.” The data were then merged and, then, with the qualitative information shown in Table 1, were added into a matrix. The matrix was calculated by the Apriori algorithm (support = 0.063, confidence = 0.25, maxlen = 2) and association rules were extracted (lift > 1) using the R package “arules.” The calculation condition was decided according to previous reports [17]. The association rule network was depicted by Gephi, an open-source software (https://gephi.org/, accessed on 29 April 2022).

5. Results

5.1. Overview

The heatmap illustrating the correlation between bacterial variances and the NMR signal is shown in Figure S1. The advanced “i-means” method was able to process unsupervised data into the supervised random forest algorithm. Figure 2 and Figure S3 show the results of mining the characteristic variance of the groups that were determined using K-means clustering based on the importance calculated by random forest. The information that was classified using i-means analysis with additional qualitative parameters is presented in Table 1, and it showed better separation than category-based separation (Figure S4). Following zero-one data conversion, we extracted 22,462 association rules (Figure 3A). In this study, we verified the K-means validity by using the random forest error rate and not internal validation metrics, such as the silhouette index. The positive judgment rate converged to the plateau by repeating i-means (Figure S5). However, the silhouette index did not correlate the error rate with the random forest (Table S2).

The axis of ordinate indicates Gini impurity (importance). The graph shows the importance ranking of bacteria. Bacterial importance has three high importance factors.

5.2. Association Rules Focused on the Bacterial Data

The relationship between bacterial/NMR classification and qualitative parameters extracted from the association network is illustrated in Figure 3B,C. Moreover, representative association rules are summarized in

Table 3. Extracted association-rule-correlated bacterial class information from i-means result to other variances.

Source	Target	Support	Confidence	Lift
Bac_class1	High_Firmicutes	0.08	0.96	6.44
	Ecology2 (Territorial-staying)	0.06	0.77	1.29
	Food4 (Crustacea-eater)	0.07	0.88	1.08
Bac_class2	High_Actinobacteria	0.10	0.32	2.42
	Low_Bacteria.Other	0.13	0.43	1.83
	High_Proteobacteria	0.23	0.77	1.44
	Season3 (September–November)	0.13	0.45	1.24
	Color5 (Black)	0.15	0.51	1.22
	Color4 (Brown)	0.08	0.28	1.17
	Body3 (Cubicshape)	0.09	0.30	1.17
	Scales3 (Shiny)	0.14	0.47	1.14
	NMR_class1	0.20	0.68	1.11
	Habitat3 (Coast)	0.14	0.46	1.07
	High_Gln.Malic.acid	0.17	0.57	1.07
	Tail1 (two-distinct)	0.12	0.41	1.05
	Length1(50–150 mm)	0.09	0.31	1.04
	High_Thr	0.24	0.82	1.03
	High_Ile	0.29	0.99	1.01
	High_Creatine.G3P.GPC	0.29	0.99	1.01
	High_Lactate	0.29	0.99	1.01
	High_Val	0.29	0.99	1.01
	High_Ile	0.29	0.99	1.01
	High_Ala	0.29	0.99	1.01
	High_Gly	0.29	0.99	1.01
	High_TMAO	0.29	0.99	1.01
	High_Creatine	0.29	0.99	1.01
	High_Leu	0.29	0.98	1.00
Bac_class3	Low_Actinobacteria	0.07	0.26	2.30
	Low_Firmicutes	0.08	0.29	2.17
	Low_Bacteroidetes	0.10	0.34	1.78
	High_Proteobacteria	0.26	0.90	1.68
	Body1 (Spindle shaped)	0.08	0.28	1.37
	Season2 (June–August)	0.07	0.26	1.18
	Ecology2 (Territorial-staying)	0.20	0.69	1.15
	Tail1 (two-distinct)	0.13	0.44	1.15
	Low_Bacteria.Other	0.08	0.27	1.14
	Length1 (50–150 mm)	0.10	0.33	1.12
	Habitat4 (Offshore)	0.09	0.32	1.10
	Food6 (Fish-eater)	0.16	0.57	1.10
	High_Gln.Malic.acid	0.17	0.58	1.08
	Color5 (Black)	0.13	0.44	1.07
	Scales3 (Shiny)	0.13	0.44	1.07
	NMR_class2	0.10	0.33	1.06
	Body2 (Compressed)	0.14	0.48	1.05
	Season1 (March–May)	0.08	0.29	1.01
	High_Acetate	0.26	0.92	1.00

,

Table 4. Extracted association-rule-correlated NMR class information from i-means result to other variances.

Source	Target	Support	Confidence	Lift
NMR_class1	Body1 (Spindle-shaped)	0.17	0.27	1.33
	Scales3 (Shiny)	0.32	0.53	1.26
	Tail1 (Two distinct)	0.29	0.48	1.24
	Low_Bacteria.Other	0.17	0.28	1.18
	Color5 (Black)	0.29	0.47	1.13
	Ecology2 (Territorial-staying)	0.41	0.67	1.12
	High_Gln.Malic.acid	0.36	0.59	1.11
	Bac_class2	0.20	0.33	1.11
	High_Proteobacteria	0.36	0.58	1.09
	Food5 (Mollusca eater)	0.17	0.28	1.08
	Length3 (200–300 mm)	0.18	0.30	1.06
	Length2 (150–200 mm)	0.20	0.32	1.05
	Color1 (Red)	0.17	0.28	1.04
	Season3 (September–November)	0.23	0.38	1.03
	High_Thr	0.50	0.82	1.03
	High_Acetate	0.57	0.94	1.02
	High_Leu	0.61	1.00	1.02
	High_Ile	0.61	1.00	1.02
	High_Creatine.G3P.GPC	0.61	1.00	1.02
	High_Lactate	0.61	1.00	1.02
	High_Val	0.61	1.00	1.02
	High_Ile	0.61	1.00	1.02
	High_Ala	0.61	1.00	1.02
	High_Gly	0.61	1.00	1.02
	High_TMAO	0.61	1.00	1.02
	High_Creatine	0.61	1.00	1.02
	Food6 (Fish eater)	0.32	0.52	1.01
	Season1 (March–May)	0.17	0.29	1.00
	Body2 (Compressed)	0.28	0.46	1.00
	Food4 (Crustacea eater)	0.50	0.82	1.00
NMR_class2	Color4 (Brown)	0.12	0.39	1.65
	Ecology3 (Rockfish)	0.10	0.32	1.54
	Scales2 (Glossy)	0.18	0.58	1.36
	Food3 (Polychaeta eater)	0.20	0.63	1.32
	Season2 (June–August)	0.09	0.27	1.26
	Length1 (50–150 mm)	0.12	0.37	1.25
	Body3 (Cubic-shaped)	0.10	0.31	1.22
	Tail2 (Others)	0.23	0.74	1.22
	Habitat4 (Offshore)	0.10	0.32	1.11
	High_Lactate	0.14	0.44	1.09
	Habitat3 (Coast)	0.15	0.46	1.08
	High_Thr	0.27	0.85	1.06
	High_Acetate	0.31	0.98	1.06
	Bac_class3	0.10	0.30	1.06
	Food4 (Crustacea eater)	0.27	0.87	1.06
	Length2 (150–200 mm)	0.10	0.32	1.05
	High_Ile	0.31	1.00	1.02
	High_Creatine.G3P.GPC	0.31	1.00	1.02
	High_Lactate	0.31	1.00	1.02
	High_Val	0.31	1.00	1.02
	High_Ile	0.31	1.00	1.02
	High_Ala	0.31	1.00	1.02
	High_Gly	0.31	1.00	1.02
	High_TMAO	0.31	1.00	1.02
	High_Creatine	0.31	1.00	1.02
	Body2 (Compressed)	0.15	0.46	1.02
	High_Leu	0.31	0.99	1.01

and

Table 5. Extracted association rules to high TMAO.

Source	Target	Support	Confidence	Lift
Body1 (Spindle-shaped)	High_TMAO	0.20	1.00	1.02
Body3 (Cubic-shaped)		0.26	1.00	1.02
Color4 (Brown)		0.24	1.00	1.02
Ecology1 (Migratory)		0.11	1.00	1.02
Ecology4 (Demersal fish)		0.08	1.00	1.02
Food1 (Plankton eater)		0.13	1.00	1.02
Habitat2 (Brackish water)		0.07	1.00	1.02
Habitat5 (Deep sea)		0.15	1.00	1.02
High_Acetate		0.92	1.00	1.02
High_Ala		0.98	1.00	1.02
High_Creatine		0.98	1.00	1.02
High_Creatine.G3GPC		0.98	1.00	1.02
High_Gln.Malic.acid		0.53	1.00	1.02
High_Gly		0.98	1.00	1.02
High_Ile		0.98	1.00	1.02
High_Ile		0.98	1.00	1.02
High_Lactate		0.41	1.00	1.02
High_Lactate		0.98	1.00	1.02
High_Leu		0.98	1.00	1.02
High_malic.acid		0.16	1.00	1.02
High_Thr		0.80	1.00	1.02
High_Val		0.98	1.00	1.02
Length2 (150–200 mm)		0.31	1.00	1.02
Low_Actinobacteria		0.11	1.00	1.02
Low_Cyanobacteria		0.09	1.00	1.02
Low_Fusobacteria		0.08	1.00	1.02
Low_Ile		0.09	1.00	1.02
Low_Planctomycetes		0.11	1.00	1.02
Low_Tenericutes		0.07	1.00	1.02
Low_Tyr		0.11	1.00	1.02
NMRclass1		0.61	1.00	1.02
NMRclass2		0.31	1.00	1.02
Scales1 (Matte)		0.17	1.00	1.02
Season2 (June–August)		0.22	1.00	1.02
Season4 (December–February)		0.15	1.00	1.02
Food6 (Fish eater)		0.51	0.99	1.01
Tail1 (Two distinct)		0.38	0.99	1.01
Season3 (September–November)		0.36	0.99	1.01
Length1 (50–150 mm)		0.30	0.99	1.01
Bac_class2		0.29	0.99	1.01
Habitat4 (Offshore)		0.29	0.99	1.01
Color1 (Red)		0.27	0.99	1.01
Food5 (Mollusca eater)		0.26	0.99	1.01
Low_Bacteria.Other		0.23	0.99	1.01
Habitat3 (Coast)		0.43	0.99	1.00
Ecology3 (Rockfish)		0.21	0.98	1.00
Scales3 (Shiny)		0.41	0.98	1.00
Length4 (Over 300 mm)		0.18	0.98	1.00

. Interestingly, all bacterial classes exhibited associations with other classes, but two NMR classes showed no association with any variables. Specifically, Bacterial Class 1 was linked with a higher ratio of Firmicutes and feeding behavior category 4. Bacterial Classes 1 and 3 both showed associations with Ecology Category 2. Bacterial Classes 2 and 3 were both found to have high proportions of Proteobacteria compared with other bacterial groups, and showed associations with four qualitative parameters (color category 5, length category 1, scale category 3 and tail category 1) and one NMR-annotated signal (Table 3). Bacterial Class 2 was associated with high levels of Actinobacteria and four qualitative parameters (body category 3, color category 4, habitat category 3 and season category 3) and many amino acid signals in NMR Class 1. Bacterial Class 3 did not exhibit a high proportion of any one bacterial group, but had a low ratio of three common bacteria (Actinobacteria, Bacteroidetes and Firmicutes) and was found to be associated with six qualitative parameters (body category 1, 2; food category 6; habitat category 4; and season 1, 2) and a high acetate NMR signal (Table 3). In contrast, NMR Classes 1 and 2 had many common NMR signals and were associated with three qualitative parameters but with no bacterial variables (Table 4).

5.3. Association Rules Focused on NMR Signals

We identified two associations for NMR Class 1, one with a high proportion of Proteobacteria and one with a lower proportion of other bacterial groups. Finally, the NMR signal network indicated that some source factors might affect the TMAO ratio in the NMR signals of intestinal content samples. Overall, the network was found to have 48 association rules consisting of seven bacterial factors including one class factor, and 18 NMR factors including two class and 23 qualitative factors (Figure 3D, Table 5). These qualitative factors had four habitats, in Categories 2–5.

6. Discussion

We have developed a two-part advanced analysis scheme, involving “i-means” and market basket analyses (Figure 1). The i-means step includes two machine learning algorithms, “K-means clustering” and “random forest.” While K-means is a well-known clustering algorithm, the appropriate number of categories is usually not known and the result depends on the initial point of the centroids, which is known as the initialization trap [33]. To overcome the problem of selecting an appropriate number of groups, we referred to previous reports [22,31,32]. To address the issue of the centroids, we checked the validity of initial centroid coordinates by comparing the random forest error rate to the K-means result. Using i-means, we were able to successfully divide a complicated dataset because the error rate from random forest improved with repeated i-means analyses (Figure S4). The plateauing of the error rate indicates that repeated i-means analysis might converge the best clusters of the dataset (Figure S5). The divided classes of NMR and Miseq datasets were able to be mined for unique variables and unique groups were identified (Figures S3 and S4). However, clustering the internal validation index, which is a well-known silhouette coefficient, did not correlate with the random forest error rate. This could mean that better clustering does not always provide meaningful separation to focus on unique characteristics of datasets (Table S2). The result was thought to be due to the data being from the field, not experimental, and they included several parameters and relevant variance but not ideal data distribution. This likely indicates that it is difficult to judge complicated data only by the clustering of the internal validation index, such as the silhouette coefficient. Additionally, qualitative data were able to be attached to classify parameters from i-means analysis and market basket analysis enabled the extraction of some association rules to the additional parameters from i-means classification (Figure 3B,C). These results demonstrate that our proposed scheme is suitable for identifying notable variables and proposing related parameters in challenging datasets.

We performed i-means analysis on the MiSeq dataset in order to divide the data into three groups because it has been previously reported that there are three human enterotypes based on gut microbiota at the phylum level, predominantly Firmicutes, Bacteroidetes, or Actinobacteria [31]. Previous studies have proposed that Proteobacteria represent the prominent microbial phyla in the fish gut [34]. The present study revealed the fish gut microbiota to be predominantly composed of Firmicutes, Actinobacteria, or Proteobacteria, rather than Bacteroidetes as in human gut enterotypes (Figures S3 and S4). Bacteroidetes were in low-to-high abundance in individual fish gut microbiomes (Figure S3A), but the composition ratio had an impact on group separation because the importance value calculated by i-means was lower than Proteobacteria, Firmicutes, and Actinobacteria (Figure S3), which may indicate that Proteobacteria and Bacteroidetes play similar roles in the gut microbiota of fish and humans, respectively. In marine biofilms, a considerable amount of carbohydrate metabolism is carried out by Bacteroidetes and Gammaproteobacteria, including Proteobacteria; thus, Proteobacteria in the fish gut might decompose consumed carbohydrates [35]. Our analysis revealed that Bacterial Class 3 was associated with a high abundance of Proteobacteria and a low abundance of Bacteroidetes, and the high-intensity NMR signal that was annotated as acetate might be derived from carbohydrate degradation (Figure 3B).

Bacterial Class 1 was found to have fewer association rules by our proposed scheme than Bacterial Classes 2 and 3. The associations were a higher proportion of Firmicutes, while food category 4 included crustacea eaters and was found to share ecology category 2 with Bacterial Class 3. Crustacea use chitin and/or chitosan to form their exoskeletons [36], and it has been reported that many species of bacteria, including bacillus which are Firmicutes, express enzymes with chitinase activity [37]. Fish with a Bacterial Class 1 gut microbiome may, therefore, digest chitin in their feed via enzymes originating from their intestinal Firmicutes. This suggestion indicates that food influences the intestinal microbiome [38,39].

We performed i-means analysis on NMR datasets in order to divide the data into four groups. This revealed that the intestinal contents of fish have many variables, but the high importance value of dividing the NMR signal clusters calculated by i-means was only TMAO, which suggested that the number of high impact factors (Gini impurity) was less than bacterial factors (Figure S3). The most important variable for dividing the NMR dataset was annotation as a TMAO signal. The group with the highest TMAO proportion was NMR Class 4 (Figure S4B). However, this group did not exhibit associations with other NMR classes, bacterial classes, or qualitative information. It is known that TMAO is involved in osmoregulation in the muscle of marine creatures [40,41]. TMAO indicated association rules to four habitat factors: brackish water, coast, offshore, and deep sea (not fresh) water (Figure 3D).

Freshwater fish have been reported to have lower TMAO concentrations in their muscle tissue than seawater fish [42,43,44]. Our scheme could be used to identify the characteristics of organisms using qualitative and non-numerical information.

7. Conclusions and Future Work

Our scheme overcomes some of the limitations of K-means clustering and random forest by combining both methods in an iterative process. Moreover, machine learning clustering methods provide a useful approach to extract association rules when combined with market basket analysis. This method may enable researchers to cluster and characterize challenging unsupervised datasets which have few correlations, making it possible to determine notable points and identify important characteristics. We could not discuss all association rules; some rules may warrant further research.

In the future, we hope to develop other combination schemes; for example, our method “i-means” would replace random forest to support vector machine or, alternatively, K-means to K-means++. Support vector machine might work as well as i-means, but K-means++ may not. The method improves the elapse time of K-means results based on the silhouette index. In this study, the silhouette index was not necessarily appropriate and it is necessary to suggest a better index to extract meaningful characteristics. We must select a better analysis scheme depending on the situation and dataset obtained. Our scheme may enable the extraction of characteristics and association rules from unsupervised or multiple non-numerical variable datasets obtained from uncontrolled and nonexperimental situations, such as environmental field data. However, it is necessary to experimentally assess the extracted rules by our scheme to clarify the meaning.