Development of machine learning models for the prediction of the skin sensitization potential of cosmetic compounds (2025)

Introduction

An important endpoint for consumer and occupational safety assessment is chemical-induced skin sensitization, which is an allergic reaction caused by repeated skin exposure to a single substance or mixture (United Nations, 2017). Continuous skin sensitization can lead to allergic contact dermatitis (ACD), a severe skin condition that affects consumer health and the experience of using cosmetic products (Kimber et al., 2002). Following the implementation of Cosmetics Regulation 1223/2009 and similar legislation in 2009, traditional animal testing for cosmetic products and ingredients has been progressively banned. As a result, there has been a growing demand for non-animal testing methods, with increasing requirements for accuracy (European Union, 2009). Based on the complexity of the biological processes involved in allergic reactions, the OECD Test Guidelines break down the allergy Adverse Outcome Pathway (AOP) into four key events (KEs): (i) covalent binding of electrophilic substances to nucleophilic sites in skin proteins initiates downstream activation (Enoch et al., 2011); (ii) keratinocyte activation (Emter, Ellis & Natsch, 2010); (iii) dendritic cell (DC) activation (Ashikaga et al., 2010); (iv) proliferation of T cells in the lymph nodes (Gerberick et al., 2007). From key events (i) to (iii), various methods such as direct peptide reactivity assay (DPRA), KeratinoSens^TM and the human cell line activation test (h-CLAT) have been developed, but their accuracy remains unsatisfactory. To improve the accuracy, the following integrated approaches to testing and assessment (IATA) methods have been established (Rovida et al., 2015).

The Genomic Allergen Rapid Detection (GARD) method employs transcriptome wide microarray analysis in a myeloid cell line to assess 200 transcriptomic biomarkers, which function as predictive signatures (Johansson et al., 2011, 2014). Microarrays, which provide limited gene-related data, use pre-designed probes to capture specific mRNA sequences (Karakach et al., 2010). With the advancement of next-generation sequencing (NGS) technology, RNA-Seq has been widely applied to protein-coding genes, non-coding RNA, and immune genes, with the cost and sequencing time rapidly decreasing (Qiao et al., 2024). The researchers used RNA-Seq to analyze the changes in the expression of the THP-1 gene and molecular mechanism under different treatment conditions (Zhang et al., 2016). Weighted gene co-expression network analysis (WGCNA), which can cluster genes with similar expression patterns and examine the relationships between these modules, is widely applied in studies investigating the association between phenotypic traits and gene expression (Langfelder & Horvath, 2008).

Machine learning (ML), a branch of artificial intelligence, comprises numerous algorithms such as support vector machine (SVM) (Cortes & Vapnik, 1995), random forest (RF) (Breiman, 2001), K-nearest neighbors (KNN) (Begum, Chakraborty & Sarkar, 2015) and linear discriminant analysis (LDA) (Goksuluk et al., 2019). With the increase in data and advances in computing power, ML has found extensive applications in RNA-Seq data analysis, including gene expression data classification (Wang et al., 2018), feature selection (Li et al., 2016), differential gene expression analysis, and the development of classification and prediction models (Johansson et al., 2014). The GARD method achieves high sensitivity and accuracy by using SVM to classify and predict tested compounds. It has been established as one of the OECD 442E standard methods (OECD, 2024). Co-culture systems offer several advantages over monoculture-based assays, including facilitating intercellular crosstalk and enabling the combination of multiple key events in a single assay (Thélu, Catoire & Kerdine-Römer, 2020). To date, numerous THP-1 co-culture methods have been employed, primarily involving the co-culture of different cell types and the co-culture of reconstructed human epidermis. Hennen developed a co-culture model comprising HaCaT keratinocytes and THP-1 cells, which has demonstrated excellent potential for identifying skin sensitizers (Hennen & Blömeke, 2017). Galbiati et al. (2020) showed that keratinocytes promote DC activation by co-culturing the epithelial-like cell line NCTC2544 with THP-1 cells. Schellenberger et al. (2019) combined RHE in coculture with THP-1 cells placed underneath RHE, applying compounds using topical exposure to avoid issues related to water solubility.

