2.2. Proteograph^TM-Enabled Proteomic Profiling of Serum Specimens

Serum from the study populations was processed using Seer Proteograph^TM, and resulting peptides from the 5-nanoparticle enrichment were analyzed through label-free dia-PASEF MS (Figure 1). No samples were removed upon analysis for outlier protein group identifications, which were defined for this study as <25% of the average number of proteins identified across all samples. At the time of writing, this study represents one of the largest MS proteomic discovery studies in serum. Across all samples, 5468 protein groups were identified in at least one subject from 24,655 total peptides, and 3628 of those protein groups were detected in at least 50% of the subjects, corresponding to 13,122 peptides (Figure 2A). The median number of protein groups detected per subject across any of the five nanoparticles was 3601, with the nanoparticle-specific enrichments ranging from 1456 to 2784 median protein groups identified (Figure 2B). Correspondingly, the median peptides per sample was 13,597 (Figure 2C). Individual nanoparticles varied from 5458 to as high as 9819 median peptides identified across all specimens. Subjects showed similar numbers of protein or peptides detected regardless of cancer status (Figure 2A–C). The percent coefficient of variation (%CV) was determined for the unnormalized MaxLFQ intensities [39] for each protein group per nanoparticle across each of the subjects separated by cancer status. The distribution of %CV, composed of both biological and technical sources of variability, is displayed as violin plot–boxplots in Figure 2D. The median protein intensity %CV was <69% for each of the nanoparticles. Peptide level %CV was equivalent to that observed for protein intensities (Figure S3). The technical variability of the study, which includes sources due to sample preparation and MS acquisition, was estimated using an on-deck process control present on each plate of 15 samples across the entirety of the study. The results from this process control showed that the median estimated technical variability was ~35% (Figure 2E). Protein groups identified in these subjects were mapped back to a deep plasma proteome dataset from Keshishian et al. [40] and showed broad-based coverage across the dynamic range of the blood proteome, with many of these detected proteins annotated in the Open Targets database for prostate cancer (Figure 2F) [41].

Open in a separate window

Figure 2

Proteomic characterization of subjects from serum prepared through Proteograph^TM. (A) The number of protein groups and peptides identified at a given detection percentage across all samples, including subjects with cancer and without cancer. Proteins and peptides identified for the entire study and at least 50% of subjects are indicated on the plot. Distribution of the number of protein groups (B) and peptides (C) identified for all subjects, including those with cancer and without cancer, from any number of NPs (left) or per NP (right). Inner boxplots show the 25, 50, and 75th percentile of data, with whiskers representing ±1.5 × interquartile range (IQR). Median protein groups or peptides are indicated above each plot. Violins encompass all observations. (D) Distribution of percent coefficient of variation (%CV) calculated from MaxLFQ values for each protein group per NP across cancer and no cancer subjects. Missing values are not included in the calculation. Boxplots display the same information as above except outliers, defined as equal to or beyond the IQR, which are not shown. Outer violins only capture data points up to the limits shown on the axis. (E) %CV for proteins enriched for the PC3 technical control across all batches of samples per NP. Parameters for the plot are the same as previously described. (F) Protein groups observed in any NP across all subjects are matched to plasma proteins characterized in Keshishian et al. [40] and their respective relative intensity and sorted rank based on published concentration. The median rank for these matched identifications is displayed by the vertical dashed line. Identified proteins matching those in the Open Targets database for prostate cancer with a globalScore of ≥0.2 are colored in green, while those with a globalScore < 0.2 or not matching any proteins in the list are in grey.

2.3. Differential Protein Expression across Cancer Subgroups and No Cancer Controls

