We downloaded the RNA sequencing data of 1,600 human pancreatic islet cells from the GEO (Transcript Expression Omnibus) database under the accession number of GSE81608 at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE81608 (Xin et al., 2016). There were 949 pancreatic islet cells from six T2D patients and 651 pancreatic islet cells from 12 non-diabetic donors. Within the 949 pancreatic islet cells from T2D patients, there were 569 alpha, 296 beta, 30 delta, and 54 PP cells. Within in the 651 pancreatic islet cells from non-diabetic donors, there were 377 alpha, 207 beta, 28 delta, and 39 PP cells. The expression levels of 39,851 genes were quantified as RPKM (Reads Per Kilo bases per Million reads). The processed gene expression profiles of these cells downloaded from https://ftp.ncbi.nlm.nih.gov/geo/series/GSE81nnn/GSE81608/suppl/GSE81608_human_islets_rpkm.txt.gz were used. Despite islet cells containing different cells, this work expects to identify the common gene signatures for T2D across multiple cell types.

In this study, we first used the Monte Carlo feature selection (MCFS) (Draminski et al., 2008) to evaluate the importance of all genes, obtaining a feature list and some informative genes expressed in diabetes. For the feature list, it was fed into the IFS (Liu and Setiono, 1998) with one classification algorithm to extract optimal genes that had a strong discriminate ability between diabetes and non-diabetes samples and construct an efficient classifier. On the other hand, repeated incremental pruning to produce error reduction (RIPPER) was employed to determine interpretable rules on gene expression patterns with informative features.

For the IFS method, one classification algorithm was necessary. In this study, we tried three classic classification algorithms: 1) support vector machine (SVM) (Cortes and Vapnik, 1995), 2) K-nearest neighbor (KNN) (Cover and Hart, 1967), and 3) random forest (RF) (Breiman, 2001). Their brief descriptions were as follows.

Classification algorithms mentioned in Section 2.3 are powerful to construct efficient classifiers. However, we cannot understand their principles because they are black-box algorithms. In this case, few clues for uncovering essential differences between T2D patients and non-diabetic donors can be obtained. In view of this, we further adopted rule learning algorithms to investigate the RNA sequencing data. Although it is generally weaker than the aforementioned algorithms, it can provide rules that clearly indicate special expression patterns on T2D patients, thereby improving our understanding on T2D. The procedures were described in the following sections.

As mentioned in Section 3.2.1, the MCFS method can select some informative features. These features are quite essential to describe the characteristics of the dataset. Here, we used these features to construct classification rules via RIPPER algorithm (Cohen, 1995). Before that, the Johnson reducer algorithm (Johnson, 1974) was applied on the informative features to select the most important features, which had the similar classification ability compared to the original informative features. The selected features were fed into the RIPPER algorithm. RIPPER, proposed by Cohen (1995), is a rule learning algorithm which is capable of handling large noisy datasets effectively. RIPPER is the improved version of IREP (Johannes and Widmer, 1994) which combines both the separate-and-conquer technique used first in the relational learner FOIL (Quinlan, 1990) and the reduced error pruning strategy proposed by Brunk and Pazzani (1991). In RIPPER, the training set is first split into growing and pruning sets. Then, repeat the rule grow phase and rule prune phase until no positive samples are left in the growing set, or the description length (DL) is 64 bits greater than the smallest DL found so far, or the error rate is greater than 50%. In the rule grow phase, one rule is generated by greedily adding conditions to the rule that achieves the highest FOIL’s information gain. In the rule prune phase, the rule is pruned using reduced error pruning. Finally, global optimization strategy is applied to further prune the rule set. The aforementioned procedures for constructing rules are also implemented in the MCFS program, that is, the set of rules is one output of the MCFS program.

