Revealing Causal Genes in Nasopharyngeal Carcinoma: An Integrated Approach of Mendelian Randomization and Single-Cell Sequencing

Transcriptomic Data Analysis

Two nasopharyngeal carcinoma (NPC) bulk transcriptomic datasets, GSE53819 and GSE12452, were retrieved from the Gene Expression Omnibus (GEO) database, comprising a total of 49 NPC tissues and 28 non-cancerous nasopharyngeal tissues. Raw expression data were downloaded from their respective platforms, and probe identifiers were mapped to standardized gene symbols using R (version 4.4.1). Expression profiles within each dataset were integrated using the merge function, and genes with missing values were removed to construct independent expression matrices for subsequent analyses.

Differential expression analysis was performed using the “limma” package (version 3.60.6). The analysis was conducted on quantified gene expression matrices that had been normalized across samples during upstream preprocessing and were represented as log₂(TPM + 1) values. As the input data consisted of log-transformed continuous expression values and no longer followed the mean–variance relationship characteristic of raw count data, the voom transformation was not applied.

Linear models were fitted using the lmFit function, and contrast matrices were constructed with makeContrasts to enable group-wise comparisons. Empirical Bayes moderation was subsequently applied using the eBayes function to improve the robustness of statistical inference. Multiple testing correction was performed using the Benjamini–Hochberg method, and genes with an adjusted P value < 0.05 were considered statistically significant. Differentially expressed genes (DEGs) were extracted using the topTable function, with entries containing missing values excluded. Volcano plots were generated using the HiOmics cloud platform (https://www.henbio.com/)

⁵

.

Mendelian Randomization Analysis

The overlapping differentially expressed genes (DEGs) were converted to Ensembl gene identifiers using the Ensembl BioMart database. In the Mendelian randomization analysis, the differential expression status of genes was not directly used as the exposure variable, but served only to screen for genes with potential biological relevance. The actual exposure variable was defined as the expression level of the corresponding gene, represented by expression quantitative trait loci (eQTLs). For each candidate gene, summary statistics of associated eQTLs were obtained from the OpenGWAS database, and single nucleotide polymorphisms (SNPs) significantly associated with gene expression were selected as instrumental variables (IVs). IVs were filtered based on the significance of SNP‑gene expression associations (p ≤ 1×10⁻⁵). Selected SNPs then underwent linkage disequilibrium clumping (kb = 10,000, R² = 0.001) to ensure independence among instruments, and palindromic variants were excluded. The strength of the instrumental variables was assessed using the F-statistic, with F > 10 considered sufficient to avoid weak instrument bias

⁶

. Finally, using gene expression as the exposure and GWAS summary statistics for nasopharyngeal carcinoma from the FinnGen database as the outcome, Mendelian randomization analysis was performed to evaluate potential causal relationships between gene expression and NPC risk.

MR analyses were conducted under three core assumptions: (1) IVs are strongly associated with the exposure (gene expression); (2) IVs are independent of confounders; and (3) IVs affect the outcome (NPC) only through the exposure. The inverse variance weighted (IVW) method served as the primary approach

⁶

. Four supplementary methods were also applied: MR-Egger, weighted median, simple mode, and weighted mode

^7,

^8,

⁹

.

Sensitivity analyses included Cochran’s Q test to assess heterogeneity and MR-Egger intercept test to evaluate horizontal pleiotropy. A p-value < 0.05 indicated significant heterogeneity or potential pleiotropy

^10,

¹¹

. All analyses were performed in R using the “TwoSampleMR” (v0.6.8) and “MR-PRESSO” (v1.0) packages.

Single-Cell RNA-Seq Processing and Analysis

Two single-cell RNA sequencing datasets of nasopharyngeal carcinoma, GSE150825 and GSE120926, were downloaded from the Gene Expression Omnibus (GEO) database. GSE150825 comprises 3 samples of nasopharyngeal lymphoid hyperplasia and 11 nasopharyngeal carcinoma samples, while GSE120926 includes 16 nasopharyngeal carcinoma samples and 8 non-malignant nasopharyngeal tissue samples. All non-malignant samples and lymphoid hyperplasia samples were uniformly classified as the non-cancerous control group.

Single-cell data analysis was performed in the R environment using the “Seurat” package (version 4.1.1). Quality control is achieved by selecting cells that contain at least 200 detectable genes and a mitochondrial gene expression ratio below 15%. Subsequently, cell cycle stages were assigned to each cell using the “CellCycleScoring” function to assess potential cell cycle effects.

During data preprocessing, the gene expression matrix of each dataset was log-normalized using the “LogNormalize” method in Seurat (scale factor = 10,000). To capture the major biological variations across cells, highly variable genes were selected using the “FindVariableFeatures” function, retaining the top 2000 most variable genes. The data were then linearly scaled based on these selected genes using the “ScaleData” procedure, followed by principal component analysis (PCA).

To mitigate batch effects between the two datasets, the Harmony algorithm was applied to integrate them within the PCA space. Using the first 15 principal components after integration, cell-cell proximity was computed with the “FindNeighbors” function to construct a k-nearest neighbor graph for downstream analysis. UMAP dimensionality reduction was then performed for visualization and unsupervised clustering.

The choice of clustering resolution was evaluated using the “Clustree” package (version 0.5.0) by comparing cluster stability and hierarchical relationships across multiple resolution parameters. A resolution of 0.2 was ultimately selected. Differentially expressed genes for each cell cluster were then identified using the “FindAllMarkers” function.

Cell type annotation was performed manually based on the top 100 marker genes per cluster, combined with known canonical cell marker genes and relevant literature. Finally, the expression patterns of key genes identified through Mendelian randomization analysis across different cell subtypes were visualized using the “DotPlot” function in Seurat, where dot color and size represent the average expression level and the proportion of cells expressing the gene, respectively. In addition, stacked bar plots illustrating cellular composition between tumor and non-tumor tissues, as well as violin plots and boxplots showing the expression of key genes in corresponding cell types, were generated using the “ggplot2” package (version 3.3.6).

Identification of Differentially Expressed Genes

Differential expression analysis was performed on the GSE53819 and GSE12452 datasets by comparing tumor tissues to non-cancerous controls, using thresholds of |log₂FC| ≥ 1 and adjusted P < 0.05. In the GSE12452 dataset, 939 differentially expressed genes (DEGs) were identified, including 575 up-regulated and 364 down-regulated genes. The GSE53819 dataset yielded 4,212 DEGs, of which 1,823 were up-regulated and 2,389 were down-regulated (Figures 1A,B). A Venn diagram generated with the R package “ggvenn” (v0.1.10) revealed 494 shared DEGs between the two datasets (Figure 1C).

Mendelian Randomization Identifies Causal Genes for NPC

GWAS data for eQTLs of 313 shared DEGs were obtained from the OpenGWAS database, and summary statistics for nasopharyngeal carcinoma (FinnGen code: C3_NASOPHARYNX_EXALLC) were retrieved from the FinnGen consortium. MR analysis was conducted to estimate the causal effects of these genes on NPC risk. The Wald ratio method was applied for exposures with one SNP, and the inverse variance weighted (IVW) method was used for exposures with two or more SNPs. Four supplementary methods—MR-Egger, weighted median, simple mode, and weighted mode—were also applied for robustness. Associations with P < 0.05 in the primary method (Wald ratio or IVW) were considered statistically significant.Wald ratio analysis identified MFSD4 as a risk gene for NPC (OR = 11.020, 95% CI: 1.637–74.150, P = 0.014).And the wide confidence interval observed for MFSD4 likely reflects the limited number of instrumental variables available for this gene and the relatively large standard error of the Wald ratio estimate. (Table1)

IVW analysis of exposures with exactly two SNPs revealed that PSPH was associated with increased NPC risk (OR = 3.942, 95% CI: 1.135–13.685, P = 0.031), whereas THBS2 was associated with reduced risk (OR = 0.211, 95% CI: 0.047–0.939, P = 0.041). No significant heterogeneity was detected for these associations (Cochran’s Q P > 0.05).(Table 2)

For exposures with more than two SNPs, IVW analysis indicated that TMEM200A was linked to decreased NPC risk (OR = 0.288, 95% CI: 0.106–0.781, P = 0.014), while VPREB3 was associated with increased risk (OR = 2.903, 95% CI: 1.192–7.068, P = 0.019). (Table 3)

Heterogeneity testing revealed no significant heterogeneity between TMEM200A and VPREB3 in relation to nasopharyngeal carcinoma (Cochran's Q p > 0.05). Furthermore, pleiotropy analysis indicated no detectable pleiotropic effect of these two genes on nasopharyngeal carcinoma (p > 0.05). (Table 4)

Single-Cell Expression Profiling of Key Genes

To evaluate batch effects and determine an appropriate clustering resolution for downstream single-cell analyses, dimensionality reduction and clustering diagnostics were performed. Prior to batch correction, principal component analysis (PCA) revealed that cells were primarily separated according to dataset and sample origin, indicating the presence of pronounced batch effects (Figure 2A). After integration using Harmony, cells from different datasets and samples were well mixed in the low-dimensional space, demonstrating effective batch correction (Figure 2B).

To further assess clustering behavior across different resolutions, clustering results were systematically examined over a range of resolution parameters using a cluster tree visualization. The clustering tree showed that clusters underwent progressive and stable refinement as resolution increased. At a resolution of 0.2, major cell populations were consistently preserved with limited cluster fragmentation, providing a balance between cluster stability and granularity. Based on this evaluation, a resolution of 0.2 was selected for subsequent clustering and downstream analyses (Figure 2C).

After quality control, 157,289 cells from two single-cell RNA-seq datasets were retained for analysis. Unsupervised clustering at a resolution of 0.2 identified 11 distinct cell populations. These clusters were annotated based on established marker genes as follows: B cells (CD79A), CD4T cells (CD3D, CD4), CD8T cells (CD8A), Myeloid cells (LYZ), plasma B cells (SDC1), epithelial cells (EPCAM), fibroblasts (COL1A1), endothelial cells (VWF), mast cells (TPSAB1), plasmacytoid dendritic cells (LILRA4), and Cycling cells (MKI67).

As shown in Figure 3A and 3C, single-cell RNA sequencing data were annotated into distinct cell populations, and the expression of canonical marker genes confirmed the identity of each annotated cell type. Compared with non-tumor tissues, nasopharyngeal carcinoma tissues exhibited a marked remodeling of immune cell composition, characterized by increased proportions of CD8⁺ T cells, CD4⁺ T cells, myeloid cells, and epithelial cells, along with a relative reduction in B cells (Figure 3B). In addition, tumor tissues showed enrichment of fibroblast and endothelial cell populations, suggesting expansion of stromal and vascular components within the NPC tumor microenvironment.

To further characterize the cellular context of MR-identified genes, their expression patterns were examined across annotated cell types and between tumor and non-tumor samples. Dot plot analysis revealed distinct cell-type-specific expression profiles of the five candidate genes (Figure 3D). PSPH was predominantly expressed in epithelial cells. TMEM200A and THBS2 expression was enriched in fibroblasts; notably, TMEM200A was also detected in CD4+ and CD8+ T cells. MFSD4 was mainly found in endothelial cells, and VPREB3 was highly expressed in B cells, with a high percentage of cells expressing it. Comparison between tumor and non-tumor samples further demonstrated differential expression of these genes at the cellular level (Figure 3E). In particular, VPREB3 exhibited higher expression and a larger fraction of expressing cells in non-tumor samples, whereas other genes showed more subtle but cell-type-dependent expression differences, indicating that MR-identified genetic signals are reflected in distinct cellular expression patterns within the tumor microenvironment.

To further characterize the expression patterns of the Mendelian randomization-identified genes under disease conditions, stratified comparisons were performed between nasopharyngeal carcinoma (NPC) samples and control samples within the primary cell types expressing these genes (Figure 4). In epithelial cells, the risk gene PSPH was significantly upregulated in NPC samples compared with controls (Wilcoxon test, P<0.0001), indicating its aberrant activation in tumor-associated epithelial cells. In fibroblasts, the protective gene THBS2 exhibited significant expression changes in NPC tissues (P<0.0001), and TMEM200A also showed a statistically significant difference in this cell population (P<0.05), suggesting pronounced molecular remodeling of stromal cells during NPC development. Among immune cells, TMEM200A displayed a modest yet significant expression difference in CD8⁺T cells (P<0.05), implying its potential involvement in regulating tumor-related immune status. Furthermore, in B-cell lineages, the risk gene VPREB3 showed significant differential expression between NPC and control samples in B cells, cycling cells, and plasma cells (all P<0.0001), demonstrating consistent disease-associated expression alterations. In contrast, MFSD4 expression in endothelial cells did not differ significantly between NPC and control samples (P>0.05), despite its clear cell type-specific expression pattern at the single-cell level. These results indicate that the expression imbalance of MR-identified pathogenic and protective genes in NPC is not uniformly driven by all cell types, but is primarily mediated by specific immune and stromal cell populations.

Limitations

Several limitations should be considered in interpreting our findings. MR estimates may be biased in the presence of horizontal pleiotropy, although sensitivity analyses in this study did not detect significant pleiotropic effects. Furthermore, while MR provides evidence for causal inference, the conclusions require validation through functional experiments and clinical studies. The GWAS and eQTL data used in this study were primarily derived from European populations, which may limit the generalizability of the results to other ethnic groups. Additionally, static eQTL data do not capture dynamic gene regulation over time or context, which may affect the interpretation of gene–disease causality. Future studies leveraging larger eQTL reference panels or tissue-specific regulatory datasets may help refine the causal effect estimate of MFSD4. And more diverse cohorts and time-resolved functional genomics data will help refine these findings.