To delineate differentially expressed genes (DEGs), RNASeq gene expression data of BRCA were retrieved from the TCGA database. Differential gene expression analysis was carried out using the R package “DESeq2”, and the Wald test was applied for significance testing. Then, the genes are ranked by adjusted p-values based on the Benjamini-Hochberg algorithm. Genes with log2 FC > 1 and adj. p < 0.05 were considered upregulated, while those with log2 FC < −1 and adj. p < 0.05 were considered downregulated genes. The functional annotation of the identified DEGs was conducted using the g:Profiler, a web-based interface (Cui et al., 2024; Zhao et al., 2022). Finally, visualization of DEGs was performed using the “ggplot2” package and presented as a volcano plot. This DEG profiling workflow and code were obtained from the methodology described by Krishnamoorthy and Karuppasamy (2024). The preprocessing pipeline for TNBC DEGs is available at the following GitHub repository (https://github.com/vitmilab/lab-projects/tree/main/DEG_Identification). The workflow of the study is depicted in Figure 1.

The experimentally verified targets of miR-30a-5p were obtained from databases such as miRTarBase, miRDB, and TargetScan. miRTarBase (https://ngdc.cncb.ac.cn/databasecommons/database/id/167) focuses exclusively on experimentally validated miRNA-target interactions supported by high-throughput evidence (Chen and Wang, 2020; Zhao et al., 2024). miRDB (https://mirdb.org/) employs the MirTarget computational tool to analyze miRNA-target interactions and their regulatory effects based on high-throughput sequencing (HTS) data (McGeary et al., 2019). TargetScan (https://www.targetscan.org/vert_80/) identifies potential miRNA targets by recognizing conserved motifs aligning the miRNA seed regions (Raudvere et al., 2019).

The common targets between “upregulated DEGs vs. miR-30a-5p targets” and “downregulated DEGs vs. miR-30a-5p targets” were identified through “Venny 2.1 (https://bioinfogp.cnb.csic.es/tools/venny/) by entering the miR-30a-5p targets and disease genes in the query box (Oliveros, 2007).

Protein-protein interactions (PPIs) among identified DEGs were assessed using the STRING database v12.0 (https://string-db.org/) (Scklarczyk et al., 2023; von Mering et al., 2004). Predicted interactions for the target genes were extracted using a medium confidence score ≥0.4, as calculated by STRING’s probabilistic scoring system (Equation 1). The confidence score shows the estimated likelihood that a predicted interaction reflects a true biological association, derived from combined evidence including experiments, co-expression, curated databases, and text mining. A score of ≥0.4 was selected to balance network connectivity and specificity, as using stricter cutoffs of ≥0.7 shattered the PPI network and removed valid interactions. S = 1 − ∏ i 1 − S i (1)

Protein-interaction network (PIN) from STRING was imported and visualized through Cytoscape 3.10.3 (https://cytoscape.org/). Topological analysis was performed using the cytoHubba plugin to identify the top 10 hub genes, through six centrality measures, including maximal centrality clique (MCC), maximum neighbour component (MNC), degree, closeness, betweenness, and stress centralities (Equations 2–7) (Chin et al., 2014). For each gene, the ranks from all six algorithms were computed to obtain the mean rank, and genes with the highest composite scores were defined as hub genes. M C C   v = ∑ C ∈ S v C − 1 ! (2)

Here, the maximal clique’s collection containing v is represented by S v , while the product of all positive integers less than C is denoted by and C − 1 . Absence of an edge between the neighbour node v , then M C C   v equals degree. M N C   v = V ( M C V (3)

In this case, the maximum connected component of the G N V is denoted by M C V and G N V is the induced subgraph of G by N V . Deg   v = N V (4)

Here, the node is represented by V and the collections of its neighbours are denoted by N V . C l o   v = ∑ w ϵ v 1 d i s t   v , w (5) B C   v = ∑ s ≠ t ≠ v ϵ c v σ s t   v σ s t   (6)

In this case, the number of shortest paths between nodes s and t is represented by σ s t . S t r   v = ∑ s ≠ t ≠ v ϵ c v σ s t v (7)

Here, the number of shortest paths from the node s to t is represented by σ s t v , which uses the node v .

Gene clusters that are tightly linked within the PIN were extracted from the Molecular Complex Detection (MCODE) plugin. Core targets identified by cytoHubba were cross-referenced with the MCODE clusters with scores >3, using default parameters, to further refine the selection of biologically significant genes (Shannon et al., 2003).

The biological significance and pathways associated with the key genes were determined using ShinyGO v0.80 (https://bioinformatics.sdstate.edu/go/). The significance was assessed using Benjamini-Hochberg false discovery rate (FDR) of 0.05 (Ge et al., 2019). Subsequently, the top 10 GO pathways involved in biological process (BP), molecular function (MF), and cellular component (CC), along with Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, were selected. Fold enrichment of each term was obtained directly from the ShinyGO enrichment output. Fold enrichment shows the overrepresented genes in a particular pathway and is calculated as observed/expected overlap k / E , where k is the number of submitted genes found in the GO/KEGG term, and E is the expected overlap, which was calculated using the formula: E = K N   x   n , where n is the number of genes submitted (10 genes), K is the number of genes annotated to that term, and N is the background universe automatically selected by the enrichment tool.

The expression levels between TNBC vs. normal breast tissues were assessed using the Gene Expression Profiling Interactive Analysis 2 (GEPIA2) server (http://gepia2.cancer-pku.cn/#index) (Tang et al., 2019), with log2 FC ≥ 1.0 and p ≤ 0.05. Furthermore, through Kaplan-Meier (KM) plotter (https://kmplot.com/analysis/), the overall survival (OS) analysis in TNBC patients was evaluated.

The mutation frequency was analyzed through the Gene Set Cancer Analysis (GSCA) database (https://guolab.wchscu.cn/GSCA/#/). This database provides a comprehensive depot of genomic alterations, including single-nucleotide variants (SNVs), in a variety of cancer forms, including breast invasive carcinoma (BRCA). SNV data corresponding to the hub genes were extracted and evaluated to determine the mutation frequency and variant types observed in BRCA samples.

The target proteins’ 3D structure was obtained from RCSB Protein Data Bank (https://www.rcsb.org/) in PDB format. AutoDock Tools (ADT) 1.5.7 was used for protein preparation, where water molecules, ligands, and heteroatoms were removed. Subsequently, polar and Kollman charges were added, and the processed structures were saved in PDBQT format. The ligand (EGCG) structure was downloaded from PubChem (CID: 65,064) in SDF format and converted to 3D using Open Babel 3.0.a1, and PDBQT files were generated using AutoDock Tools. A grid box was defined for each target protein to encompass the known active/binding site. Molecular docking was performed using AutoDock Vina v1.2.x. The output PDBQT files were converted and visualized using PyMOL to analyse interactions and binding conformations. The binding affinities (kcal/mol) of the top-ranked poses were recorded. Binding energy represents the predicted free energy of ligand-protein binding Δ G b i n d i n g , where more negative values indicate stronger and more favorable interactions. Using the PLIP server, the interactions (hydrogen bonding, π–π stacking, and hydrophobic interactions) were evaluated.

GROMACS 2020.3 package with the CHARMM36 force field was used for MD simulations to ascertain the structural stability and integrity of the four gene-ligand complexes (Van Der Spoel et al., 2005). Ligand parameter and topology files were generated through the CHARMM General Force Field (CGenF). Each protein-ligand complex was placed in a dodecahedron box using a simple point charge (SPC) water model, and the system’s neutrality was ensured by introducing three chlorine counterions. The system’s energy minimization was carried out using the steepest descent algorithm to eliminate unfavourable weak Van der Waals contacts. The algorithms, like the Particle Mesh Ewald (PME) method and the linear constraint solver (LINCS), were applied to constrain the electrostatic and hydrogen bond interactions. The system was equilibrated using canonical NVT (number of particles, volume, and temperature) and isobaric NPT (number of particles, pressure, and temperature) ensembles to stabilize its thermodynamic conditions. Berendsen thermostat used to heat the system to 300 K gradually, with a time lapse of 0.1 ps and a pressure of 1 bar. Then, the MD simulation was carried out for 100 ns with a 2-fs integration time step, and the trajectory data (RMSD, RMSF, Rg, and SASA) were analyzed based on the Formula (8-10). SASA is calculated from the atomic coordinates of a native protein numerically (Ghahremanian et al., 2022). R M S D = 1 N ∑ i = 1 N x 1 m − x 1 1 2 + y 1 m − y 1 1 2 + z 1 m − z 1 1 2   (8)

Where x m , y m , z m are the initial coordinates and x 1 , y 1 , z 1 are the trajectory coordinates at frame t, and N is the number of atoms. R M S F = 1 T ∑ i = 1 T x i − x ¯ 2 (9)

Where N is the trajectory frame numbers and x ¯ is the time-averaged position. R g = 1 N ∑ i = 1 N r i − r c e n t e r 2 (10)

Where r i is the coordinates of the atom i , and r c e n t e r is the center of mass, and N is the number of protein atoms.

TNBC-specific transcriptomic data were obtained from the TCGA-BRCA dataset. The dataset was preprocessed, and the DEGs were identified using the “DESeq2 package” in R, ensuing an initial finding of 43,881 DEGs. Subsequently, the data were filtered based on log2 FC and adj. p value, resulting in a 5008 upregulated and 2,905 downregulated refined gene list within the TNBC dataset.

Functional annotation of these DEGs was performed using g:Profiler to traverse through their biological importance. Expression of DEGs was visually interpreted through a volcano plot (Figure 2a), markedly showing upregulated genes in red and downregulated genes in blue, giving a clear depiction of the expression landscape.

An integrated target gene set of 2,115 genes was compiled through the integrative curation of 2,704 hsa-miR-30a-5p-associated genes, sourced from three well-established miRNA target prediction databases: miRTarBase (734 genes), TargetScan (431 genes), and miRDB (1,539 genes) (Figure 2b). Redundancy was systematically eliminated to ensure maximal data fidelity in the gene set.

The convergence between miR-30a-5p target genes and TNBC-related genes was examined using InteractiVenn (Figures 2c,d). This integrative analysis revealed 393 overlapping genes, comprising 219 upregulated and 174 downregulated genes. All these genes are differentially expressed in TNBC and predicted to contain validated or putative miR-30a-5p binding. These represent functionally relevant targets of miR-30a-5p within the TNBC.

The PPI analysis of the DEGs from the STRING database with a confidence score of ≥0.4 generated 392 nodes and 622 edges. The network exhibited an average node degree of 3.17 and an average local clustering coefficient of 0.351. Enrichment analysis yielded a highly significant PPI enrichment p-value of 2.49e-13 (p < 0.05), indicating that the observed interactions are unlikely to be random and reflect biologically meaningful associations.

The resulting PIN from STRING was imported into Cytoscape for further topological assessments (Figure 3). The network topology showed a clustering coefficient of 0.207 and a characteristic path length of 4.187, indicating that the network holds sturdy connections, toughness, and flexibility.

CytoHubba plugin comprised the top 10 hub genes involved in TNBC survival and progression based on the ranking obtained from the six topological parameters. Furthermore, the MCODE plugin detected significant gene clusters within the PPI network. The clusters identified through MCODE were cross-referenced with the results from cytoHubba. Notably, RRM2, KIF11, ANLN, CDC20, CCNA1, AGO2, YWHAZ, DTL, SKP2, and PCNA consistently appeared across multiple cytoHubba rankings (Table 1) and were predominantly located within the MCODE clusters with a score >3 (Figure 4).

GO enrichment analysis was conducted using ShinyGO v0.80 to gain deeper insights into the biological significance of the ten key target genes regulated by miR-30a-5p in TNBC (Table 2). The GO BP analysis revealed significant enrichment of 437 pathways, with top-ranking terms including regulation of synapse maturation, mitotic cell cycle phase transition, regulation, processes, cell cycle, and its regulation. In the CC category, 67 enriched pathways were identified, with the hub genes predominantly localized to functionally critical structures and protein complexes. These included numerous complexes such as cyclin-dependent protein kinase (CDK) holoenzyme, cullin-RING ubiquitin ligase, transferase, intracellular protein-containing complex, microtubule cytoskeleton, nucleoplasm, and nuclear lumen. Within the MF, 80 significantly enriched terms were observed. The most prominent among these were purine-specific mismatch base pair DNA N-glycosylase activity, ribonucleotide-diphosphate reductase activity thioredoxin disulfide as acceptor, DNA polymerase processivity factor activity, ribonucleotide-diphosphate reductase activity, and enzyme activity associated with nucleotide metabolism and DNA replication. Primarily, the GO pathway enrichment analysis indicated that the genes are linked with cell cycle regulation, emphasizing the importance of hub genes in cell cycle regulation.

From the KEGG analysis, out of 31 distinct pathways, the most significantly enriched ten pathways are detailed in Table 3. Hub genes are prominently involved in the cell cycle pathway, followed by viral carcinogenesis, mismatch repair, hepatitis B, and oocyte meiosis pathways. The dominant involvement of these genes in cell cycle-related pathways strongly suggests that miR-30a-5p may exert its therapeutic potential in TNBC by modulating core signaling cascades. The fold enrichment and total number of genes were reported in Tables 2, 3 for Gene Ontology and KEGG pathways, respectively.

The expression analysis of the ten hub genes was analyzed using GEPIA2 to assess their expression levels between normal and tumor samples. Expression levels were represented as -log₂(TPM+1) to normalize the data. Among these, seven genes, namely, RRM2, KIF11, ANLN, CDC20, YWHAZ, DTL, and PCNA, exhibited significantly elevated expression in tumor samples (p ≤ 0.05). Whereas CCNA1, AGO2, and SKP2 did not show statistically significant differences, indicating relatively stable expression between normal and tumor tissues (Figure 5).

The prognostic significance of seven key upregulated hub genes, namely, RRM2, KIF11, ANLN, CDC20, YWHAZ, DTL, and PCNA, was evaluated for overall survival in TNBC patients using Kaplan-Meier survival analysis. From the above genes, CDC20 exhibited the strongest prognostic effect with a hazard ratio of 1.71 (confidence interval (CI): 1.41–2.07; p = 2.7e-08). This was followed by RRM2 (HR = 1.57, CI: 1.3–1.89, p = 3e-06), YWHAZ (HR = 1.53, p = 9.6e-06), ANLN (HR = 1.53, p = 0.0018), and DTL (HR = 1.52, p = 1.2e-05). Additionally, PCNA (HR = 1.31, p = 0.0046) and KIF11 (HR = 1.26, p = 0.014) also exhibited significant associations with reduced survival probability (Figure 6).

Mutation analysis revealed that KIF11 and ANLN had increased mutations (20% each), followed by CDC20 (17%) and YWHAZ (14%). Missense mutation was found to be the abundant variant classification, with SNP and C>T as the variant type and SNV category, respectively. CDC20 and YWHAZ carried missense mutations, whereas ANLN contained both missense and nonsense mutations. In contrast, frame-shift mutations were observed in KIF11, alongside missense and nonsense mutations (Figures 7a,b). Investigation of 29 BRCA samples showed that 100% of the samples displayed changes in at least one prognostic gene, with 28% of the cases in KIF11 (Figure 7c). Moreover, comparable transitions and transversions with C>T and C>G being the most occurring conversions (Figure 7d). The increased burden of mutation diversity among the four prognostic genes suggests their role in disease progression and therapeutic targeting.

EGCG exhibited the strongest affinity for YWHAZ (−10.1 kcal/mol), forming eight hydrogen (H₂) bonds in the positions of Asp 191, Trp 234, Ser 275, Trp 276, Arg 316, Trp 363, Ser 404, and Ser 448. The next most favorable interaction was with ANLN (−8.9 kcal/mol), forming seven H₂ bonds, followed by KIF11 (−8.5 kcal/mol) stabilized by eleven polar contacts. CDC20 showed the weakest binding of −7.0 kcal/mol with only four H₂ bonds. Overall, the binding hierarchy was YWHAZ > ANLN > KIF11> CDC20, suggesting that YWHAZ have the most stable and complementary binding pocket for EGCG (Figure 8). Table 4 gives a summary of docking results with the number of hydrogen bonds.

The RMSD parameter (Figure 9a) examines the structural stability across all four protein-EGCG docked complexes. Notably, YWHAZ displayed good structural stability with slight deviation (⁓0.12 ± 0.02 nm), suggesting that the native scaffold was preserved by EGCG binding without persuading extensive conformational alterations. ANLN, in contrast, showed the highest fluctuations (⁓0.34 ± 0.03 nm), while CDC20 (⁓0.20 ± 0.03 nm) and KIF11 (⁓0.24 ± 0.03 nm) occupied intermediate positions. From the RMSD profile, it was observed that YWHAZ has a structurally rigid environment for EGCG binding, whereas CDC20 and KIF11 hold on to substantial backbone flexibility, leading to transitory shifts.

Consistent with the RMSD hierarchy, the solvent-accessible surface area (SASA) values (Figure 9b) revealed YWHAZ as the most compact complex (⁓95 ± 4 nm²), followed by ANLN (⁓120 ± 5 nm²), while CDC20 (⁓170 ± 5 nm²) and KIF11 (⁓205 ± 4 nm²) displayed substantially larger solvent exposure. The lower the SASA, the greater the compactness. YWHAZ shows a conformationally restrained state having low SASA, aligning with RMSD data. But CDC20 and KIF11 display greater flexibility and more solvent-exposed binding surfaces due to their increased SASA.

The Radius of Gyration (Rg) remained stable after rapid equilibrium within the first few ns (Figure 9c). Time-averaged Rg values strengthened the compactness hierarchy: YWHAZ (1.64 ± 0.02 nm) and ANLN (1.65 ± 0.02 nm) maintained tightly packed structures, whereas CDC20 (1.91 ± 0.02 nm) and KIF11 (1.93 ± 0.02 nm) were more extended. Low RMSD, SASA, and rigid Rg for YWHAZ highlight the increased stability, while it was vice-versa for CDC20 and KIF11 with their dynamic scaffolds.

RMSF demonstrated the structural grading of YWHAZ with the slightest mobility (baseline (⁓0.03–0.06 nm, peaks ⁓0.10–0.12 nm) (Figure 10d), reflecting a uniformly firm scaffold. ANLN showed intermediate flexibility (⁓0.12–0.18 nm at termini and loop regions) (Figure 10c), while CDC20 and KIF11 displayed distinct peaks (⁓0.20–0.27 nm) (Figures 10a,b), consistent with higher global flexibility.