The analysis flowchart for the study is presented in Figure 1. A two-phase analytical approach was employed in this study, merging P/TWAS for association and MR for causation. Additionally, to guarantee the validity and reliability of the findings, a discovery-confirmatory framework was implemented in both phases.

The GWAS summary data from the Breast Cancer Association Consortium (BCAC), which specifically focused on individuals of European descent (https://bcac.ccge.medschl.cam.ac.uk/), was utilized in our study. This dataset was comprised of 122,977 breast cancer cases and 105,974 controls. The same analytical approach was also applied to ER positive and negative breast cancer. The ER-positive subtype was found to consist of 69,501 cases and 105,974 controls, while the ER-negative subtype included 21,468 cases and 105,974 controls.

Cis-pQTL data for European Americans’ (EA) plasma proteins were obtained from the ARIC cohort (nilanjanchatterjeelab.org/pwas/), generated using PLINK2 software (Purcell et al., 2007). The SeqID file names correspond to the SOMAmers (Slow Off-rate Modified Aptamers), which are utilized for measuring protein levels in biological samples by leveraging their enhanced affinity and specificity for target proteins (Rohloff et al., 2014). For external validation, cis-pQTL data from 4 extensive plasma protein cohorts of European descent were used. Additionally, our study also explored expression quantitative trait loci (eQTLs), which influence gene expression at the transcriptome level (Zhu et al., 2016). We extracted eQTL data using the SMR toolkit (Wu et al., 2021), a tool specifically designed for genetic epidemiological research, from two major sources: the Genotype-Tissue Expression Project (GTEx Consortium, 2020) and the eQTLGen consortium (Võsa et al., 2021). Detailed descriptions of each dataset are provided in Table 1.

FUSION (Boston, MA, United States) (Gusev et al., 2016), which is a software to establish associations between functional phenotype and GWAS phenotype, was used to conduct P/TWAS analysis. In our study, FUSION was implemented to identify associations between protein/gene expression levels and Breast Cancer susceptibility. Methodologically, FUSION takes two inputs: 1) Precomputed functional weights, and 2) GWAS summary statistics unified to a reference SNP panel. In PWAS, precomputed functional weights of plasma proteins were obtained from the ARIC study (Zhang et al., 2022b), and the reference SNP panel was derived from the European descent of the 1000G project (http://www.internationalgenome.org/faq/how-do-i-cite-1000-genomes-project). The primary outputs of FUSION are the Z-score and p-value, wherein Z-score quantifies the strength and direction of the associations between plasma proteins and breast cancer, while the p-value elucidates the statistical significance of this association. To enhance our findings at transcriptomic level, we incorporated TWAS for both whole blood and breast mammary tissues. The precomputed functional weights for TWAS, provided by Junghyun Jung from the Mancuso lab (http://gusevlab.org/projects/fusion/). A false discovery rate (FDR, Benjamini–Hochberg) threshold of 0.05 was applied to determine the statistical significance of the results.

Bayesian colocalization analysis (Giambartolomei et al., 2014) was utilized to evaluate the probability that the same genetic variant affects both plasma protein and breast cancer. The default parameters set by the analysis were followed, including p1 = 10e−4 (the probability of a variant being a significant pQTL), p2 = 10e−4 (the probability of a variant associated with breast cancer), and p12 = 10e−5 (the probability of a variant being significant in both protein/gene and GWAS). This analysis involved five predefined hypotheses: H0, indicating no association with either trait; H1, signifying association with trait1 only; H2, implying association with trait2 only; H3, representing associations with both traits due to different SNPs; and H4, indicating association with both traits due to a common SNP. A posterior probability of H4 (PPH4) exceeding 0.8, or in some cases 0.7, is generally interpreted as strong evidence of the same genetic variant being implicated in both traits (Giambartolomei et al., 2014).

Recent studies have investigated the causal associations between plasma proteins and diseases like colorectal cancer (Sun et al., 2023) and inflammatory bowel disease (Chen et al., 2023), utilizing a scoring system that integrates p-value and PPH4. Building on this approach, our research employs P/TWAS and Bayesian Colocalization analysis to systematically categorize the degrees of association between proteins. The scoring system was as follows: a significant adjusted p-value was awarded 1 point, and a PPH4 > 0.75 also earned 1 point. Based on the cumulative scores, associations were categorized as follows: a score between 1 and 2 indicated a “Weak” association, 3 to 4 suggested a “Moderate” association, and 5 to 6 signified a “Strong” association.

In the causal analysis, we primarily conducted further analysis on proteins with strong and moderate associations. MR analysis were based on 3 essential assumptions for genetic instrumental variables: relevance, independence, and exclusion-restriction (Davies et al., 2018). We implemented a stringent selection process for SNPs to be used as instrumental variables, requiring a p < 5e-8, or p < 5e-6 in cases when SNP was absent. Clump was applied in accordance with the default parameters. The Wald Ratio (WR) method was employed when a single SNP was used as the instrumental variable, whereas the inverse-variance weighted (IVW) method was predominant when the instrumental variables involved multiple SNPs (Burgess et al., 2019). To reinforce the robustness of our findings, we conducted several sensitivity analyses. The Steiger filtering test (Deng et al., 2022) was utilized to eliminate the possibility of reverse causal associations. Additionally, heterogeneity and pleiotropy sensitivity analyses were conducted for proteins that met the criteria (Bowden et al., 2015; Greco et al., 2015). Furthermore, to improve the reliability and applicability of our results, external validation was carried out on pQTL data derived from 4 extensive plasma protein cohorts in European populations.

To delve deeper into the intricate relationships and biological functions of significant proteins identified in our PWAS, gene ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses was performed. Given the emerging role of plasma proteins as potential therapeutic targets (Sun et al., 2018), we matched P/TWAS-MR significant proteins with the druggable genome database (Finan et al., 2017), which categorizes 4,479 genes into three druggability tiers: Tier 1 includes approved drugs and candidates in clinical trials, Tier 2 encompasses targets of biologically active molecules and those similar to approved drug targets, and Tier 3 comprises genes for secreted or extracellular proteins and other key druggable gene family members. Additionally, the significant proteins were annotated using the Therapeutic Target Database (http://db.idrblab.net/ttd/) (Zhou et al., 2022).

In this study, data analysis was executed using R software (version 4.3.1). The P/TWAS analysis followed the analytical process previously described. The Benjamini–Hochberg method was employed for multiple testing correction, with adjusted p-values <0.05 considered statistically significant. Causations were investigated using the “TwoSampleMR” package, while Bayesian colocalization analysis was carried out using the “COLOC” package. The “ClusterProfiler” package (Wu et al., 2021) was utilized for functional enrichment analysis. Data visualization was achieved through the “Forestploter” and “ggplot2” packages, and data cleaning was performed using the “tidyverse” package.

In our study, a total of 25 plasma proteins were significantly associated with breast cancer (Table 2; Figure 2A, and Supplementary Table S1). Of these proteins, 14 showed a Z-score greater than 0, denoting a positive association with breast cancer. Conversely, the remaining 11 proteins suggested an inverse association with the disease. When duplicate SOMAmers are present, we select the protein corresponding to the smallest p-value for subsequent analysis, such as RSPO3 (Supplementary Table S1). In ER subtypes analysis, 16 proteins were found to be significantly associated with ER-positive breast cancer and 6 with ER-negative breast cancer (Supplementary Table S2, S3). The PWAS Manhattan plot illustrates the distribution of significant genes across different chromosomes and their respective p-value (Figure 3A, Supplementary Figure S2A, B).

For the 25 proteins identified by PWAS, 12 showed significant associations in the whole blood TWAS analysis (P.adj < 0.05) (Supplementary Table S4). While, in the breast mammary tissue TWAS, 10 of these proteins were further validated (P.adj < 0.05) (Table 2, Supplementary Table S5). Among the 16 significant proteins in ER-positive breast cancer, with 7 were confirmed in both whole blood and breast mammary tissue analyses. Meanwhile, in ER-negative breast cancer, 2 out of the 6 significant proteins were validated (Figures 2B,C, Supplementary Table S6). It should be noted that MST1 exhibited contradictory associations in PWAS (Z = 4.194, P.adj = 0.004) and TWAS (Z = −2.547, P.adj = 0.014). This pattern was also observed in ER-positive and ER-negative subtypes. Due to the complex nature and potential biological implications of MST1’s contrasting results, we did not conduct further analysis on this protein.

Among 25 significant proteins, 9 exhibited strong genetic colocalization evidence. Additionally, 4 proteins–SNUPN (PPH4 = 93.9%), CSK (PPH4 = 84.3%), CTSF (PPH4 = 94.4%), and PARK7 (PPH4 = 96.5%)–also demonstrated the same strong genetic evidence at the whole blood transcriptomic level. Remarkably, CSK (PPH4 = 86.3%) and CTSF (PPH4 = 93.9%) were further validated in the breast mammary tissue transcriptomic level (Table 2). In ER-positive breast cancer, 5 proteins showed strong evidence of genetic colocalization. Notably, 2 of these proteins, CSK (PPH4 = 85.8%, 86.5%) and GDI2 (PPH4 = 97%, 97.2%), demonstrated the same strong genetic colocalization evidence in both whole blood and breast mammary tissues. In the ER-negative breast cancer, PEX14 showed strong genetic colocalization evidence in protein (PPH4 = 99.9%) and breast mammary tissue (PPH4 = 88.8%), but this pattern was not replicated at the whole blood transcriptomic level (PP4 = 8%, Supplementary Table S6).

In breast cancer, 25 proteins were classified: 5 as “Strong” association (red), 6 as “Moderate” association (blue), and the remaining as “Weak” association (grey) (Figure 2A). Among the “Strong” Tiers, PEX14 (Z = 4.839) and CTSF (Z = 3.681) had a positive association with breast cancer. Whereas, SNUPN (Z = −4.413), CSK (Z=−4.417), and PARK7 (Z = −3.648) showed negative associations (Table 2).

In ER-positive breast cancer, 16 proteins were classified: 3 proteins showed “Strong” association (red), with GDI2 (Z = −3.652) newly identified and negatively associated with ER-positive breast cancer. Additionally, CSK and CTSF followed the same trends with the findings from breast cancer. Besides, 4 proteins were “Moderate” (blue), and 9 proteins were “Weak” associations (grey) (Figure 2B). In ER-negative breast cancer, 6 proteins were classified: PEX14 and MST1 showed “Strong” associations. Notably, PEX14 not only showed the same trend as observed in breast cancer (Z = 4.839, p = 0.0004) but also exhibited a notably stronger effect (Z = 5.929, p = 2.02E-6). MST1 was not further analyzed due to inconsistent trends in P/TWAS. The other 4 proteins were categorized as “Weak” association (grey) (Figure 2C, Supplementary Table S6).

It is crucial to highlight that, although PGD and TLR1 were significant across all three outcomes in PWAS analyses (Figure 3A and Supplementary Figure S2A, B), their absence from the corresponding TWAS analysis relegated them to the “Weak” association. Moreover, the results of these two proteins were not sufficiently reliable in MR Analysis (Supplementary Figure S2C, Supplementary Table S7).

Upon determining the strength of associations, we supplemented the causations with MR analysis (Supplementary Table S8). We primarily focused on the causal effects of “Strong” associated proteins. Among the 5 “Strong” associated proteins, PEX14 was found to have a positive causation at the proteomic (OR = 1.201, p = 0.017) and transcriptomic level (OR = 1.17, p < 0.001). CTSF demonstrated a positive causation in the ARIC cohort (OR = 1.114, p < 0.001), and the consistent trends were also external validated in INTERVAL cohort (OR = 1.144, p < 0.001) (Sun et al., 2018) and AGES-Reykjavik cohort (OR = 1.159, p < 0.001) (Gudjonsson et al., 2022) (Figure 3B). The remaining 3 proteins, SNUPN (OR = 0.905, p < 0.001), CSK (OR = 0.962, p = 0.038), and PARK7 (OR = 0.954, p < 0.001), all exhibited negative causations with breast cancer. External validations from the deCODE cohort further confirmed the causations for SNUPN (OR = 0.797, p < 0.001) and PARK7 (OR = 0.844, p = 0.017). However, CSK’s causation at the whole blood transcriptomic level was somewhat unsignificant (OR = 0.84, p = 0.129) (Figure 3B, Supplementary Table S9).

In ER-positive breast cancer, CSK (OR = 0.955, p = 0.038) and CTSF (OR = 1.125, p < 0.001) maintained the same causal trends as observed in breast cancer (Supplementary Table S10). Additionally, GDI2 was identified as a newly negatively significant protein (OR = 0.92, p < 0.001). However, its causal effect was not significant at the transcriptomic level (OR = 1.001, p < 0.981, Figure 4A). In ER-negative breast cancer, PEX14 stood out as the sole “Strong” protein. Notably, its causal effect in this subtype (OR = 1.645, p < 0.001, Figure 4B) was further pronounced compared to breast cancer (OR = 1.201, p = 0.017). Meanwhile, we expanded our MR analyses to include “Moderate” proteins. The results revealed that their causal effects were generally less consistent and of reduced significance compared to those of the “Strong” proteins (Supplementary Figure S1, Supplementary Table S11).

Considering that the pQTLs of most plasma proteins was a single SNP, conducting sensitivity analyses for heterogeneity and pleiotropy is typically not required. As result, in ER-positive breast cancer, BTN3A3, EMILIN3, FOLR3, and NTN4 showed heterogeneity, while in ER-negative cases, this was not observed (Supplementary Table S10). BTN3A3 in ER-positive breast cancer also displayed pleiotropy. The Steiger filtering test confirmed that MR effects were due to plasma proteins affecting breast cancer outcomes (Supplementary Table S8, S10). Importantly, our “Strong” proteins exhibited neither heterogeneity nor pleiotropy.

Furthermore, the plasma proteins identified by PWAS were subjected to Gene Ontology (GO) cluster analysis. This analysis revealed a predominant association with biological processes related to oxidative stress, such as “reactive oxygen species metabolic” and “response to reactive oxygen species” terms. Additionally, for cellular components, we observed a significant enrichment in the “collagen-containing extracellular matrix” term (Supplementary Figure S3). Besides, KEGG pathway enrichment did not reveal any significantly enriched pathways (Supplementary Table S12).

Plasma proteins are not only crucial as diagnostic biomarkers but also serve as potential drug targets. In our study, we evaluated the significant proteins for their potential as drug targets. By aligning our findings with the druggable genome database (Finan et al., 2017), we determined that 16 of the 25 proteins have druggable targets. These include 3 proteins in Tier 1; 3 in Tier 2, and 10 in Tier 3 (Supplementary Table S13, Left column). Furthermore, we compared our results with the Therapeutic Target Database (Zhou et al., 2022), 11 of these 16 proteins were identified as targets of existing or potential drugs. This group comprised 3 Successful targets, 3 Patented-recorded Targets, 1 in clinical trials, and 4 documented in literature (Supplementary Table S13 Right column). Among the “Strong” proteins, CSK and CTSF were found to be drug targets with patent records, categorized under Tiers 1 and 2 respectively. CTSF has been documented to be used in the treatment of bone cancer and chronic obstructive pulmonary disease (Li et al., 2017) (Table 3). However, the remaining “Strong” proteins have not yet been reported.