Tear Proteomic Predictive Biomarker Model for Ocular Graft Versus Host Disease Classification.

Purpose: Diagnosis of ocular graft-versus-host disease (oGVHD) is hampered by a lack of clinically-validated biomarkers. This study aims to predict disease severity on the basis of tear protein expression in mild oGVHD.
Methods: Forty-nine patients with and without chronic oGVHD after AHCT were recruited to a cross-sectional observational study. Patients were stratified using NIH guidelines for oGVHD severity: NIH 0 (none; n = 14), NIH 1 (mild; n = 9), NIH 2 (moderate; n = 16), and NIH 3 (severe; n = 10). The proteomic profile of tears was analyzed using liquid chromatography-tandem mass spectrometry. Random forest and penalized logistic regression were used to generate classification and prediction models to stratify patients according to disease severity.
Results: Mass spectrometry detected 785 proteins across all samples. A random forest model used to classify patients by disease grade achieved F1-measure values for correct classification of 0.95 (NIH 0), 0.8 (NIH 1), 0.74 (NIH 2), and 0.83 (NIH 3). A penalized logistic regression model was generated by comparing patients without oGVHD and those with mild oGVHD and applied to identify potential biomarkers present early in disease. A panel of 13 discriminant markers achieved significant diagnostic accuracy in identifying patients with moderate-to-severe disease. Conclusions: Our work demonstrates the utility of tear protein biomarkers in classifying oGVHD severity and adds further evidence indicating ocular surface inflammation as a main driver of oGVHD clinical phenotype. Translational Relevance: Expression levels of a 13-marker tear protein panel in AHCT patients with mild oGVHD may predict development of more severe oGVHD clinical phenotypes. Copyright 2020 The Authors.

Introduction

Graft versus host disease (GVHD) is a common complication that occurs in patients after allogenic hematopoietic cell transplantation (AHCT), which is performed increasingly to treat a variety of hematologic disorders. Ocular involvement of GVHD occurs in approximately 50% of AHCT recipients. Continuous improvements in long-term survivorship of posttransplant patients has led to an increase in the number of individuals living with ocular GVHD (oGVHD). Manifestations of chronic oGVHD include extensive inflammation and fibrosis at the ocular surface, which can lead to severe ocular surface dryness (sicca). Inadequately-treated sicca can result in corneal keratopathy, vascularization, ulceration, and eventual loss of vision. Subjective symptoms, which include ocular surface pain and photophobia, significantly reduce patient quality of life (QOL).

Currently, diagnosis of oGVHD depends largely on observation of symptoms and clinical signs, many of which have low specificity for the disease. Success of early detection depends on the skill of the attending physician, as well as the regular engagement of the patient with the ophthalmology team. The burden of treatment for other, more life-limiting, complications of their condition often precludes patients from attending specialist ophthalmology clinics. Despite this, it is widely accepted that late intervention with corticosteroids or other anti-inflammatory therapies substantially reduces the efficacy of treatment of oGVHD. Therefore there remains an unmet need for the identification of objective diagnostic or prognostic biomarkers that could be used to determine risk of disease progression. Because of the possibility for noninvasive and repeated sampling, the use of tear proteins as biomarkers of ocular surface disease is an area of growing research interest.,

Biomarker studies, particularly those using omics technologies, can generate large amounts of data per patient sample. Machine learning (ML) techniques, which lie at the interface of mathematics, statistics, and computer science, can be applied to make predictions about future datasets based on models created using known data. Supervised ML models are trained to generate predictions about future data through provision of input data that is labeled with classifier information. In recent years, ML and deep learning approaches have been applied to diverse problems in medicine to aid physicians in diagnosis of patients.– In ophthalmology, such methods have been applied to OCT and color fundus images and tear proteomic data to classify diabetic retinopathy patients and to color fundus photographs to classify glaucoma patients. Biomarker studies in conjunction with ML represent an important potential future avenue for disease classification.

To date, several studies have proposed panels of biomarkers for the diagnosis of oGVHD. These groups have focused solely on measuring inflammatory marker expression at the protein level in tears and as mRNA in conjunctival cells.– We previously described the tear proteomic profile of AHCT recipients with moderate-to-severe oGVHD compared with those without oGVHD, in a cross-sectional analysis. In that study, we used an unbiased mass spectrometry discovery-proteomics workflow to uncover the protein content of tear samples from oGVHD patients, versus oGVHD-negative controls. In the current study, we report tear proteomic analysis from patients with oGVHD through all severity grades, rated according to the NIH disease burden scale. Using the large number of tear protein expression data obtained from each patient, the study aimed to evaluate whether ML methods can effectively classify disease severity based on a predictive tear protein expression signature.

Methods

Patients