In this study, we used six measurements to evaluate the performance of all classifiers under 10-fold cross-validation (Kohavi, 1995; Li Z. et al., 2022; Ding et al., 2022; Tang and Chen, 2022), including sensitivity (SN) (same as recall), specificity (SP), accuracy (ACC), Matthew correlation coefficient (MCC), precision, and F1-measure (Matthews, 1975; Zhao et al., 2018; Zhao et al., 2019; Jia et al., 2020; Liang et al., 2020; Zhang et al., 2021a; Zhang et al., 2021c; Pan et al., 2021). Their formulations are written as follows: S N = R e c a l l = T P T P + F N , (1) S P = T N T N + F P , (2) A C C = T P + T N T P + T N + F P + F N , (3) M C C = T P × T N − F P × F N ( T P + F P ) ( T P + F N ) ( T N + F P ) ( T N + F N ) , (4) P r e c i s i o n = T P T P + F P , (5) F 1 − m e a s u r e = 2 × R e c a l l × P r e c i s i o n R e c a l l + P r e c i s i o n , (6) where TP represents the number of truly positive samples, FP represents the number of false-positive samples, TN represents the number of truly negative samples, and FN represents the number of false-negative samples. Among these six measurements, we selected F1-measure as the key one, whereas others were provided for reference.

