Transcriptome-Conditioned Personalized De Novo Drug Generation for AML Using Metaheuristic Assembly and Target-Driven Filtering

Raw exon-level expression data was first mapped to gene-level profiles to provide a biologically interpretable format. Using the UNC v2 exon-to-hg19 probe mapping reference, exon identifiers were aligned with their corresponding Human Genome Organisation (HUGO) gene symbols. To derive a single expression value for each gene, after computed the arithmetic mean of all exon probes associated with that specific gene locus. This step condensed the high-dimensional exon dataset into a gene-level expression matrix: $M_{\mathrm{genes}\times\mathrm{samples}}$ , ensuring that transcriptional activity was represented per functional gene unit.

It was essential to exclude non-coding RNA and pseudogenes. The metohed is utilized the ”mygene” Python library to query the type_of_gene field for all gene symbols in the aggregated dataset against the human reference database. A strict inclusion filter was applied to retain only genes classified as ”protein-coding.” This filtration step ensures that all subsequent candidate biomarkers represent physically translatable structures capable of forming ligand-binding pockets.

To account for the heteroscedasticity inherent in RNA-seq data and to focus the analysis on biologically distinct signals, the filtered dataset underwent normalization and feature selection:

The bimodal distribution is characteristic of high-quality data: the left peak represents low-abundance background noise, while the right peak represents actively expressed genes. The lack of deviant curves confirms the absence of major batch effects. Low-Expression Filtering is used to reduce technical noise, genes with a mean log-transformed expression value below a threshold of 0.5 across the cohort were discarded. Moreover, the HVG Selection is a primary step that prioritized genes contributing the most to patient heterogeneity. Variance was calculated for each gene across all samples. The genes were ranked by variance, and the top N=2,000 Highly Variable Genes (HVGs) were selected. These HVGs serve as the input for the Weighted Gene Co-expression Network Analysis (WGCNA) in Phase 2, ensuring the network is built on drivers of molecular variation rather than housekeeping genes.

The heatmap in the figure 4 validates the selection process, revealing distinct gene modules. The ”checkerboard” pattern indicates that these genes effectively capture the molecular heterogeneity necessary to distinguish potential subtypes within the dataset. Following HVG selection, a validation process occurred in the data structure using linear and non-linear dimensionality reduction techniques to ensure the dataset contained discernable biological signals as.

The pairwise Pearson correlation matrix reveals low global correlation blue areas, confirming biological heterogeneity. However, distinct blocks of high correlation by red diagonal identify cliques of samples with similar transcriptional profiles, suggesting the presence of molecular subtypes.

PCA at Figure 6A reveals that the top two components explain $\sim$20% of the variance, indicating a complex, high-dimensional dataset. Both t-SNE and UMAP at Figure 6A & 6B show a largely continuous distribution of samples rather than discrete, separated clusters, though a small, distinct island of samples is visible in the non-linear projections, potentially representing a rare subtype.

The Silhouette analysis at figure 7B shows generally low scores ($<$0.08), consistent with the continuous nature of the data seen in PCA. While K-means K=3 partitions the PCA space as shown in figure 7C the boundaries are not strict. The Hierarchical Dendrogram as shown in figure 7A identifies roughly five major clades, suggesting that while the data is continuous, it possesses a hierarchical taxonomy of subtypes that can be exploited for the AML drug target discovery.

The objective of this analysis was to transform raw gene expression profiles into a systems-level understanding of gene–gene interactions in AML, thereby enabling the identification of AML-specific co-expression modules, biomarker candidates, and putative therapeutic targets. Co-expression network analysis was selected because it captures higher-order transcriptional organization and coordinated gene activity that cannot be detected through single-gene differential analysis alone.

A weighted gene co-expression network was constructed using the 2000 HVGs selected in Phase 1. First, a symmetric gene–gene Pearson correlation matrix was computed across 173 AML samples, capturing linear co-expression patterns:

• •

  High positive correlation indicates genes that exhibit coordinated behavior across AML patients, often reflecting shared pathway involvement.

• •

  Low or near-zero correlation suggests independent or uncoordinated transcriptional activity.

To enhance meaningful biological signals while suppressing noise, the correlation matrix was transformed into a weighted adjacency matrix using a soft-thresholding power function, as shown in Equation (1).

$$A_{ij}=\lvert r_{ij}\rvert^{\beta}$$

