From chemical compounds to herbal interventions: a transfer learning-based perturbational transcriptome prediction framework for cancer drug discovery

Chemical compound perturbation data

Compound perturbation transcriptomic data were obtained from the GEO database and derived from CMAP project datasets (GSE70138 and GSE92742)³⁴. Specifically, Level 2 (GEX) L1000 processed gene expression profiles were used from these datasets, which are deconvoluted absolute expression values for 978 landmark genes generated by the L1000 high-throughput gene expression assay platform. These datasets include transcriptomic profiles collected from multiple cell lines subjected to various chemical compounds at different dosages and treatment durations with paired pre- and post-perturbation conditions.

The datasets were first integrated during preprocessing, ensuring the presence of corresponding control groups for each experimental condition. All chemical compounds were required to have valid SMILES representations in public databases, such as PubChem, to guarantee data quality and comparability³⁵. A total of 956,245 experimental and 47,653 control samples were retained after filtering, covering 18,973 chemical compounds and 82 cell lines. Each transcriptomic profile consisted of expression measurements for 978 genes, broadly representing cellular gene expression states. To enable alignment with the herbal perturbation transcriptomic data, only the intersection of genes shared between the two datasets was retained, resulting in a final set of 977 genes for downstream analyses.

Herbal perturbation data

Herbal perturbation transcriptomic data were obtained from the KORE-Map database and comprised multiple GEO datasets, including GSE244687, GSE244707, GSE244694, and GSE245912, with an additional dataset (GSE273185) included. All these datasets provide batch-corrected quasi-count expression matrices generated from RNA-sequencing of herbal medicine-treated human cell lines. The GSE273185 dataset was incorporated to expand the herbal perturbation transcriptomic dataset and the same technical standards for transcriptome sequencing and analysis as the KORE-Map datasets were adopted to ensure data compatibility, thus improving the robustness of subsequent integrated analyses. Each sub-dataset in the KORE-Map database corresponds to one unique cell line and includes perturbation profiles of multiple herbal agents on that cell line. In contrast, the GSE273185 dataset contains one unique herbal agent and includes perturbation profiles across multiple cell lines treated with this herb. Accordingly, the subsequent out-of-distribution generalization experiments for herbal agents and cell lines involve validation on entirely independent datasets. After integration, a total of 378 samples were collected, covering 4 tumor cell lines (HT29, SW1783, HepG2, and A549) and 11 herbal agents with exactly 9 samples per herb-cell line combination, ensuring a fully balanced experimental design (Table 1). HT29, SW1783, HepG2, and A549 are derived from a human colon adenocarcinoma, a human brain glioblastoma, a human hepatocellular carcinoma, and a human lung adenocarcinoma, respectively. Each combination of herbal agent and cell line was tested at doses of 500, 100, and 20 μg/mL with a fixed treatment duration of 24 h. No missing values were apparent in the dose and time metadata. Prior to model input, these experimental conditions were all subjected to standardization. To enable alignment with the chemical compound perturbation data, only genes shared between the herbal and compound transcriptomic datasets were retained, resulting in a final set of 977 genes for subsequent analyses.

Feature extraction for chemical compounds and herbal agents

SMILES representations of the compounds were retrieved from the PubChem database for chemical compound perturbation data. Molecular descriptors were then converted into molecular graphs using the RDKit toolkit. The InfoGraph algorithm³⁶ was subsequently applied to transform each molecular graph into a 300-dimensional vector representation, which effectively captured compound-level features.

Feature extraction for herbal agents followed a similar strategy. Because herbal agents are composed of multiple chemical constituents and the biological effects typically arise from the synergistic actions of these constituents, the constituent compounds for each herbal agent were first obtained from the HERB 2.0 database³⁷. Notably, only the chemical composition information of herbal agents were retrieved from the database for feature construction with no additional data (including target, pathway, or transcriptomic profile information) involved in this process. Molecular features of these constituent compounds were extracted using the InfoGraph algorithm based on graph neural networks. A 300-dimensional feature vector for each herbal agent was then generated by aggregating the feature vectors of the constituent compounds via max-pooling, in which the maximum value across all constituent compounds was retained for each feature dimension. Let X_herb ​ ∈ ℝ^d denote a herbal agent feature composed of n constituent compounds and let c_i ​ ∈ ℝ^d be the d-dimensional feature vector of the i-th compound (I = 1,2,…,n). The feature vector, X_herb, was computed via max-pooling as follows:

