We selected 11 ART medications commonly used in combination regimens, of which 10 are currently prescribed and one is no longer commonly used in the U.S. but at one time was the core drug for the single most common ART regimen, to experimentally evaluate their effect on the three cell lines (Figure 2). The dose selection was based on the conversion from human equivalent dose used in clinical settings to dosages applicable to in vitro settings, using the body surface area normalization method (Reagan-Shaw et al., 2008). The 11 medications were applied to the cell lines in clinically relevant conditions, with drug dosages titrated and tested for viability (Cell Titer Fluor, Promega, cat. no. G6081), cytotoxicity (LDH-Glo Cytotoxicity Assay, Promega, cat. no. J2381) and apoptosis/necrosis (RealTime-Glo Annexin V Apoptosis and Necrosis Assay, Promega, cat. no. JA1012), following manufacturer’s instructions. The final concentration was selected based on both the in vitro testing results and by referring to dosages previously used in the same cell line settings (see Supplementary Methods). The selection of the three cell lines was based on their relevancy to widely observed adverse events. The human monocytic cell line THP-1 serves as an in vitro model for foam cell formation, which is involved in atherosclerosis at all stages of development; Human Kidney-2 (HK-2) is a model for kidney disease; and hepatoblastoma-G2 (HepG2) is a model for liver disease. The human cell lines used in this study were obtained from the American Type Culture Collection (ATCC) and cultured in 384-well plates according ATCC’s instructions. Briefly, THP-1 cells (10,000 cells/well; ATCC, cat. no. TIB-202) were cultured in RPMI-40 medium (Gibco, cat. no. 11875-093) supplemented with 10% fetal bovine serum (FBSGibco, cat. no. 10438026) and 1% penicillin/streptomycin (P/S, Corning, cat. no. 30-002) at 37 °C, 5% CO₂. Macrophage differentiation was induced in vitro with phorbol-12-myristate-13-acetate (PMA, 50 ng/mL; Sigma) for 72 h at 37 °C in 5% CO₂, as described (Tsuchiya et al., 1982; Ibanez et al., 2012; Talukdar et al., 2016). HK-2 cells (5,000 cells/well; ATCC, cat. no. CRL-2190) were cultured in Keratinocyte Serum Free Medium (K-SFM), Kit (Gibco, cat. no. 17005-042) at 37 °C, 5% CO₂. HepG2 cells (5,000 cells/well; ATCC, HB-8065) were cultured in Eagle’s Minimum Essential Medium (EMEM; ATCC, cat. no. 30-2003), supplemented with 10% FBS and 1% P/S at 37 °C, 5% CO₂. After 24 h of culture for HepG2 and HK-2 cells, and 72 h for THP-1 macrophages, the cells were placed in their respective culture media supplemented with 2% FBS for 16 h. Subsequently, they were treated for 6 or 24 h with either vehicle (DMSO) or with ART medications at a concentration of 0.1, or 0.5 µM (Supplementary Table S1). For all treatments, three biological replicates were used for each cell line per condition. After treatment, cells were lysed using the TCL Buffer (Qiagen, cat. no. 1031576) and incubated for 30 min at 37 °C. The lysates were then stored at −80 °C. This study does not involve human subjects as defined by federal guidelines, and therefore did not require Institutional Review Board (IRB) review.

Gene expression was quantified using the L1000 assay, a high-throughput profiling method previously described (Subramanian et al., 2017) (Figure 2). This approach directly measures the expression of 978 landmark genes. The resulting dataset was then normalized relative to the expression of 80 invariant landmark control genes. In this assay, the expression of additional genes is inferred via a linear model trained on thousands of Affymetrix gene expression datasets from the gene expression omnibus (GEO). In this study, we utilized the L1000 data with the 12,328 measured and inferred mRNAs.