(1)

With $\beta=6$ selected to approximate a scale-free topology based on Fig. 8A & 8B. This value represents the inflection point that maximizes the scale-free topology fit ($R^{2}>0.85$) while maintaining sufficient mean connectivity to ensure the network remains biologically interconnected and not overly sparse. This step mirrors the core principle of WGCNA, ensuring that strong gene–gene relationships contribute disproportionately to module formation. Self-loops were removed by setting all diagonal entries to zero.

Each module was summarized by its corresponding module of eigengene, defined as the first principal component (PC1) of the module’s normalized expression matrix. For modules containing at least two genes, eigengenes were computed using PCA and served as representative expression signatures that capture dominant transcriptional variation within each module.

To interrogate the systemic relationships between the eight identified co-expression modules, the constructed Module Eigengene Network was implemented. For each module, a ”Module Eigengene” (ME) defined as the first principal component of the module’s expression matrix was calculated to represent its summary expression profile. Then computed the Pearson correlation matrix between all pairwise MEs. Intramodular hub genes were then identified by quantifying intramodular connectivity, where genes with the highest summed edge weights within a module were classified as hubs. Hub genes are particularly important because, with a total of 8 modules that shown at figure 9,

While the preceding co-expression network analysis considered the broader landscape of highly variable genes, this step focused on narrowing the gene set to a manageable subset of highly informative candidates. To achieve this, a ranking of WGCNA selection was employed, leveraging the assumption that genes exhibiting the largest variance across AML patients are most likely to reflect meaningful biological differences rather than stochastic noise. Hemoglobin-related genes were excluded from consideration, as their high expression is often non-specific and can obscure the detection of true AML-relevant signals. Following this filtering process, the top 20 genes were retained as candidate biomarkers, representing a balance between variability and biological relevance.

The top 20 prioritized biomarkers reveal a distinct targetable landscape characterized by the upregulation of myeloid-lineage surface receptors, metabolic transporters, and secreted inflammatory mediators. The analysis identified three primary functional clusters suitable for therapeutic intervention. First, Solute Carrier (SLC) Transporters such as SLC7A7, SLC11A1, SLC15A3 and the metabolic enzyme HK3 represent high-value targets for small-molecule inhibition. The upregulation of the cationic amino acid transporter SLC7A7 suggests a metabolic dependency on exogenous arginine uptake to fuel leukemic blast proliferation, a vulnerability that can be exploited by starvation strategies or high-affinity competitive inhibitors. Similarly, HK3 overexpression underscores the ’Warburg effect’ in AML, offering a defined enzymatic pocket for direct pharmacological blockade. Second, the immune-modulatory surface receptors SIGLEC9, LILRA5/6, and CD300E highlight mechanisms of immune evasion. SIGLEC9, an inhibitory checkpoint receptor, recruits phosphatases to dampen NK and T-cell cytotoxicity; blocking its ligand-binding interface could restore anti-leukemic immunity. These surface-accessible proteins are optimal candidates for both small-molecule modulation (targeting the sialic acid binding cleft) and antibody-drug conjugate (ADC) delivery. Finally, the secretome is dominated by the S100 Alarmin Family such as S100A8, S100A9, S100A12, and FGL2. These factors condition the bone marrow niche to support chemoresistance and suppress T-cell function. While traditionally targeted by neutralizing antibodies, the structural analysis suggests that their oligomerization interfaces and receptor-binding domains possess druggable hotspots suitable for novel protein-protein interaction (PPI) inhibitors generated by the proposed de novo pipeline.

Hierarchical clustering was then applied to the top 20 genes to assess co-expression patterns and potential subgroupings among patients. A clustered heatmap was generated, which visualized both gene-gene relationships and patient clustering simultaneously. This approach allowed the identification of co-expressed gene clusters, reflecting putative functional modules within this reduced biomarker set.

To elucidate the systemic interrelationships among the top 20 biomarkers, using a constructed undirected weighted co-expression network based on pairwise Pearson correlations. Only gene pairs exceeding a correlation threshold of 0.7 were connected, with edge weights reflecting co-expression strength. This network topology revealed hub genes with highly connected nodes likely to act as central drivers of the identified AML transcriptional program. Additionally, the network structure uncovered distinct modular subgroups, suggesting shared regulatory mechanisms or functional pathways suitable for coordinated therapeutic targeting. The heatmap, combined with hierarchical clustering, reveals sample-level patternsand identifies potential biomarker subgroups. Dendrograms, computed using Euclidean distance and average linkage, offer structural insight into how expression profiles cluster together, indicating whether biomarkers define coherent biological modules.