This cross-sectional study was performed at the Department of Ophthalmology, University Hospital Basel, Basel, Switzerland. The study was approved by the local ethical review board (Ethikkommission Nordwest- und Zentralschweiz) and adhered to the tenets of the Declaration of Helsinki. All participants gave informed consent.

Participants were recruited during outpatient clinics at the Department of Ophthalmology, University Hospital Basel. The study enrolled 49 participants who had previously undergone AHCT for hematological disorders. All patients were examined by a single attending physician (DG). After ocular examination, patients were classified by oGVHD status according to the National Institutes of Health (NIH) consensus criteria: grade 0 (n = 14), grade 1 (n = 9), grade 2 (n = 16), grade 3 (n = 10). In brief, the NIH grading system classifies oGVHD based on severity of impairment of activities of daily living (ADL), extent of therapy required and degree of vision loss arising from sicca (0 = normal, no symptoms; 1 = mild, no effect on ADL, hydrating drops <3 times/day; 2 = moderate, some effect on ADL, hydrating drops >3 times/day or punctal plugs, no vision loss; 3 = severe, significant effect on ADL, loss of vision caused by KCS). In this study, all data and analyses are presented for one eye per patient. Eye selection was not randomized; sample choice was based on the overall clinical profile observed by taking into account the results of the Oxford score and Schirmer and Tear Break Up Time diagnostic tests. For the patients with oGVHD (NIH 1-3) the worse eye was chosen. For patients without oGVHD (NIH 0), the best eye was chosen.

Clinical Examination and Sample Collection

Patients underwent a full clinical ocular surface examination (best corrected visual acuity, slit lamp examination of ocular surface, Oxford score grading, tear break up time), including clinical and therapeutic history. Patients provided information on quality of life using the ocular surface disease index (OSDI) questionnaire (Allergan Inc, Irvine, CA, USA). Tear samples were collected using Schirmer strips (Schirmer-Plus; Gecis, Neung-sur-Beuvron, France). To reduce the risk of contamination, gloves were worn throughout the procedure. No topical medication was used before the Schirmer test was performed. Strips were inserted at the lateral third of the lower eyelid according to standard procedure, without the use of topical anesthesia, and collected after five minutes. Tear volume was noted, and strips were immediately frozen on dry ice and stored at −80°C until analysis.

Proteomic Analysis of Tear Samples

Mass spectrometric analysis of tear proteins was performed as described in detail previously. In brief, tryptic digests were prepared from small pieces of Schirmer test strips and acidified with trifluoroacetic acid (TFA; Applied Biosystems, Rotkreuz, Switzerland) at a 1% final concentration. Peptides were desalted on C18 SepPak cartridges (Waters, Dättwil, Switzerland) and washed with a solution of 0.1% TFA. Forty-microliter aliquots of sample were prepared for analysis by liquid chromatography tandem mass spectrometry (LC-MS/MS) using an Orbitrap FT hybrid system (Thermo Fisher Scientific, Waltham, MA, USA). Full details of the injection and flow parameters are provided by Gerber-Hollbach et al.

Tryptic digests from a group of healthy volunteers without ocular surface disease were pooled to provide material for use as an internal standard in the analysis. The LC-MS/MS data were searched with proteomics software (Proteome Discoverer 1.4; Thermo Fisher Scientific) set to Mascot search engine with 10 ppm precursor ion tolerance, whereas the fragment ions were set to 0.6 Da tolerance. The following modifications were used during the search: carbamidomethyl-cysteine was set as a fixed and oxidized methionine as a variable modification. Match-searches were carried out at a 1% false-discovery rate. To analyze equal amounts of peptides, pairwise digests (patient sample versus pooled control) were run. An initial run was carried out on the mass analyzer instrument (Thermo Fisher Scientific), and the resulting chromatogram was integrated. To compensate for individual differences in peptide amount obtained from the Schirmer strips, the integrated chromatogram from patient and internal standard were adjusted to ensure equal peptide material was injected. Three technical replicates were run for each sample. The search files from Proteome Discoverer were loaded into Scaffold (Proteome Software, Portland, OR, USA) and the proteins were quantitated by spectral counting. The peptide false-discovery rate was set to 1%, whereas the protein threshold was 95%. For protein quantification, the number of usable peptides was set to 3. Protein expression was calculated as ratios versus the internal standard, and these data were used in downstream analysis.

Statistical Analyses

Demographic and Clinical Data

Statistical analysis was performed by a biostatistician (A.S.). Demographic and clinical data were summarized as mean ± standard deviation or as percentages, where appropriate. When data were normally distributed, ordinary one-way ANOVA with Tukey's post-hoc test was used to detect significant differences between groups. In cases of nonparametric data, the Kruskal-Wallis test was used to compare between groups.

Development and Validation of Predictive Tear Protein Models