Some rules can be extracted via the Johnson reducer and RIPPER algorithms, which involved several features (genes), called rule genes, in the following text. We performed Gene Ontology (GO) enrichment analysis using R package topGO (http://bioconductor.org/packages/release/bioc/html/topGO.html, v.2.24.0) on these rule genes. The genes of interest were set as rule genes, and the gene background was set as all the available genes. The p-value threshold was set at 0.001.

Although we can construct efficient classifiers to identify cells of T2D patients through three classification algorithms, these classifiers were absolute black-box algorithms, which prevented us from uncovering the essential differences between cells of T2D patients and non-diabetic donors. As mentioned in Section 2.4, rule learning algorithms were employed.

According to the output of the MCFS program, 235 features were selected as informative features. To test the utility of the classification rules yielded by Johnson reducer and RIPPER algorithms, we performed the 10-fold cross-validations three times, obtaining the F1-measure of 0.910, which was lower than that of the optimal SVM classifier but higher than that of the optimal KNN and RF classifiers. The SN was 0.898, SP was 0.891, ACC was 0.895, MCC was 0.784, and precision was 0.923. Although such performance was lower than that of the optimal SVM classifier, the RIPPER algorithm can construct a group of rules, which made the classification procedure completely open and provided more insights. Thus, the Johnson reducer and RIPPER algorithms were applied to all samples, producing nine different classification rules, as listed in Table 2. These rules are able to accurately screen patients with T2D from non-diabetic population. Although these rules were mainly for non-diabetes, based on the aforementioned evaluation results (SP = 0.891), it was believed that these rules were statistically shown to cover almost all possible non-diabetes samples. Thus, investigation on these rules can also figure out the characteristics of T2D patients in an opposite aspect.

The MCFS method can directly output some informative features. These features can capture essential information of the dataset. Here, as mentioned in Section 3.3, 235 features were selected as informative features. We can directly use them to construct classifiers with different classification algorithms. These classifiers were also evaluated by 10-fold cross-validation. The main measurement, F1-measure, of these classifiers is listed in Table 3. For KNN, F1-measure varied between 0.839 and 0.849. The F1-measure of RF changed between 0.889 and 0.897. Also, SVM provided the F1-measure varying between 0.631 and 0.886. Compared with the F1-measure yielded by the optimal classifier based on the corresponding classification algorithm, the classifier using informative features generated a lower F1-measure, suggesting that such a classifier was inferior to the optimal classifier. The employment of the IFS method can help us construct more efficient classifiers.

Because the optimal SVM classifier provided the best performance, which used top 745 features (genes), we focused on these 745 genes. However, it is impossible to analyze them one by one. Here, only top seven genes were analyzed, which are listed in Table 1.

The first predicted gene, WDR45-like pseudogene (100128906), is the pseudogene of gene WDR45. According to recent publications, it encodes a functional lncRNA associated with the regulation of WDR45 (Tsuyuki et al., 2014; Lebovitz et al., 2015). WDR45 has been functionally related to autophagy (Lebovitz et al., 2015). Considering that abnormal autophagy has been well known to contribute to the pathogenesis of T2D (Lee, 2014), it is reasonable to speculate that the expression level of WDR45 and its upstream regulator (i.e., our predicted gene LOC100128906) may have quite different expressions in diabetic pancreatic islets cells compared to normal cells.

The next identified gene is MTND4P24 (100873254), which is shown to have quite different expression levels in diabetic and normal tissues containing multiple cell subtypes. As an lncRNA-encoding pseudogene, the expression level of such a gene is able to reflect the regulatory ability of lncRNAs on its target gene, MT-ND4 (Torrell et al., 2013; Mella et al., 2016). Recent publications also confirmed that the expression level of the target gene MT-ND4 is functionally related to cellular insulin sensitivity in rat models (Houstek et al., 2012). Therefore, as one regulator of MT-ND4’s expression, the expression pattern of MTND4P24 may involve in the pathogenic insulin sensitivity decreasing in type 2 diabetic cells. Similarly, a homolog of MTND4P24 and MTND2P28 (100652939) has also been predicted to have different expression levels in multiple cell subtypes from pathogenic or normal pancreatic islets. Considering its similar regulatory mechanisms and the biological function of MTND2, it is also quite convincing to regard such a gene as a potential distinctive standard for diabetic and non-diabetic cells (Mathews et al., 2005).

The predicted gene, RPS14P1 (100271063), is also a pseudogene, contributing to the regulation of ribosomal protein S14’s expression (Aubert et al., 1992). Meanwhile, the function of ribosomal protein S14 is widely reported to participate in p53-dependent cell-cycle arrest by interacting with MDM2 (Zhou et al., 2013), which is abnormally activated during the pathogenesis of diabetes (Golubnitschaja et al., 2006; Garufi et al., 2017). Thus, it is a reasonable assumption that ribosomal protein S14 together with RPS14P1 has different expression levels in normal and diabetic cells.

Apart from such predicted pseudogenes, we also identified some functional lncRNAs that may have different expression patterns in normal and diabetic cells. LINC00486 (285045) is a predicted lncRNA that contributes to the distinction of normal and diabetic cells. According to recent publications, various functional lncRNAs (Liu et al., 2014; Pullen and Rutter, 2014), including LINC00486, have been confirmed to contribute to the initiation and progression of T2D (Pullen and Rutter, 2014).

The following predicted gene, named ZBTB8OSP2 (729898), is a pseudogene and has been reported to contribute to anti-saccade response and eating disorders (Cornelis et al., 2014; Broer and van Duijn, 2015). As a transcriptional regulator for ZBTB8, such genes may indirectly contribute to a specific complication of T2D, the refractory diabetes insipidus, especially in adolescent male patients (Soto et al., 2014). Therefore, we can infer that such genes together with their downstream binding targets may have respective specific expression patterns in normal and diabetic cells.

The next predicted gene is THRAP3P1 (391524), the pseudogene of THRAP3. The post-transcriptional regulatory target of THRAP3 has been confirmed to dock on phosphoserine 273 of PPAR-gamma and further contribute to the pathogenic programming of diabetic genes, inducing insulin resistance (Choi et al., 2014). Therefore, to accomplish the regulatory role, such gene has a high expression level in normal cells compared to diabetic cells.

As we have discussed earlier, we identified multiple pseudogenes associated with T2D. Pseudogenes are nonfunctional segments with similar or reverse sequences of actual coding genes. The biological functions of pseudogenes are still unclear. It has only been speculated that pseudogenes participate in the post-transcriptional regulation via generating siRNAs, piRNAs, microRNAs, or other small RNAs (Guo et al., 2009). Although pseudogenes cannot generate protein products, the regulatory effects of such group of genes may still be significant under physical and pathological conditions (Tay et al., 2014). For transcriptomics analyses, especially for single-cell transcriptomics analyses, multiple pseudogenes have been identified as candidate biomarkers for different systematic diseases in studies specifically focusing on pseudogenes’ effects (Kalyana-Sundaram et al., 2012; Poliseno et al., 2015). For most previous studies, the pseudogenes were removed in the data preprocessing. Therefore, most previous studies have not identified a lot of pseudogenes as potential candidate biomarkers for diabetes. In our study, we did not filter out the pseudogenes and for the first time confirmed that pseudogenes with potential transcriptomic regulatory effects may further contribute to the regulation of specific diseases via regulating the biological functions of their respective recognized protein-coding genes.

Here, in this study from other perspective of view, we applied several machine learning algorithms to identify new potential biomarkers for T2D patients. Multiple previous publications have already identified a group of T2D biomarkers such as HbA1c, advanced glycation end-products (AGEs), and pigment epithelial-derived factor (PEDF) (Lyons and Basu, 2012). Also, for the publication from which we retrieved the single-cell sequencing data, unique biomarkers like LINC00486, ZNF445, and SYBU have also been identified for T2D (Xin et al., 2016). Compared with these prediction results, first, we identified a group of confirmed biomarkers like LINC00486, validating the efficacy and accuracy of our results. Second, we identified a group of new biomarkers like MTND4P24 and THRAP3P1. Although such genes have been shown to be functionally correlated with T2D, previous studies have not identified such genes as potential biomarkers of T2D. There are two major advantages of our studies compared to previous studies, which may lead us to find novel biomarkers:

1) First, compared with previous studies, we used the single-cell level data with the gene expression profiling of different cells and not just an averaged comprehensive value for each patient. Therefore, we can identify potential biomarkers that are missing due to the averaging procedures.