where X_herbj is the j-th dimension of the feature vector for the herbal agent and c_i,j is the j-th dimension of the feature vector for th i-th compound. This approach effectively preserves the most salient bioactivity signals from individual compounds, better capturing the dominant action patterns of herbal agents and providing a meaningful basis for downstream biological mechanism analysis.

Model architecture

The compound perturbation prediction algorithm developed in this study was based on a generative model with an encoder-decoder architecture and incorporated a self-attention mechanism to capture compound features and the effects on transcriptomic responses. The model encoded compound features, dosage, and treatment duration together with features derived from control samples to generate predicted perturbational transcriptomic profiles. The model consisted of three main components:

(1) Compound encoder. The compound encoder adopted a three-layer multilayer perceptron (MLP) architecture³⁸. The input consisted of 300-dimensional compound feature vectors, which were sequentially transformed through hidden layers of 256, 128, and 32 dimensions, yielding a 32-dimensional latent representation. This module was designed to extract representative latent features from compound structural information for subsequent perturbation generation.

(2) Control encoder. The control encoder had a structure similar to that of the compound encoder and employed a separate MLP to extract features from the control transcriptomic data. The input consisted of 977-dimensional gene expression profiles from control samples, which were mapped through layers of 512, 256, 128, and 64 dimensions to generate a 64-dimensional latent representation. This module was intended to capture baseline transcriptional states and model the relationship between control and perturbed conditions.

(3) Perturbation decoder. The decoder concatenated the latent representations from the compound and control encoders with compound dosage and treatment duration as inputs. These concatenated features were processed through multiple fully connected layers with dimensions of 98, 128, 256, 512, and 977 to generate a 977-dimensional predicted perturbational transcriptomic profile. In addition, a self-attention mechanism was embedded within the decoder to dynamically adjust the contributions of individual input features. This design enabled the model to focus on features that had a greater influence on transcriptomic changes and improved predictive performance. By capturing the complex dependencies among input features, the self-attention mechanism further enhanced the modeling of intricate interactions between compounds and transcriptomic responses.

Model training and validation

All input data were standardized during model training and validation to account for differences in value ranges among transcriptomic profiles, dosage, treatment duration, and compound features. Transcriptomic data were first log-transformed to reduce distributional skewness and outlier impacts. Per-gene standardization was subsequently performed using StandardScaler³⁹, which independently scales the value of each gene to a zero mean and unit variance. This preprocessing ensured comparability across all feature dimensions, alleviated transcriptomic data biases, improved model training stability and convergence, and justified the use of the mean squared error (MSE), R², and Pearson correlation coefficient as an evaluation metric.

Because multiple control samples were available for each cell line, principal component analysis (PCA) was applied to reduce the dimensionality of the 977 landmark gene expression profiles to 50 principal components (PCs) when selecting control samples. The similarity between perturbation and control samples was measured using Euclidean distance in the PCA space and the nearest control sample (1-nearest neighbor) was matched for each perturbation sample within the same cell line. This approach enabled the identification of control samples most like each perturbation condition, thereby enhancing the ability of the model to learn differences between perturbed and control states. Importantly, PCA was only used during the data preprocessing stage for control matching and the subsequent training/validation split was performed by grouping unique experimental conditions to ensure that samples from the same condition did not appear in both the training and validation sets, thus preventing any potential data leakage across splits.

The dataset was split into training and validation sets in an 8:2 ratio. To avoid data leakage and rigorously evaluate generalization of the model to unseen experimental conditions, a strict multi-dimensional stratified split strategy was adopted. Samples were grouped into unique experimental conditions based on four core attributes: cell line; perturbagen (herbal agents or compounds); dosage; and treatment duration. All samples belonging to the same experimental condition were assigned exclusively to either the training or validation set, but not both, thereby eliminating data leakage.

Model training was performed with a fixed batch size of 16,384 and a global random seed of 0 for the Python random module, NumPy, and PyTorch to ensure full reproducibility, and the MSE was used as the loss function. Parameter optimization was performed using the AdamW optimizer⁴⁰ with an initial learning rate of 1e-3 and a weight decay of 1e-4, and a CosineAnnealingLR learning rate scheduler (with T_max = 50 and eta_min = 1e-5) was applied to dynamically adjust the learning rate and mitigate overfitting. An early stopping strategy was also implemented with a patience of 20 epochs and a min_delta of 1e-4, whereby training was automatically terminated when the validation loss failed to show a significant improvement over the predefined number of epochs. In terms of model architecture, a single-head scaled dot-product self-attention module followed by LayerNorm was integrated into the decoder to refine latent features and stabilize model training.