We used two methods to generate predictive models for oGVHD severity based on tear protein expression. In the first approach, the scikit-learn Python library (version 0.21.1) was used to train a random forest (RF) model with 20 estimators and a maximum depth of 10. Missing values (below level of detection) were imputed by providing a value that was half of the lowest detected value. Importantly, in this approach, no data were removed from the dataset to be analyzed. To prevent overfitting, the unfiltered data were divided into two parts: 80% used for model training and 20% used for model testing. The training set was further divided using stratified 10-fold cross validation, with 20% kept for validation. Classification performance was reported as a confusion matrix and the individual feature importance was reported for the top 100 proteins.

In the second method, any variables with more than 20% missing values (below level of detection) were removed before analysis. Missing values in the remaining dataset were imputed using an algorithm for left censored data, according to the method of Wei et al. Additionally, data were log-transformed for further analysis. To compare protein expressions among all severity grades (NIH 0-3), nonparametric analysis (Kruskal-Wallis Test) was performed, and corresponding P-values were reported. P-values were FDR-controlled using the BH procedure. Results were sorted by adjusted P-values to give a top list of the proteins with statistically significant differential expression and visualized using boxplots.

To create a predictive model for oGVHD severity based on tear protein expression, penalized logistic regression was performed using the package “Glmnet” in the statistical software R. Details are described elsewhere. The model was optimized using tear proteome data from patients with mild disease (group 2, NIH 1) versus those with no disease (group 1, NIH 0). To prevent overfitting, regression was internally optimized using 20 iterations of five-fold cross-validation. Optimization was done using the AUC of the corresponding ROC curves. Results were reported indicating median and quantiles of the AUCs. A feature selection based on “Glmnet” was presented reporting odds ratios (OR). Sensitivity and specificity were not reported because they are not appropriate for small sample sizes. Subsequently, the model was validated on group 1 (NIH 0) versus groups 3 and 4 (NIH 2 and 3), presenting predicted probabilities with corresponding boxplots and descriptive statistics. Internal cross-validation and prediction was done using the package “caret” within R.

To quantify correlation between the expression levels of the selected proposed biomarker proteins and clinical parameters, Spearman's ranked correlation coefficient was calculated. A P-value < 0.05 was considered significant. All evaluations were done using R software version 3.5.2 (2018, https://www.r-project.org/).

Functional analysis was carried out to determine the biological significance of the features selected by the model. The list of selected protein IDs were converted to gene IDs using the www.uniprot.org database conversion tool (release February 2019). Gene Ontology (GO) term enrichment was assessed using the GO and panther databases (www.pantherdb.org, release October 2019) with the human reference genome as the background dataset. GO terms for biological process were ranked according to fold enrichment.

Results

Demographics of Study Population

The demographic data are outlined in Table 1. All participants had previously undergone AHCT for various hematological disorders. Group 1 included 14 AHCT recipients (nine men) with no clinical evidence of oGVHD, group 2 included nine patients (seven men) with mild oGVHD, group 3 included 16 patients (nine men) with moderate oGVHD, and group 4 included 10 patients (nine men) with severe oGVHD after AHCT, according to the NIH grading criteria. Patient groups were age-matched. The mean ages at the time of inclusion in the study were as follows: group 1 (56.1, range 38–73 years), group 2 (48.4, range 25–69 years), group 3 (52.6, range 28-74 years), and group 4 (52.6, range 24–69 years). There were no significant differences in time elapsed between AHCT procedure and study inclusion. Mean elapsed times in months were as follows: group 1 (82.7, range 2–265 months), group 2 (66.1, range 20–170 months), group 3 (88.9, range 16–206 months), and group 4 (90.6, range 27–199 months). All participants, with the exception of two in group 1, received peripheral blood stem cell transplantation.

Table 1.

Demographic Data of Patients Included in the Study

	Group 1 (NIH 0)	Group 2 (NIH 1)	Group 3 (NIH 2)	Group 4 (NIH 3)