Initial data quality assessment was performed using the L1000 preprocessing pipeline. To ensure internal consistency, touchstone perturbations were included on the plates, utilizing drugs with well-characterized, robust, and highly reproducible effects, such as Vorinostat (Supplementary Methods; Supplementary Figures S1–S3). To validate the quality of L1000 gene expression profiles with the use of external data, we conducted a quality control (QC) and filtering process by comparing them to publicly available ARCHS4 (Lachmann et al., 2018) reference profiles for the same cell lines HK2, HepG2, and THP-1. This step aimed to identify and exclude L1000 samples that diverged from the expected gene expression patterns of these cell lines, preserving only reliable profiles for further analysis. Using the archs4py Python package, we retrieved ARCHS4 gene expression profiles via metadata searches with the terms “HK2,” “HepG2,” and “THP-1“ to match samples to each cell line. For every L1000 profile, we computed its correlation with the corresponding ARCHS4 profiles by focusing on the 11,888 overlapping genes and aligning plates to their respective cell lines. We then calculated the average correlation of each L1000 profile with all matching ARCHS4 samples. Next, a Gaussian mixture model (GMM) with two components (k = 2) was applied to separate the profiles into clusters of higher and lower average correlation to the ARCHS4 references. We discarded all samples in the low-correlation cluster, retaining only those with strong alignment to the ARCHS4 profiles (Supplementary Methods; Supplementary Figure S4). This process ensures that only high-quality L1000 data, consistent with established expression patterns, are used in subsequent analyses (Supplementary Figure S5).

To generate replicate-consensus signatures for drug perturbations, we first computed the modulated z-score (MODZ) for each plate after removing gene expression profiles that failed to meet previously established QC standards following the protocol used in the original L1000 publication (Keenan et al., 2018; Subramanian et al., 2017). Each drug combination was tested in triplicate. We began by determining the z-score for each drug perturbation on the plate. Next, we constructed a pairwise Spearman correlation (Wissler, 1905) matrix among all replicates, setting the diagonal self-correlations to NaN to exclude them. Next, we determined the average correlation of each replicate with the others. For the final consensus signature, we summed the z-score signatures, weighting each by its average correlation. Replicates with a negative average correlation were assigned a weight of zero, ensuring that only those with positive correlations contributed to the consensus signature. This method prioritizes reliable replicates in forming the consensus profile.

For each of the three gene expression matrices extracted from the ARCHS4 database (THP-1, HepG2, and HK-2), we first filtered genes based on expression levels, requiring a gene to have a count of at least 20 in 10% or more of the samples. This step eliminates genes with low expression and pseudogenes from subsequent analysis. After filtering, the resulting datasets included 16,395 genes across 3,439 samples for HepG2, 21,303 genes across 6,431 samples for HK-2, and 17,619 genes across 2,606 samples for THP-1. We then normalized the gene count data for library size using the Trimmed Mean of M-values (TMM) method (Baguio, 2008). Following normalization, we computed pairwise Pearson correlations for all gene pairs. To construct the gene neighborhood network, we connected genes if they ranked among the top 100 correlated genes for each other. Finally, we calculated the normalized enrichment score (NES) for each MODZ gene expression signature using the gene neighborhood network as the gene set library with blitzGSEA (Lachmann et al., 2022).

Genome-wide association study (GWAS) gene sets were obtained from the GWAS Catalog (Sollis et al., 2023) (2019) library via the Enrichr database (Kuleshov et al., 2016). Prior to performing enrichment analysis, we excluded ribosomal genes from the gene expression profiles. These genes exhibit strong correlations with one another, which can lead to their overrepresentation in enrichment results. This step is necessary because enrichment analyses typically assume statistical independence among genes, an assumption that ribosomal genes violate due to their high intercorrelations. In this study, we sought to determine whether gene sets associated with specific biological themes exhibited coordinated up- or downregulation of genes sets associated with known side effects of ART therapy. To this end, we curated groups of gene sets corresponding to common adverse events associated with ART (Supplementary Table S2). To investigate the collective behavior of gene sets, we first performed gene set enrichment analysis (GSEA) using the blitzGSEA method (Lachmann et al., 2022) on a total of 996 GWAS-derived gene sets. This initial step generated a ranking for each gene set, indicating the significance of up- or downregulation of its constituent genes within drug response signatures. Next, we organized these gene sets into biologically related categories, specifically focusing on low-density lipoprotein cholesterol (LDL) levels, body fat composition, and glomerular filtration rate; each category comprised multiple gene sets. To assess whether these thematic groups exhibited coordinated regulation, we applied a secondary enrichment analysis using blitzGSEA on the rankings of the gene sets within each category. This approach determined whether the gene sets within a group were significantly enriched, indicating that the genes in these sets were consistently up- or downregulated in the context of the drug signatures, which indicates a coordinated regulation. Furthermore, to pinpoint potential driver genes underlying these biological effects, we examined the genes most frequently appearing in the leading edges of the enriched gene sets, as identified by the blitzGSEA analysis. This dual strategy allowed us to assess thematic gene set regulation but also identify key genes that contribute to the observed phenotype.

