Classifying Oryza sativa accessions into Indica and Japonica using logistic regression model with phenotypic data

Phenotypic data

A set of phenotypic data and subpopulation information was used, which was originally generated and analyzed by Zhao et al. (2011). The data set consisted of 413 accessions originating from 82 countries, obtained at http://ricediversity.org/data/sets/44kgwas. Out of 32 phenotypic variables, 24 variables were selected because they were quantitative and not confined to a certain geography. The selected variables were divided into morphology (culm habit, flag leaf length, flag leaf width), yield components (panicle number per plant, plant height, panicle length, primary panicle branch number, seed number per panicle, florets per panicle, panicle fertility), seed morphology (seed length, seed width, seed volume, seed surface area, brown rice seed length, brown rice seed width, brown rice volume, seed length/width ratio, brown rice length/width ratio), stress tolerance (straighthead susceptibility, blast resistance) and quality (amylose content, alkali spreading value, protein content). Every accession in the data set belonged to one of the following subpopulation groups: admixed (62), aromatic (14), aus (57), indica (87), temperate japonica (96) and tropical japonica (97). In this study, temperate japonica and tropical japonica were combined into japonica.

Stepwise variable selection and parameter estimation

The LRM formula used in this study can be denoted as follows: (1)${\displaystyle P\left(japonica|x_1,\dots ,x_n\right)=\frac{1}{1+e^{-\left({\beta }_0+{\beta }_1{x}_1{\cdots +\beta }_n{x}_n\right)}}}$

where P(japonica|x₁, …, x_n) is the probability of an accession being japonica (>0.5) or indica (<0.5) given the predictor variables, x₁,…, x_n; β₀ is the constant term; β₁,…, β_n are the parameters for the predictor variables, x₁, …, x_n, respectively.

In Eq. (1), the response variable is binary between indica and japonica, and the predictor variables (phenotypic variables) are quantitative. To identify a set of predictor variables that can maximize the indica/japonica separation power, n sets were prepared from 1D to nD, in which the nD is a set containing all possible combinations of n different phenotypic variables (e.g., 1D contains every single phenotypic variable, 2D contains all possible pairs of phenotypic variables, and so forth). The n was increased by one until the maximum P(japonica|x₁, …, x_n) was reached, in which every single selection from the nD set was subject to the following steps:

1. Calculate a set of parameters by fitting the LRM (Eq. (1)) with the phenotypic variables in a selection.

2. Applying the phenotypic variables in the selection used in Step 1 to Eq. (1) with the parameters estimated in Step 1. The resulting value must be between 0 and 1, from which the indica or japonica can be determined at 0.5.

In order to estimate the indica/japonica separation power, a receiver operating characteristic (ROC) curve was implemented using an R package called pROC (Robin et al., 2011). In an ROC space, the x- and y- axes ranging between 0 and 1 represent the false positive rate (FPR or specificity) and true positive rate (TPR or sensitivity), respectively. Thus, an area under an ROC curve (AUC) can range between 0 and 1. The closer the AUC is to 1, the higher the indica/japonica separation power. By referring to the resulting AUCs, a set of phenotypic variables that maximized the indica/japonica separation power was identified, and the resulting parameters were used to define the customized LRM.

Applications of the customized LRM

1. Estimating the indica/japonica classification accuracy: given the 280 indica/japonica accessions, the customized LRM along with the related phenotypic variables were taken into 100 iterations of ten-fold cross-validations, from which the 100 values were obtained. The resulting values were averaged into the indica/japonica classification accuracy.

2. Estimating the variable interaction: as the LRM uses multiple phenotypic variables, some portion of indica/japonica classification power might be derived from interactions between variables. The interaction magnitude for every variable with all other variables was calculated using the H-statistic (Frieman & Popescu, 2008). The resulting H-statistic can range between 0 and 1 with no interaction resulting in 0 and full interaction resulting in 1. The H-statistic was calculated using an R package, iml (Molnar, Casalicchio & Bischl, 2018).

3. Classifying accessions in each minor subpopulation group into indica and japonica: the customized LRM was applied to each minor subpopulation group (aromatic, aus, admixed) to divide accessions into indica and japonica. This examination aimed to observe how accessions in each minor subpopulation group are phenotypically divided from the indica/japonica perspective.

Dendrogram-based indica/japonica classification

To draw dendrograms, the genomic data set was obtained at http://ricediversity.org/data/sets/44kgwas (Zhao et al., 2011). The genomic data set consisted of the accessions in the phenotypic data, genotyped with 36,901 SNPs. Two dendrograms were drawn; one included the 280 indica/japonica accessions, and the other included all 413 accessions. Then, the dendrogram-based indica/japonica classifications were compared with the LRM-based indica/japonica classifications. Each dendrogam was graphed based on a genetic distance matrix in which the genetic distances between two different accessions were calculated using the following equation: (2)${\displaystyle GeneticdistancebetweenAandB=2-{IBS}_{A,B}}$