In this present study, we focused on the evaluation of RNA-Seq analysis of THP-1 in co-cultured system to predict sensitizing and non-sensitizing compounds. The co-cultured combines two different cell types (HaCaT transfected with a luciferase plasmid driven by the AKR1C2 promoter and THP-1) that address KE2 and KE3 of AOP. The overall experimental design, as depicted in Fig. 1. In the first step, we co-cultured HaCaT and THP-1 cells; compounds were added, and cytotoxicity was assessed. In a second step, we collected the THP-1 cells, extracted the RNA and performed transcriptome sequencing. Finally, data processing, machine learning model construction and prediction were performed. Based on the results of this study, we established a robust experimental design that compared different data processing methods, feature selection, and various machine learning models to evaluate their impact on prediction accuracy. These results provide valuable insights for future applications.

Materials and Methods

Results

Discussion

In this study, eight machine learning models were utilized to predict the allergenic potential of chemical compounds. Predicting the potential for skin sensitization is crucial for the risk assessment of the use of these compounds. The RNA-seq data utilized in this study were obtained from THP-1 cells derived from HaCaT/THP-1 co-culture system. In recent years, RNA-Seq technology has matured significantly and has become a routine method for detecting gene expression changes (Wang, Gerstein & Snyder, 2009). We utilized RNA-Seq data as input for our machine learning models to explore new perspectives that could be used in allergy prediction.

Among common sensitizers, pre-haptens and pro-haptens (Table 1) require spontaneous air oxidation or enzyme-mediated activation to convert into sensitizers, as keratinocytes cells containing key metabolic enzymes such as cytochrome P450 (Eskes et al., 2019; Swanson, 2004). In the co-culture system involving THP-1 and HaCaT/keratinocytes cells (Hennen & Blömeke, 2017), which represent KE2 and KE3 in AOP for skin sensitization. Eskes et al. demonstrated that HaCaT cells co-cultured with THP-1 cells have the potential to increase the response to pro-haptens chemicals (Eskes et al., 2019). In contrast, Sawada et al. (2022) utilized keratinized normal human epidermal keratinocytes (NHEK) and THP-1 cells. We linked the AKR1C2 promoter to luciferase plasmid, electroporated plasmid into HaCaT cells, and selected the cells, which were used to represent KE2. The HaCaT cells can also be used in OECD 442D luciferase assays to assist in verifying skin sensitization. Most skin sensitization assays involving THP-1 co-culture systems use flow cytometry to detect CD86 and CD54, which rely on single or few markers for classification (Nukada et al., 2011). Therefore, we used RNA-Seq for comprehensive analysis of co-cultured THP-1 transcript expression, which can help identify better biomarkers and provide deeper insights into the underlying mechanisms. Simultaneously, to broaden the model’s applicability and mitigate bias due to the limited sample size, we conducted experience on a more diverse selection of chemical compounds including varied molecular structures, antigenic properties, and potency categories.

Microarrays based on hybridization of fluorescently labeled cDNA to pre-designed probes on a chip have been widely used for clustering and classification of continuous expression data (Stears, Martinsky & Schena, 2003). The initial GARD workflow used Affymetrix’s RMA algorithm for normalization (Johansson et al., 2013). Subsequently, the advanced GARD workflow converted microarray data to the NanoString nCounter platform to facilitate better assessment of biomarker characteristics (Forreryd et al., 2016). Robinson, Wang & Storey (2015) discovered that PCR quantitative validation of microarrays may exhibit greater systematic bias than RNA-seq, particularly for genes with low expression intensity. Due to superior performance, integrative RNA-seq analysis and machine learning methods are applied to explore the key genes, disease classification (Bostanci et al., 2023), and potential mechanisms.

In the context of machine learning analysis of RNA-Seq data, feature extraction is a critical step. When selecting feature extraction methods, factors such as dataset size and complexity, computational efficiency, and diversity of results are important concerns. Methods such as recursive feature elimination for gene selection (Escanilla et al., 2018) and PCA for dimensions reduction to extract important features are viable options. In this study, we compared three methods: the most basic variance (var100), DE analysis, and WGCNA. Var is the most direct and simplest method, focusing purely on the data perspective. DE analysis is a conventional approach for group comparisons. Figure 3A shows the PCA distribution where the samples are not completely separated from the perspective of PC1 and PC2, which is generally consistent with the result of Henrik. Figures 3B and 3C presented heatmap and volcano plot related to the distribution of differentially expressed genes. Figures 3D and 3E show the enrichment of differential genes in pathways such as cell chemotaxis, granulocyte migration, and MAPK signaling pathway, which are all related to DC maturation and allergy (Miyazawa et al., 2007).

