Comprehensive searches of Gene Expression Omnibus (GEO) (14) and ArrayEpress (AE) (15) databases were performed on October 5, 2021. The search terms were as follows: (“Inflammatory Bowel Disease”) OR (“Crohn’s disease”) OR (“Ulcerative colitis”) AND (“high throughput sequencing”[Platform Technology Type]) AND “homo sapiens”[Organism] AND (“RNA-seq” OR “RNAseq” OR “RNA sequencing”). This systematic review was performed according to the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) guidelines (16). The flow chart for the selection of studies is shown in Figure 2 .

Only studies meeting specific inclusion criteria were considered for the mega-analysis. These criteria included the availability of RNA-seq raw data, transcriptomes obtained from whole tissue samples (punch biopsy), direct sequencing without any further manipulations (e.g. tissue culture or isolated cells), the use of the illumina sequencing platform, and the availability of metadata. Studies based on microarray platform, single-cell RNA-seq and small RNA-seq were excluded from the analysis. Additionally, uninflamed biopsies were excluded due to the limited size of the retrieved cohort.

We uniformly processed the RNA sequencing data from the selected eligible studies. Raw data were downloaded using the SRA-Toolkit (https://hpc.nih.gov/apps/sratoolkit.html). Briefly, we utilized the fasterq-dump tool to download FASTQ files, followed by adapter trimming using Trim Galore (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) and quality checking with FASTQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Next, we aligned the high-quality reads to the reference genome (GRCh38) using HISAT2 (17). The number of reads mapped to each annotated genes was then counted using featureCounts software (18). For the statistical analysis we retained only samples with a library size greater than 10 million counts, and genes with less than 1 count per million reads (CPM) across more than 90 samples were filtered out using the ‘filterByExpr’ function implemented in edgR package (19).

Normalization and differential gene expression were conducted using the edgeR (19) and limma (20) R packages. In brief, raw read counts were normalized using the TMM method and then a voom transformation was applied to approximate a normal distribution, resulted in a dataset with logCPM values. Prior to initiating the differential expression analysis, we performed Multidimensional Scaling (MDS) analysis to represent the data in a reduced-dimensional space which enable the visualization of batch effects associated with the dataset. We found the “study source” factor to have the largest contribution to the variation, with a notable correlation to the tissue source (small intestine vs. large intestine) ( Supplementary Figure 1 ). Following the consideration of study batch effects, we applied linear modelling to detect differential expression. Statistically significant differentially expressed genes (DEGs) were defined as those having fold change > 2 and False Discovery Rate (FDR) < 0.05.

To distinguish between UC and CD, we utilized samples from the large and small intestines, specifically from the Rectum (UC= 233, non-IBD control= 46) and from the Ileum (CD= 155, non-IBD control= 24), respectively. We added an interaction term to the linear model, in order to detect which genes respond differently to the disease in UC compared to CD. To model the tissue effect we added a comparison between the small intestine and the large intestine among the non-IBD control samples. The study feature was added as covariate to the model to account for its effect.

Linear equation used to distinguish between each disease type and control samples as well as analysis details are described in Supplementary Materials under methods section.

To identify canonical pathways and upstream regulators, we used the Ingenuity Pathway Analysis (IPA) (21) based on list of DEGs. Briefly, we performed the core analysis to identify significant canonical pathways that were enriched among the DEGs and to predict significant activated upstream regulators. Negative z score value implies an overall pathway’s inhibition and a positive z-score value suggests an overall pathway’s activation. Pathways and upstream regulators having p-value < 0.05 and absolute z‐score ≥ 2 were considered statistically significant.

We employed supervised machine learning (ML) techniques on the gene expression data for classifying various diagnoses: UC (n= 233, Rectum) vs. non-IBD control (n= 46, Rectum), CD (n= 164, Small intestine) vs. non-IBD control (n= 33, Small intestine), and IBD (388 cases, including 233 UC rectum and 155 CD ileum) vs. non-IBD control (70 samples, including 46 rectum and 24 ileum). These diagnoses were treated as the output variables. Each diagnosis was subjected to separate analysis using various ML algorithms to develop a predictive model for classification.

To ensure an unbiased evaluation of the model’s performance, we split our data into training (80%) and test (20%) sets using the train_test_split function from the Scikit-learn (22) library in Python (23). After applying feature selection, we trained and tested different classifier models, followed by validation using the GSE193677 dataset as an external and independent data source (as detailed below).

1.Feature selection

To reduce the number of features and select input variables that have the strongest relationship with the output variables we utilized the SelectKBest algorithm from the Scikit-learn library. This algorithm ranks the features based on their individual ANOVA F-values and selects the top k features with the highest F-values as the most relevant for prediction. This approach is suitable for categorical output variables, as in our study. Evaluation of the selected features was performed using a Logistic Regression classifier with different random states and cross-validation. This allowed us to assess the effectiveness of the selected k features in accurately predicting the target variable. By iterating through different k values and random states, we were able to identify the combination that yielded the highest accuracy, indicating the optimal number and significance of features. The k values and random state are described in Supplementary Table 1 .

2.Classifiers

In our analysis, we utilized six classifier models from the Scikit-learn (22) and xgboost (24) packages: SVM, KNN, Random Forest, Extra Trees, XGBoost and Naive Bayes. To optimize their performance, we implemented a nested loop to iterate over different hyperparameter values for each algorithm. These hyperparameters were tuned to enhance predictive capabilities and mitigate overfitting for each model ( Supplementary Table 1 ).

To distinguish between UC and CD, we first run DE analysis as described above then, we selected genes meeting two criteria: 1) significant interaction (fold change > 2 and FDR < 0.05) and 2) non-significance in tissue comparison (defined as fold change < 1.3 and FDR > 0.5). This resulted in 365 genes. These genes were subsequently used as features in a ML analysis. The expression matrix used as input to ML analysis contained normalized expression values that were corrected for tissue and study batch effect using removeBatchEffect() function from the limma package (20). The ML analysis, included 233 UC-inflamed Rectum samples and 155 CD-inflamed Ileum samples for training the models. For validation, we utilized data from the GSE193677 dataset, which included 60 UC-inflamed Rectum samples, 68 CD-inflamed Rectum samples and 60 CD inflamed ileum. Feature selection and ML analysis was done as described above.