To comprehensively evaluate model performance, the coefficient of determination (R²) and the Pearson correlation coefficient were used in addition to the MSE. Specifically, MSE quantifies the absolute discrepancy between predicted and observed values, R² reflects the proportion of variance explained by the model, and the Pearson correlation coefficient measures the linear concordance between predictions and ground truth. All metrics were calculated on the standardized gene expression scale consistent with model training. No inverse standardization was performed before MSE calculation. This ensures training-evaluation consistency and avoids bias from differences in original gene expression magnitudes. All genes underwent identical standardization, unifying expression values to the same scale and ensuring the comparability of MSE values.

Model architecture

A herbal perturbation prediction algorithm was developed based on the compound perturbation prediction model and incorporated into a transfer learning strategy. Given the limited availability of transcriptomic data on herbal perturbations, transfer learning enabled the effective reuse of knowledge learned from compound perturbation datasets, thereby reducing reliance on herbal data. The model architecture was identical to that used for compound perturbation prediction, retaining the encoder-decoder generative framework and self-attention mechanism to model herbal features and the effects on transcriptomic responses.

Model training and validation

The same standardization procedures used for the compound perturbation data were applied during training and validation to ensure feature comparability and improve training stability and predictive accuracy under transfer learning. The data were split into training and validation sets owing to the limited size of the herbal dataset using a 7:3 ratio to reduce evaluation variance and enhance reliability. The same multi-dimensional stratified train-validation splitting strategy used for the compound perturbation dataset was also applied to the herbal perturbation dataset to ensure consistency and avoid data leakage.

Given the scarcity of herbal data, one matched control sample was randomly selected from all available control samples for each experimental instance during training. This strategy mitigated the bias introduced by fixed control selection and increased dataset diversity across training iterations.

Model training followed the same protocol as that used for compound perturbation prediction, including the AdamW optimizer, MSE loss function, learning rate scheduler, early stopping strategy, and identical evaluation metrics, with the only adjustment being the batch size set to 64 to match the data scale of herbal perturbation experiments. The key distinction lay in the model initialization. Specifically, pretrained parameters obtained from the compound perturbation prediction model were used as initial weights. Through transfer learning, knowledge learned from large-scale compound perturbation data provided effective parameter initialization for herbal transcriptome prediction, thereby accelerating convergence and improving predictive accuracy, and enabling the model to achieve relatively strong performance even under conditions of limited herbal perturbation data.

Model training performance

The herbal perturbation prediction model, referred to as Model_H, was evaluated based on 20 repeated experiments. The average performance metrics were an MSE of 0.1395 (0.1311, 0.1479), R² of 0.8561 (0.8484, 0.8638), and a Pearson correlation coefficient of 0.9258 (0.9218, 0.9298; Figure 3). These results indicated that Model_C achieved strong predictive performance on herbal perturbation data after transfer learning with accurate fitting to experimental observations and a high correlation with true transcriptomic profiles.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3

Performance of the herbal perturbation prediction model. (A) Loss curves during training and validation. (B) Scatter plot of predicted vs. observed transcriptomic values. Figure created using Python 3.10.14.

These findings suggested that Model_H effectively leveraged knowledge transferred from the compound perturbation model despite the limited size of the herbal dataset, thereby improving prediction accuracy and alleviating data scarcity.

Baseline model performance comparison

Comparisons were made against three baseline models to further evaluate the performance of Model_H: logistic regression (LR_H)⁴¹; random forest regression (RF_H)⁴²; and support vector regression (SVR_H)⁴³. Each model was trained and evaluated based on >20 repeated experiments. Boxplots of the three evaluation metrics (MSE, R², and Pearson correlation coefficient) are shown in Figure 4A; the mean values are summarized in Table 2. The results demonstrated that Model_H consistently outperformed all baseline models across all evaluation metrics, indicating a clear performance advantage.

• Download figure
• Open in new tab
• Download powerpoint

Figure 4

Performance comparison of herbal perturbational transcriptome prediction models. (A) Baseline model performance comparison. Model denotes the proposed method, LR denotes the logistic regression model, RF denotes the random forest regression model, and SVR denotes the support vector regression model. Each model was evaluated over 20 repeated experiments. (B) Transfer learning performance comparison. Transfer learning was applied using different proportions of herbal datasets, with each experiment repeated 20 times. LR, logistic regression; MSE, mean squared error; RF, random forest regression; SVR, support vector regression. Figure created using Python 3.10.14.