Characteristic	n = 14	n = 9	n = 16	n = 10
oGVHD severity	None	Mild	Moderate	Severe
Age, yrs
Mean ± SD	56.1 ± 9.6	48.4 ± 15.4	52.6 ± 14.0	52.6 ± 15.2
Range	38–73	25–69	28–74	24–69
Female, n (%)	5 (36)	2 (22)	7 (44)	1 (10)
Primary diagnosis, n (%)
NHL	0 (0)	1 (11)	4 (25)	1 (10)
MM	2 (14)	0 (0)	0 (0)	2 (20)
ALL	0 (0)	2 (22)	3 (19)	2 (20)
AML	8 (57)	3 (33)	7 (44)	3 (30)
MDS	2 (14)	0 (0)	2 (13)	0 (0)
CML	1 (7)	0 (0)	0 (0)	1 (10)
CLL	0 (0)	1 (11)	0 (0)	0 (0)
Other	1 (7)	2 (22)	0 (0)	1 (10)
HCT source, n (%)
PB	12 (86)	9 (100)	16 (100)	10 (100)
BM	2 (14)	0 (0)	0 (0)	0 (0)
No. of HLA-mismatch, n (%)
1	1 (7)	1 (11)	2 (13)	1 (10)
2	0 (0)	1 (11)	2 (13)	0 (0)
Unrelated donor, n (%)	3 (21)	5 (56)	8 (50)	5 (50)
Time elapsed between AHCT and ocular exam (mo)
Mean ± SD	82.7 ± 85.7	66.1 ± 50.8	88.9 ± 62.6	90.6 ± 68.7
Range	2–265	20–170	16–206	27–199

No statistically significant differences were detected in age or time elapsed since AHCT between the groups.

Therapy

Because of the chronic nature of oGVHD, individuals included in this study were in treatment for their ocular surface disease. An outline of therapies used, showing numbers and percentages of participants using each therapy, is provided in Supplementary Table S1. Two patients (14%) in group 1 had punctal plugs historically placed, and, additionally, two patients (14%) used artificial tears infrequently to supplement ocular surface hydration but had no clinical features to warrant classification as having oGVHD. In group 2, six patients (67%) had punctal plugs in place. Artificial hydration was used infrequently (< 5x/day) by 4 patients (44%) and more often (<10×/day) by three patients (33%). One patient (11%) used topical steroids. In groups 3 and 4, more than 60% of patients in each group used topical steroids daily. One patient (6%) in group 3 and four patients (40%) in group 4 used cyclosporine. Between 80% to 100% of groups 3 and 4 patients had punctal plugs in place. The majority (69%) of patients in group 3 required frequent supplementary hydration (> 5×/day). In group 4, 70% of patients used artificial tears very frequently (more than 10 times daily). Additionally, autologous serum was used by two patients in group 3 (13%) and one patient in group 4 (10%). Full details of systemic immunosuppressive therapies in use are provided in Supplementary Table S2.

Ocular Parameters of Study Population

The results of the ocular examinations are outlined in Table 2. Extent of ocular surface keratopathy was assessed using fluorescein staining and scored according to the Oxford scale. All participants in group 1 scored 0 or 1 (mean 0.1 ± 0.4). Participants in groups 2, 3, and 4 had scores ranging from 1 to 5. Mean Oxford scores were 1.4 ± 0.5 (group 2), 3.2 ± 1.2 (group 3), and 4 ± 0.7 (group 4). Mean tear production, as assessed by Schirmer test I, was within normal range for group 1 (19.3 ± 8.9 mm) and significantly progressively depressed in all oGVHD groups. Mean Schirmer test readings were 8.9 ± 8.3 mm, 2.8 ± 3.6 mm, and 2.3 ± 2.9 mm for groups 2, 3, and 4, respectively.

Table 2.

Ocular Characteristics of Patients Included in the Study

	Group 1 (NIH 0)	Group 2 (NIH 1)	Group 3 (NIH 2)	Group 4 (NIH 3)
Diagnostic Criteria	n = 14	n = 9	n = 16	n = 10
Oxford Score
Mean ± SD	0.1 ± 0.4	1.4 ± 0.5	3.2 ± 1.2	4 ± 0.7
Range	0–1	1–2	1–5	3–5
Sign. vs. NIH 0		P > 0.05	P < 0.0001	P < 0.0001
Schirmer Test I (mm)
Mean ± SD	19.3 ± 8.9	8.9 ± 8.3	2.8 ± 3.6	2.3 ± 2.9
Range	7–35	2–25	0–11	0–10
Sign. vs. NIH 0		P < 0.05	P < 0.0001	P < 0.0001
TBUT, s
Mean ± SD	4.2 ± 2.3	2.9 ± 1.6	2.1 ± 0.8	1.9 ± 0.6
Range	2–10	1–6	1–3	1–3
Sign. vs. NIH 0		P > 0.05	P < 0.01	P < 0.01
OSDI score
Mean ± SD	6.5 ± 5.6	19.8 ± 12.4	37.8 ± 17.4	67.1 ± 12.8
Range	0–20.8	10.4–43.8	11.4–68.2	43.2–81.8
Sign. vs. NIH 0		P < 0.05	P < 0.0001	P < 0.0001

When data were normally distributed, significant difference (sign.) was tested using one-way ANOVA with Tukey's post hoc test; when data were not normally distributed Kruskal-Wallis test was used.