Such advantages explained why we identified novel protein biomarkers to distinguish T2D patients from normal controls. As we have discussed earlier, some identified biomarkers have been functionally correlated with T2D, implying that it is reasonable to regard such genes/transcripts as potential biomarkers for T2D.

We also screened out a group of functional quantitative rules of the gene expression pattern to distinguish non-diabetic cells from diabetic ones with more interpretability, which are listed in Table 2. Many qualitative rules can be validated according to the gene expression level in existing databases and recent reports on gene expression trends, which support the efficacy and accuracy of the rules. The detailed analysis of each expression rule is widely discussed as follows:

The first rule (rule1) involved four genes including LOC100128906 [(100128906), RPL3P4 (326307), PSPHP1 8781], and PTCHD1 (100873065). As mentioned earlier, gene LOC100128906 has been reported to have quite different transcriptomics patterns between normal and diabetic cells, inhibiting autophagy (Lebovitz et al., 2015). As the antagonistic gene of diabetes-associated autophagy, such genes are reasonable to have high expression in normal cells compared to diabetic cells. As for gene RPL3P4, the regulatory target of such pseudogene, RPL3 has been reported to have a quite low expression level in diabetic cells compared to normal cells (Tsai et al., 1994), corresponding with this rule. As for PSPHP1 (8781), it has been shown to be associated with the macrophage-related inflammation processes (Walker et al., 2015). Considering that during the initiation and progression of diabetes, regional and systematic inflammation have been widely observed (Donath et al., 2003; Lontchi-Yimagou et al., 2013), it is reasonable to predict such genes as quantitative parameters for the distinction of non-diabetes and diabetes. As for PTCHD1, although no direct evidence confirms its contribution on diabetes, it has been confirmed that such a gene is associated with the eye and ear complications of diabetes (Gambin et al., 2017), consistent with this rule.

As for the second rule (rule2), four genes were involved including MICOS10P3 (100462954), CTBP1 (1487), RPL3P4 (326307), and MTND4P24 (100873254). Few publications have reported the biological contribution of MICOS10P3; therefore, it is hard to interpret such gene’s contribution on T2D. As for gene CTBP1, it has been reported to participate in the abnormal phosphorylation processes (Kim et al., 2013) and shows a quite high expression level in diabetic cells compared to normal controls. As for gene RPL3P4, the regulatory target of such a pseudogene, RPL3 has been reported to have a quite low expression level in diabetic cells compared to normal cells (Tsai et al., 1994), corresponding with such a rule. MTND4P24 and its homolog, MTND5P11, have been confirmed to regulate a group of functional mitochondrial-encoded NADH ubiquinone oxidoreductase. According to recent publications, during the pathogenesis of diabetes, MT-ND4 has a quite low-expression pattern and on the contrary, MT-ND5 has a relevantly higher expression level, corresponding with the prediction expression level of their agonists individually (Elango et al., 2014; Urbanova et al., 2017).