UniProt accessions were used to fetch predicted 3D structures from the AlphaFold3 Protein Structure Database [26]. HTTP requests were used to download each structure, and all successfully retrieved PDB files were stored in a structured directory for later inspection. VCAN biomarker logged and flagged for manual examination because it doesn’t have an available AlphaFold pdb strucutre. For those with existing renderings in the project directory, all images were assembled into a unified figure, ensuring that structural data could be directly visualized downstream.

Finally, all computed metrics values, composite scores, functional annotations, PCA coordinates, clustering results, correlation matrices, and structural availability information were merged into a single organized Excel workbook. This consolidation step ensures that researchers can navigate seamlessly from statistical ranking to biological interpretation and structural modeling, forming a complete biomarker discovery framework applicable to AML.

To identify and characterize druggable pockets within the target protein, we used the DOGSiteScorer engine provided through the ProteinsPlus API [27]. The workflow consists of four major stages: structure upload, job execution, pocket detection, and feature extraction/scoring. The protein target structure in PDB format was submitted to the DOGSite endpoint of the ProteinsPlus API. Upon submission, the service returned a unique job_id corresponding to a server-side processing task. This asynchronous workflow ensures that binding pocket detection is completed prior to downstream analysis.

DOGSite internally performs grid-based geometric and physicochemical analysis on the protein surface to identify concave regions that exhibit cavity-like properties. The algorithm examines:

• •

  3D surface topology.

• •

  Local curvature.

• •

  Enclosed volume regions.

• •

  Spatial arrangement of residues and atoms.

After detection, pocket-level physicochemical and geometric features are retrieved with each entry including the following quantitative descriptors, which are central for ligandability assessment.

These features collectively represent the binding potential of the pocket as shown in table IV. To prioritize the most promising pockets, a composite scoring function was applied:

	Score	$\displaystyle=$	$\displaystyle 0.3\times\text{Volume}+0.2\times\text{Depth}+0.2\times\text{Enclosure}\times 100$		(2)
			$\displaystyle+0.1\times\text{Hydrophobicity}\times 100$
			$\displaystyle+0.1\times\text{Aromaticity}+0.1\times(\text{Donors}+\text{Acceptors})$

This scoring formula balances geometric enclosure (important for ligand binding), cavity size, and chemical interaction potential. After scoring:

• •

  Scores were normalized to a 0-1 scale.

• •

  Pockets were sorted in descending order.

• •

  Enclosed volume regions.

• •

  The top-scoring pocket was selected as the most druggable region for downstream docking and molecular modeling (taking the top 3 pockets).

The structural analysis of the biomarker A0AV96 revealed multiple potential ligand-binding sites distributed across its protein surface. Using the DOGSiteScorer engine, the 3D structure of A0AV96 was systematically scanned to identify surface cavities that exhibit favorable geometric and physicochemical characteristics for ligand interaction. The algorithm detected six distinct binding pockets, each varying in depth, enclosure, hydrophobicity, and aromatic contribution. These pockets represent the most probable druggable and functionally relevant regions on the biomarker. Visual inspection of the pocket architecture shows a heterogeneous distribution of cavity sizes, including both deep, well-enclosed pockets and more solvent-exposed shallow regions, suggesting the biomarker may support diverse modes of ligand engagement.

Fragment Spatial Alignment Using Rodrigues 3D Rotation.

The generation process utilizes a genetic algorithm with a population size of $N=60$ over $G=40$ generations. The algorithm follows a reaction-first strategy, prioritizing chemical validity over random connectivity. The initial population is generated by randomly selecting fragments and spatially aligning them to the identified binding hotspots. Rodrigues’ rotation formula is employed to align the principal axis of a fragment vector $\vec{v}_{\mathrm{frag}}$ with the corresponding target hotspot vector $\vec{v}_{\mathrm{target}}$:

$$R=I+\left(\sin{\theta}\right)K+\left(1-\cos{\theta}\right)K^{2}$$

(3)

Where K is the cross-product matrix of the unit vectors. This matrix ensures that fragments are correctly oriented toward the interaction sites before linkage. For each placed fragment, the following computations were performed: its 3D center of mass, its principal axis (defined as the farthest-atom vector), and the target direction derived from the hotspot point.