Tear film stability was evaluated by length of tear breakup time (TBUT). In all groups, TBUT was shortened on average, when five seconds is used as the standard cutoff for healthy individuals. Group 1 individuals had mean TBUT of 4.2 ± 2.3 seconds, falling just outside the normal range. Group 2 had a mean TBUT of 2.9 ± 1.6 seconds (P > 0.05 compared with group 1). Groups 3 and 4 had significantly shortened TBUTs, measured at 2.1 ± 0.8 seconds and 1.9 ± 0.6 seconds, respectively (P < 0.01 compared with group 1).

Subjective negative impact on QOL, measured by the symptom-based OSDI score, was observed in all oGVHD patients. Compared with a mean OSDI score of 6.5 ± 5.6 in patients without oGVHD, patients in group 2 (mean 19.8 ± 12.4, P < 0.05), group 3 (mean 37.8 ± 17.4, P < 0.0001), and group 4 (mean 67.1 ± 12.8, P < 0.0001), all had significantly more symptomatic complaints.

Tear Protein Expression in Ocular GVHD Samples

Relative quantitative protein analysis from human tears detected 785 proteins in the patient cohort studied. Using the raw dataset a RF classifier model was trained to classify patients according to disease severity. The results of this classification are illustrated by means of a confusion matrix (Fig. 1A) and summary table (Table 3). The F1-measure values achieved by the classifier were 0.95 (group 1), 0.8 (group 2), 0.74 (group 3), and 0.83 (group 4). Ranking of selected features importance for the top 100 proteins did not reveal any features with distinctly high contribution to patient classification (Fig. 1B; Supplementary Dataset S1). The top five features (MDI > 0.025) were proline-rich protein 4 (PRR4), polymeric immunoglobulin receptor (PIGR), immunoglobulin J chain (JCHAIN), proline-rich protein 1 (OPRPN), and deleted in malignant brain tumors 1 protein (DMBT1).

Figure 1.

A. Confusion matrix illustrating performance of a random forest model in classifying oGVHD disease severity. Columns represent predicted disease severity, whereas rows represent actual disease severity. B. Feature importance computed for the top 100 proteins, with the red bar indicating feature importance and the black line showing intertree variability.

Table 3.

Performance of Random Forest Classifier in Predicting Disease Severity in oGVHD Patients

Disease Severity	Precision	Recall	F1-Measure
Group 1 (no oGVHD)	0.90	1.00	0.95
Group 2 (mild oGVHD)	0.80	0.80	0.80
Group 3 (moderate oGVHD)	0.78	0.70	0.74
Group 4 (severe oGVHD)	0.83	0.83	0.83

Because of the known possibility of overfitting in datasets with large numbers of features compared with the number of individual samples, a more strict statistical analysis was subsequently performed. After a data cleaning step, in which proteins not detected in at least 80% of individuals were removed, a final list of 46 robustly detected proteins was generated (Supplementary Dataset S2). In this set of 46 proteins, there was no significant difference in the median variance of protein expression between control and patient samples (P value = 0.7129) when taking the multiple replicates measured for each individual into account. Using nonparametric analysis, a top list of 12 proteins significantly differentially expressed with advancing disease severity was generated (FDR-corrected P value < 0.05). Six proteins—lactotransferrin (LTF), lysozyme C (LYZ), PIGR, JCHAIN, prolactin inducible protein (PIP) and immunoglobulin heavy constant alpha (IGHA1)—were progressively downregulated with increasing disease severity (Table 4; Fig. 2). Four proteins—actin (ACTB), annexin A2 (ANXA2), glutathione S-transferase P (GSTP1), and phosphoglycerate mutase 1 (PGAM1)—were progressively upregulated with greater disease severity (Table 3; Figs. 3A–3D). Two proteins—keratin type II cytoskeletal 6A (KRT6A) and pyruvate kinase PKM (PKM)—were upregulated in tears from moderate-to-severe disease, but not in patients with mild disease (Table 4; Figs. 3E, 3F).

Table 4.

Tear Proteins With Significant Differential Expression Among Patients With oGVHD of All Severity Grades

Gene ID	Accession ID	Protein Name	Profile	P Value
LTF	P02788	Lactotransferrin	↓	<0.0001
LYZ	P61626	Lysozyme C	↓	<0.0001
PIGR	P01833	Polymeric immunoglobulin receptor	↓	<0.0001
JCHAIN	P01591	Immunoglobulin J chain	↓	0.0009
PIP	P12273	Prolactin-inducible protein	↓	0.0015
IGHA1	P01876	Immunoglobulin heavy constant α1	↓	0.0020
ACTB	P60709	Actin, cytoplasmic 1	↑	0.0009
ANXA2	P07355	Annexin A2	↑	0.0104
GSTP1	P09211	Glutathione S-transferase P	↑	0.0205
PGAM1	P18669	Phosphoglycerate mutase 1	↑	0.0205
KRT6A	P02538	Keratin type II cytoskeletal 6A	↓↑	0.0107
PKM	P14618	Pyruvate kinase PKM	↓↑	0.0184

