Mendelian Randomization Transcriptomics and Network Pharmacology for Identification of Osteoarthritis Drug Targets

Summary

We present a methodological protocol that integrates Mendelian randomization, transcriptomic analysis, and network pharmacology to systematically identify candidate drug targets for osteoarthritis. This approach enables the prioritization of targets through genetic, expression, and interaction evidence, using publicly available data and computational tools.

Abstract

This article delineates an integrated methodology that combines Mendelian randomization (MR), transcriptomic analysis, and network pharmacology to identify and prioritize potential therapeutic targets for osteoarthritis (OA). It is designed to guide researchers in implementing this multimodal pipeline to investigate drug repurpose and the development of novel therapeutic interventions for OA. The methodological framework comprises five sequential stages: first, MR analysis pipeline, employing two-sample MR to identify putative causal plasma proteins associated with OA, followed by Steiger filtering and phenome-wide association scanning to assess causal directionality and potential off-target effects; second, transcriptomic sequencing workflow, integrating RNA-seq data to identify and refine candidate protein targets; third, integration strategy, merging MR and transcriptomic results to prioritize candidate proteins; fourth, network pharmacology and molecular docking procedures, involving the construction of protein-protein interaction networks, functional enrichment analysis, and molecular docking to explore ligand-target interactions; and fifth, intended application, focusing on the prioritization of candidate compounds and natural products with potential therapeutic relevance. By organizing the workflow into distinct analytical phases, this framework provides a reproducible approach for transitioning from genetic and transcriptomic discovery to computational drug-target evaluation, without presenting specific experimental outcomes. The methodology facilitates systematic and hypothesis-driven investigation into OA therapeutics using publicly accessible datasets and computational tools.

Introduction

Osteoarthritis (OA) is characterized by the progressive deterioration of cartilage, changes in subchondral bone, and excessive bone growth¹. It is now a major cause of disability and economic burden worldwide². Current treatments primarily involve steroidal or nonsteroidal anti-inflammatory drugs (NSAIDs) to alleviate pain and inflammation. However, NSAIDs often lead to undesirable side effects in the gastrointestinal and cardiovascular systems³. Joint replacement surgery offers an alternative, but approximately one-third of patients do not experience significant pain relief or functional improvement⁴. As our understanding of OA's molecular mechanisms improves, several prospective therapeutic targets have emerged, enabling more precise and personalized treatments⁵. Given the critical roles of human proteins in biological processes, they represent key candidates for the development of OA drugs. Nelson et al. have suggested that selecting genetically supported targets may help address the high failure rate of drugs in clinical development⁶.

Mendelian randomization (MR) is a genetic approach for inferring causal relationships between exposures and outcomes by using single nucleotide polymorphisms (SNPs) from genome-wide association studies (GWAS) as instrumental variables^7,8. Compared with observational studies, MR leverages the random assortment of genotypes at conception to reduce confounding and reverse causation⁹. Sensitivity analyses help account for polygenic pleiotropy, providing inferences that approximate causality. Due to these strengths, MR is increasingly used to identify therapeutic targets for diseases such as multiple sclerosis and rheumatoid arthritis¹⁰. However, MR-based predictions still require experimental validation, which can be facilitated through transcriptomic analysis to prioritize key disease-related genes.

Network pharmacology, an interdisciplinary approach combining systems biology and computational tools¹¹, proposes that diseases often result from disruptions in biological networks rather than single-gene mutations¹². This method enables the analysis of multi-component drugs and the prediction of their therapeutic targets at a systems level by mining drug-disease databases¹³. Molecular docking, a key component of this framework, simulates binding affinities between drugs and targets to assess druggability. Despite these advances, few studies have integrated MR with network pharmacology for OA.

This article presents a protocol that integrates MR, transcriptomic analysis, and network pharmacology to identify plasma proteins as viable therapeutic targets for OA. The aim is not only to identify causal proteins but also to prioritize candidate therapeutic compounds through computational analysis.

Protocol

We obtained ethical approval and informed consent from the Biomedical Research Ethics Committee of the First Affiliated Hospital of Nanchang University. Ethics Number: (2025)CDYFYYLK(08-007).

MR analysis