Following normalization to adjust for technical variability, a linear mixed model was used to determine differential protein expression between subjects with cancer vs. no cancer (case definition 1) and GG ≥ 2 vs. GG 1 plus no cancer (case definition 2). The model also adjusted for batch effects and for the extracellular content of a subject’s serum (Section 4.9). For case definition 1, there were 130 protein groups with statistically significant differential expression (FDR < 0.05), including 58 upregulated and 72 downregulated (Figure 3A). Case definition 2 showed fewer differentially abundant protein groups, 38, with 14 and 24 of those up- and downregulated, respectively. Only 10 common protein groups were differentially abundant in both case definitions, but agreement in fold change for all proteins that were statistically significant across either comparison showed high correlation, R = 0.81 (Figure 3B). Functional enrichment analysis via 1D annotation enrichment was performed on the genes mapped to the protein groups in the dataset (Table S3) [42]. Selected biological processes meeting a significance threshold of FDR < 0.05 are shown in Figure 3C for upregulated enriched pathways across various cancer subtype comparisons and Figure 3D for those downregulated. Interaction networks from STRING [43] with highlighted proteins annotated for prostate cancer in the Open Targets database are displayed for proteins with upregulated (Figure 3E) and downregulated (Figure 3F) abundances. Enriched upregulated processes are observed for cytoskeletal organization and its regulation, which has known involvement in prostate cancer progression and proliferation specifically through its role in the epithelial-to-mesenchymal transition [44] and several metabolic pathways, such as glycolysis and pyruvate metabolism. Processes associated with gene regulation and translation were downregulated in multiple cancer subtypes relative to the no cancer subjects. Proteins annotated for these processes are indicated by different node shapes in the interactomes.

Open in a separate window

Figure 3

Differential protein expression for case definitions 1 and 2 and resulting functional enrichment analysis. (A) Volcano plots showing differential protein abundance from a linear mixed model for case definition 1 (left) and case definition 2 (right). Linear coefficients represent the log2 fold change in protein intensity between the respective indicated groups. (B) Scatterplot comparing the log2 fold change linear coefficients of protein groups with statistically significant (adjusted p-value < 0.05) differential abundance in case definition 1 vs. case definition 2. The dashed line in black, y = x, represents R = 1 correlation between the log2 fold changes in protein abundance and the solid line in teal is the resulting linear regression fit to the displayed data. The Pearson correlation coefficient, R, between these fold changes is displayed in the plot. (C) Upregulated and (D) downregulated processes according to 1D annotation enrichment analysis for different combinations of cancer and no cancer comparisons. (E) Interactome of the upregulated proteins and the downregulated proteins (F) for case definition 1. Edges represent interactions with a minimum STRING overall confidence score of 0.4.

2.4. Protein Multimarker (MM) Classifier for Prostate Cancer Diagnosis

A custom machine learning pipeline was developed to build a multimarker signature from the protein intensities for each nanoparticle in the training set. After data preprocessing, protein–nanoparticle features were filtered based on ranked p-values from a Wilcoxon rank sum test and pairwise correlation and then combined into a random forest classifier (Section 4.10.2). We focused our classifier building on the subset of subjects with PSA 4–10 ng/mL, as these are the intended population for PSA reflex tests. Data from subjects meeting this criterion in the training set were used to build a classifier separately for cancer vs. no cancer (case definition 1) and GG ≥ 2 vs. both GG 1 and no cancer (case definition 2). Optimization of the parameters for data preprocessing and feature selection was performed exclusively in this training set (Figures S8 and S9). This model was evaluated against those generated from AutoML, an automated machine learning pipeline that uses a Bayesian framework for model selection and hyperparameter tuning. AutoML was used separately on the entirety of the training data and to replace random forest as the classifier within our custom pipeline. This was performed to assess our custom pipeline for bias and overfitting (Figures S10 and S11, and Table S4). The protein multimarker signature for our optimized classifier achieved an average AUC of 0.64 (95% CI: 0.58–0.70) for case definition 1 using cross-validation compared to 0.60 (95% CI: 0.54–0.67) for the Prostate Cancer Prevention Trial (PCPT) risk calculator version 1.0 [35], which served as our clinical benchmark in the same set of patients. PCPT utilizes clinical factors, encompassing subject age, race, family history, DRE, and prior biopsy results in addition to PSA measurement to determine an assessment of risk for prostate cancer. For case definition 2, we observed an average AUC of 0.62 (95% CI: 0.56–0.68) with the protein feature classifier versus 0.60 (95% CI: 0.52–0.67) with the PCPT risk score. Full details of classifier performance are summarized in Table 2.

2.5. Validation of the Multimarker Signature

