Integrated Single-Cell RNA Sequencing and Proteome-Wide Mendelian Randomization Identifies Therapeutic Targets for Psoriasis and Associated Complications

Introduction

Psoriasis is a chronic inflammatory autoimmune disease that primarily affects the skin tissues.¹ In individuals without psoriasis, skin cells are produced and shed every three to four weeks. However, in psoriasis patients, this process is significantly sped up, and skin cells can cycle every three to seven days. This rapid turnover leads to the buildup of skin cells on the surface, forming the characteristic thick, red, scaly patches known as plaques.² Generally, it is estimated that psoriasis affects about 2–3% of the world’s population, although these figures can differ between regions and populations.³ In the United States, for example, psoriasis affects about 2–2.6% of the population, according to the National Psoriasis Foundation.⁴ It is less common in African and Asian populations and more common in Western countries. In some northern European and Scandinavian populations, the prevalence can be higher, reaching up to 8%.^5,6

The exact cause of psoriasis is not fully understood, but it is believed to involve the immune system, genetics, and environmental factors. Treatments range from topical ointments and creams for mild to moderate cases, to phototherapy (light therapy), and systemic medications, including biologic drugs for more severe cases. The condition can also have a psychological impact, as visible skin conditions often do, potentially leading to decreased self-esteem, depression, and social isolation.^7,8 Newer biologic drugs target-specific pathways and molecules associated with the inflammatory process that drives psoriasis, such as interleukins (IL) and tumor necrosis factors (TNF), including IL-23 inhibitor, IL-17 inhibitor and Phosphodiesterase 4 (PDE4) inhibitor. However, it is important to note that less than half of the patients see a significant improvement in their condition when using inhibitors of TNF-a or IL-17A.^9–11 This highlights the fact that, since psoriasis is currently incurable, most patients need continuous treatment with medications that may have side effects to manage their symptoms. Although IL23 and IL17 are traditional markers for psoriasis, however, these traditional markers only focus on one cellular inflammatory pathway, and not meet the unmet needs of personalized therapy strategy. More importantly, some patients respond poor and have side effects to these biologics. Therefore, it is critical to find more novel markers of psoriasis, to prevent disease progression, reduce mortality, and improve the quality of life for those patients.

Currently, mendelian randomization (MR) study has gained traction as a powerful method for linking existing drugs and identifying new targets for therapy.¹² Genome-wide association studies (GWAS) have pinpointed specific single nucleotide polymorphisms (SNPs) situated on chromosomes that are instrumental in controlling protein expression. These SNPs correspond to the levels of protein expression and are known as protein quantitative trait loci (pQTLs).¹³ MR leverages these pQTLs as instrumental variables, allowing researchers to probe into the potential cause-and-effect links between exposures and outcomes, which is invaluable for identifying new drug targets and biomarkers.^13,14 Combining the single cell RNA-seq with MR analysis can illustrate the cellular origin of different eQTLs or pQTLs, and related functional pathways.

MR offers a more reliable assessment of causality compared to traditional observational studies by mitigating the influence of confounding variables.¹⁵ Furthermore, the integration of phenome-wide association studies (PheWAS) enables the prediction of potential adverse reactions related to these targets.¹⁶ Plasma proteins, which are vital in a myriad of biological functions including signaling, transport, growth, repair, and immune responses, often become dysregulated in various diseases. This dysregulation makes them prime candidates for drug development. In pursuit of new treatments for psoriasis, we embarked on a comprehensive proteome-wide MR study to identify promising therapeutic targets.¹⁷

Our approach began with a two-sample MR analysis to assess the causal influence of plasma proteins on psoriasis. We then carried out colocalization analyses to confirm the credibility of our results. Finally, we assessed the potential side effects of the identified proteins that could be targeted for psoriasis treatment by utilizing PheWAS.

Methods

Approval of Research Ethics and Study Blueprint

Our study harnessed comprehensive GWAS summary statistics from original investigations, which were conducted with participants’ informed consent. The ethics requirement has been approved by the local ethical committee in Shanghai Pudong New area People’s Hospital.

Characterization of Plasma Protein QTLs