Transfer learning performance evaluation

Models were trained using different proportions of the herbal dataset to assess the effectiveness of transfer learning and compared to models trained from scratch using herbal data alone. When the proportion of the herbal dataset was relatively large, the performance of the transfer learning model was comparable to that of the model trained from scratch (Figure 4B). However, as the number of herbal datasets decreased, the transfer learning model exhibited a clear performance advantage. This advantage was particularly pronounced when the number of herbal samples was very small, demonstrating the effectiveness of transfer learning under data-scarce conditions. Furthermore, the results showed that when the herbal dataset was limited, the box width of the boxplot for the transfer learning model was narrower than the model trained from scratch, indicating superior stability. In summary, transfer learning effectively leveraged knowledge obtained from compound data to improve the prediction accuracy of the model when herbal data were scarce.

In addition, when the proportion of the herbal dataset was zero (i.e., when Model_C was directly applied for prediction), the MSE and R² of the model were close to the MSE and R² of a random model. This result suggested that the transfer model could not adequately adapt to the characteristics of herbal medicines in the complete absence of herbal data. However, the Pearson correlation coefficient remained relatively high, indicating that the transfer model could partially capture the trends in gene expression changes induced by herbal interventions. This phenomenon highlights the specificity of herbal data and implies that the performance of transfer learning in practical applications improves rapidly in the early stages as the volume of herbal data increases.

Out-of-distribution (OOD) generalization performance

Systematic OOD generalization experiments were performed to rigorously evaluate the generalization ability of the model on completely unseen data⁴⁴ in two scenarios with independent validation datasets adopted for core external OOD validation: unseen-herb testing; and unseen-cell-line testing (Figure 5A–D). The fully independent external dataset (GSE273185), which only contained Aconiti Lateralis Radix Praeparata was used as the dedicated independent validation dataset for the OOD test of this unseen herbal agent. Strict leave-one-herb-out OOD cross-validation was performed for the remaining 10 herbal agents based on the KORE-Map and GSE273185 dataset. Each cell line was validated using a complete, independent dataset (or combined independent datasets) as the exclusive independent validation set for the unseen-cell-line OOD testing with no cross-over of datasets between model training and testing. The independent validation datasets for each cell line were as follows: A549 (GSE244687); HepG2 (GSE244707); HT29 (GSE244694 + GSE273185); and SW1783 (GSE245912 + GSE273185).

• Download figure
• Open in new tab
• Download powerpoint

• Download figure
• Open in new tab
• Download powerpoint

• Download figure
• Open in new tab
• Download powerpoint

Figure 5

Out-of-distribution (OOD) generalization performance of the herbal perturbational transcriptome prediction model. (A) Performance metrics from batch-wise testing on unseen herbal agents, repeated 100 times. (B) Performance metrics from batch-wise testing on unseen cell lines, repeated 100 times. (C) Performance metrics from leave-one-herb-out testing, repeated 20 times. (D) Performance metrics from leave-one-cell-line-out testing, repeated 20 times. (E) t-SNE dimensionality reduction and clustering plots of the herbal perturbational transcriptome, annotated by herbal category, dose, and cell line category, respectively. MSE, mean squared error; OOD, out-of-distribution; t-SNE, t-distributed stochastic neighbor embedding. Figure created using Python 3.10.14.

In the batch-wise unseen-herb tests the herbal dataset was split into training and validation sets at a 7:3 ratio and testing was repeated 100 times. The results indicated that Model_H exhibited favorable generalization performance when applied to unseen herbal agents with average MSE, R², and Pearson correlation values of 0.1824 (0.1725, 0.1923), 0.8131 (0.8045, 0.8217), and 0.9043 (0.8997, 0.9089), respectively. These findings demonstrated the ability of the model to generalize to herbal agents not observed during training. One herbal agent was used as the validation set in the subsequent leave-one-herb-out tests, while the remaining agents used for training and evaluation were repeated 20 times. Prediction performance remained satisfactory for most herbal agents with relatively favorable MSE, R², and Pearson correlation values, indicating effective learning of transcriptomic features. However, the model exhibited relatively weak performance on Aconiti Lateralis Radix Preparata. The analysis revealed that this herbal medicine was derived from an additional dataset outside the KORE-Map database and potential batch effects may have reduced the performance of the model compared to other datasets. Nevertheless, the Pearson correlation coefficient for this herbal medicine still reached approximately 0.76, indicating that the model can accurately predict the trends in transcriptomic changes induced by herbal interventions.