We next sought to validate the protein multimarker signature from our custom pipeline in the withheld validation set. The performance of the protein feature classifier for case definition 1 was 0.48 (95% CI: 0.39–0.58) for AUC when evaluated in the reserved validation set compared to the PCPT risk score of 0.70 (95% CI: 0.61–0.79), with a one-sided p-value of 1.0 (Table 2 and Table S5). Similarly, the protein feature classifier for case definition 2 was unable to surpass the performance of the clinical benchmark, reaching an AUC of 0.51 (95% CI: 0.40–0.62) versus 0.68 (95% CI: 0.58–0.78) for PCPT; the one-sided p-value was 0.99.

2.6. Bootstrapped Bias-Corrected Cross-Validation

Given the protein multimarker signature in the held-out validation set did not surpass the clinical benchmark, we combined samples from the training and validation sets to utilize the entirety of the study population and maximize the amount of data available to train the protein classifier. Along with considering only individuals in the PSA 4–10 ng/mL range, we also developed models using all subjects. To reduce bias and overestimation of out-of-sample performance from the selection of model hyperparameters outside of a fully nested cross-validation scheme, we utilized bootstrapped bias-corrected cross-validation (BBC-CV) [45] to select optimal model configurations, estimate out-of-sample performance with bootstrapping, and provide a measure of the algorithm’s ability in the selection, fitting, and tuning of the model (Figure S13). Details of the BBC-CV employed are further described in Section 4.11.2. As displayed in Figure 4, the performance of the protein multimarker signature in subjects restricted to PSA 4–10 ng/mL for case definition 1 as estimated using the BBC-CV was 0.59 (95% CI: 0.52–0.65) for AUC. For the PCPT risk score, the AUC was 0.63 (95% CI: 0.56–0.70). In case definition 2, the protein multimarker signature achieved an AUC of 0.53 (95% CI: 0.48–0.58) compared to 0.62 (95% CI: 0.54–0.69) for PCPT. Similarly, a classifier generated for all subjects approached but did not exceed the diagnostic performance of the PCPT benchmark (Table 2 and Figure S12).

Open in a separate window

Figure 4

Classifier performance of the protein multimarker signature. ROC curves are shown for the protein-level data and the PCPT risk score in subjects with PSA 4–10 ng/mL for cancer vs. no cancer, defined here as case definition 1 (A) and GG ≥ 2 vs. both GG 1 and no cancer, case definition 2 (B). Average AUC for the protein signature and the PCPT risk score are displayed for each respective plot, with the 95% confidence interval indicated in brackets.

4.2. Study Population

We used serum samples from two cohorts that collected blood/serum samples from men who subsequently underwent prostate biopsies. Both cohorts used sample handling and storage protocols from the Early Detection Research Network (EDRN). The first cohort, the Finasteride Challenge Study (FCS) [37], consisted of participants in a trial investigating the effect of a 3-month dosage of Finasteride on biomarkers predicting the outcome of a prostate biopsy. For the purposes of this study, baseline serum samples were used to exclude Finasteride contamination. Between 2011 and 2016, FCS enrolled men from this cohort 55 years or older without prior history of prostate cancer who had an estimated risk of prostate cancer between 20 and 60% based upon the Prostate Cancer Prevention Trial (PCPT) risk calculator version 1.0 [35]. At 90 ± 14 days, subjects underwent a 12-core prostate biopsy. Each core was fixed in formalin, processed for pathologic assessment, and individually read for the presence of prostate cancer. After biopsy, subjects were followed clinically using the local standard of care.

We also drew serum specimens from UT Health San Antonio’s SABOR (San Antonio Center of Biomarkers for Risk for Prostate Cancer) cohort [38]. SABOR is a multi-ethnic, multi-racial cohort of 4170 men from San Antonio and South Texas with no prior diagnosis of prostate cancer who agreed to undergo routine prostate cancer screening at the time of enrollment. Between 2000 and 2011, subjects were seen annually for PSA testing and DRE, with an average of approximately 6 visits per participant. Since 2011, subjects were followed with PSA only, and those with PSA less than 1.0 ng/mL were seen biennially. Biopsies were performed and evaluated according to community standards. Included in our sample were all men in the SABOR cohort who had a prostate biopsy during the course of their participation, excluding the small subgroup of men who requested a biopsy without clinical indication. For men in our study, we used the serum sample collected most proximally (within six months of) and prior to the first biopsy in the sample. If the serum corresponding to the first biopsy was not available, we used the serum corresponding to the second biopsy.

4.3. Cancer Outcomes

