Whole blood transcriptional profiles of 120 patients with SLE and 58 healthy individuals (4) were analyzed. Disease activity at the time of blood sampling was assessed by the modified Systemic Lupus Erythematosus Disease Activity Index 2000 (SLEDAI-2K), after exclusion of the serologic features (anti-dsDNA and complement levels) (clinical SLEDAI) (5). Remission was defined as a clinical SLEDAI-2K = 0 and daily prednisolone dose of ≤5 mg (6, 7). Active disease was defined as a clinical SLEDAI-2K ≥4. Irreversible organ damage was assessed using the SLICC damage index (SDI) (8). Lupus Severity Index was calculated for each patient (9).

We employed CoCena² (construction of co-expression network analysis-automated, https://github.com/UlasThomas/CoCena2), using the 10,000 most variable genes as input, to determine modules of co-expressed transcripts. Next, agglomerative hierarchical clustering of patients, based on their group fold changes (GFC) for each cluster of co-expressed genes, defined the disease molecular endotypes. Functional enrichment analysis was performed using clusterProfilerR package (10).

DEseq2 was used to identify differentially expressed genes (DEGs) specific for each patient’s endotype (11). We obtained gene expression signatures of drugs that are incorporated in the treatment recommendations for SLE (12), or have failed to reach SLE clinical trials endpoints and are included in the following iLINCS sublibraries: i) iLINCS chemical perturbagens (LINCSCP); ii) iLINCS targeted proteomics signatures (LINCSTP); iii) Disease-related signatures (GDS); iv) Connectivity Map signatures (CMAP); v) DrugMatrix signatures (DM); vi) Transcriptional signatures from EBI Expression Atlas (EBI); vii) Cancer therapeutics response signatures (CTRS); and viii) Pharmacogenomics transcriptional signatures (PG). These were downloaded using the iLINCS API (https://github.com/ucbd2k/ilincsAPI/blob/master/usingIlincsApis.Rmd). Statistically significant DEGs from each drug signature were ordered by decreasing fold change magnitude. The top 300 DEGs were selected and upregulated/downregulated genes were identified. Gene set enrichment analysis (GSEA) was performed using fgsea R package (13). To determine the optimal number of drug clusters for k-means clustering, the elbow method was applied.

Drug repurposing results were prioritized using the bioinformatic tool CoDReS (Computational Drug repositioning score) (14), which enables the exploration of compound drugability, based on an algorithm that combines functional and structural scores. The functional score quantifies the pharmacodynamic potential of a compound by assessing its association to SLE hallmarks. This potential includes the binding affinity to SLE molecular targets (enzyme, receptor, transcription factor, etc.), as well as the overlap of its genomic targets with genes implicated in the pathogenesis of the disease. The structural score pertains to the pharmacokinetic properties of compounds and contains information related to the hydrophilic-lipophilic balance, solubility, permeability, as well as oral bioavailability of a drug candidate, based on the “Lipinski rules of 5” (15) and “Veber’s rule” (16).

Applying the CoCena² pipeline, we identified nine modules of co-expressed transcripts illustrated with different colors in Figure S1 . Hierarchical clustering of samples according to each module’s group fold changes (GFC) reassigned patients into five groups (G1 to G5) ( Figures 1A, B ). To define disease-driving molecular mechanisms, we investigated the CoCena²-derived modules enrichment in each patient group ( Figures 1C, D ). Interestingly, groups displayed distinct enrichment patterns, each exhibiting unique major module predominance. Platelet activation and hemostasis were identified as two group 1 specific signals (G1, “Hemostasis” group), overrepresented in the orchid module. Detailed functional enrichment analysis of the dark-grey module revealed that autophagy-associated signatures were prominently enriched in patient group 2 (G2, “Autophagy” group). Macroautophagy disturbances in G2 are accompanied by deregulation of pathways involved in neutrophil activation and toll-like receptor (TLR) cascade. Combined enrichment of the pink module, linked to aberrancies of mRNA splicing and mRNA surveillance mechanisms, and the dark-orange module, implicated among others in mitochondrial dysfunction, efficiently distinguished group 3 (G3, “Metabolism” group). Heightened expression of the indian-red module, which predominantly consists of genes implicated in neutrophil activation and degranulation, defines group 4 (G4, “Neutrophil” group). Enrichment of the dark-green module, which comprises genes (such as CD38, BLNK, IGHA1, TNFRSF17, CD22, CD79A, MS4A1, IGHD) linked to B-cell and plasmablast-mediated responses, was indicative of group 5 (G5, “B cell” group). Interestingly, G5 displays a concurrent increased expression of the steel-blue module, which is associated with type I interferon (IFN) signaling.

To evaluate the clinical implications of molecular endotype characterization, we next assessed each group’s clinical features, including demographics, clinical manifestations, serologic features, and administered treatments. The “Neutrophil” group (G4, n= 11, 9.1% of the total cohort) almost uniformly encompassed patients with active lupus nephritis (n=9/11) ( Figure 2 ). Patients of this cluster also exhibited high serologic and clinical activity; the majority were treated with cyclophosphamide at the time of blood sampling ( Figures 2 , S2 , S3 ). The “B-cell” group (G5, n=18, 15% of the total cohort) was characterized by high prevalence of constitutional symptoms. Although statistical significance was not reached, a tendency to a higher frequency of hematological and neurological manifestations was apparent in this cluster. Mucocutaneous and musculoskeletal manifestations were most common in the “Metabolism” group (G3, n=30, 25% of the total cohort), occurring in 63% and 50% of patients, respectively, while a history of neuropsychiatric SLE (NPSLE) was reported in 27%. Interestingly, the clinically heterogenous “Hemostasis” group (G1, 24.2% of the total cohort) was characterized by high frequency of male patients, while Disease Modifying Anti-Rheumatic Drugs (DMARDs) were the most commonly used therapy. Finally, the “Autophagy” group (G2, n=32, 26.7% of the total cohort) consisted of patients with mild to moderate SLE. Accordingly, photosensitivity and malar rash were found in 59,3% and 81,2% of the patients of G2, respectively.

To explore personalized therapeutic solutions, we identified compounds tailored to each group’s molecular fingerprint. This was achieved through leveraging our CoCena² based co-expression analysis, to establish an in silico, signature-based, drug prediction pipeline. As group-specific signatures, we employed the DEGs resulting from the comparison of each SLE endotype with a pool of 58 healthy controls.

To this end, we initially collected the transcriptional profiles corresponding to cellular responses against drugs that are either currently used in clinical practice, are or have failed in SLE clinical trials and are listed in the iLINCS prediction databases ( Table S1 ). Our query returned 3,900 drug signatures ( Table S2 ). Using SLE group-specific transcriptional profiles as input, we performed GSEA against the datasets of the top upregulated and top downregulated DEGs for each drug signature and normalized enrichment scores (NES) were defined. Next, we calculated the difference (ΔNES) of the NES from the downregulated gene set and the NES from the upregulated gene set for each drug signature per SLE cluster ( Table S3 ). Accordingly, a positive ΔNES indicated compounds that were predicted to reverse the group-specific transcriptomic aberrancies. To determine endotype-specific drug candidates, we applied k-means clustering, in order to group drug signatures according to ΔNES ( Figures S4 , S5 ). Drug signatures with the highest ΔNES within each drug cluster induce cellular transcriptional alterations which most efficiently counteract group-specific SLE signatures.

In G5, the top signatures were linked to azathioprine (ΔNES=2.76) and ixazomib (ΔNES=2.67) ( Figure 3A ), whereas in G2 to the proteasome inhibitor bortezomib (ΔNES=2.84). Signatures related to the SYK kinase inhibitor tamatinib (ΔNES=2.81) were top ranked in G3 subgroup ( Figure 3B ). In G4 group, signatures related to bortezomib occurred in high frequency (76%) among the top 50 signatures, starting with a ΔNES score 2.54 and, together with the calcineurin inhibitor cyclosporine (ΔNES score 2.49), might represent alternative G4-specific therapeutic options ( Figure 3C ). Finally, in both groups 4 and 5, signatures related to vitamin D derivatives (such as seacalcitol) prevailed, with a ΔNES score 3.04 and 2.97, respectively.

Since the majority of the G4 patients were treated with cyclophosphamide at sampling, an agent that could drastically alter the whole blood transcriptional landscape, we divided G4 into two subgroups; one treated with cyclophosphamide (G4A, n=6/11) and a “cyclophosphamide-free” subgroup (G4B, n=5/11) and we applied the drug prediction pipeline. In accordance with our initial findings, bortezomib was overrepresented among the top 10 signatures in both subgroups ( Tables S4 , S5 ).

Finally, we sought to propose new SLE therapeutic agents. To this end, we used a drug repurposing pipeline identifying patented compounds with potentially unrecognized efficacy in SLE. Using the iLINCS and CLUE platforms, we identified novel compounds that could reverse the previously defined SLE group-specific signatures. To sort out the top perturbagens derived from the iLINCS platform, we applied a negative concordance score cut-off of ≤ -0.5. Regarding the CLUE based analysis, only compounds exhibiting an inhibitory score of ≤ -50 were selected. Lastly, group-specific perturbagens were determined, as shown in the Venn diagram ( Figures 4A, B ). To enhance the performance of our approach, group-specific compounds were ranked, according to their druggability (“druggability prediction”). For this purpose, we used the bioinformatic tool CoDReS (Computational Drug Repositioning Score) (14). Uploading the iLINCS- and CLUE-derived compound lists (which were related exclusively to each SLE endotype) to the CoDReS platform resulted in the re-ranking of the repurposed drugs, according to their biological and pharmaceutical potential ( Tables S6 – S15 ).