We extracted plasma protein quantitative trait loci (pQTLs) data from the UK Biobank-PPP, which encompasses proteomic data from 54,219 individuals (https://www.synapse.org/#!Synapse:syn51365303). A meticulous mapping was performed for 2,940 proteins.¹⁸ pQTLs are typically situated near their respective gene, known as “cis-pQTLs”, suggesting a regulatory effect via the proximate gene.¹⁹ In contrast, “trans-pQTLs” are located further away or on different chromosomes, presumably exerting effects through other genes. The distinction between cis and trans pQTLs often hinges on set distance thresholds, commonly 500 kb or 1,000 kb.^20–22

For our proteome-wide MR analysis targeting drug discovery, we adopted pQTLs as instrumental variables with strict inclusion criteria: 1) SNP proximity within ±1 Mb of the gene region (cis-pQTLs); 2) a genome-wide significance level of P<5x10⁻⁸ for SNPs strongly associated with plasma proteins; 3) an LD threshold of r² < 0.001 and a genetic distance of 10,000 kb to ensure SNP independence; 4) an F-value over 10 to minimize weak instrument bias.¹²

Genetic Association Data for Psoriasis

Genetic data for psoriasis was sourced from the Finnish database’s R10 release (https://r10.finngen.fi/). The dataset includes 6408 patients and 397564 controls.²³

Conducting MR Analysis

We implemented a two-sample MR analysis, considering plasma proteins as the exposure and psoriasis as the outcome. The selection of pQTLs adhered to previously specified criteria. We utilized the “TwoSampleMR” R package (version 0.6.0) for analysis, applying the Wald ratio for single SNP scenarios and the inverse variance-weighted (IVW) method for multiple SNPs.¹² To account for multiple testing, we used the False Discovery Rate (FDR) method to adjust P-values, with a Pfdr < 0.05 denoting statistical significance. The SNPs of druggable genes are obtained from the literature review.²⁴

Analysis of Colocalization

Colocalization analysis seeks to ascertain shared genetic variants influencing both exposure and outcome, further corroborating MR findings. For proteins with positive MR outcomes, we scrutinized SNPs within a ±1MB window of the gene (cis-pQTLs) for their colocalization with psoriasis.²⁵ Five hypotheses are tested in colocalization analysis, determining whether SNPs are related to the protein, the disease, both, or neither. We concentrated on genes with a combined posterior probability (PPH3+PPH4) ≥ 0.8 due to the limited power in colocalization analyses.²⁶

Phenome-Wide Association Analysis

PheWAS, the inverse of GWAS, is employed to discern associations between SNPs or phenotypes and a broad spectrum of phenotypic traits.¹⁶ This methodology is especially useful for identifying potential side effects of drug targets.²⁶ In our study, we focused on plasma proteins with positive MR results as the exposure, maintaining the same criteria for instrumental variables. The outcome entailed phenotypic data from the Finnish database R10 version, which includes 2,272 phenotypes across 46 categories, used for a comprehensive phenome-wide MR analysis. A Pfdr <0.05 was indicative of statistical significance.

GEO Data Download

We downloaded the datasets (GSE54456, GSE121212) from NCBI website. GSE54456 has 174 blood samples including 92 psoriasis patients and 82 healthy controls.²⁷ The GSE121212 dataset has 28 patients with skin lesion, and 38 healthy controls.²⁸

Identifying Clinical Relevance of Key Genes

To determine the potential impact of key genes in clinical settings, we utilized the “pROC” and “xgboost” software packages to develop a machine learning model that assesses the predictive value of our diagnostic model. The outcomes of this analysis were presented in the form of ROC (Receiver Operating Characteristic) and PRA (Precision-Recall Analysis) curves. XGBoost, an acronym for eXtreme Gradient Boosting, is a sophisticated version of the gradient boosting technique.²⁹ It has been notably effective in various supervised learning scenarios, often surpassing other machine learning approaches in competitions and practical applications. Its effectiveness stems from both its ability to produce accurate and reliable predictive models and its efficiency in doing so.³⁰

The ROC curve is a visual tool used to assess the effectiveness of a classification model by plotting the True Positive Rate (TPR, also referred to as recall or sensitivity) against the False-Positive Rate (FPR) at different threshold settings. This curve can sometimes give a misleading impression of a model’s performance, especially when there is a significant imbalance in the number of positive and negative cases.³¹ The PR curve, on the other hand, illustrates Precision (the proportion of true positives among all positive predictions) against Recall for various threshold levels. This curve is particularly useful for evaluating models in situations where there is a class imbalance, as it focuses on the model’s ability to correctly predict the minority class.

GSVA Analyses

Utilizing the “c2.cp.all.v7.0. symbols” gene set, we employed the GSVA R package to calculate the metabolic scores for relevant pathways based on the gene expression matrix of each sample, using the ssGSEA method.³² The pathways that showed significant differences between the two clusters were displayed in a heatmap, with statistical significance defined as adj.P.val<0.05. Subsequently, we used Pearson analysis to identify the relationships between these significant metabolic pathways and risk genes.

Immune Infiltration Analysis

The MCP counter employs linear support vector regression to deconvolute the transcriptome expression matrix, thereby estimating the composition and abundance of immune cells in a mixed cell population. We uploaded the gene expression matrix data to IOBR, excluded samples with p > 0.05, and derived an immune cell infiltration matrix. We then used Pearson analysis to explore the connections between significant immune cells and risk genes.

Single Cell Sequencing Data Obtained and Processing

The single-cell data were sourced from the NCBI website under the accession number GEO230842. We performed data quality control, focusing on cells with less than 10% mitochondrial genes, and selected cells with a gene count ranging from 200 to 10,000 that were expressed in at least three cells. The number of highly variable genes was set at 2000, aligning with previous research.^33–35 The six samples were integrated using SCT correction, followed by dimensionality reduction through the uMAP method. Single cell-RNA sequencing was utilized to identify key genes and determine their cellular origins.³⁶ The Cellchat R package (version 1.4.0) was employed to investigate interactions between different cell clusters. We compared the extent of glycolysis in each specific cell using the addmodulescore method.

Results

The schematic plot of the project is shown in Figure 1.

Figure 1 Study design overview.

Plasma Proteins Related to Psoriasis

Upon strict adherence to the instrumental variables screening criteria in this study, eventually a total of 1280 plasma proteins were included in the MR analysis. It is worth noting that, among the subset of 1280 plasma proteins, the MR analysis, based on IVW or Wald ratio results (Pfdr<0.05), revealed positive causal associations with Psoriasis for six plasma proteins—APOF (Apolipoprotein F), ATP6V1G2 (V-type proton ATPase subunit G 2), GRELD1 (Protein disulfide isomerase CRELD1), IFNLR1 (Interferon lambda receptor 1), PRSS8 (Prostasin 88), and TNF (Tumor necrosis factor). The odds ratios and 95% confidence intervals for these associations are 1.69 (95% CI:1.40–2.03), 5.98 (95% CI:2.70–13.25), 1.14 (95% CI:1.08–1.21), 1.21 (95% CI: 1.12–1.31), 1.74 (95% CI:1.35–2.24) and 2.14 (95% CI:1.59–2.89) respectively. For comprehensive details, please refer to Figures 2 and 3.

Figure 2 Volcano plot of MR results: Causal relationship between plasma proteins: APOF, TNF, FNLR1, CRELD1, PRSS8, ATP6V1G2 and psoriasis onset with OR>0.

Figure 3 Forest plot of the MR results: Effects of 6 plasma proteins on psoriasis. Abbreviations: CI, confidence interval; OR, odds ratio.

Sensitivity Analysis for Plasma Proteins Associated with Psoriasis

For these 6 plasma proteins, we conducted gene colocalization analysis within a range of ±1MB upstream and downstream of their respective genes to explore potential associations with psoriasis. The results indicate that all these potential plasma proteins may share a causal variant in this region (PPH3+PPH4>0.8). For detailed information, please refer to Figure S1. This suggests that all these 6 plasma proteins may serve as potential targets for treating psoriasis.

Phenome-Wide Associations Analysis for 6 Plasma Proteins Linked to Psoriasis

To evaluate the potential beneficial or deleterious effects of the 6 plasma proteins associated with psoriasis on other phenotypes, we conducted a phenome-wide association analysis, screening 2408 phenotypes across 45 categories from the Finnish database (version R10). We observed significant causal associations between APOF and 6 phenotypes (Pfdr<0.05), ATP6V1G2 and 15 phenotypes (Pfdr<0.05), CRELD1 and 5 phenotypes (Pfdr<0.05), IFNLR1 and 8 phenotypes (Pfdr<0.05), PRSS8 and 17 phenotypes (Pfdr<0.05), and TNF and 95 phenotypes (Pfdr<0.05). For detailed data, please see Figure 4. These phenotypes may be therapeutic objects or deleterious effects for the target protein.

Figure 4 Manhattan plot of result of PheWAS analysis of associations between 6 plasma proteins and other disease outcomes.

Associations Between 6 Plasma Proteins Linked to Psoriasis Complication

As psoriasis patients can develop several complication, we further analyzed the relationship between pQTL of each protein and the related complications, and found all six proteins were related to psoriatic vulgaris, and psoriatic arthropathies except PRSS8. However, only PRSS8 is correlated with pustular psoriasis and only APOF and ATP6V1G2 are associated with Guttate psoriasis (Figure 5). These results suggest that as for the choice of potential drug inhibitor, the complication relation of each protein should be taken into consideration.

Figure 5 Heatmap result of correlational analysis of between 6 psoriasis-associated plasma proteins and disease complications. The * indicates the correlation has significant meanings (p<0.05).

Validate the Role of Plasma Proteins in Psoriasis

We then used the GEO dataset to validate the role of these protein coding genes in psoriasis, we used five genes (except Apof, as it is not found in the chipset data) to construct the ROC and PRC curve in the training dataset (GSE 54456) which included 92 psoriatic skin and 82 normal human skin tissue and testing dataset (GSE 121212) which includes 28 lesion skin and 38 healthy skin, showing that they have a good diagnostic effect in the training data with ROC at 0.989 and PRC at 0.991; while, with ROC at 0.916 and PRC at 0.83 in the testing dataset (Figure S2).

The Immune Infiltration and Metabolic Pathway in Hub Genes

Then we found both PRSS8 and TNF are positively related to several immune cells: T cells follicular helper, CD4 memory T cells and CD8 T cells. Regarding the metabolic change, we found both PRSS8 and TNF are associated with Glycolysis, Oxidative phosphorylation (Figure S3).

We finally applied the single-cell RNA-sequencing data to validate our results. First, we saw that the glycolysis pathway is more enriched in the psoriasis group compared to healthy controls (Figure S4A). Then, both TNF and PRSS8 are positively associated with the glycolysis in the keratinocytes (Figure S4B and C).

The Druggable Gene for Psoriasis

Applying MR method, we found 341 druggable gene in psoriasis. We found most genes are related to the immune system, sphingolipid pathway and cholesterol metabolism (Figure S5A). Then, we investigated the targets and drugs and found most targets are related to the clinical medicine. Among then, simvastatin is related to TNF (Figure S5B).

We next applied the molecular docking for both simvastatin and TNF-alpha and found the binding score is around −10.0, which indicated that they have a strong binding (Figure S6).

Discussion

Identifying new therapeutic targets is therefore crucial for enhancing treatment outcomes. In this study, the plasma proteins that were positively correlated with psoriasis were APOF, ATP6V1G2, GRELD1, IFNLR1, PRSS8, and TNF. TNF is a vital cytokine that plays a central role in the human immune system and is implicated in the development of various inflammatory and autoimmune diseases, including psoriasis.^9–11 The role of TNF in psoriasis has been extensively researched, and medications that inhibit TNF, such as adalimumab, golimumab, infliximab, certolizumab pegol, and etanercept, have demonstrated considerable efficacy in treating psoriasis. This reinforces the validity of our findings that TNF is a reliable target for psoriasis treatment. However, the use of TNF inhibitors is not without challenges, such as varying effectiveness, cost concerns, and potential adverse effects.^9–11 The PheWAS findings from this study seem to support the possibility of more serious adverse effects, including thyroiditis, follicular lymphoma, and skin malignancies. Hence, it is essential to be vigilant about these potential adverse reactions when using TNF inhibitors and to conduct regular assessments to ensure timely detection and treatment. It is also important to note that observational studies have not generally indicated an increased risk of malignancy with the use of tumor necrosis factor inhibitors.³⁷ Most importantly, we showed that TNF can combine to simvastatin, and TNF is an important target for druggable gene in psoriasis, which indicated that simvastatin might be useful to treat psoriasis. Clinically, TNF inhibitors are effective for many immune-mediated diseases but carry specific risks, most notably increased risk of serious infections, infusion/injection reactions, and rare immune³⁸ or neurological complications;^39,40 While, statin application has very few side effects of these.

Apolipoprotein F (APOF) is a component of high-density lipoprotein (HDL) particles and has been studied for its role in lipid metabolism.⁴¹ It is known to modulate the function of lecithin-cholesterol acyltransferase (LCAT), an enzyme that is important for the maturation of HDL particles.⁴¹ Dyslipidemia, characterized by abnormal lipid levels, is a common comorbidity in patients with psoriasis and is part of the metabolic syndrome that is more prevalent in these patients. The inflammatory state of psoriasis is believed to contribute to the development of dyslipidemia, and conversely, dyslipidemia may exacerbate the inflammatory processes that drive psoriatic lesions.¹

ATP6V1G2 is a subunit of the vacuolar ATPase (V-ATPase) proton pump, which is involved in acidifying a variety of intracellular compartments in eukaryotic cells.⁴² V-ATPase activity can be implicated in immune system function and inflammation, which are both key components in the pathogenesis of psoriasis. For example, the acidification of intracellular compartments is essential for the processing and presentation of antigens by immune cells, and dysregulation of this process could potentially affect immune responses.⁴³

Cysteine-rich with EGF-like domains 1 (CRELD1) is a protein that in humans is encoded by the CRELD1 gene. The role of specific proteins like CRELD1 in psoriasis would likely relate to the pathways that regulate inflammation and cell proliferation.

In the context of psoriasis, which is a chronic inflammatory skin condition characterized by dysregulated immune responses and keratinocyte hyperproliferation, IFNLR1 might be involved in antiviral response and immune regulation. In psoriasis, this pathway could influence the disease course by affecting how the skin’s immune system responds to viral infections, which can sometimes trigger or exacerbate psoriatic lesions.

PRSS8, which stands for Serine Protease 8, also known as Prostasin, is a serine protease that is found on the surface of epithelial cells.⁴⁴ In the context of skin biology, PRSS8 is involved in the regulation of the epidermal barrier function and is considered important for maintaining skin homeostasis.⁴⁵

We further used the GEO dataset combing single-cell RNA sequencing and bulk RNA sequencing to validate our results. We showed that all hub genes coding these plasma proteins (except APOF) have a good prognostic effect in psoriasis with both ROC value more than 0.90. Then, we showed that the hub gene PRSS8 and TNF are associated with the immune cell infiltration and glycolysis level. This was also further validated in the single-cell RNA data. These indicate that both PRSS8 and TNF could be potential targets for treating psoriasis.

Study Design and Methodological Strengths: The study is underpinned by a robust Mendelian Randomization (MR) design, which serves to mitigate the influence of confounding factors and bias stemming from reverse causality. This methodology is further enhanced by the incorporation of cis-acting protein Quantitative Trait Loci (cis-pQTLs), which offer a higher level of evidence compared to trans-acting pQTLs and expression Quantitative Trait Loci (eQTLs). Gene colocalization analyses were employed to augment the statistical robustness of the findings, thereby bolstering the credibility of the study’s conclusions. Additionally, a comprehensive phenotypic association analysis was conducted to delve into the potential side effects of drug targets with greater thoroughness. We advocate for the adoption of the Phenome-Wide Association Study (PheWAS) approach by the scientific community to further elucidate the potential side effects of drug targets, thereby enriching the collective knowledge in this domain.

Limitations of the Study

The study’s generalizability may be constrained by the homogeneity of the GWAS participants, all of whom were of European descent. Despite the UK Biobank Pharma Proteomics Project (UKBPPP) database containing data on 2,940 plasma proteins, only 1,908 were included in the MR analysis due to the constraints of the available instrumental variables. The power of the colocalization analyses was limited, necessitating a focus on cases where the joint posterior association probability (PPH3+PPH4) exceeded 0.8. The feasibility of incorporating cellular and animal experimental data was limited by practical conditions; however, we plan to consider integrating such experiments into our research agenda or collaborating with other laboratories to advance the field. While the study utilized an MR approach to screen for causal associations between plasma proteins and psoriasis, it did not apply a genetics-led drug target prioritization method (Priority index, PI) to prioritize less-explored targets. However, molecular docking predictions in therapeutic validation are limited due to the rigid docking methods; however, sometimes they are flexible docking.

Conclusion

Our research has delineated the causal associations between six plasma proteins—APOF, ATP6V1G2, IFNLR1, CRELD1, PRSS8, and TNF—and psoriasis, thereby identifying novel therapeutic targets for this condition. We anticipate further exploration of these drug targets, their potential therapeutic applications, and their implications for the treatment of dermatological disorders in subsequent studies.