In both cohorts, cancer cases are defined based on having a diagnosis of prostate cancer within 18 months of the index biopsy evaluation. The Gleason group (GG) was obtained based on the most recent biopsy within 18 months of the blood draw. In the FCS, we had data on GG on prostatectomy if the patient had one within 18 months of the blood draw (n = 20). For our multiple marker signature and differential expression analyses, we considered two case definitions: (1) cancer versus no cancer and (2) clinically significant cancer (GG 2 or higher) versus low-grade (GG 1) or no cancer.

4.4. Preparation and Batch Randomization

Subject blood samples were drawn at UTHSA in outpatient clinics for the FCS cohort and in the community before transport to UTHSA on ice for the first five years of SABOR before transitioning for follow-up at UTHSA clinics by year 5 and onwards. The maximum time between draw and serum preparation was 4 h. All subject serum specimens were prepared at a common site (UTHSA) by trained operators using unified EDRN protocols [36]. All serum specimens selected for this study had a maximum of one freeze–thaw cycle during storage at UTHSA. Serum aliquots and associated clinical metadata were transferred from UTHSA to OHSU via the CEDAR (Cancer Early Detection Advanced Research Center) Specimen Biorepository with pre-assignments to the order and batching of the Proteograph^TM and mass spectrometry processing. Those handling the specimen or sample derivatives were blinded to the case/control (cancer/no cancer) status at all steps in the process.

Samples were processed in batches of 16 (15 experimental samples, 1 processing control, described in more detail below). To minimize systematic bias from batch effects, we used a stratified randomization procedure to assign case and control samples to each batch, ensuring that approximately equal numbers of cases/controls appeared in each block. Cases and controls within a block were randomly assigned. Prior to finalization of the randomization list, we examined the distribution of age at biopsy, PSA, cohort, and race/ethnicity across the batches and rerandomized if there were obvious imbalances in these variables across the batches.

We also used stratified randomization to assign 70% of specimens to training and 30% to internal validation sets. The randomization was stratified by case/control status, race/ethnicity, age, and PSA to ensure proportionally equal representation of these variables in training and validation sets. The training set was intended to be used for development of a multimarker panel. The validation was intended to be used as an independent test set to assess the performance of the multimarker panel.

4.5. Automated Blood Serum Sample Processing with Proteograph^TM Workflow

Serum specimens (N = 922 from both cohorts) and control samples were prepared on two SP100 automation instruments with a Proteograph^TM Assay Kit included in the Proteograph Product Suit (Seer, Inc., Redwood City, CA, USA) using five distinctly functionalized nanoparticles (NPs) as previously described [22,29], with specific details as follows. Each Proteograph plate accommodated 16 samples and 3 fixed, on-deck control samples (Process, Digest, and MPE Controls). Specific to this study, 1 of the 16 sample slots on each plate was reserved for analysis of PC3, a pooled plasma sample that served as an additional end-to-end workflow technical control. Upon thawing of serum specimens in advance of Proteograph preparation, qualitative observations regarding specimen color and turbidity were recorded. Immediately after SP100 on-deck dispensing of serum specimens for nanoparticle enrichment, residual serum from the original 250 μL input was manually collected and stored at −80 °C. Following completion of SP100 Automation Instrument on-deck peptide cleanup and elution, purified peptides were quantitated using the Pierce Quantitative Fluorometric Peptide Assay (Thermo Fisher Scientific, Waltham, MA, USA) using an on-deck SP100 routine and a standard plate reader instrument (Spark20M, Tecan, Männedorf, Switzerland). Samples with outlier peptide recovery, defined as greater than 2-fold or less than 0.5-fold the ratio of sample peptide mass to median peptide mass for all samples on the plate for the corresponding nanoparticle, were flagged for possible preparative issues and subsequent re-processing with a second aliquot of serum. To further assess preparation quality control, 65 μL of the 3 on-deck control samples was dried to completion in LowBind Eppendorf tubes and analyzed through mass spectrometry (described below). Elutions from all 80 sample NP enrichments and the remainder of on-deck control samples for each plate were transferred to a 96-well shallow plate, dried to completion overnight in a vacuum concentrator, and then frozen at −80 °C. This sample preparation routine was employed for individual analysis of all 922 patient specimens as well as for preparation of pooled patient specimen serum samples, either from residual collection or from supplemental patient serum aliquots when available, to support generation of a study-specific spectral library (described in more detail below).