The ARTexpress Signature Portal was developed in Python 3 using a FastAPI backend. The backend serves signature information based on gene-centric or signature-centric queries. The gene search function lists the gene activity of a gene of interest for all drug perturbations and cell lines. The signature-centric search lets a user specify the cell line, time point, and drug perturbation and returns the corresponding gene activity profile. ARTexpress supports enrichment analysis using single gene sets or paired up/down gene sets. Enrichment analysis of drug signatures is computed using blitzGSEA (Lachmann et al., 2022). For paired gene sets, enrichment is calculated separately for the up and down gene sets, and NES for both gene sets are multiplied to derive a composite enrichment score. The absolute composite enrichment score is reported. Based on directionality, ARTexpress labels drug signatures as either mimickers, reversers, or ambiguous. A signature is a mimicker of a paired gene set query if the “up genes” are upregulated in the signature and the “down genes” are downregulated. A signature is labeled as a reverser when the “up genes” are downregulated, and the “down genes” are upregulated. Otherwise, a signature is labeled as ambiguous; for example, if both the “up genes” and “down genes” are upregulated. For a single gene set enrichment, a signature is marked as a mimicker if the genes in the gene set are upregulated, and as a reverser when the genes are downregulated. Modulated z-score signatures, gene activity profiles, and tissue-specific gene networks are available for download (Supplementary Methods). ARTexpress is openly accessible from: https://maayanlab.cloud/artexpress. The source code of ARTexpress is available at: https://github.com/maayanlab/artsexpress.

A Uniform Manifold Approximation and Projection (Becht et al., 2018) (UMAP) plot of L1000 drug treatment cell line-specific activity profiles for THP-1, HepG2, and HK-2 cell lines at 6-h and 24-h time points showed that drug signatures were clustered by cell type, with no significant differences observed between the two time points (Supplementary Figure S6A). However, the availability of high-quality drug perturbation signatures was not uniform across cell lines and time points, resulting in missing drug perturbations due to removal during the QC step (Supplementary Figure S6B). For example, after filtering gene expression profiles as described in the Methods, certain drug effects, like those with Elvitegravir (EVG), only passed QC in THP-1 cell lines at the 24-h mark, showing low similarity to reference cell line profiles and being reported in one specific condition only. Signature strength analysis revealed that ART drug perturbations elicited less pronounced transcriptional responses compared to the touchstone perturbations, leading to reduced reproducibility. On average, touchstone perturbations exhibited significantly higher replicate correlation (Supplementary Figures S1–S3, S5). Some plates, such as those using HepG2 at 24-h, showed low concordance with other HepG2 cell line data collected from ARCHS4, resulting in nearly half the samples being removed in the QC step. In contrast, samples on other plates, such as HK-2 at 6-h, mostly passed QC (Supplementary Figure S4).

To investigate whether the transcriptional signatures of single drugs are preserved in the combination ART signatures (Figure 3), we analyzed the response of ART combinations compared to the single drugs within those combinations. We found some of the ART combinations that induced responses similar to those observed for single drugs included in the regimen. For example, the transcriptional response induced by Abacavir (ABC) closely resembled the response by ABC/Dolutegravir/Emtricitabine (ABC/DTG/FTC; Figure 3A), indicating that the effect of ABC dominates that of DTG and FTC in the combination. Similarly, the perturbation signature of Bictegravir/FTC/Tenofovir Alafenamide Fumarate (BIC/FTC/TAF) resembled that of BIC but inversely related to that induced by TAF (Figure 3B). Moreover, the transcriptional signature induced by the combination Darunavir/FTC/Ritonavir/Tenofovir Disoproxil Fumarate (DRV/FTC/RTV/TDF) aligned with the RTV signature while being opposite of the DRV signature (Figure 3C). Additionally, the signature induced by DTG/FTC/TAF was opposite to that of TAF (Figure 3D). In other instances, ART combination signatures were distinct from the individual perturbations within the regimen. For instance, Cobicistat/EVG/FTC/TAF (COBI/EVG/FTC/TAF) and Efavirez/FTC/TDF (EFV/FTC/TDF; Figures 3E,F) displayed transcriptional patterns that differed from those induced by the single drugs, suggesting that the combination ART regimens induce unique gene expression profiles.