Data retrieval
The plasma pQTL data were obtained from the study by Zheng et al.¹⁴, which integrated five GWAS datasets^{15,16,17,18,19}, and from the study by Ferkingstad et al. The inclusion criteria for the data were as follows: (i) genome-wide significant associations (p < 5 × 10⁻⁸); and (ii) plasma proteins as potential therapeutic targets for OA. The study design is summarized in Figure 1. First, we identified candidate therapeutic targets using GWAS data from the IEU OpenGWAS and plasma pQTL data from the studies by Zheng¹⁴ and Ferkingstad²⁰(Supplemental Table S1 and Supplemental Table S2). Steiger filtering and phenotype scanning were then conducted to validate the robustness of the results. The IEU OpenGWAS (https://gwas.mrcieu.ac.uk/) was used to obtain summary statistics for hip or knee OA (n = 417,596), knee OA (n = 403,124), and hip OA (n = 393,873)²¹.

SNP filtering commands

SNPs with genome-wide significance (p < 5 × 10⁻⁸) were subjected to a clumping process (r² < 0.001, F-statistics > 10, window size = 10,000 kb) prior to MR analysis.

MR analysis

To investigate potential drug targets, MR analysis was performed using plasma proteins as exposures and OA as the outcome, implemented via the "TwoSampleMR" package in R (v4.3.1). When only a single pQTL was available for a protein, the Wald ratio was used; otherwise, inverse variance weighted MR (MR-IVW) was applied, followed by heterogeneity and pleiotropy assessments. Bonferroni correction was used to account for multiple testing, with a threshold of p < 5.63 × 10⁻⁵ for prioritizing proteins.

Steiger filtering and phenotype scanning

To assess reverse causality, we conducted Steiger filtering. A result of “TRUE” with p < 0.05 indicated no reverse causality. Phenotype scanning was conducted using LDtrait (https://ldlink.nih.gov/?tab=ldtrait#home-tab)²² to examine associations of pQTLs with other traits. The thresholds were R² = 0.1 and a ±500,000 base pair window. Pleiotropic effects were assigned to pQTLs meeting both of the following: (i) genome-wide significant association (p < 5 × 10⁻⁸), and (ii) association with known OA risk factors.

Phenome-wide association study

To account for gene pleiotropy and off-target effects, we conducted a phenome-wide association study (PheWAS) using the AstraZeneca PheWAS Portal (https://azphewas.com/), which contains 15,500 binary phenotypes and 1,500 continuous phenotypes from ~450,000 UK Biobank participants²³. Thresholds were set to default values to minimize false positives.

Protein - protein interaction (PPI) network

To visualize interactions among potential protein targets identified by MR, we used GeneMANIA (https://genemania.org/) for protein-protein interaction analysis and result visualization²⁴.

Enrichment analysis

To investigate biological relevance, we conducted enrichment analysis using bioinformatics tools from https://www.bioinformatics.com.cn for data analysis and visualization.

Transcriptomic workflow

Total RNA was extracted using the RNA extraction reagent kit following the manufacturer's guidelines. RNA quality was assessed using an automated RNA quality assessment system; only samples with RIN ≥7.0 were used. Quality was confirmed by RNase-free agarose gel electrophoresis (1.5% gel). Eukaryotic mRNA was enriched using Oligo(dT) beads; prokaryotic mRNA was enriched using the RNA elimination Magnetic Kit. mRNA was fragmented (200-700 nt) and converted to cDNA using the RNA Library Prep Kit. The cDNA library was end-repaired, A-tailed, ligated to adapters, purified using DNA-purifying magnetic beads (1.0×), and PCR-amplified. Sequencing was performed on a high-throughput next-generation sequencing platform. Differentially expressed genes were defined by log₂FC > 1 and adjusted p < 0.05.

Network pharmacology

To identify potential drugs for target proteins, we used BATMAN-TCM (http://bionet.ncpsb.org.cn/batman-tcm/index.php)²⁵. A score cutoff of 0.74 (LR = 32.5) was used to select known and predicted compounds. Herbal components were retrieved from TCMSP (https://old.tcmsp-e.com/index.php) and filtered with OB > 30% and DL > 0.18²⁶.

Molecular docking

Molecular docking was used to assess binding interactions. Protein structures were retrieved from the PDB (https://www.rcsb.org/). UCSF Chimera was used to preprocess structures by removing ligands and solvents. AutoDock Tools was used to calculate Gasteiger charges and define box centers and sizes. Drug structures were obtained from PubChem (https://pubchem.ncbi.nlm.nih.gov/) and preprocessed similarly. Docking was performed using AutoDock Vina. Box dimensions varied by target. Binding affinities were calculated, and results were visualized in UCSF Chimera.

Results

Screening the proteome for osteoarthritis causal proteins

Using plasma pQTL data from the study by Zheng et al., Mendelian randomization (MR) analysis identified six protein-knee or hip OA associations, one protein-knee OA association, and four protein-hip OA associations that met the Bonferroni significance threshold (P < 5.63 × 10⁻⁵) (Table 1 and Figure 2A

Discussion

This is a two-sample MR analysis method that has been employed to predict potential proteins related to osteoarthritis (OA), supplemented by sequencing validation of blood samples from OA patients and network pharmacology to predict drug targets for OA. We utilized plasma proteomics data and GWAS data for OA to conduct the two-sample MR analysis, identifying 19 proteins associated with OA, including EFEMP1, PTHLH, PCSK1, ECM1, ISLR2, SEMA3G, SPOCK2, LRIG3, TMEM190, CRYZ, KNG1, OMG, DDX19A, MAPK3, PTPN9, USP8, DDX19B, PRO...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Agilent 2100 Bioanalyzer	Agilent Technologies	G2939AA	Automated RNA quality assessment system
AutoDock Vina	The Scripps Research Institute	Vina 1.1.2	Open-source command-line software for computational molecular docking
Illumina NovaSeq 6000 platform	Illumina	SY-415-1001	High-throughput next-generation sequencing platform
NEBNext Ultra RNA Library Prep Kit	New England Biolabs	NEB #E7770	Fragment and transcribe enriched mRNA into cDNA
Ribo-Zero Magnetic Kit	Epicentre	MRZH11124	Eliminate rRNA
Trizol reagent kit	Invitrogen	AM9738	Extract total RNA
UCSF Chimera	University of California, San Francisco	Chimera-alpha-win64	Molecular visualization and structure analysis software.