As an independent validation dataset we used the transcriptomic data from the Argmann et al. study (8), available on GEO database (GSE193677). The raw data was downloaded and underwent reprocessing steps, including alignment to the GRCh38 reference genome and voom transformation. The dataset downloaded consisted of 240 samples. For our analysis, we included 60 samples from patients diagnosed with UC, showing severe to moderate inflammation, collected from the rectum, and 60 samples from patients with CD, exhibiting similar inflammation levels, obtained from the ileum. Additionally, we incorporated 120 samples from non-IBD control participants, 60 samples from the rectum and 60 samples from the ileum. We maintained an equal representation of male and female participants in all groups. The dataset was used to validate the six different classifier models. The accuracy and area under the ROC curve (AUC) were calculated to assess models’ performance.

The blood samples dataset was obtained from the Argmann et al. study (8) available on GEO database (GSE186507) which comprises blood RNA-seq data collected from whole blood samples during the participant’s endoscopy visits, along with comprehensive clinical, and endoscopic evaluations. We downloaded and reanalyzed (as described above) a total of 111 samples collected from patients diagnosed with both UC and CD. These samples were from individuals exhibiting clinical symptoms along with either severe or inactive endoscopic states. Additionally, we included samples from patients without clinical symptoms and with an inactive endoscopic state, as well as samples from non-IBD controls.