4.6. Data-Independent Acquisition (DIA) Liquid Chromatography-Coupled Mass Spectrometry (LC-MS)

In real-time, during sample preparations using Proteograph^TM, the three standard on-deck controls were subjected to label-free dia-PASEF (data-independent acquisition parallel accumulation serial fragmentation)-mode liquid chromatography–mass spectrometry (LC-MS) in positive mode using a timsTOF Pro MS system (Bruker, Billerica, MA, USA), CaptiveSpray source (Bruker), and an inline-coupled nanoElute liquid chromatography instrument (Bruker). This served to conduct a proteomics assessment of the automated Proteograph workflow with minimal delay to avoid the further loss of subsequent samples due to problems with sample processing that were not obvious solely from peptide recovery. The timsTOF m/z and tims dimensions were linearly calibrated before each batch of samples or every 4–5 days for continual acquisitions using 3 selected ions from the Agilent ESI LC-MS tuning mix [m/z, 1/K0: (622.0290, 0.9915), (922.0098, 1.1986), (1221.9906, 1.3934). Dried control samples were resuspended manually in 2% ACN and 0.1% FA containing 5 pmol/μL of PepCal standard peptide mixture (Promega, Madison, WI, USA) for a final peptide concentration of 200 ng/μL and loaded into autosampler vials. For dia-PASEF, the MS scan range was 100–1700 m/z, the ion mobility window was 1/K0 0.6–1.6, the ramp and accumulation time was 100 ms, the source voltage was 1700 V, the dry gas was 3.0-L/min, and the dry temperature was 200 °C. The mobilogram shape and ion mobility window (1/K0 0.6–1.6) were tailored to Proteograph enrichments of serum. MS-coupled nanoflow reversed-phase chromatography at a flow rate of 0.5 μL/min was accomplished on the nanoElute instrument using mobile phase Solvents A (water, 0.1% formic acid) and B (acetonitrile, 0.1% formic acid), an inline PepMap C100 C18 cartridge trap (5 mm/0.3 mm/5 µm particle, Thermo Fisher Scientific), a C18 analytic column (15 cm × 150 µm × 1.5 µm particle size, PepSep, Marslev, Denmark) housed in a Column Toaster (Bruker) at 50 °C, and a 20 μm Classic Emitter (Bruker) with a 43 min total time acquisition (0 min, 5%B/28 min, 35%B/35 min, 40%B/36.67 min, 95%B/43 min, 95%). Baseline proteomic outputs for protein group and peptide identifications were initially set with 8 injections of the 3 Proteograph technical control types run on healthy serum controls. In concert with measurements of peptide recovery, peptide and protein ID rates of the 3 controls based on library-free database searches using DIA-NN [74] (details below in Section 4.8) were monitored on PAS and observed for drops below 2 × sd (standard deviation) from the mean (Figure S2). LC-MS performance of the instrument measuring these control samples was assessed every set of two plates with a K562 system suitability sample for numbers of protein groups and peptides. Deviations of either of these metrics below 2 × sd triggered an assessment of the instrument. Additionally, reconstructed peaks for PepCal spike-in peptides were also monitored qualitatively with Panorama AutoQC 21.1.0.158 [75].

For interrogation of Proteograph enrichments from cohort specimens and the PC3 control, DIA LC-MS was performed on individual nanoparticle enrichments from each specimen with separate injections. These acquisitions were performed using label-free dia-PASEF-mode LC-MS using two additional timsTOF Pro MS platforms, each coupled to an UltiMate 3000 nanoLC instrument (Thermo Fisher Scientific). Dried peptides from Proteograph enrichments were reconstituted in 2% ACN, 0.1% FA containing 5 pmol/μL of PepCal standard (Promega) to a final peptide concentration of 200 ng/μL using an on-deck Proteograph SP100 routine. Peptide (200 ng) from each NP enrichment was analyzed per injection. Nanoparticles 1–3 were analyzed on one timsTOF and 4–5 on the other to minimize technical MS acquisition bias across specimens. A standard dia-PASEF method was used with settings identical to those indicated above, except the ion mobility range was set to 1/K0 0.6–1.4 and a smaller mobilogram was used for faster cycle time. MS-coupled reversed-phase chromatography was accomplished using mobile phase A (water, 0.1% formic acid) and B (acetonitrile, 0.1% formic acid), an inline cartridge C18 trap (Thermo Fisher Scientific), a micropillar array C18 column (50 cm uPAC C18, PharmaFluidics, Gent, Belgium), and a 20 µm Classic Emitter (Bruker) with a 33 min total time method. Quality control of MS acquisitions ($>$4500 injections for the entire cohort analyzed on a per-nanoparticle basis) included assessment of peptide and protein identification rates from library-free DIA-NN database searches, determination of retention times and chromatographic profiles for reconstructed PepCal peaks using Skyline 21.1 and Panorama AutoQC 21.1.0.158 [75] and Compass DataAnalysis (Bruker), as well as monitoring of total intensity chromatograms and LC backpressures. A small subset of samples with issues flagged by the above quality control checks was prepared again from remaining serum aliquots and re-analyzed. For one sample, subject 476, acquisitions for NP1 and NP4 were unable to be completed and were not included in subsequent data analysis.

4.7. Spectral Library Generation

A subset of specimens from the two-cohort training set was selected as input material for spectral library assembly. First, training set serum samples exhibiting distinguishing color and/or turbidity (n = 18, the “qualitatively different” or “QD” samples) as detected through qualitative observation upon thawing were selected for spectral library generation. A second subset of the remaining training set specimens (“qualitatively normal” or “QN” samples) comprising the majority of the specimens were chosen as follows. During progression of DIA LC-MS acquisition for all cohort specimens, emerging raw data were monitored in real-time using library-free DIA-NN through the Seer cloud-based Proteograph^TM Analysis Suite (PAS, Seer, Inc.) to determine aggregate unique protein and peptide identifications and the contribution provided by the incremental specimen dataset addition. Saturation of new, unique peptide identifications (18,000–19,000 unique peptides) on resulting accumulation plots (Figure S1A) indicated that the subset of protein and peptide identifications observed through 450 samples would be approximately representative of the total proteomic diversity in the remainder of the cohort. Several approaches were tested for sample selection in the generation of a spectral library, including random sample selection and using maximal distance between selected samples from principal component analysis (PCA). A custom “MaxID” algorithm was developed, which identifies the samples that will provide the maximum number of possible peptide identifications based on the set of interrogated samples (deposited in GitHub on 8 July 2024 at https://github.com/CEDAR-Proteomics/Orion, Figure S1B). The same method can also be used to model the total number of peptides identified from a given number of samples using the MaxID method, allowing for the determination of the desired sample size used as input for the spectral library (Figure S1C). A targeted subset of 75 QN specimens was prioritized for spectral library inputs using this algorithm to reasonably represent overall cohort proteomic diversity.

Several approaches were tested for spectral library generation using these selected samples, and an approach using deep offline fractionation provided the deepest library. In the primary approach, pools of peptides for the QD and QN samples were generated with the five nanoparticle enrichments for each kept separate (10 pools total). To build these pools, both (a) peptides left over from initial Proteograph preparations and (b) supplementing peptides from second-round Proteograph preparations performed on additional serum aliquots of the same specimens provided by UTHSA were leveraged. For fractionation of these 10 resulting pools, offline, basic reversed-phase high-pressure liquid chromatography was performed by the University of Reno Las Vegas Proteomics Core Facility. Pooled samples submitted were dried down completely and reconstituted in 6 μL of Solvent A (10 mM of ammonium hydroxide in water, pH 10). Basic pH reversed-phase fractionation of 5 μL of this resuspension was performed on an UltiMate 3000 LC instrument (Thermo Fisher Scientific) fitted with an Acquity UPLC BEH C18 column (100 mm × 1.0 mm × 1.7 μm particle size, Waters, Milford, MA, USA) using a flow rate of 75 μL/min with an increasing gradient of Solvent B (10 mM of ammonium hydroxide in acetonitrile, pH 10; 2–55% over 47 min, total time 62 min). Fractions were collected every 50 s for 62 min and concatenated to 24 aggregate fractions for each peptide pool. All concatenated fractions were returned to OHSU, dried to completion, and resuspended in 7 μL of MS-grade water, 0.1% formic acid.

Each reconstituted fraction was analyzed through LC-MS using data-dependent acquisition parallel accumulation serial fragmentation (DDA-PASEF) [31] on a timsTOF Pro1 MS system (Bruker) using a CaptiveSpray source (Bruker) and inline-coupled UltiMate 3000 liquid chromatography instrument equipped for nanoflow (Thermo Fisher Scientific). For DDA-PASEF, the MS scan range was 100–1700 m/z, the ion mobility window was 1/K0 0.6–1.4, the ramp and accumulation time was 166 ms, the number of PASEF ramps = 8, the source voltage was 1700V, the dry gas was 3.0-L/min, and the dry temperature was 200 °C. MS-coupled nanoflow reversed-phase chromatography at a flow rate of 0.5 μL/min was accomplished on the nanoELUTE instrument using mobile phase Solvents A (water, 0.1% formic acid) and B (acetonitrile, 0.1% formic acid), an inline PepMap C100 C18 cartridge trap (5 mm/0.3 mm/5 μm particle, Thermo Fisher Scientific), a C18 analytic column (25 cm × 150 µm × 1.5 μm particle size, PepSep) housed in a column heater (MonoSLEEVE, Analytical Sales and Services, Flanders, NJ, USA) at 50 °C, and 20 µm Classic Emitter (Bruker) with a 120 min total time acquisition (0 min, 4%B/92 min, 30%B/102 min, 35%B/104 min, 95%B).

Serum from the same sample sets as used in the primary spectral library strategy detailed above were also pooled together in an alternative approach and processed with Proteograph. From these NP enrichments of pooled serum, peptides were analyzed through DDA-PASEF, as described above but using a 240 min total acquisition time. Library depth was compared to the first strategy, with fractionation providing the highest number of proteins and peptides. Ultimately, these were combined into one spectral library from both approaches for use in this study.

4.8. Raw MS Data Processing

Raw DIA MS (.d) data for each of the subjects and controls were searched in PAS with DIA-NN 1.8 against the canonical Uniprot human protein sequence fasta database with isoforms included (version from 2 March 2022). For DIA-NN library-free, default settings were used, except MBR was not enabled. The MS1 and MS2 Mass Accuracy were each set to 10 ppm. Precursor FDR (Q Value), Global Precursor FDR (Global Q Value), and Global Protein Group FDR (Global PG Q Value) were all 0.01. N-term M excision and C carbamidomethylation were set to True, while M oxidation and Acetyl (N-term) were False. For searches with the study-specific spectral library, the same parameters were used as described above, and Reannotate set to True.

To generate the study-specific spectral library, the MSFragger-DIA-NN workflow for spectral library generation was used [76]. FragPipe (version 18.0), MSFragger (version 3.5) [77], Philosopher (version 4.4.0) [78], and EasyPQP (version 0.1.30) were used to search the raw DDA-PASEF data files using the same fasta database employed above and to create a spectral library from the results. Reverse decoys were appended to the database using Fragpipe.

4.9. Differential Expression Analysis and Functional Annotation

Data processing and analysis were performed in R Statistical Software (version 4.2). Protein group intensities per NP quantified through MaxLFQ were log2-transformed and fold-change-normalized to mitigate technical variance across sample NP MS runs. For normalization, we initially identified protein groups present in at least 80% of the samples. Next, we computed the fold change per protein group per NP by subtracting the median intensity across samples from the intensity. Subsequently, the shift of fold change per sample per NP was determined as the fold change minus the median fold change across samples. Finally, we normalized the log2-transformed intensities of all proteins by subtracting the fold change shifts.

Protein–NP pairs absent in 60% or more of the samples were excluded (40% minimum sample detection), resulting in 3578 quantified protein groups eligible for differential analysis. Missing protein intensities were imputed for the NP with the highest number of non-missing intensity measurements for a given protein. The imputation process assumes that data are not missing at random, and that missing values are likely to be low values (below the limit of detection of the mass spectrometer). To that end, missing values were sampled from a normal distribution with parameters Normal $\left(I_0–1.8\sigma ,0.3\sigma \right)$, where $I_0$ represents the mean and $\sigma$ denotes the standard deviation of protein log2 intensities in the entire dataset [79].

A linear mixed model was employed to identify proteins exhibiting significant up- or downregulation in cancer compared to non-cancer (case definition 1) or in high-grade compared to low-grade cancer (case definition 2). Additionally, a model incorporating a three-level version of the PC status variable (non-cancer, low-grade cancer, and high-grade cancer, with non-cancer serving as the reference group) was fitted. For brevity, we present the model with single case status coding for case definition 1, denoted as $PCStatus$. The lme4 (version 1.1-35.1) package was utilized to fit the model separately for each protein, pooling information across NPs.

Each model adjusted fixed effects, including age at serum collection, year of serum collection, outlier status of samples based on appearance, and sample composition. Sample composition was quantified as a continuous covariate that accounts for the variability of blood plasma subfractions enrichment in the NP protein corona of a given Proteograph^TM sample. This variability stems from various factors, such as sample collection methods and age-related changes in blood composition. The variable was calculated as the median log2 fold change (relative to the median protein log2 intensity across the cohort) of proteins annotated as belonging to the “Extracellular” GO Cellular Compartment in each Proteograph sample.

For each protein, given plate i, NP j, and subject k, the linear mixed effects model of the normalized log-intensity values $Y_{ijk}$ was as follows (omitting additional fixed effects for brevity):

where $u_i$, $v_j$, $w_{ik}$, and $z_{k\left(i\right)}$ represent random effects for the plate, NP, plate $\times$ NP, and subject nested within the plate, respectively. For proteins with only a single NP having non-missing intensity values, the model excluded random effects $v_j$, $w_{ik}$, and $z_{k\left(i\right)}.$ This structure effectively captures the technical variance attributable to plate effects and accounts for the correlation between NP–protein pairs within individuals. The coefficient corresponding to the $PCStatus$ variable represents the mean increase or decrease in fold-change-normalized log2 intensity by case status, adjusting for other covariates.

p-values associated with the estimated $PCStatus$ coefficient, derived from the linear models for each protein group, were adjusted for multiple testing using the Benjamini–Hochberg method. Proteins with a corrected p-value below 0.05 were subjected to functional enrichment analysis using 1D annotation enrichment analysis [42]. These differentially abundant networks were further assembled into interaction networks with STRING 12.0 (accessed on 26 May 2024) [43] and further visualized in Cytoscape 3.10.2 [80].

4.10. Protein Multimarker Classifier

4.10.2. Custom Pipeline

The steps in the custom pipeline are illustrated in Scheme 1. Initially, features with more than 60% missingness were excluded. Subsequently, the data were fold-change-normalized using the previously described method. Next, the fold-change-normalized log2 intensities were adjusted for batch effects and sample composition by fitting the linear mixed model

$$Y_{ik}={\beta }_0+u_i+SampleCompositio{n}_{ik}+{ϵ}_{ik},$$

(2)

where $Y_{ik}$ represents the fold-change-normalized protein–NP intensity for plate $i$ and subject $k$, $u_i$ represents random effects for the plate, and $SampleCompositio{n}_{ik}$ is the sample composition covariate computed as previously described. The residuals obtained from this model replaced the log2 max-LFQ intensity as features in the subsequent stages of the pipeline. Finally, missing values in transformed intensities were imputed with zero imputation, under the assumption that data are not missing at random and that missing values are not present in the sample.

Open in a separate window

Scheme 1

Overview of the steps of the machine learning pipeline employed.

Following pre-processing, features undergo ranking based on Wilcoxon tests of differential expression in cases versus controls. Subsequently, pairwise correlations are computed among the top 1000 features. For a given feature, denoted as F, other features are identified whose correlation with F surpasses a predefined cutoff value, denoted as p. If any of these correlated features exhibit higher coverage than feature F, it is excluded from the set. The algorithm iterates through all features, retaining those that pass the filtering process. Post correlation filtering, the remaining top k features, identified through the Wilcoxon tests, are preserved. These remaining features are then utilized in the classifier.

4.11. Validation Analysis

We express our sincere appreciation to James McGann and Khatereh Motamedchaboki for their help in coordinating this collaborative inter-institutional study and to Max Mahoney and Paul Pease for Proteograph^TM technical support. We are grateful to OHSU CEDAR for providing funding to support this study and to the M. J. Murdock Charitable Trust in helping support the purchase of our Bruker timsTOF mass spectrometry platform at OHSU.