Prognostic Significance of R-Loops in Lung Adenocarcinoma: Implications for Immune Response and Drug Sensitivity

Introduction

LUAD is the most common histological subtype of lung cancer and remains a leading cause of cancer-related mortality worldwide.^1,2 Despite advances in diagnosis and treatment, the prognosis of LUAD patients remains poor, largely due to the high heterogeneity of the disease and its complex tumor microenvironment.^3,4 Consequently, there is an urgent need to identify novel biomarkers and therapeutic targets to improve prognostic predictions and optimize personalized treatment strategies.

R-loops are triple-stranded structures formed during transcription, consisting of an RNA-DNA hybrid and a displaced single-stranded DNA region.⁵ These structures are known to play critical roles in maintaining genomic stability, regulating transcription, and influencing chromatin dynamics.^6,7 While these canonical roles are well established, recent findings suggest that R-loops may also contribute more broadly to cancer progression. Increasing evidence indicates that R-loops influence not only genome integrity but also the tumor microenvironment. For instance, R-loop–derived cytoplasmic RNA–DNA hybrids can activate the cGAS–STING pathway and type I interferon signaling, affecting antigen presentation and thus modulating tumor immunity.^8,9 Single-cell transcriptomic studies have revealed that R-loop patterns are associated with metabolic reprogramming and T cell exhaustion in cancer, characterized by elevated glycolytic activity and mitochondrial dysfunction.¹⁰ These analyses also highlight substantial heterogeneity in R-loop regulatory gene expression across immune and metabolic compartments within the tumor microenvironment, underscoring their potential as functional biomarkers and therapeutic targets.

R-loops are commonly detected using high-throughput techniques such as DNA-RNA immunoprecipitation sequencing (DRIP-seq) and R-ChIP, which enable genome-wide profiling of their distribution and abundance. Several genes involved in R-loop formation or resolution—such as DDX5, EIF4A3, and RNASEH2A—have been implicated in transcriptional regulation, DNA damage response, and tumor progression; however, their specific roles in LUAD remain largely undefined.^11–13 While aberrant R-loop accumulation has been associated with various diseases,¹⁴ both their mechanistic roles in LUAD and their value as prognostic or predictive biomarkers have not been systematically explored.

Recent advances in high-throughput sequencing technologies and machine learning have provided new approaches for investigating cancer mechanisms.^15,16 Zhang et al utilized single-cell technologies to investigate the connections between R-loops and immune evasion and metabolic reprogramming, highlighting their impact on cancer progression.¹⁰ Although a previous study has explored the relationship between R-loops and the prognosis of LUAD,¹⁷ it was likely limited by a small sample size and lacked a comprehensive analysis of pathway enrichment, immune characteristics, and drug sensitivity associated with R-loops. Therefore, there is an urgent need for a more detailed investigation into a prognostic model based on R-loops.

In this study, we conducted an integrated analysis of R-loops genes in LUAD to identify key regulatory genes. By developing a novel R-loops prognostic risk score model (R-loops Score), we aimed to assess its predictive value for overall survival and explore its associations with immune characteristics and drug sensitivity. Additionally, we investigated the underlying biological pathways and clinical implications of the R-loops Score. To further validate our findings, qPCR and Western blot analyses consistently confirmed EIF3B as a key gene, supporting its potential as a prognostic biomarker. This comprehensive analysis seeks to provide new insights into the role of R-loops in LUAD progression and offers potential guidance for personalized treatment strategies.

Methods

Data Source