Significance was tested using nonparametric Kruskal-Wallis test.

Figure 2.

Relative expression of proteins with significantly downregulated expression with advancing disease severity. P values were calculated with Kruskal-Wallis test.

Figure 3.

Relative expression of proteins with significantly upregulated expression with advancing disease severity. Expression levels of KRT6A (E) and PKM (F) were significantly higher only in patients with moderate-to-severe disease. P values were calculated with Kruskal-Wallis test.

We aimed to identify a protein expression signature present early in disease that could also stratify patients with moderate-to-severe disease, using multiple penalized logistic regression. A feature selection of tear proteins capable of discriminating between patients without oGVHD (group 1) and those with mild oGVHD (group 2) was performed. Calculated ORs for each of the panel of 13 selected proteins are shown in Table 5. GO analysis revealed enrichment of biological process terms related to tissue homeostasis and immune regulation in the set of selected features (Supplementary Dataset S3). As a means of testing the capacity of the model to predict oGVHD in patients with moderate-to-severe disease, the remainder of the patient cohort (groups 3 and 4) was used as a validation dataset. The predicted probabilities of the classified patients are shown in Figure 4. Analysis of correlation between the suggested biomarkers and the clinical measures of disease severity using Spearman's rank order correlation test. Of the 13 proteins, three were detected as having significant negative or positive correlations with various clinical measures (Table 6, Supplementary Table S3). Both LYZ and PIGR showed significant negative correlations with NIH score, corneal fluorescein staining (Oxford score) and clinical symptoms (OSDI) (P < 0.0001). Both proteins were positively correlated with TBUT and Schirmer test results (P < 0.0001). PGAM1 showed significant positive correlation with NIH score, Oxford score and OSDI (P < 0.01). Negative correlation was detected between PGAM1 and TBUT (P < 0.01) and Schirmer test (P < 0.0001).

Table 5.

Tear Proteins Selected With Multiple Logistic Regression (glmnet) as Discriminating Features of oGVHD, Shown With Corresponding Odds Ratios

Gene ID	Accession ID	Protein Name	Odds Ratio
PGAM1	P18669	Phosphoglycerate mutase 1	3.8574
KRT9	P35527	Keratin type I, cytoskeletal 9	2.3164
KRT1	P04264	Keratin type 2, cytoskeletal 1	1.6653
FABP5	Q01469	Fatty acid binding protein, epidermal	1.6161
PFN1	P07737	Profilin-1	0.8869
IGKC	P01834	Immunoglobulin κ constant	0.6005
DCD	P81605	Dermicidin	0.3012
S100A4	P26447	Protein S100-A4	0.2725
LYZ	P61626	Lysozyme C	0.2417
PIGR	P01833	Polymeric immunoglobulin receptor	0.1130
GAPDH	P04406	Glyceraldehyde 3 phosphate dehydrogenase	0.1130
ALB	P02768	Serum albumin	0.0240
GSN	P06396	Gelsolin	0.0036

Figure 4.

Predicted probabilities of oGVHD group calculated on the basis of selected features.

Table 6.

Correlation of Expression of Selected Proteins With Clinical Parameters

Gene ID		NIH Score	Oxford Score	TBUT	Schirmer Test	OSDI
PGAM1
	rs	0.453	0.404	−0.476	−0.517	0.382
	P value	0.001	0.004	0.001	< 0.0001	0.007
LYZ
	rs	−0.757	−0.665	0.594	0.67	−0.723
	P value	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001
PIGR
	rs	−0.807	−0.729	0.605	0.607	−0.772
	P value	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001

Correlations were determined for the entire study population using Spearman's rank coefficient (rs). Table shows only selected biomarkers shown to be significantly positively or negatively correlated with clinical measures.

Discussion

The clinical management of ocular GVHD following AHCT is still hindered by a lack of robust and specific early clinical markers of the disease. This can lead to a delay in commencement of treatment and less favorable outcomes. The objective of this study was to investigate whether ML techniques could be applied to tear proteomic expression data to classify patients with ocular GVHD after AHCT. More specifically, we analyzed the ability of a predictive model to identify changes in protein expression in early (mild) disease that could subsequently be used to correctly identify patients with moderate-to-severe disease.

Several studies have reported the expression of inflammation-related genes and proteins in the conjunctival cells and tear film, respectively, of oGVHD patients.–,– Cocho et al. used a support vector machine (SVM) model to demonstrate the predictive power of a panel consisting of four genes—EGFR, IL-6, IL-9 and NAMPT—expressed in conjunctival cells to discriminate oGVHD patients from controls. Subsequently, this group showed, using a logistic regression model adjusted for age and sex, that tear concentrations of IL-8 and CXCL10 could predict oGVHD status. To our knowledge, we were the first group to use untargeted proteomic methods to investigate oGVHD tear protein expression. In this way, we have broadened the number of investigated targets beyond inflammatory pathways alone.