where the IBS_A,B is the IBS (identical by state) coefficient between A and B, which can range between 0 and 2.

The IBS matrix was computed using Numericware i (Kim & Beavis, 2017) which is freely available at https://figshare.com/articles/Numericware_i/3496787.

Data availability

With the exception of calculating the IBS matrix, all other computations were conducted using R (R Core Team, 2016). The data set and R scripts used in this study are freely available at https://github.com/bongsongkim/logit.regression.rice.

Stepwise variable selection

Figure 1 and Table 1 suggest that more phenotypic variables led to stronger indica/japonica separation power of the LRM, and the fully accurate indica/japonica separation power (AUC = 1) was achieved with a 7D selection. Table 2 summarizes the phenotypic variables yielding the maximum indica/japonica separation power in each set (1D to 7D). The set of phenotypic variables yielding AUC = 1 comprises panicle number per plant, seed number per panicle, florets per panicle, panicle fertility, straighthead susceptibility, blast resistance and protein content.

Estimating the indica/japonica classification accuracy

Given the seven phenotypic variables yielding AUC = 1, the indica/japonica classification accuracy was estimated using the 280 indica/japonica accessions, for which 100 iterations of ten-fold cross-validations were implemented. As a result, the indica/japonica classification accuracy of 100% was obtained. This indicates that the seven phenotypes were certainly impactful for the fully accurate indica/japonica classification. Assuming that some portion of classification power might be derived from interactions between variables, the H-statistic was calculated for the purpose of estimating how much each variable generates the classification power in collaboration with other variables. Table 3 summarizes the resulting H-statistic values, indicating nearly no interactions between variables.

Dendrogram-based indica/japonica classification

Figure 2 shows the dendrogram-based indica/japonica classification with the 280 indica/japonica accessions, in which the indica- and japonica-varietal clades were accurately divided, and the japonica accessions were further accurately divided into temperate japonica and tropical japonica (Fig. S1). This result shows that the dendrogram-based indica/japonica classification accuracy was 100%.

Applying the LRM to each minor subpopulation group

The LRM customized with parameters yielding AUC = 1 was used to investigate how it classifies the accessions in each minor subpopulation group (admixed, aromatic, aus) into indica and japonica. The customized LRM was as follows: (3)${\displaystyle P\left(japonica|x_1,\dots ,x_7\right)=\frac{1}{1+e^{-\left(84446.0-5832.5x_1+35800.4x_2-38008.7x_3-50943.2x_4-491.6x_5+353.7x_6-214.6x_7\right)}}}$

where P(japonica|x₁, …, x₇) is the probability of an accession being japonica (>0.5) or indica (<0.5) given the predictor variables, x₁,…, x₇; x₁ = the predictor variable for panicle number per plant; x₂ = the predictor variable for seed number per panicle; x₃ = the predictor variable for florets per panicle; x₄ = the predictor variable for panicle fertility; x₅ = the predictor variable for straighthead susceptibility; x₆ = the predictor variable for blast resistance; x₇ = the predictor variable for protein content.

The absolute values for parameters in Eq. (3) in descending order are 50943.2 (-) for panicle fertility, 38008.7 (-) for florets per panicle, 35800.4 (+) for seed number per panicle, 5832.5 (-) for panicle number per plant, 491.6 (-) for straighthead susceptibility, 353.7 (+) for blast resistance and 214.6 (-) for protein content. These suggest that, when it comes to the indica/japonica separation power, the panicle fertility is most impactful, followed by florets per panicle, seed number per panicle, panicle number per plant, straighthead susceptibility, blast resistance and protein content.

Because of missing phenotypic records, the size of each minor subpopulation group was reduced from 62 to 52 for the admixed group, 14 to 12 for the aromatic group and 57 to 26 for the aus group. Figure 3 shows that Eq. (3) split each subpopulation group into indica and japonica in ratios of 5:47 for the admixed group, 5:7 for the aromatic group and 18:8 for the aus group, respectively. Meanwhile, the dendrogram drawn with the whole accessions (413) shows two major clades (upper and lower), in which temperate japonica, tropical japonica and aromatic formed accurately separate groups within the upper clade (hereafter called japonica-varietal clade), while indica and aus formed accurately separate groups within the lower clade (hereafter called indica-varietal clade). The admixed accessions were spread across all subpopulation groups (Fig. S2). Figure 4 shows three Venn diagrams, each of which represents comparison between the LRM-based and dendrogram-based classifications for each minor subpopulation group; the agreements were 92.3% (48/52) in the admixed group, 58.3% (7/12) in the aromatic group and 69.2% (18/26) in the aus group.