The drug repurposing analysis was conducted using the Connectivity Map (CMap) (25) through the web application CLUE (https://clue.io) (26). CMap is a collection of genome-wide transcriptional expression data responses from human cell lines that have been treated with chemical compounds known as perturbagens. These perturbagens induce changes in cellular gene expression patterns. Our aim was to identify perturbagens that induced opposing expression patterns, to the top significant 300 upregulated and downregulated genes within the IBD vs. non-IBD control DEGs, and within the comparison between each of the diseases and non-IBD control samples separately. Querying these expression patterns with the CMap reference, particularly the Gene Expression (L1000) Touchstone dataset, generated a heatmap displaying the connectivity score (tau) of 2837 small-molecule compounds as perturbagens. Our focus was on compounds demonstrating the highest anti-correlation and CMap median tau score below −95 in the HT-29 cell line column. These compounds are candidates for drug repurposing. We chose to look at HT-29 cells, since this cell line, originating from colon adenocarcinoma, serves as a relevant in-vitro model for studying processes like absorption, transport, and secretion by intestinal cells (27–29).

Statistical analyses were conducted using the R statistical framework (v.4.2.2) (30). Boxplots illustrating expression levels were generated using the ggplot2 R package (v.3.4.0) (31). Unless indicated otherwise, comparisons between different groups were calculated using t test or ANOVA, with a significance threshold set at p < 0.05. The statistical analysis involved training and evaluating six machine learning algorithms was done using Python (v.3.10.11) (23). Model performance was assessed using the accuracy_score function to compute classification accuracy, while Receiver Operating Characteristic (ROC) curves were utilized to evaluate the model’s ability to distinguish between classes, with the AUC (Area Under the Curve) serving as a measure of classifier performance. The scikit-learn (v.1.2.0) (22) library was employed for computation, and matplotlib (v.3.6.2) (32) for visualization.

Based on our inclusion criteria, 11 datasets were selected and subsequently were included in the mega-analysis. The datasets included clinical data and biopsy samples from the intestines of a total of 697 participants, comprising 569 IBD patients and 128 non-IBD controls ( Figure 2 ). Distribution of expression data pre and post normalization is given in ( Supplementary Figure 2 ).

The 11 eligible studies included in the analysis are described in Table 1 . Within our UC cohort, which is derived from eight of these studies, we obtained a total of 386 inflamed biopsies specifically from the large intestine. For the CD cohort, which is derived from four studies, we identified 183 inflamed biopsies primarily from the small intestine. Among the included studies, five provided non-IBD control samples, mostly from subjects suspected of having IBD but with radiographic, endoscopic, and histologic findings. The remaining samples obtained from normal bowel sections, located more than 10 cm away from tumors in patients undergoing bowel resection for sporadic colon cancer. Overall, the control group included 95 samples from the large intestine and 33 samples from the small intestine.

We identified 4,027 DEGs (2396 upregulated and 1631 downregulated genes, Figure 3A ) when comparing gene expression from inflamed biopsies of the UC group and the control group, both taken from the large intestine. Given the involvement of CD in both the small and large intestines, we merged data of the biopsies obtained from inflamed regions of both intestinal segments. We then compared this combined dataset to non-IBD control biopsies collected also from both the small and large intestines. By conducting this comparative analysis, we sought to gain insights into the gene expression patterns associated with CD across both tissues. We found 1,008 DEGs including 629 upregulated genes and 379 downregulated genes in inflamed samples vs. non-IBD controls ( Figure 3B ). Finally, we chose to investigate IBD as a unified group when compared to the control group. In the IBD vs. control group we identified 2099 DEGs including 1387 upregulated and 712 downregulated. ( Figure 3C ). The top ten up- and down-regulated genes are displayed in Figure 3D . Several genes consistently appear across all comparisons. Notably, a lncRNA ENSG00000254645, stands out as it ranks among the top ten downregulated genes in the IBD vs. Control comparison. This transcript has previously been associated with IBD (33). In our analysis, this transcript demonstrates significant downregulation in both the CD vs. Control and UC vs. Control comparisons ( Supplementary Figure 3 ).

We utilized IPA to identify molecular pathways enriched within the list of DEGs obtained by comparing inflamed to non-IBD control tissues in the different comparisons ( Figures 4A–C ). As expected, activated pathways in all comparisons, were enriched in immune response pathways, including pathogen induced cytokine storm signaling, phagosome formation pathway and wound healing signaling. IL-17 was also upregulated in all comparisons, although did not rank in the top five upregulated pathways in CD and IBD comparisons. The downregulated pathways observed in all comparisons encompassed various metabolism and catabolism processes. For example; serotonin degradation, nicotine degradation III, thyroid hormone metabolism II (via conjugation and/or degradation), and melatonin degradation I. Along with the metabolism pathways there are also downregulated nuclear receptor (NR) pathways which related directly to gut microbiota and play protective roles in intestinal epithelial integrity (34). Notably, LXR/RXR Activation signaling, which is a NR pathway, was downregulated only in CD and IBD comparisons, but not in UC. IPA Comparison analyses revealed activated and inhibited upstream regulators ( Figure 4D ). Upon the common activated regulators we found pro-inflammatory cytokines (IL1B, TNF), and the NFkB complex which trigger immune responses. Upon the common inhibited regulators we identified alpha catenin and the IL10RA which their inhibition is known to increase intestinal permeability (35) and promotes inflammation (36), respectively. Overall, these results show that similar pathways and upstream-regulators are involved in both diseases and suggest that similar molecular mechanisms operate in both.

A machine learning analysis aimed at distinguishing between UC and CD identified ten genes ( Supplementary Table 2 ). To verify the accuracy of the trained model, we employed two approaches. Initially, we utilized independent UC and CD samples collected from the same tissues used in the training set (rectum and ileum, respectively). This yielded an impressive AUC of 0.97 in ROC analysis (see Supplementary Figure 4A ). However, validation failed when utilizing external data from the rectum, which represents the same tissue for both diseases, resulting in AUC scores hovering around 0.6 (see Supplementary Figure 4B ). Acknowledging the limitations posed by tissue effects in discriminating between the diseases, we subsequently, developed diagnostic models to identify genes capable of effectively discriminating between inflamed samples of UC, CD and IBD each compared separately against non-IBD controls. Six ML algorithms were employed with hyperparameter tuning (see methods). We found ten genes distinguishing between UC and non-IBD samples, another 9 genes for the comparison between CD and non-IBD samples, and 34 genes for the differentiation between IBD and non-IBD samples (the genes are indicated in Figure 5A ). The most accurate models were generated for the classification of UC from non-IBD control samples, followed by the IBD classification, and finally CD classification. Specifically, the highest accuracy scores were achieved with the validation data using KNN and Random Forest. In the UC analysis, Random Forest achieved 0.947 accuracy, KNN reached 0.933. For IBD, both models achieved 0.9 accuracy. In CD, KNN had the highest accuracy of 0.9, followed by Naive Bayes at 0.89. ROC curves were generated for all models and analyses on the validation data ( Figures 5B–D ), illustrating their performance in terms of the true positive rate versus the false positive rate. The curves demonstrated good discriminatory ability, with a notable separation between the diseased and non-IBD control samples in all comparisons. Overall, the results indicate that the classification models achieved a robust and reliable performance in distinguishing between different disease groups and non-IBD controls.

Among the genes identified through our ML analyses, we recognized known immune related genes and diagnostic markers. Some of the genes were common and intersected between the different models, with AQP9 present across all of them, FCGR3A is present in both CD and IBD analyses, while SAA1, SAA2, and KCNDA3 intersect between UC and IBD analyses ( Figure 5A ). Of particular interest, IBD analysis revealed three long non-coding RNA (lncRNA) genes, two are novel transcripts, ENSG00000285744 and ENSG00000287626, along with known lncRNA MIR4435-2HG. To verify the differential expression of the six intersected genes along with the lncRNAs, box plots were generated using the normalized data and significant expression level differences were observed in both our mega-analysis dataset and the GSE193677 validation dataset ( Figures 6A, B ), with lower expression of these genes in non-IBD control samples compared to CD and UC. Notably, all genes obtained from the ML analysis showed significant expression levels differences ( Supplementary Figure 5 ).

To evaluate the potential of the ML-selected genes as non-invasive biomarkers, we assessed their expression in blood samples from IBD patients ( Supplementary Figure 6 ). Out of the 34 genes, 29 were available in the blood expression dataset. We observed a significant upregulation of 12 genes (p-value < 0.01) whose expression was higher in individuals with severe endoscopic manifestations of IBD when compared to the control group. The expression of six of these genes, namely ADGRG3, KCNJ15, AQP9, MIR4435-2HG, S100A8, and S100A9, which showed the most significant differences based on their p-values, are depicted in Figure 7 .

Finally, we utilized the Broad Institute Connectivity Map (CMap) (25) to predict potential novel therapeutic compounds based on differential gene expression results. For each the three DEG comparisons described above, we run a separate CMap analysis. In each analysis we used the top 300 DEGs, comprising 150 upregulated and 150 downregulated genes as a query to the CMap system. In addition we run another query that contained the 34 genes selected in IBD ML analysis, 30 of them were recognized by the CMap tool. These four sets of genes were employed as queries in seeking potential connections or similarities in gene expression patterns with drug-induced gene signatures present in the CMap database. It’s worth noting that the CMap system employs different cell lines for its analyses. Given the intestinal nature of our diseases of interest, we placed a specific emphasis on results derived from the HT-29 cell line.

After the tool analyzed the valid genes, it revealed 22 compounds with median_tau_score values lower than -95. This indicates a strong and significantly opposing biological effect directed towards the gene signatures associated with the various diseases (as depicted in Figure 8 ). The compounds and various families identified in the CMap suggest potential implications for IBD treatment. Notably, Antimycin-A exhibits the highest connectivity score (-99.53), along with Histone deacetylase (HDAC) inhibitors (Pyroxamide, tacedinaline, and trichostatin-a),Monoamine oxidase inhibitor (Procarbazine), BCR-ABL kinase inhibitors (Dasatinib) and RAF inhibitors (Vemurafenib, AZ-628, PLX-4720) mTOR inhibitors (WYE-354) and SYK inhibitors (Fostamatinib) MEK inhibitors (PD-0325901, MEK1-2-inhibitor) Cannabinoid receptor antagonist (O-1918), cyclooxygenase inhibitor (Valdecoxib) and Histamine receptor antagonist (Loratadine).