The current study enrolled 35 patients with oGVHD of various degrees of severity, graded according to the NIH criteria and evidenced by the progressively worsening ocular surface clinical phenotype (Table 2). Fourteen individuals without oGVHD after AHCT were recruited as controls. The patients included in this study were not treatment-naïve, with the majority of patients in the oGVHD groups receiving some form of immunosuppression (Supplementary Table S2). However, our extensive clinical examination demonstrated that the patients included in the oGVHD groups experienced significant negative signs and symptoms of their disease, in spite of their treatment regimen (Table 2).

In this study, unbiased discovery mass spectrometry was used to report the proteomic profile of human oGVHD tears. We used relative quantification of tear proteins versus an internal standard consisting of tears from a pool of healthy volunteers. Initial analysis revealed 785 unique proteins detected across the entire dataset. The overarching aim of this study was to apply modelling techniques to classify oGVHD patients by disease severity based on their tear protein expression profile. Previous studies have demonstrated the better diagnostic capacity of panels of molecules, compared to individual biomarkers.

Random forest is a popular machine learning classification tool used to deal with large datasets. A RF model trained on 80% of the data and tested on the remaining 20% demonstrated moderate ability to classify patients according to disease severity (Fig. 1A, Table 3). However, a ranking analysis of the importance of individual proteins in decision-making using RF did not reveal any features or set of features with a distinct contribution to the model (Supplementary Dataset S1). The top five proteins (using a cutoff of MDI > 0.025)—PRR4, PIGR, JCHAIN, OPRPN and DMBT1—have all previously been demonstrated as differentially expressed in various forms of ocular surface dryness. The use of RF models for classification problems may have good performance, but interpretation is complicated because they generate no clear feature selection. For this reason, we performed further analysis using penalized logistic regression (glmnet).

The initial RF analyses were carried out on the raw dataset, with missing values handled by providing a standard value that was equal to half the limit of detection observed. However, to allow for more consistent analysis and reduce bias introduced by the imputed values, a data-cleaning step was used, whereby only proteins detected in 80% of patients were further investigated. This resulted in a concise list of 46 tear proteins (Supplementary Dataset S2), which was used for further analysis.

Six proteins were detected as progressively downregulated with increasing disease severity (Fig. 2). In agreement with other proteomic profiles of dry eye syndromes,– expression of the lacrimal gland defense proteins LYZ, LTF and PIP was depressed in oGVHD patients, indicating substantial lacrimal impairment. Furthermore, we detected significant downregulation of the secretory immunity proteins PIGR and JCHAIN. Secretory immunity is an important defense mechanism in all mucosal epithelia, including at the ocular surface, and the molecular machinery for the production of its main effector, secretory immunoglobulin A (sIgA), has been shown in the lacrimal gland and conjunctival epithelium. PIGR is responsible for the transcytosis of immunoglobulin A (IgA) to the apical surface of epithelial cells where it can be cleaved and secreted as a pathogen binding complex known as secretory component (SC). Reduction of PIGR expression in tears has been shown in a rabbit model of Sjögren's syndrome (SS) and in human SS tears. JCHAIN is an adaptor molecule that binds IgA dimers and facilitates their transport across the apical epithelial cells of the conjunctiva and lacrimal gland. Reduced expression of PIGR and JCHAIN has previously been shown in patients with moderate to severe dry eye, and implies extensive damage to the lacrimal and conjunctival epithelial cells. In turn, the reduced SC availability may leave the ocular surface more vulnerable to infection in oGVHD patients. We also detected reduced expression of IGHA1, which is produced by B lymphocytes and involved in humoral immunity. Previous tear proteomics studies have shown downregulation of IGHA1 in aqueous-deficient dry eye (ADDE) and SS patients.

In this study, four proteins were detected as continuously upregulated with progressive oGVHD severity—namely ACTB, ANXA2, GSTP1 and PGAM1. Two proteins—KRT6A and PKM—were upregulated in moderate and severe oGVHD (Fig. 3). ACTB is a cytoskeletal protein involved in cellular structure and motility. The presence of extracellular ACTB implies widespread cellular damage; it is released from dying cells into the extracellular space, where it causes significant further toxicity. Interestingly, Perumal et al. reported ACTB as being downregulated in human ADDE tears, whereas, in agreement with our findings, other groups have demonstrated increased ACTB in tears with advancing age and in SS patients. Plasma gelsolin (GSN) is an actin-scavenging protein with a role in promoting corneal re-epithelialization. GSN has been shown to have increased expression in ADDE tears. Similarly, we detected increased GSN in tears in moderate and severe oGVHD (Supplementary Dataset S2).