In the third rule (rule3), apart from genes we have discussed earlier, the gene EIF5AL1 (143244) has also been predicted to have a higher expression pattern in normal cells but not in diabetic cells. Considering the abnormal endocrine stress responses of diabetic cells (Siddiqui et al., 2015), the lower expression level of EIF5AL1 may also contribute to the identification of diabetic cells. FXYD2 (486) has been shown to contribute to the pathogenesis of diabetes (Ding et al., 2019). Another specific gene in rule3 is the homolog of RPL3P4, RPL9P7, which may also participate in the regulation of the pathogenesis of T2D with similar expression patterns to RPL3P4.

From the fourth to eighth rules, most of the involved genes occurred in the top three rules or were the top T2D-associated genes just with different combination patterns. Specific genes, like RPL9P9 (388147) and RPL36AP21 (100271332) for rule4, MT-TS1 (4574) for rule5, RPL26P30 (653147) and RPL6P20 (285900) for rule6, RPL35P1 (440737) for rule7, and RPS3AP26 (644972) for rule8, have been identified in our quantitative rules. As we can see from such typical rule associating biomarkers, most of the genes are ribosome-associated genes like RPL3P4 (326307) as discussed earlier. Although no direct evidence confirmed the associations between such genes and T2D, it is still reasonable to speculate that such genes may play an irreplaceable role in the identification of T2D. As for MT-TS1, such genes have already been reported as potential biomarkers for T2D (Mannino and Sesti, 2012), corresponding with our prediction.

There are two potential applications for identified T2D-associated genes: 1) potential biomarkers for T2D diagnosis and monitoring; 2) potential drug target for T2D therapy.

For the identified T2D-associated genes, considering that such genes are identified from pancreatic tissues, they can reflect the original tissue alterations during T2D initiation and progression. Therefore, such genes can be used as biomarkers for direct pancreatic islet biopsy examinations. Apart from that, the candidate genes as potential drug targets can also be manually regulated to prevent the initiation and progression of T2D. Using high-throughput drug screening, antibodies or chemicals specifically targeting the candidate genes can be identified and developed as potential target drugs for T2D.

For the quantitative T2D-associated rules, although we have already identified a group of genes associated with T2D, it is still quite difficult to diagnose T2D. With specific quantitative rules, the identification of T2D patients can be more accurate and efficient. Also, the rules can also be summarized as clinical guidelines for T2D diagnosis using pancreatic tissue single-cell sequencing techniques.

As listed in Table 2, we identified quantitative rules associated with T2D. The GO enrichment analyses on rule genes were conducted. Table 4 lists the enriched GO terms of these rule genes. It was indicated that most rules are shown to be associated with ribosome-associated biological processes. According to recent publications, ribosome-associated biological processes have been widely shown to be associated with the pathogenesis of T2D. In 2019, in a metabolic study on pancreatic tissues, ribosome-associated genes have been shown to participate in the ERK/hnRNPK/DDX3X pathway in pancreatic islet cells and further regulated the initiation and progression of T2D (Good et al., 2019), consistent with our results. Apart from that, in 2020, DIMT1, as a regulator of ribosomal biogenesis has been shown to participate in the physical biological processes of pancreatic tissue, further validating our results.

In this study, for the first time, we adopted several machine learning algorithms to identify disease-specific biomarkers using the mixed single-cell sequencing data. Such analyses may not only identify biomarkers from the single-cell level, getting rid of the bias generated by the averaged transcriptomics using the bulk sequencing method, but also overcome the sample size restriction of traditional single-cell analysis. Compared with traditional single-cell analysis, we did not focus on the classification of different cell subtypes but just the patients and control subjects, improving the analysis accuracy. However, there still remain three major limitations of current analyses on pancreatic single-cell sequencing data:

1) First, the dataset we used is still a relatively small dataset, with around 20 subjects. A larger single-cell sequencing dataset may improve the efficacy and accuracy of our results.