We collected 803, 612, and 469 R-loop-associated proteins from three independent studies.^18–20 From the R-loopBase database (https://rloopbase.nju.edu.cn/index.jsp), we downloaded the file “gene_list_of_R-loop_regulators.gz” and extracted 1185 human R-loop regulators. These datasets were integrated, and all protein or gene identifiers were standardized to gene symbols using R. After removing duplicates, we compiled a comprehensive, non-redundant list of 1,771 R-loop-associated genes.

Gene expression data and clinical information for TCGA-LUAD were obtained from the UCSC Xena database (https://xenabrowser.net/datapages/). Additionally, we retrieved expression and clinical data for four Affymetrix microarray cohorts (GSE50081, GSE31210, and GSE29013) from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/).^21–24

For TCGA-LUAD, gene expression data in FPKM format were downloaded and converted to TPM format using R.²⁵ For the Affymetrix microarray datasets, we downloaded the raw CEL files for each dataset. Background correction, normalization, and summarization were performed using the Robust Multi-array Average (RMA) method with the R affy and simpleaffy packages.²⁶

Construction of R-Loops Prognostic Risk Score Model

Differentially expressed genes (DEGs) were identified using the edgeR R package and intersected with the R-loop-associated gene set. Subsequently, univariate Cox regression, LASSO (Least Absolute Shrinkage and Selection Operator) regression, and multivariate Cox regression analyses were performed to identify nine core genes. These nine genes were used to construct the R-loops prognostic risk score model. The scoring formula for the model is as follows:

Expr represents the gene expression level, and Coef represents the weight coefficients of the genes.

Differential Gene Expression Analysis and Functional Enrichment Analysis

Differential gene expression analysis was conducted using the edgeR R package,²⁷ with thresholds set at p < 0.05, FDR < 0.1, and |log2FC| > 1. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the clusterProfiler R package.^28,29 Gene Set Enrichment Analysis (GSEA) was conducted using the GSVA R package.³⁰

Immune Infiltration Analysis

Immune infiltration analysis was conducted to estimate immune cell composition from bulk RNA-seq data. The Immunedeconv R package was employed, integrating multiple immune infiltration methods, including MCP-counter, quanTIseq, EPIC, ESTIMATE, ssGSEA, xCell, TIMER, and CIBERSORT.^31–38

Immune Therapy Response and Drug Sensitivity

The response to immunotherapy was evaluated using the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu/), and TIDE scores were calculated for each TCGA-LUAD patient. Patients were then stratified into high and low TIDE score groups based on the median value, and differences in R-loops Scores between the two groups were assessed using the Wilcoxon rank-sum test. Patients were categorized as responders (True) or non-responders (False) according to TIDE predictions. The association between R-loops Score groups (high vs low, based on median) and immunotherapy response status was evaluated using the Chi-square test.

Drug sensitivity analysis was performed by predicting the IC50 values of chemotherapeutic agents in TCGA-LUAD patients using the oncoPredict R package,³⁹ based on cell line IC50 data from the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org). Patients were then stratified into high and low R-loops Score groups according to the median score. The predicted IC50 values for each drug were compared between the two groups using the Wilcoxon rank-sum test.

RNA Extraction and Quantification of Key Gene Expression

Total RNA was extracted from normal human bronchial epithelial cells (Beas-2B) and three LUAD cell lines (PC9, A549, and HCC827). All the cell lines were obtained from CHINA CENTER FOR TYPE CULTURE COLLECTION(CCTCC). Additionally, fresh tumor and adjacent non-tumor tissues were collected from LUAD patients undergoing surgical treatment at Nanfang Hospital, Southern Medical University. The mRNA expression levels of key genes in these cell lines, normal cells, and tissue samples were analyzed using qPCR, with primer sequences detailed in Table S1. The study protocol was approved by the Institutional Review Board of Nanfang Hospital, Southern Medical University (Ethics Approval No.: NFEC-202408-K24). Patient-derived lung cancer samples and paraffin-embedded tissues were obtained from the Department of Thoracic Surgery at Nanfang Hospital, with all patients providing written informed consent.

Western Blot

Total cellular proteins were extracted using RIPA lysis buffer (Biotech, China) on ice for 10 minutes. Protein concentrations were determined with a BCA assay kit (Biotech, China) according to the manufacturer’s instructions. Equal amounts of protein samples were resolved by SDS-PAGE with 8% or 10% gels, followed by electrotransfer onto PVDF membranes. The membranes were blocked with 5% non-fat dry milk in TBST at room temperature for 1 hour and then incubated overnight at 4 °C with the appropriate primary antibodies. After washing, membranes were incubated with HRP-conjugated secondary antibodies at ambient temperature for 1 hour. Finally, protein signals were detected using NcmECL SuperUltra chemiluminescence reagent (NCM Biotech, P10060).

Statistical Analysis

Statistical analyses were performed using R software (version 4.3.2). Kaplan-Meier survival curves were generated with the survival package, and the Log rank test was applied to determine statistical significance. Univariate and multivariate Cox proportional-hazards regression analyses were conducted using the survival package to assess the relationship between gene expression, clinical variables, and patient survival. Prognostic gene selection was performed through LASSO regression using the glmnet package. Time-dependent ROC analysis was carried out using the timeROC package to evaluate the predictive performance of the risk model, with results presented as the area under the curve (AUC).

Statistical tests, including the Chi-square test, Wilcoxon rank-sum test, Kruskal–Wallis test, and Pearson correlation analysis, were conducted with the stats package. A p-value < 0.05 was considered statistically significant. The following symbols were used to indicate significance levels: * indicates p < 0.05, ** indicates p < 0.01, and *** indicates p < 0.001.