Upregulation of the calcium-binding protein ANXA2 has been shown previously in tears and conjunctival cells of patients with various forms of DED. ANXA2 is a pleiotropic molecule with various groups suggesting functions in corneal wound healing, macrophage activation, and regulation of secretions. The role of extracellular ANXA2 in oGVHD needs to be further investigated. We detected upregulation of GSTP1 in patients with moderate and severe oGVHD. GSTP1 is an enzyme involved in the detoxification of reactive oxygen species and its expression has been shown in corneal epithelial cells. Increased expression of GSTP1 has been shown in tears of ADDE patients., PGAM1 and PKM are enzymes involved in glycolysis. Upregulation of glycolytic pathways conceivably may be related to increased cellular stress and subsequent activation of repair pathways in corneal epithelial cells, as has been shown previously shown in epidermal keratinocytes.

We sought to build a predictive model based on tear protein expression in mild disease that could identify patients with more severe forms of oGVHD. To achieve this, we used a form of penalized logistic regression (glmnet)., Glmnet fits a generalized model by applying maximum penalized likelihood which shrinks unimportant features in class prediction to zero. It is suited for problems with smaller numbers of samples as it is less prone to overfitting. Glmnet is readily interpretable, since it generates predicated probabilities for class membership in the form of odds ratios.

In this study, penalized logistic regression selected 13 proteins expressed in mild oGVHD tears as features that could correctly classify patients with moderate and severe oGVHD (Table 4). Interestingly, of these, only three proteins—PGAM1, PIGR and LYZ—were detected as significantly differentially expressed with advancing disease severity. The main benefit of using an unbiased proteomics approach is the ability to identify possible novel players in disease progression, beyond already-known pathways. However, analysis of gene ontology (GO) annotation in the marker panel (Supplementary Dataset S3) primarily showed enrichment of terms related to tissue homeostasis and immune response, corroborating previous oGVHD gene and protein biomarker studies.,, We performed correlation analysis to determine the relationship between the 13 chosen biomarker proteins and the clinical findings (Table 6). We detected significant positive or negative correlations for three of the proteins—namely LYZ, PIGR, and PGAM1. Intuitively, reductions in the lacrimal gland and defense proteins LYZ and PIGR were associated with higher symptom (OSDI) and corneal fluorescein staining (Oxford) scores. Positive correlations were detected between these proteins and Schirmer and TBUT scores, indicating pathological tear film dynamics. Levels of both proteins were strongly inversely related to overall NIH score. Regarding PGAM1, weaker positive and negative correlations with clinical measures were detected. Levels of PGAM1 were positively associated with epithelial pathology (Oxford score) and symptomatic complaints (OSDI). Levels of this protein were negatively correlated with Schirmer test and TBUT scores.

This study has a number of limitations, primarily the use of immunosuppressive therapies to control disease in oGVHD patients, particularly those with moderate-to-severe disease. However, we contend that because these patients still have persistent symptoms and clinical signs, their inclusion is important to give an accurate picture of the day-to-day clinical reality of oGVHD. Perhaps related to this unbiased patient selection, we have identified some discrepancies between our results and those of other groups that have studied tear proteomics in untreated dry eye syndromes. Furthermore, here we report the proteomic profile of a relatively small number of oGVHD patients, when patients are subdivided by disease severity group (minimum group size: n = 9 in group 2). Therefore it would be interesting to further validate our findings in larger numbers of patients and external datasets. More extensive evaluation of the specificity of our model could be achieved by applying our findings to proteomic datasets obtained from other ocular surface disorders such as MGD, SS or Stevens Johnson syndrome. Mass spectrometry is generally considered unsuitable for widespread clinical diagnostic use. Follow-up validation of our chosen markers with enzyme-linked immunosorbent assay or multiplex array will go some way to providing clinical translatability. Finally, it is important to acknowledge that the protein changes detected here can also be seen in other forms of dry eye, because functional analysis indicates them as markers of ocular surface and LFU damage. However, the goal of this article was to identify molecules that are changed in a subset of patients with mild oGVHD following AHCT. Our objective was to assess whether these proteins could predict more severe disease in this patient subgroup.

In summary, we have compared two approaches to extract predictive biomarkers from tear proteomic profiles of patients with oGVHD of worsening severity grades after AHCT. Using strict filtering of raw proteomics data and application of a glmnet model, we have identified a 13-marker panel of proteins expressed in mild disease that is capable of identifying more severe forms of oGVHD. Our work has provided further evidence to support previous studies demonstrating ocular surface inflammation as a main driver of oGVHD clinical findings.