Genetic Operators Offspring ligands are generated through probabilistic application of three operators:

1. 1.

   Crossover (Recombination): Two parent molecules are combined using either:

   • •

     Reaction Linking: Application of SMARTS-based reaction templates to ensure synthetic feasibility.

   • •

     BRICS Crossover: Decomposition of parents into BRICS (Biosynthetic Responsible Interaction of Chemical Scaffolds) fragments, followed by shuffling and re-merging.

2. 2.

   Mutation:

   • •

     Fragment Injection: Addition of a new fragment from the library linked via reaction rules.

   • •

     Conformer Relaxation: Stochastic UFF energy minimization to relieve internal strain.

   • •

     Scaffold Hopping: Replacement of a structural subunit with a bioisostere from the library.

3. 3.

   Safety Filters: All operations utilize a safe merge protocol, which validates valence constraints and atomic connectivity before accepting a new structure.

To ensure the generation of high-quality drug candidates, employed a composite fitness function F(m) that balances physicochemical properties, geometric fit, novelty, and conformational energy.

$$F\left(m\right)=w_{p}\cdot S_{proxy}+w_{f}\cdot S_{fit}+w_{n}\cdot S_{novelty}+w_{s}\cdot S_{strain}-\lambda\cdot P_{SA}$$

(4)

Where the weights were empirically set to: $w_{p}=0.45$, $w_{f}=0.35$, $w_{n}=0.15$, and $w_{s}=0.2$

Proxy Score $S_{proxy}$: A weighted average of QED (Quantitative Estimate of Drug-likeness), LogP (Lipophilicity), Molecular Weight normalization, and Rotatable Bonds.

Pocket Fit $S_{fit}$ is calculated using:

$$S_{fit}=0.6\cdot F_{inside}+0.4\cdot\left(1-\frac{{\overline{d}}_{hot}}{{\overline{d}}_{hot}+4.0}\right)$$

(5)

Where $F_{inside}$ the fraction of heavy atoms within the pocket radius, and $bard_{hot}$ is the mean distance to the nearest hotspot.

• •

  Novelty $S_{novelty}$: Calculated as $1-\max{\left(T_{c}\right)}$ where $T_{c}$ Tanimoto similarity to known reference ligands using Morgan Fingerprints of 2048 bits.

• •

  Strain Energy $S_{strain}$: A penalty derived from the per-atom UFF energy after constrained embedding, ensuring molecules do not adopt high-energy, unrealistic conformations; scaled 0–1.

• •

  SA Penalty $P_{SA}$: A synthetic accessibility penalty derived from fragment complexity and ring topology.

The generative process operates through evolutionary iterations. In each generation, numerous fragment combinations are sampled, spatially guided, chemically merged, and evaluated. Only the highest-fitness molecules are retained for subsequent generations, enabling progressive refinement toward ligands with strong drug-likeness, AML relevance, optimal pocket compatibility, and stable 3D conformations. After the predefined number of generations, top-scoring candidates from all cycles are aggregated and reported as the final set of de novo AML-conditioned ligands.

The selected ligands demonstrated an impressive balance between high binding affinity and desirable pharmaceutical attributes, supporting the efficacy of the multi-objective optimization process.

• •

  Pocket Fit Score: TAll ten prioritized ligands achieved high Pocket Fit Scores, ranging from 0.706 to 0.725. Ligand L4 exhibited the highest score of 0.725, while L1 and L7 also demonstrated strong predicted binding with scores of 0.706 and 0.711, respectively.

• •

  Drug-likeness (QED and Proxy Score): The Quantitative Estimate of Drug-likeness (QED) scores for the top candidates spanned from 0.48 to 0.62. Ligand L8 and L9 showed the highest QED values of 0.62 and 0.61, respectively, indicating a favorable balance of desirable physicochemical properties for drug development.

• •

  Synthetic Feasibility (SA Score): The synthetic accessibility of the ligands was assessed using the SA Score, where lower values indicate easier synthesis. The scores were excellent, ranging from 0.83 to 2.01. Ligand L5 achieved the lowest SA Score of 0.83, classifying it as very easy to synthesize.

• •

  Lipophilicity (LogP): The calculated partition coefficient (LogP) for the ligands varied from 0.56 to 2.99, placing them well within the range considered suitable for oral bioavailability and membrane permeability.