Results

Construction of R-Loops Prognostic Risk Score Model (R-Loops Score)

The study design is illustrated in a flowchart (Figure 1). We collected 803, 612, and 469 R-loops-associated proteins from three independent studies^18–20 and supplemented this collection with 1185 R-loop regulators(human) obtained from the R-loopBase database (https://rloopbase.nju.edu.cn/index.jsp). After removing duplicates, we compiled a non-redundant list of 1771 R-loops genes. Additionally, differential expression analysis was performed on the TCGA-LUAD dataset, resulting in the identification of 4863 differentially expressed genes (DEGs) (Figure 2a). By intersecting the DEGs with the R-loops gene set, we identified 197 R-loops genes that were differentially expressed in LUAD (Figure 2b).

Figure 1 Workflow for the construction and validation of the R-loops prognostic risk score model (R-loops Score). Note: Red text highlights key methods and validation steps.

Figure 2 Identification of hub R-loops genes. (a) Volcano plot showing the differentially expressed genes (DEGs) between cancer and normal tissues. Red points represent upregulated genes, blue points represent downregulated genes, and gray points represent non-significant genes. The horizontal dashed line represents the FDR threshold, and the vertical dashed lines indicate log2 fold-change thresholds. (b) Venn diagram illustrating the overlap between DEGs (red circle) and R-loop-associated genes (blue circle). A total of 197 genes were identified as both DEGs and R-loop-associated. (c) Forest plot presenting the hazard ratios (HRs) and 95% confidence intervals (CIs) for significant genes derived from univariate Cox regression analysis. Genes with significant HRs (p < 0.05) are highlighted, indicating their potential prognostic value. (d) LASSO regression coefficient profiles of the candidate genes. Each curve corresponds to a gene, and the horizontal axis represents the log-transformed penalty parameter (λ). (e) Partial likelihood deviance plot for LASSO regression. The optimal λ value was determined using cross-validation, as indicated by the dashed line. (f) Forest plot summarizing the results of multivariate Cox regression analysis. Significant genes with adjusted HRs and 95% CIs are shown, highlighting their independent prognostic potential. *p-value < 0.05.

To assess the prognostic significance of these genes, univariate Cox regression analysis was conducted, identifying 37 R-loop-associated genes significantly linked to LUAD prognosis (p < 0.0005, Figure 2c). These genes were further refined using LASSO regression and multivariate Cox regression, resulting in a final set of 9 hub R-loops genes (Figure 2d–f). Using these 9 genes, we developed a prognostic risk assessment model, referred to as the “R-loops Score”. The R-loops Score was calculated using the following formula (Expr represents the gene expression level, and Coef represents the weight coefficients of the genes):

Evaluation and Application of the R-Loops Prognostic Risk Score Model (R-Loops Score) in TCGA Cohort

The prognostic significance of the R-loops-based risk model was evaluated through multiple analyses. Patients were stratified into high-risk and low-risk groups based on the median R-loops Score. Kaplan-Meier survival analysis (Figure 3a) demonstrated that patients in the high-risk group had significantly shorter overall survival (OS) compared to those in the low-risk group (p < 0.0001). Time-dependent ROC curves (Figure 3b) showed that the model had good predictive performance, with AUC values of 0.704, 0.689, and 0.641 for 1-year, 2-year, and 5-year OS, respectively.

Figure 3 Prognostic performance of the R-loops prognostic risk score model (R-loops Score) in predicting overall survival. (a) Kaplan-Meier survival curves for patients in the high R-loops Score (blue) and low R-loops Score (Orange) groups. The number of patients at risk at each time point is displayed below the plot, and the p-value indicates significant differences between the two groups. (b) Time-dependent receiver operating characteristic (ROC) curves for predicting 1-year, 2-year, and 5-year OS. The area under the curve (AUC) values for each time point are provided. (c) Forest plot of multivariate Cox regression analysis for overall survival. Variables include age, gender, stage, and risk group. Hazard ratios (HRs) and 95% confidence intervals (CIs) are shown. (d) Nomogram integrating age, stage, and risk score for predicting 1-year, 2-year, and 5-year OS. The total points correspond to the estimated survival probabilities. (e) Calibration plot assessing the agreement between predicted and observed OS at 1 year, 2 years, and 5 years. The diagonal dashed line represents perfect calibration, while colored lines indicate model performance. (f) ROC curves of the nomogram for predicting OS at 1 year, 2 years, and 5 years. The AUC values reflect the model’s predictive accuracy over time.

Univariate Cox regression analysis across clinical subgroups (Figure 3c) revealed that the R-loops Score was an independent prognostic factor. Patients in the high-risk group exhibited a hazard ratio (HR) of 2.07 (95% CI: 1.51–2.82), while advanced stages (Stage III and IV) were also strongly associated with worse prognosis (HR = 2.98 and 3.52, respectively).

To provide individualized survival predictions, a nomogram was constructed (Figure 3d) incorporating the R-loops Score, age, and clinical stage. Calibration curves (Figure 3e) showed strong concordance between predicted and observed OS for 1-year, 2-year, and 5-year outcomes. Additionally, ROC curves for independent validation (Figure 3f) confirmed the robustness of the model, with AUC values of 0.732, 0.713, and 0.719 for 1-year, 2-year, and 5-year OS, respectively.

These findings suggest that the R-loops-based risk model has the potential to provide valuable insights for prognostic assessment in LUAD patients; however, further validation is necessary to confirm its clinical applicability.

Validation of the R-Loops Prognostic Risk Score Model (R-Loops Score) in GEO Cohorts

To assess the generalizability of the R-loops prognostic risk score model (R-loops Score), we applied it to three independent GEO cohorts: GSE31210, GSE72094, and GSE50081. Patients in each cohort were stratified into high-risk and low-risk groups based on the median R-loops Score.

In the GSE31210 cohort, Kaplan-Meier survival analysis (Figure 4a) demonstrated that patients in the high-risk group had significantly poorer overall survival (OS) compared to those in the low-risk group (p < 0.0001). The time-dependent ROC curves (Figure 4b) showed that the R-loops Score achieved AUC values of 0.779, 0.767, and 0.736 for 1-year, 2-year, and 5-year OS, respectively, indicating good predictive performance. Similarly, in the GSE72094 cohort, Kaplan-Meier survival analysis (Figure 4c) revealed a significant survival difference between the high-risk and low-risk groups (p = 0.00022). The ROC curves (Figure 4d) demonstrated AUC values of 0.693, 0.671, and 0.709 for 1-year, 2-year, and 5-year OS, respectively. In the GSE50081 cohort, Kaplan-Meier analysis (Figure 4e) also showed a significant difference in survival between the two risk groups (p = 0.011). The ROC analysis (Figure 4f) yielded AUC values of 0.744, 0.652, and 0.639 for 1-year, 2-year, and 5-year OS, respectively.

Figure 4 Validation of the prognostic value of the R-loop score in multiple independent datasets. (a) Kaplan–Meier survival curves comparing overall survival (OS) between high and low R-loop score groups in the GSE31210 dataset. (b) Time-dependent receiver operating characteristic (ROC) curves for 1-, 2-, and 5-year OS prediction in GSE31210. (c) Kaplan–Meier survival curves comparing OS between R-loop score groups in the GSE72094 dataset. (d) Time-dependent ROC curves for 1-, 2-, and 5-year OS prediction in GSE72094. (e) Kaplan–Meier survival curves comparing OS between R-loop score groups in the GSE50081 dataset. (f) Time-dependent ROC curves for 1-, 2-, and 5-year OS prediction in GSE50081.

These results suggest the potential robustness and generalizability of the R-loops Score in predicting overall survival across multiple independent GEO cohorts, highlighting its promise as a reliable prognostic tool.

Association Between the R-Loops Score and Clinical Characteristics

To investigate the association between the R-loops Score and various clinical characteristics, we performed statistical analyses on the TCGA-LUAD dataset.

We found no significant correlation between the R-loops Score and either age (p = 0.13, Figure 5a) or smoking history (p = 0.46, Figure 5b). However, the R-loops Score was significantly associated with gender, with male patients exhibiting higher scores than female patients (p = 0.032, Figure 5c). In terms of clinical staging, higher R-loops Scores were observed in patients with more advanced stages. As shown in Figure 5d–f, patients in stages III and IV had higher scores compared to those in stages I and II. Similarly, patients with higher T stages (T3 and T4) had significantly elevated R-loops Scores compared to those in T1 and T2. Additionally, patients with lymph node metastasis (N1-N3) exhibited higher R-loops Scores compared to those without lymph node involvement (N0). However, no significant difference in R-loops Scores was observed between patients with and without distant metastasis (pM0 vs pM1, p = 0.15, Figure 5g).

Figure 5 Analysis of R-loops Score across clinical and pathological features in patients. (a) Correlation between R-loops Score and patient age. The scatter plot displays the distribution of R-loops Score by age, with a Pearson correlation coefficient (R) of −0.068 and a p-value of 0.13, indicating no significant correlation. (b) Comparison of R-loops Score between patients with and without a smoking history using the Wilcoxon rank-sum test (p = 0.46), showing no significant difference. (c) R-loops Score between male and female patients compared by the Wilcoxon rank-sum test (p = 0.032), showing a statistically significant difference. (d) Distribution of R-loops Score across cancer stages (I–IV) analyzed using the Kruskal–Wallis test (p = 3e-09), revealing significant differences among stages. (e) Comparison of R-loops Score across pathological tumor stages (T1-T4) using the Kruskal–Wallis test (p = 2e-07), showing significant differences. (f) Distribution of R-loops Score among lymph node metastasis (pN0-N3) analyzed using the Kruskal–Wallis test (p = 6e-07), showing significant differences. (g) Comparison of R-loops Score between patients with and without distant metastasis (pM0 vs pM1) using the Wilcoxon rank-sum test (p = 0.15), showing no significant difference.

Enrichment Analysis of High and Low R-Loops Score Groups with Multiple Datasets

To investigate the biological pathways associated with the high and low R-loops Score groups, we conducted GO and KEGG analyses across four datasets (TCGA-LUAD, GSE50081, GSE72094, and GSE31210). By intersecting the enriched pathways from these datasets, we identified 111 shared pathways (Figure 6a). Among these, Figure 6b highlights pathways with a Count > 30. The enrichment analysis revealed that the high R-loops Score group were significantly enriched in pathways and biological processes associated with the cell cycle, chromosome segregation, and organelle assembly.

Figure 6 Functional and pathway enrichment analysis of genes in R-loops Score model. (a) Venn diagram illustrating the overlap of significantly enriched pathways identified in LUAD datasets (TCGA-LUAD, GSE50081, GSE72094, and GSE31210). (b) Enrichment analysis of the overlapping pathways from (a), including Gene Ontology (GO) categories—biological process (BP), cellular component (CC), and molecular function (MF)—and KEGG pathways. Bars represent the top enriched GO terms and KEGG pathways, with their lengths corresponding to -log10 transformed p-values, indicating statistical significance. (c) Venn diagram showing the overlap of hallmark pathways identified by GSEA analysis across different datasets. The intersections highlight the shared hallmark pathways among these datasets. (d) Gene Set Enrichment Analysis (GSEA) results for hallmark pathways shared in (c). The lower right subpanel displays the hallmark gene sets with significant enrichment, along with their p-values and adjusted p-values.

Key biological processes (BP) included negative and positive regulation of the cell cycle, mitotic cell cycle phase transition, chromosome segregation, sister chromatid segregation, and non-membrane-bounded organelle assembly. Interestingly, the enrichment of non-membrane-bounded organelle assembly suggests that R-loops Score may be associated with the formation and regulation of membraneless organelles, which are driven by phase separation processes. For cellular components (CC), enriched pathways included the chromosomal region, midbody, spindle, and microtubule-associated complexes, emphasizing their involvement in mitotic and chromosomal dynamics. Molecular functions (MF) such as microtubule binding and microtubule motor activity were also significantly enriched, indicating their roles in structural organization and motor processes essential for cell division. KEGG pathway enrichment highlighted “Motor proteins” and “Cell cycle” underscoring the critical roles of these genes in cellular division and chromosomal segregation.

In addition, we employed the GSEA method to analyze hallmark pathways across the four datasets. The analysis revealed 12 pathways with consistent and statistically significant enrichment across all datasets (Figure 6c). The detailed results for these pathways are presented in Figure 6d, including key cancer-related pathways such as DNA repair, E2F targets, epithelial-mesenchymal transition (EMT), G2M checkpoint, and mitotic spindle, highlighting the association between high R-loops Score and tumor progression.

Correlation Analysis Between R-Loops Score and Immune Characteristics

To investigate the relationship between R-loops Score and immune-related characteristics, we performed a series of correlation analyses. Scatter plots revealed a significant positive correlation between R-loops Score and ImmuneScore (R = 0.16, p = 0.00025) as well as ESTIMATEScore (R = 0.12, p = 0.0072). However, no significant correlation was observed between R-loops Score and StromalScore (R = 0.06, p = 0.18) (Figure 7a–c). A heatmap analysis demonstrated the association between R-loops Score and immune checkpoint molecules (Figure 7d). Notably, R-loops Score showed a weak but significant positive correlation with CD274 (R = 0.11, p = 0.017) and CTLA4 (R = 0.1, p = 0.019), as illustrated in Figure 7e and f.

Figure 7 Correlation between R-loops Score, immune microenvironment, and immune checkpoint. (a–c) Correlation of R-loops Score with ImmuneScore (a), StromalScore (b), and ESTIMATEScore (c) based on ESTIMATE analysis. (d) Heatmap showing the correlation between R-loops Score and the expression of immune checkpoint molecules. (e and f) Scatter plots showing the correlation of R-loops Score with CD274 (PD-L1, (e)) and CTLA4 (f), with significant positive correlations (CD274: R = 0.11, p = 0.017; CTLA4: R = 0.1, p = 0.019). (g) Venn diagram illustrating the overlap of significantly enriched immune-related pathways identified across LUAD datasets (TCGA-LUAD, GSE31210, GSE72094, and GSE50081). The central intersection represents common pathways shared across datasets. (h) Heatmap of correlations between R-loops Score and immune cell components shared in (g). (i) Correlation heatmap showing the relationships between the 9 R-loops genes and immune cell components estimated by CIBERSORT. *p < 0.05, **p < 0.01, ***p < 0.001.

Subsequently, we performed immune infiltration analysis on TCGA-LUAD, GSE31210, GSE72094, and GSE50081 datasets. By identifying immune cell components significantly associated with R-loops Score in each dataset and taking their intersection, we found seven immune components that were consistently and significantly correlated with R-loops Score across all four datasets (Figure 7g). A heatmap provided detailed correlation information for these shared immune components (Figure 7h).

Finally, we analyzed the correlation between the 9 genes included in the R-loops Score and immune cell components estimated using cibersort algorithm. As shown in Figure 7i, the results highlight the relationships between these genes and various immune cell components, providing deeper insights into the interplay between the R-loops Score and immune regulation.

Correlation Between R-Loops Score and Drug Treatment Response

To evaluate the relationship between R-loops Score and response to immunotherapy or drug treatments, we performed Tumor Immune Dysfunction and Exclusion (TIDE) analysis and drug sensitivity assessments.

The results revealed that tumors with high R-loops Scores exhibited significantly higher TIDE scores compared to those with low R-loops Scores (Wilcoxon test, p = 3.2e-11, Figure 8a), suggesting a stronger capability for immune evasion. Furthermore, the distribution of immune responders revealed a lower proportion of immune-responsive patients in the high R-loops Score group (chi-squared test, p = 2.6e-06, Figure 8b), which is consistent with the previous conclusion that tumors with High R-loops Scores exhibit stronger immune evasion.

Figure 8 Association between R-loops Score and therapeutic response. (a) Violin plot showing a significant difference in TIDE (Tumor Immune Dysfunction and Exclusion) scores between low and high R-loops Score groups (Wilcoxon test, p = 3.2e-11), indicating immune dysfunction in high R-loops Score tumors. (b) Proportion of responders and non-responders to immunotherapy in low and high R-loops Score groups. A significant association was observed (chi-square test, p = 2.6e-06). (c–h) Violin plots showing sensitivity (as measured by drug response scores) to six chemotherapeutic agents between low and high R-loops Score groups. Significant differences are observed for cisplatin ((c), p = 2e-04), paclitaxel ((d), p = 1.1e-15), docetaxel ((e), p < 2.2e-16), vincristine ((f), p = 3.4e-07), gefitinib ((g), p = 2.5e-05), and erlotinib ((h), p = 3.6e-06), with higher sensitivity associated with specific R-loops Score groups.

Next, we assessed the association of R-loops Score with sensitivity to various chemotherapeutic agents. Tumors with high R-loops Scores demonstrated significantly sensitivity to cisplatin (Wilcoxon test, p = 2e-04), paclitaxel (p = 1.1e-15), docetaxel (p < 2.2e-16), vincristine (p = 3.4e-07), gefitinib (p = 2.5e-05), and erlotinib (p = 3.6e-06) compared to tumors with low R-loops Scores (Figure 8c–h). This trend was generally observed across stages I–III when stratified by clinical stage, as shown in Figure S1, whereas the differences were attenuated in stage IV, potentially due to increased tumor heterogeneity.

These findings suggest that patients with low R-loops Scores may be more responsive to immunotherapy but relatively resistant to traditional chemotherapy and targeted therapies, which may help inform clinical treatment strategies.

Evaluation of Gene Expression in the R-Loops Prognostic Risk Score Model

We employed qPCR to assess the mRNA expression levels of multiple genes in LUAD cell lines (A549, PC9, H23) and the normal human bronchial epithelial cell line (Beas-2B, referred to as 2B) (Figure 9a). The results showed that, compared to the normal cell line (2B), IGF2BP1, SFXN1, LDHA, FKBP4, CCT6A, FSCN1, and EIF3B were upregulated in multiple lung cancer cell lines, whereas WDHD1 and ALDOA were significantly downregulated across all LUAD cell lines. Further analysis of these genes in LUAD patient tumor tissues and paired normal tissues (Figure 9b) revealed that, apart from ALDOA and EIF3B, the expression differences of the other genes between tumor and normal tissues were not statistically significant (p > 0.05). Specifically, ALDOA was significantly downregulated in tumor tissues compared to normal tissues (p < 0.05), while EIF3B was significantly upregulated in tumor tissues (p < 0.05).

Figure 9 The qPCR and Western blot analysis of gene expression in LUAD cell lines and patient tissues. (a) Relative mRNA expression levels of selected genes in LUAD cell lines (A549, PC9, H23) and the normal human bronchial epithelial cell line (Beas-2B, referred to as 2B). Gene expression levels were normalized to an internal control, and statistical comparisons were performed relative to the 2B cell line. (b) Relative mRNA expression levels of the same genes in LUAD tumor tissues and paired adjacent normal tissues. Statistical significance was determined using Student’s t-test. (c) Western blot analysis of EIF3B and ALDOA protein expression in LUAD cell lines and the 2B normal cell line. β-tubulin was used as a loading control. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001. Data are presented as mean ± SEM. Abbreviation: ns, not significant.

To validate gene expression at the protein level, we performed Western blot (WB) analysis in LUAD cell lines and the 2B normal cell line (Figure 9c). The results showed that EIF3B protein levels were consistently upregulated in cancer cell lines, which aligns with the transcriptomic and qPCR findings. In contrast, ALDOA protein levels were elevated in LUAD cell lines, showing an opposite trend compared to its mRNA expression, suggesting possible post-transcriptional regulation or differential protein stability. In addition, siRNA-mediated knockdown of EIF3B or ALDOA resulted in no marked changes in drug sensitivity (Figure S2), supporting the view that drug response is likely driven by the integrated activity of multiple genes rather than any single gene alone.

The consistent upregulation of EIF3B at both mRNA and protein levels highlights its potential as a robust biomarker in LUAD.

Discussion

In this study, we developed and validated a novel R-loops Score model based on 9 hub R-loop-associated genes to predict prognosis, immune characteristics, and drug response in LUAD. Our findings demonstrate the complex roles of R-loops in cancer progression, highlighting their potential as both prognostic biomarkers and therapeutic targets.

Patients with high R-loop scores exhibited significantly shorter overall survival (OS) across multiple cohorts, further validating the robustness of the model. By integrating the R-loop score with clinical variables into a nomogram, we achieved highly accurate survival predictions, highlighting its clinical relevance. These results align with the pivotal roles of R-loops in maintaining genomic stability and regulating transcriptional dynamics.^40,41

Enrichment analysis revealed pathways such as cell cycle, mitotic cell cycle phase transition, chromosome segregation, and sister chromatid segregation, all of which are strongly associated with genomic instability and align with established findings. Interestingly, non-membrane-bounded organelle assembly, a process linked to phase separation, was also enriched. These findings suggest that R-loops not only contribute to mitotic progression and cellular organization but may also play a role in phase separation-mediated biological processes, such as the assembly of membraneless organelles.^42,43

Our analysis uncovered a positive correlation between the R-loops Score and immune evasion, including higher TIDE scores and increased expression of immune checkpoint molecules. Tumors with high R-loops Scores exhibited a lower proportion of immune-responsive patients, suggesting that R-loops may facilitate immune evasion and reduce responsiveness to immunotherapy. Conversely, tumors with high R-loops Scores demonstrated increased sensitivity to multiple chemotherapeutic and targeted agents, including cisplatin, paclitaxel, and EGFR inhibitors. This dual observation highlights the potential of R-loops as predictive markers for both immunotherapy resistance and chemotherapy responsiveness, emphasizing the importance of tailoring treatment strategies based on R-loops profiling.

To further validate the nine hub genes included in the R-loops Score model, we performed qPCR analysis in LUAD cell lines and patient tissues. Our results demonstrated that most of these genes (IGF2BP1, SFXN1, LDHA, FKBP4, CCT6A, FSCN1, and EIF3B) were significantly upregulated in LUAD cell lines compared to normal human bronchial epithelial cells, while WDHD1 and ALDOA were consistently downregulated. However, when analyzing their expression in LUAD tumor tissues and paired adjacent normal tissues, only EIF3B and ALDOA exhibited significant differential expression, with EIF3B upregulated and ALDOA downregulated in tumor tissues. These findings suggest that while the model’s genes are relevant in LUAD, their expression patterns in cell lines may not fully reflect clinical conditions, likely due to tumor microenvironmental influences and transcriptional regulatory differences in vitro. The differential expression of EIF3B and ALDOA in patient tissues highlights them as potential biomarkers for LUAD and warrants further investigation into their functional roles. To further validate their expression at the protein level, Western blot analysis was performed in LUAD cell lines. EIF3B protein levels were consistently upregulated in cancer cells, corroborating the transcriptomic and qPCR findings. In contrast, ALDOA protein levels were elevated in LUAD cell lines despite being downregulated at the mRNA level, suggesting a possible post-transcriptional regulatory mechanism. These findings underscore the robustness of EIF3B as a candidate biomarker and highlight the importance of integrating multi-omics approaches when evaluating gene function and clinical relevance.

To better understand the biological relevance of the R-loops Score model, we further explored the functional roles of the nine hub genes. EIF3B, a translation initiation factor, promotes cell proliferation and cisplatin resistance in LUAD, and its overexpression correlates with poor clinical outcomes.⁴⁴ ALDOA, a glycolytic enzyme, promotes proliferation and cisplatin resistance in lung cancer cells.⁴⁵ IGF2BP1, an m⁶A reader, stabilizes oncogenic mRNAs, driving tumor progression, metastasis, and therapy resistance in LUAD.⁴⁶ WDHD1 (AND-1), involved in DNA replication and repair, is essential for NSCLC growth and its inhibition increases radiosensitivity.⁴⁷ FKBP4 participates in protein folding and chromatin regulation; while detailed studies in LUAD are limited, it interfaces with genome integrity mechanisms linked to R-loops. LDHA, a key enzyme in glycolysis, contributes to stemness, immune suppression, and drug resistance under hypoxic conditions in lung tumors.⁴⁸ FSCN1, an actin-bundling protein, is upregulated in NSCLC and promotes cell migration, invasion, and metastasis.⁴⁹ CCT6A, part of the chaperonin family, is co-amplified with EGFR, and mediates LUAD progression via metabolic reprogramming.⁵⁰ SFXN1, a mitochondrial serine transporter, has prognostic relevance in LUAD and plays roles in mitochondrial function and tumor cell survival.⁵¹ Collectively, these genes participate in key oncogenic pathways, including genomic stability maintenance, metabolic reprogramming, and post-transcriptional regulation, many of which are intimately connected to R-loop biology. Their functional provides a mechanistic basis for the prognostic and therapeutic predictive value of the R-loops Score model in LUAD.

Despite these promising insights, several limitations should be addressed. Firstly, this study utilized data from publicly available datasets, which may introduce sampling bias. In certain cohorts, the limited patient sample size could restrict the generalizability of our findings to larger and more diverse populations. Moreover, factors such as racial differences, clinical characteristics, and treatment histories could influence the results and should be considered. Secondly, while our study offers strong statistical evidence linking R-loops to immune regulation and chemotherapy sensitivity, it lacks direct experimental validation. Specifically, the absence of functional studies, such as in vitro or in vivo experiments, limits our ability to establish a causal relationship between R-loops-associated genes and the observed clinical outcomes. Conducting targeted experiments would provide critical mechanistic insights and significantly enhance the translational potential of our findings. Lastly, the absence of direct data on immunotherapy responses in LUAD limits our ability to fully explore the relationship between R-loops and treatment outcomes. This limitation highlights the necessity for future studies that integrate comprehensive clinical and treatment response data to provide a more holistic understanding.

Conclusion

In this study, we developed and validated a novel R-loops Score model based on nine genes to predict prognosis, immune landscape, and drug sensitivity in LUAD. This score demonstrated robust prognostic value across multiple cohorts and was associated with key biological processes, including genomic instability, immune evasion, and therapeutic responsiveness. Our findings highlight the potential of R-loops not only as mechanistic drivers of tumor progression but also as clinical biomarkers. These results provide a foundation for future translational research and underscore the promise of R-loop–targeted strategies in precision oncology.