In contrast, model performance declined when the same evaluation strategy was applied to unseen cell lines, particularly with respect to R², which showed a substantial deviation from 1. This result indicated that accurate prediction became more challenging when cell lines were not observed during training. Nevertheless, the Pearson correlation coefficients remained moderate, suggesting that the model was still able to capture overall expression trends, although prediction accuracy at the gene expression level was constrained by cell line-specific biological differences. These limitations were likely due to pronounced biological heterogeneity across cell lines, leading to substantial divergence in transcriptomic profiles and increased difficulty in generalization. Similar trends were observed in the leave-one-cell-line-out tests, in which one cell line was used as the validation set and the remaining cell lines were used for training; evaluation was repeated 20 times. Model performance for unseen cell lines remained suboptimal in these tests. Although the Pearson correlation indicated partial capture of overall expression trends, prediction accuracy was still affected by strong cell line-specific biological characteristics.

To further investigate these observations, t-distributed stochastic neighbor embedding (t-SNE) was applied to the herbal perturbation transcriptomic data (Figure 5E). The results showed that herbal perturbation profiles from different cell lines were widely separated in the low-dimensional space, indicating substantial differences in perturbation responses across cell lines. In contrast, perturbation profiles remained relatively close within the same cell line, even across different herbal agents and dosages. These findings suggested that cell line identity had a critical role in perturbation response patterns and should be explicitly accounted for in perturbational transcriptome prediction models.

Overall, the model demonstrated stable and accurate generalization to unseen herbal agents, producing reliable predictions for herbal samples that were not included during training. However, performance variability was observed for unseen cell lines and for samples with strong cell line-specific characteristics. These results indicated that, the model could effectively capture latent patterns in herbal perturbation data under data-limited conditions, whereas cell line specificity and sample size remain key factors influencing generalization performance.

Case studies

Antitumor and cancer-suppressing herbal medicine validation

Six categories of antitumor herbal medicines reported in the existing literature were selected as the test set in this study to verify the reliability of the herbal perturbational transcriptome prediction model in dissecting the mechanisms of action of antitumor herbal medicines⁵¹. These categories included primary therapeutic drugs, broad-spectrum antitumor drugs, antitumor drugs for liver cancer and hepatic system tumors, antitumor drugs for lung cancer and pulmonary system tumors, antitumor drugs for colorectal cancer, and antitumor drugs for brain tumors. Paired tests were performed for each category of herbal medicines and the corresponding tumor cell lines (e.g., liver cancer herbal medicines with the HepG2 cell line and brain tumor herbal medicines with the SW1783 cell line) using the prediction model.

First, the model was used to predict the transcriptomic perturbation effects of the herbal medicines on the target cell lines with each experiment repeated 20 times. DEGs following herbal intervention were then identified. GO and KEGG functional enrichment analyses were subsequently performed on the identified DEGs using Fisher’s exact test. The genome-wide human gene set was used as the background for enrichment analysis to capture the general functional bias of DEGs in the context of the entire human genome. The Benjamini-Hochberg method was applied for multiple testing correction of enrichment P-values. Tumor-associated pathways were filtered using core pathway keywords related to “tumorigenesis-proliferation-metastasis-drug resistance” (e.g., “cell proliferation,” “apoptosis,” “PI3K-Akt,” and “metastasis”), and the number of significantly enriched tumor pathways for each herbal medicine was counted. The final results are presented in Figure 6A.

• Download figure
• Open in new tab
• Download powerpoint

Figure 6

Case validation results. (A) Number of cancer-related pathway enrichments for antitumor and cancer-suppressing herbal medicines across four cell lines. (B) Pathway enrichment results predicted from the perturbational transcriptome of each single herbal medicine in Jingling Oral Liquid. (C) Heatmap showing the occurrence frequencies of differentially expressed genes in the predicted pathways of each single herbal medicine in Jingling Oral Liquid. Figure created using R 4.4.1.

The results showed that 60% of the herbal medicines enriched >30 pathways, reflecting the “multi-target and multi-pathway” intervention mode of herbal medicines. Moreover, primary therapeutic and broad-spectrum antitumor drugs enriched a considerable number of tumor-associated pathways across different cell lines, embodying the principle of “treating different diseases with the same method” in TCM. These findings indicated that the herbal perturbational transcriptome prediction model could predict the interactions between antitumor herbal medicines and the corresponding tumor cell lines, providing a computational prediction tool for association analysis of the “cell line-action pathway-clinical indication” relationship of antitumor herbal medicines.