Focusing on specific genes targeted by combination ART regimens across different cells, we found that ACTG1, or actin gamma 1 — a gene that encodes gamma actin, a protein crucial for the localization of the HIV-1 reverse transcription complex—was consistently downregulated by the four combination ART regimens, ABC/DTG/FTC (Figure 3A), BIC/FTC/TAF (Figure 3B), DTG|FTC|TAF (Figure 3C), and EFV/FTC/TDF (Figure 3F). Additionally, the host-HIV-associated gene, PPP2R1A—which encodes the alpha (α) subunit of Protein Phosphatase 2A (PP2A) — was downregulated with COBI/EVG/FTC/TAF (Figure 3D) and DTG/FTC/TAF (Figure 3C). Among the ten most differentially expressed genes in the ART regimens (Figure 3), we identified genes associated with metabolic and cellular processes components, including genes involved in RNA methylation (PRMT1), transcription (DDX54, CNBP), and glycolysis (ENO1, ALDOA, and PKM). Despite these findings, most of the top target genes were specific to treatment regimen. Notably, EEF1G and EEF1B2 exhibited ∼4-fold downregulation in response to DTG/FTC/TAF regimen (Figure 3D). This significant downregulation was unique to the combination signature and was not observed when DTG, FTC, and TAF were used individually.

Next, we investigated whether genes exhibiting the greatest transcriptional response across combination ART regimens influenced cellular metabolic pathways relevant to viral replication and persistence. By examining genes encoding known host-HIV interacting proteins (Chatr-aryamontri et al., 2009) within the ART transcriptional signatures, we discovered that the gene expression signatures of ART combinations, averaged across the cell types, were significantly enriched for genes encoding human proteins known to interact with HIV viral proteins (Figure 4; Supplementary Table S3). This enrichment was more pronounced for combination ART regimens compared to the individual drugs that comprise these combinations (t-test on NES values of combination vs. single drugs, p = 0.045; Supplementary Table S3), consistent with a stronger potential effect on HIV suppression.

In addition to ACTG1, which was the top downregulated gene by four combination drug regimens yet upregulated by TAF alone, four additional HIV-associated genes (YBX1, TUBB4B, ARF1, and PPP2R1A) were downregulated by at least two different ART regimens. Conversely, five single drugs, namely EFV, FTC, TDF, COBI, and RTV, were found to upregulate ORM2 (Orosomucoid 2) encoding alpha-1-acid glycoprotein 2 (AGP2). Notably, AGP2 is a key acute phase plasma protein that binds HIV envelope glycoprotein and CCR5 which plays a role in HIV infection (Rabehi et al., 1995; Atemezem et al., 2001). Certain single ART drugs such as EFV, DRV, DTG, and TAF each upregulated two or more host-HIV-interacting genes, whereas other single drugs regulated just one. Overall, across different cells, combination ART regimens appeared to target distinct HIV-associated transcriptional targets compared to individual drugs and were more likely to downregulate the expression of genes associated with the viral life cycle.

Although single ART medications are sometimes effective in suppressing viral load temporarily, they carry the potential to cause severe adverse events (Panel on Clinical Practices for Treatment of HIV Infection, 2001), which are often associated with a specific drug or a drug class. Therefore, to evaluate whether combination ART regimens have safer molecular profiles, we cross-referenced tissue-specific transcriptional signatures of drug response to single drugs and combination ART regimens with gene sets associated with lipid profiles, body weight-related traits, and renal function, the most common adverse events of ART, using the GWAS catalog (Sollis et al., 2023) (see Methods and Supplementary Table S2).

The ARTexpress signature portal enables users to browse all pre- and post-treatment gene activity data across three cell lines (HepG2, HK-2, THP-1) and 44 ART regimens (11 single agents, 20 two-way, 11 three-way, and 2 four-way combinations). Users can also browse by regimen and cell line, select specific time points, and view the correspondent transcriptional signatures summarized as modulated z-score gene activity profiles. The portal ARTexpress also performs enrichment analysis for single gene sets and paired up/down gene sets using blitzGSEA, reporting NES, composite scores for paired sets, and directionality labels (mimicker, reverser, ambiguous). All signatures, gene activity profiles, and tissue-specific gene networks can be downloaded for offline analysis.