The set of top ligands represents diverse chemical scaffolds, supporting the claim of successful exploration of the chemical space.

• •

  Pocket Fit Score: TAll ten prioritized ligands achieved high Pocket Fit Scores, ranging from 0.706 to 0.725. Ligand L4 exhibited the highest score of 0.725, while L1 and L7 also demonstrated strong predicted binding with scores of 0.706 and 0.711, respectively.

• •

  Aromatics and Small Molecules: Ligands L4 and L6 are relatively small, simple halogenated aromatic compounds, while L9 is a simple N-methylbenzamide.

• •

  Heterocycles: Structures L1 and L7 contain pyrimidine-aniline scaffolds, known for their privileged structures in kinase inhibition

• •

  Aliphatic Structures: Ligands such as L5 and L8 incorporate saturated rings and flexible aliphatic chains, respectively, contributing to the overall structural breadth of the generated library.

The computational predictions from ADMET Lab 3.0 provide an initial, high-throughput screen for potential liabilities and favorable pharmacokinetic properties.

Histograms at figure 24; showing the distribution of the MWs, the octanol-water Partition Coefficient $\log P$ and TPSA for the compound set. Vertical lines denote established thresholds for favorable oral bioavailability. The distribution of MW shows where most compounds lie relative to the Lipinski Rule of Five (Ro5) limit of 500 Da. Compounds with MW greater than 500 Da may face challenges with passive permeability and oral absorption. The plot assesses the lipophilicity of the set. Compounds are ideally balanced, with $\log P$ typically between 1 and 5. Highly lipophilic compounds log(p) ¿ 5 often lead to poor solubility and non-specific binding, while highly polar ones $\log P$ ¡ 0 struggle with membrane permeability. The histogram indicates the central tendency and spreads of lipophilicity. Finally, the TPSA relates to the size and polarity of a molecule, directly impacting passive transport. Most compounds should fall below the oral limit $140\,\text{\AA }^{2}$. For CNS-targeting drugs, an even stricter limit is often applied $<90\,\text{\AA }^{2}$. The distribution here reveals how many candidates meet the criteria for general vs. CNS penetration.

At figure 25; Bar chart showing the number of molecules out of N=10 predicted to be inhibitors probability $>0.5$ of the three most clinically relevant cytochrome P450 (CYP) isoforms. Drug-Drug Interaction (DDI) Risk: Inhibition of CYP enzymes, particularly the highly prevalent CYP3A4, is a major cause of clinical failure due to DDI risk. While in the Lead Prioritization; Candidates showing high inhibition counts for any major isoform, especially CYP3A4, should be flagged as high-risk and are generally de-prioritized or require significant structural modification to remove the metabolic liability. Compounds with zero predicted inhibition across all isoforms are preferred for minimizing DDI risk.

The overall in silico ADMET and physicochemical profiling of the ligand candidates revealed an exceptionally promising drug-like profile. All compounds in the set successfully adhere to major drug-likeness criteria, recording zero violations across the Lipinski, Pfizer, GSK, and Golden Triangle rules, suggesting strong prospects for good oral bioavailability and membrane permeability. Physicochemically, the set’s distribution across MW, log(P), and TPSA confirms that the molecules reside within the favorable chemical space for oral drug candidates. Furthermore, the library exhibits minimal predicted liabilities concerning xenobiotic metabolism, with zero predicted inhibition probability ¿ 0.5 reported for the three major Cytochrome P450 isoforms CYP1A2, CYP2C9, and CYP3A4, which significantly minimizes the risk of adverse DDIs. The property map provides a clear visualization of the compounds’ log (P) versus TPSA space, which, when coupled with the predicted Blood-Brain Barrier (BBB) scores, enables efficient prioritization of leads for either CNS or non-CNS indications. Collectively, these results indicate that the ligand set possesses excellent pharmacokinetic and safety attributes, positioning them as high-quality starting points for further development.

Molecular docking was conducted using two independent platforms: the widely-used SwissDock server and the internal CB.DOCK2 [30] ; both are using AutoDock Vina docking environment, to ensure the robustness and convergence of the predicted binding pose as well as auto blind docking. The active site was defined based on conserved catalytic residues identified from sequence homology and crystallographic data and predicted based on the top consensus binding pocket. For both methods, the lowest-energy binding mode was selected for detailed structural analysis. To ensure a comprehensive search of the conformational space and to avoid entrapment in local energy minima, a high sampling exhaustivity of 35 was applied. This rigorous sampling protocol allowed for an intensive exploration of the ligand’s degrees of freedom within the protein’s binding architecture, ensuring that the identified global minimum is robust and reproducible.