Many algorithms have been proposed to classify based on continuous log intensities obtained from microarray experiments. It is reported that RNA sequencing data, derived from genome alignment raw counts, are excessively dispersed, with the count variance exceeding the mean (Goksuluk et al., 2019). Therefore, using appropriate data processing can not only enhance accuracy but also enable analysis in conjunction with microarray algorithms. In this study, we applied machine learning algorithms to the classification of RNA-Seq data using the MLSeq (2.22.0) R package. Normalization methods containing deseq median ratio (deseq) and trimmed mean of M means (TMM) were used to correct systematic variations biases in RNA-Seq data. After completing the normalization process, NBLDA and PLDA can be used directly for classification. Witten assumed that RNA-Seq counts follow a Poisson distribution and consequently proposed sparse PLDA as a classification algorithm for RNA-Seq data (Witten, 2011). However, PLDA performs deteriorates when the counts are over-dispersed. The results from Tables 1 and 2 show that among the three feature extraction methods, only WGCNA achieved 90% performance with PLDA, while the others performed relatively poorly. Dong used edgeR to perform the gene selection for NBLDA classifier and He also noted the use of network information of pathways for better testing. In our results, we found that the WGCNA feature processing method combined with the DESeq normalization method achieved relatively good results (NBLDA).

Common transformation techniques included variance stabilizing transformation (vst) (differential expression), logarithm of counts per million reads (logcpm) (Robinson, McCarthy & Smyth, 2010) and variance modelling at observational level (voom) (Law et al., 2014) that can be used to make discrete RNA-seq data hierarchically more similar to microarray data. VoomNSC is a sparse classifier that models the relationship between the mean and the variance using the voom method. It incorporates the precision weights from voom into the NSC classifier through weighted statistics (voomDDA). Dudoit extended diagonal discriminant classifiers to voomDLDA and voomDQDA, which use all features for classification (Dudoit, Fridlyand & Speed, 2002). In the research results, voomDQDA and voomDLDA demonstrated better training performance compared to voomNSC, regardless of the feature selection method or data processing approach used. Notably, voomDQDA consistently achieved accuracy rates above 85%. We believe this advantage stems from two main factors: the robustness of the voom method and the strengths of diagonal algorithms. Traditional machine learning methods such as KNN, RF, and SVM were also used for testing. Both GARD and the improved GARD utilized SVM as the classifier, demonstrating SVM excellent performance. The effect of different kernels (Radial and Poly) on test and train accuracy was also compared. The results showed that the SVM accuracy was below 90% when using the vars feature selection combined with deseq-vst data processing. Additionally, the Poly kernel required a longer model construction time. In Forreryd’s (2016) study, the RF model achieved an accuracy of 100%. However, in our study, RF reached 100% accuracy only when using deseq-vst combined with WGCNA for feature selection. Balanced accuracy (Table S3) can balance sensitivity and specificity, providing a more equitable view of model performance in the presence of class imbalance. Although we have employed various models and rigorous validation techniques to mitigate the limitations posed by the sample size and to ensure that the model can capture relevant patterns, further testing of more molecules is still needed to ensure broader ranges.

Conclusions

Our research presents an integrated approach combining co-culture systems, RNA-Seq, and machine learning techniques to provide a more comprehensive assessment of allergenic compounds. This multifaceted strategy offers a holistic perspective, high-throughput data, and advanced analysis capabilities. However, we acknowledge several challenges, including sample size limitations, model accuracy, and model interpretability. To date, there has been no alternative testing method that can sub-categorize skin sensitizers into UN GHS subcategories 1A and 1B. By pursuing these future directions, we aim to enhance the robustness, accuracy, and applicability of our approach, ultimately contributing to safer consumer products and more effective risk assessment strategies.

Supplemental Information