The protein is displayed in a gray solvent-accessible surface or ribbon view, highlighting the deep embedding of the ligand within the catalytic cleft. This orientation corresponds to Model 1, representing the global minimum energy. Moreover, figure 26 demonstrates the high degree of geometric complementarity between the ligand scaffold and the protein’s interior contours. The distribution of polar and non-polar regions suggests an optimized fit that accounts for the stable binding energies observed across the top-ranked docking clusters.

Docking Results and Thermodynamic Stability

The docking simulations successfully identified several high-affinity clusters. The use of blind docking via CB-DOCK2 allowed for the identification of the most probable binding cavities, while SwissDock provided detailed energetic scoring for the ligand poses. The top-ranked models, their calculated affinities, and the specific grid sizes ;docking dimensions are summarized below:

The calculated binding free energy $\Delta$ G for the top-ranked pose is -6.571 kcal/mol by SwissDock and -6.800 kcal/mol using CB-Dock2. The energy magnitude has a values in the range of -6.0 to -7.0 kcal/mol are indicative of spontaneous binding with moderate-to-high affinity, typically corresponding to an inhibition constant of $K_{i}$ in the low micromolar range. Moreover, the boltzmann distribution with the narrow energy gap between Model 1 and Model 10 a difference of only $~0.4kcal/mol$ suggests a ”flat” energy landscape where the ligand can adopt several slightly different but closely related orientations, all of which are energetically favourable. The docking results provide strong computational evidence that Ligand L1 is a viable candidate for the A08A96 biomarker which is one of the candidates; and emphasizing on the metaheuristic model to generate highly promising ligands. The consensus in binding energy, the high sampling exhaustivity, and the favorable thermodynamic scores confirm that L1 possesses the structural requirements to form a stable inhibitory complex with the target protein.

Without experimental testing, critical biological aspects such as cellular uptake, metabolic stability, toxicity, immunogenicity, and off-target effects remain unknown. Molecular docking and ADMET predictions provide useful early insights, but their reliance on simplified force fields, static protein models, and approximations of solvation and dynamics may lead to false positives or overestimated binding affinities. Additionally, dependence on AlphaFold3-predicted structures introduces uncertainty, as even high-confidence models may deviate from native conformations, particularly in flexible regions crucial for ligand binding.

Another limitation stems from using bulk RNA-seq data from the TCGA-LAML cohort. Bulk transcriptomics averages signals across heterogeneous cell populations, potentially masking the molecular profiles of rare leukemic stem cells and therapy-resistant subclones that drive relapses. As a result, target selection may be biased toward dominant clones while overlooking clinically relevant vulnerabilities in minor subpopulations.

Future work will prioritize experimental validation and iterative optimization of top-ranked ligands. Selected candidates will be synthesized and tested in vitro for binding and potency, followed by in vivo evaluation in patient-derived xenograft models to assess efficacy, pharmacokinetics, pharmacodynamics, and preliminary toxicity. These results will feed into a closed-loop Design-Make-Test-Analyze cycle, enabling continuous refinement of the computational model and medicinal chemistry optimization.

To address tumor heterogeneity, future iterations will integrate single-cell RNA sequencing and spatial transcriptomic data. This will allow identification of molecular vulnerabilities in leukemic stem cells and resistant subclones, support subclone-specific targeting and guide the design of combination or multi-target therapies. Incorporating profiles from non-malignant bone marrow cells, such as immune and stromal populations, will enable microenvironment-aware drug design by highlighting targets involved in pro-leukemic interactions and immunomodulatory pathways.

The pipeline’s modular, disease-agnostic design allows its extension beyond AML to other hematologic malignancies, solid tumors with high heterogeneity, and diseases with limited treatment options. Integrating additional data modalities, including proteomics, epigenomics, and metabolomics, will provide a more comprehensive representation of disease biology and reveal previously inaccessible druggable targets. Collectively, these features position the pipeline as a flexible, generalizable computational platform for precision drug discovery across oncology and other therapeutic areas with unmet clinical needs.