Cross Trial Meta and Time Series Pharmacodynamic Analyses of Apremilast Effects

Linking molecular disease biomarkers, drug-induced pharmacodynamic effects and response status in human trials can provide number of values to patient beneﬁts: 1) advancement in understanding the mechanism-of-action of current therapy; 2) back-translation utilities to fast-track preclinical development of next-generation therapeutics. Both opportunities are predicated on proactive generation of molecular biomarker datasets that capture longitudinal trajectories in patients before and after pharmacological intervention. Here, we present the largest plasma proteomic biomarker datasets available to date from placebo-controlled Phase III clinical trials of the phosphodiesterase type 4 inhibitor apremilast in psoriasis (PSOR), psoriatic arthritis (PsA), and ankylosing spondylitis (AS), and the corresponding meta-analyses. Using approximately 150 plasma analytes tracked across three time points from 526 subjects, we found robust biomarkers of PSOR severity, IL-17A and KLK-7, which were also reduced in apremilast responders. Furthermore, when compared to placebo, we observed a shared apremilast-speciﬁc pharmacodynamic pattern across all three diseases, with IL-17A universally represented in this pattern. Other notables in this pattern included MDC and HE4, common to the AS and PSOR cohorts, and IL-1RII in PSOR and PSA. Taken together, these ﬁndings provide foundational surrogate biomarkers surrogates to accelerate preclinical efforts including compound differentiation, synergy, and repurposing therapeutic strategies. pro-inﬂammatory genes, including inducible nitric oxide synthase (iNOS), IL-12/IL-23 subunit p40, IL-23 subunit p19, IL-17A, IL-22, and IL-8. In a placebo-controlled double-blinded phase III PSOR trial (ESTEEM2), treatment with apremilast 30 mg BID resulted in signiﬁcant reductions of plasma IL-17A, IL-17F, IL-22, and TNF- α protein levels 15 . In a placebo-controlled double-blinded phase III study in PsA, apremilast treatment was associated with decreased plasma levels of IL-6, IL-8, MCP-1, MIP-1 β , TNF- α , ferritin, and a small increase in von Willebrand factor (vWF) plasma protein levels. Among these, the changes in TNF- α and vWF were signiﬁcantly associated with achieving an American College of Rheumatology (ACR) 20 clinical response. After 40 weeks of apremilast treatment, signiﬁcant decreases in plasma IL-6, IL-23, IL-17, and ferritin, and signiﬁcant increases in IL-10 and IL-1RA were observed 11 .


Introduction
Psoriasis (PSOR) is an autoimmune disease and manifests as thickened erythematous patches of skin covered with silvery scales that can arise across the entire body but often at predisposed areas such as the extensor aspects of elbows and knees, nails, scalp, palms, soles, and intertriginous areas. The psoriatic lesions are frequently associated with pain and pruritus 1,2 . Histologically, the skin lesions are characterized by acanthosis from rapid keratinocyte proliferation, hypogranulosis, parakeratosis from aberrant keratinocyte differentiation, erythema caused by dilation of blood vessels in the papillary dermis, and a dense inflammatory infiltrate consisting of dendritic cells and T cells and dendritic cells. Psoriatic arthritis (PsA) is a chronic inflammatory arthritis present in up to 42% of individuals with PSOR, and is classified as one of the seronegative spondyloarthropathies, based on the presence of spinal involvement in about 50% of patients 3 . In the majority of patients, PSOR precedes PsA by several years, especially in patients with nail involvement 4 . PsA is manifested as one of 5 subtypes: distal interphalangeal joint involvement only, asymmetrical oligoarthritis, symmetrical polyarthritis, spondylitis, and arthritis mutilans 5 . Ankylosing spondylitis (AS) is a systemic chronic, inflammatory disorder primarily affecting the spine and characterized by sacroiliitis and lower back pain that can result in destruction and fusion of spinal vertebrae 6 . PSOR, PsA, and AS are pathologically related conditions sharing common clinical manifestations and immunological drivers characterized by dendritic cells and helper T cells producing Th1 and Th17 cytokines including IL-23, TNF-α, IFN-γ, IL-17, and IL-22 [7][8][9] .
Taken together, these results suggest that partial inhibition of several proinflammatory mediators, including TNF-α and IL-17 may play a key role in the mechanism of action of apremilast in the treatment of psoriatic diseases. Expansion of larger protein biomarker sets that link disease and drug effect will help organizing next-generation effective therapies.

Results
The abundance of 155 plasma proteins was measured by Myriad RBM sanwich immunoassays and analyzed in subjects with AS, PsA, or psoriasis PSOR who participated in 3 independent Phase III apremilast clinical trials. Subjects were randomly selected to have received apremilast (20mg or 30 mg arms combined together in the current study for downstream analysis) or placebo during the first 16 weeks of respective trials. For the present analyses, subjects were considered to be responders at Week 16 if they met the following criteria based on disease severity scores (Table 1). In the AS study, clinical responses as measured by a 20% improvement in the Assessment of SpondyloArthritis international Society (ASAS20) response rates were similar in the apremilast treatment groups compared with placebo 16 . In the PSOR trial at Week 16, significantly more apremilast-treated subjects achieved a 75% improvement in Psoriasis Area and Severity Index (PASI) score (PASI 75) (28.8%), PASI50 (55.5%) and static Physician's Global Assessment (sPGA) score of 0 or 1 (20.4%) vs. placebo (5.8%, 19.7%, 4.4%, respectively; P < 0.001) 17 . In the PsA study at Week 16, clinical response was measured by a 20% improvement in the modified American College of Rheumatology (ACR20) response criteria (p<0.001) 18 . At Week 24, plasma levels of IL-8, TNF-α, IL-6, MIP-1β, MCP-1 and ferritin were significantly reduced from baseline in the apremilast 20mg BID or 30mg BID treated subjects versus placebo, and the ACR20 response correlated with change in TNF-α level with both apremilast doses 11 . For the current analyses, analytes were measured at baseline, Week 4 and Week 16. Expressions for 121 analytes in AS, 122 in PSOR, and 155 in PsA passed QC and were used for subsequent analyses.

Systematic discovery of disease biomarkers
We performed systematic correlations of analyte expression with disease scores across 3 diseases in these trials. In order to take advantage of disease score changes captured in the trials, we explored different ways to slice the data including separating and combining measurements from all time points, with and without apremilast treatments. In addition, delta-delta correlations between changes of analyte level and severity scores were also calculated. Univariate regression analysis controlling for age and gender were performed. Best regression models were selected based on F-statistics (p<0.05 with FDR correction) and analytes with significant correlations by the model (p<0.05 with FDR correction) were reported.
In AS subjects, the analytes IL-6 and LRG1 significantly correlated with Ankylosing Spondylitis Disease Activity Score (ASDAS) as reported in Table 2. Statistical significance of LRG1, IL-6 and AGP-1 was confirmed in delta-delta correlations between changes in analytes and changes in ASDAS (Supplementary Tables 1-3). In PSOR subjects, significant correlations were reported between IL-17A and KLK-7 with PASI total scores (Table 3) and confirmed in delta-delta correlations . No significant correlations were found between the Disease Activity Score of 28 joints (DAS28) score of PSA subjects and analytes measured on the DiscoveryMAP immunoassay panel. Next, we explored the possibility of additive effects by evaluating the strength of combining individually significant biomarkers found in the univariate analyses. We first explored the dependencies amongst analytes by lasso regularization (Supplementary Tables 7-8). Analytes shown to be independent in the lasso models were included for subsequent multivariate analyses. We showed that additive models significantly captured more variability than univariate models in PSOR ( Figure 1A, Supplementary Tables 9-10). For potential clinical utilities, variables in the additive models were subsequently combined into one score to monitor apremilast effect. In PSOR subjects, the combined score of IL-17A and KLK-7 was significantly downregulated after apremilast treatment ( Figure 1B). Similarly in AS subjects, the combined score of LRG-1 and IL-6 was a stronger biomarker than either alone (Supplementary Tables 11-12). Interestingly, the combined expression-based score in AS showed that apremilast treatment had little effect on the biomarkers, offering a plausible molecular link consistent with an overall lack of clinical response (Supplementary Figure 1).

Apremilast effects at each time point
We next examined the apremilast effect at each snapshot in each trial, and by combining all trials at the same time points. Differential expression analyses were applied to compare placebo and apremilast, as well as between apremilast responders and non-responders. Very few analytes were robustly different between placebo and apremilast when comparing at weeks 4 and 16 separately (Supplementary Figure 2). Meta-analyses pooling from both time points increased the detection power within each disease. In psoriasis subjects, MDC, IL-16, AXL, CA-15-3, FasL, IGFBP-2, Maspin and SP-D were significantly downregulated by apremilast. Of note, AXL and its upstream stimulus, IL-16 are known to be involved in Th1 19,20 , Treg 21 , and Th17 22 regulation. MDC and IL-16 were also recently reported among biomarkers of atopic dermatitis severity 23 , a disease with overlapping clinical symptoms as psoriasis. In AS subjects, IP-10, IL-17A, RAGE were downregulated while EPO was upregulated in the apremilast arm. EPO is known for its inflammation reducing properties 24 , and together with IP-10, is a product of the JAK-STAT signaling pathway (KEGG:hsa04630) 25 . In PsA subjects, IL-1RII, BAFF, SOD-1 and YKL-40 were downregulated while CgA and MPIF-1 were upregulated in the apremilast arm. It is understood that IL-1RI, IL-1RII and IL-1B form complexes that amplify innate responses and promote Th17 differentiation [26][27][28] . Overall, we observed biomarkers of inflammation downregulation by apremilast that represented Th17, Th1 and JAK-STAT pathways.
We did not observe overlap of differentially abundant proteins between diseases in these snapshot statistical models. Metaanalyses of pooled effect across diseases increased detection power, however, and identified a set of analytes that underwent significant changes induced by apremilast (Supplementary Figure 2). Among these analytes, several were found to be associated with disease progression in the previous section: IL-17A as confirmation to existing knowledge, and LRG-1 as a new biomarker that links between apremilast impact and disease progression in AS patients. TNFα, IL-6, IL-8, CCL4 were confirmed as apremilast pharmacodynamic markers in the shorter week 4 -week 16 time frame as in previous studies at weeks 24 and 40 11 . Outside of our own analyses, literature search revealed strong disease relevancy of the other pharmacodynamic biomarkers reported here from the meta-analyses. Notably, there were therapeutic targets (TNFα, IL-2A, IL-6, IL-8, IL-12A, IL-12B, IL-18) and reported biomarkers specific to these 3 diseases: CRP, VEGF, MMP3, OPG in all three diseases [29][30][31] ; Osteocalcin in AS, PSA 30,32 ; CCL3, CCL4, YKL-40 in PSA and PSOR 29, 33-36 ; TGFβ and TIMP-1. Fibrinogen, was also shown to be upregulated and SOD-1 to be downregulated in PSOR 33 .
We applied the same meta-analysis approach to distinguish the apremilast response effect. PGI was found to be upregulated in apremilast PSOR responders; KLK-7 and angiogenin were upregulated in PSA responders (Supplementary Figure 3). We were unable,however, to construct predictive models using baseline measurements in any of the 3 diseases (Subblementary Figure 4).

Apremilast pharmacodynamic effect by time-series analyses
We reasoned that using static models to identify significant apremilast effect at each time point was a high bar to investigate underlying biology in inflammatory diseases, which are known to involve overlapping and compensatory pathways. As an alternative, we sought to capture dynamic patterns along the longitudinal axis to obtain an overall picture of the processes via time-series analyses. We adapted an algorithmic workflow (see Methods) by applying the CoGAPS algorithm 37 , a recently developed method based on Markov chain Monte Carlo non-negative matrix factorization, to identify time-dynamic patterns. When jointly applied with PatternMarkers 38 , the association of analytes with each pattern were probabilistically determined by trajectory distance calculation. We applied the workflow to identify patterns associated with both apremilast pharmacodynamic (apremilast vs. placebo comparisons) and Week 16 response status (apremilast responder vs. non-responder comparisons).
In apremilast vs. placebo comparisons, we observed a pharmacodynamic pattern that was consistent across diseases ( Figure  2A). The pattern captured a "bump" or upregulation of a group of analytes in the placebo group at Week 4 before a return to baseline level at Week 16. By comparison, apremilast induced a sustained downregulation of these analytes at Week4 and Week 16. Although a majority of the analytes associated with this pattern vary by disease, we observed IL-17A as a common analyte across all diseases and consistent with existing knowledge about apremilast effect. In addition to IL-17A, IL-1RII was present in this pattern in both PSOR and PsA. MDC and HE4 were present in this pattern in both AS and PSOR.
By nature of the algorithm, analytes differentially regulated by apremilast treatment in the snap-shot based meta-analysis were captured by various patterns. Of interest, LRG1, Hemopexin, AGP-1 and SAA, observed in an AS pattern, were all shown in the first part of our study to be ASDAS disease score biomarkers. In this pattern, subjects in the apremilast arm had higher biomarker levels, although the clinical study was randomized. This difference was not directly observed by the ASDAS score, however, which may suggest the utility of biomarkers to capture subtle signals that are not fully quantifiable by clinical manifestations.
Several important disease-specific patterns were captured in apremilast non-responder vs. responder comparisons. In psoriasis subjects, we observed a pattern represented by IL-17A and KLK-7 ( Figure 2B, middle panel, pink pattern) that were downregulated longitudinally in both the non-responder and responder groups. Despite similarities in downward direction, the mean starting level was lower and the downward magnitude larger in the responder group. The earlier discovery that both IL-17A and KLK-7 are strong biomarkers for disease severity provides rationale to explain the differential pharmacological benefits observed. When combined as an expression-based score, we observed by conventional between-group statistics a descending downward trend of the combined biomarkers between placebo, non-responders, and responders ( Figure 1B).
The most frequent pattern for ankylosing spondylitis subjects represented by Prostate Specific Antigen, free (PSA-f), displayed downregulation at Week 4 in the responder group. Overall, AS is strongly associated with male gender, as reported in a recent global survey among 1548 axial and peripheral spondylorthritis patients which found that 64.9% were men 39 . The ASAS20 response rate for apremilast 30 mg BID was 36% in males vs 30% in females. ASAS20 response rates were not significantly greater, however, in the apremilast-treated comparing to placebo-treated subjects in either sex (33% in males vs 28% in females). We were therefore unable to interpret the observed pharmacodynamic significance of PSA-f in the responder group.

Discussion
Disease biomarker discoveries along with drug-induced pharmacodynamics studies are necessary to bring appropriate treatment for patients. The benefits of such linkage are twofold: patient segmentation for precision medicine and back-translation into next-generation of therapeutics. Large consortia are emerging to facilitate disease subsets by molecular biomarkers 40 toward these goals. We advocate that complementary biomarker datasets in well-controlled clinical trials focusing on treatment-induced pharmacodynamics and responses in the same diseases are necessities to power our continued innovations in medicine of autoimmune diseases. Availabilities of these datasets and subsequent analyses rely on joint efforts from pharmaceutical companies and technology partners in well-powered clinical trial settings.
In this study we demonstrate the values of linking disease-centric and drug-centric biomarker studies through meta-analyses in 3 apremilast Phase III trials of psoriasis, ankylosing spondylitis, and psoriatic arthritis. Plasma analytes were profiled in participating subjects in the trial, and several computational models were applied to explore implications to disease severity, apremilast effect and clinical response status.
In psoriasis subjects, we identified IL-17A and KLK-7 as robust positive biomarkers to PASI total scores that were significantly downregulated by apremilast in all subjects. The magnitude of downregulation was trending more in responders over non-responders, making the two biomarkers uniquely positioned to link disease-and apremilast-centric biomarkers. IL-17A was a previously known driver for PSOR and its downregulation shown to be a mechanism of action for the efficacy of apremilast. This study confirmed these findings and added further evidence that inadequate downregulation of IL-17A is a key differentiation factor between clinical non-responders to apremilast compared to responder counterparts. KLK-7 displayed a similar pattern to IL-17A, prompting consideration of its causal role in psoriasis manifestation and apremilast efficacy. Kallikrein-related peptidase 7 (KLK-7) is a serine protease involved in skin barrier functions. A systematic literature survey revealed a general association to primarily skin-related diseases (Figure 3), including key publications on its dysregulation in atopic dermatitis 41,42 and psoriasis 43 . To our knowledge, our profiling effort is the largest study to confirm misregulation of KLK-7 in psoriasis subjects. For validation, we surveyed KLK-7 mRNA expression patterns in human psoriasis subjects in a series of public gene expression datasets with and without treatment (Figure 4). Upregulation of KLK-7 mRNA expression patterns was observed across multiple studies, confirming our proteomic measures in psoriasis subjects. KLK-7 upregulation was consistent in both peripheral blood and lesional skin biopsies. Interestingly, like our findings for apremilast, two other known therapeutics for psoriasis, broadalumab (anti-IL-17) and etanercept (anti-TNF), downregulate KLK-7 as well. Taken together, these findings support KLK-7 as a therapeutic biomarker with a potential causal mechanism for disease etiology. As a primer for back-translation, we propose that the combined expression-based score of IL-17A and KLK-7 protein levels may serve as a preclinical surrogate to estimate the effectiveness of next-generation psoriasis therapy ( Figure 1B).
Bridging IL-17A and KLK-7 in disease-centric and treatment-centric biomarker analyses appears to apply exclusively to psoriasis in our study. In the AS cohort, we did not discover the same for robust disease-centric biomarkers LRG-1 and IL-6. LRG-1 was upregulated by apremilast without differentiation impact on efficacy. As a result, we were not able to interpret its impact by apremilast as on-target and causal.
Although clear individual biomarker signals exist (IL-17A, KLK-7, IL-6, LRG-1), synchronized dynamics from a group of analytes may better capture the underlying biological complexity. Clinical trials present unique opportunities whereby these coordinated dynamics are induced and subsequently detected in the presence of perturbational conditions. There were two types of perturbation in our study design: longitudinal and pharmacological. By utilizing both types of perturbation, we were able to utilize time-series analyses to capture subgroups of apremilast-specific pharmacodynamics. One universal pharmacodynamic pattern showed sustained apremilast-induced downregulation across all three diseases and was commonly represented by IL-17A. Other than IL-17A, analytes included in this pattern varied by disease: MDC and HE4 were common in AS and psoriasis, while IL-1RII was common in PSOR and PSA.
Another PDE4 inhibitor, roflumilast, was tested in vitro using sputum cells from patients with chronic obstructive pulmonary disease, but no effect on spontaneous MDC (CCL22) production was found 44 . There are no reports of PDE4 inhibitor effects on HE4, which has recently been proposed as a biomarker of ovarian cancer 45 . IL-1RII has previously been reported to be overexpressed in the basal compartment of PSOR lesional epidermis 46 . There are no previous reports of apremilast or any other PDE4 inhibitor affecting IL-1RII expression. Therefore, the findings here that apremilast can reduce levels of MDC, HE4, and IL-1RII are novel and exemplify the learnings efforts of an unbiased meta-analysis. We speculated that these subtle differences may shed light into shared biology between these related diseases.
In time-series analyses of apremilast responders vs. non-responders in the psoriasis cohort, we noted that most patterns displayed similar apremilast effects in these two groups at weeks 4 and 16, but with different starting points at baseline ( Figure  2B). These observations prompted subsequent explorations of baseline predictors for eventual response status. The application of machine learning methods did not yield robust predictive models (Supplementary Figure 4), likely due to the combination of uneven distribution of the groups, and small samples after slicing for time points/treatment arms. In meta-analyses, that pool data from both time points to increase power to detect differential protein expression within each disease, MDC, IL-16 and AXL were significantly downregulated by apremilast in the PSOR study. MDC is discussed above. Expression of AXL (AXL Receptor Tyrosine Kinase) has been shown previously to be reduced upon cAMP treatment in glioblastoma cells 47 , although its relevance in psoriasis is unclear. IL-16 is a product of M1 macrophages and its expression is induced by chronic IL-17 exposure 48,49 . Apremilast reduces IL-17A levels, thus it is therefore not surprising that other IL-17-dependent cytokines such as IL-16 would be decreased as a result.
Interestingly, IL-10, IL-17A, and RAGE are downregulated by apremilast in AS subjects. Although apremilast has been reported to enhance IL-10 production by monocytes, when IL-10 is produced by T cells apremilast has an inhibitory effect 50,51 . This suggests that the IL-10 measured in the plasma of AS patients is mainly derived from T cells rather than monocytes. RAGE (S100A12) is a Th17 pathway-dependent protein produced by keratinocytes 52 , and we have recently reported that apremilast significantly reduces RAGE gene expression in primary human epidermal keratinocytes exposed to IL-17 in vitro 53 . Therefore changes in IL-10, IL-17A, and RAGE observed in apremilast-treated AS patients are consistent with previously observed preclinical mechanistic studies of apremilast. Such confirmation highlights the biological plausibility of the biomarker changes observed with the current analysis.
Using 185,360 plasma biomarker measurements collected from 3 independent apremilast Phase III trials in ankylosing spondylitis, psoriasis and psoriatic arthritis, we identify IL-17A and KLK-7 as disease-, apremilast-, and efficacy-relevant biomarkers with insights into causal mechanisms in psoriasis patients. In addition, IL-6 and LRG-1 are identified as diseasecentric biomarkers for ankylosing spondylitis disease severity. Finally, time-series analyses identify IL-17A, KLK-7, MDC and AXL as cross-disease apremilast pharmacodynamic biomarkers. These findings suggest high potential for use of these biomarkers as surrogates in preclinical development of next-generation therapeutics including combination strategies and compound differentiation. We demonstrate value and advocate the incorporation of exploratory molecular biomarker profiling into future clinical trials.

Protein quantification
Protein profiling was performed by Myriad RBM that offered absolute protein quantitation in a CLIA certified laboratory with their MAP immunoassay panels for 150 analytes and Simoa ultrasensitive immunoassays for 5 analytes (Supplementary  Table 13). For each trial, plasma samples were taken from each patient at each time point, and assayed with this platform. Pre-processing was performed on the raw data by removing individuals and proteins that had more than 50% missing values or values under the lower limit of quantitation (LLOQ). Remaining missing values were set to the protein average, and remaining

5/13
. CC-BY-NC-ND 4.0 International license certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was not this version posted May 30, 2019. . https://doi.org/10.1101/652875 doi: bioRxiv preprint <LLOQ values were set to 1 2 the protein LLOQ. Values were then log2 scaled and quantile normalized. In each trial, roughly 15% of the proteins were measured as <LLOQ in nearly all patients, and these proteins were largely consistent between the trials. This suggests possible mis-calibration of certain analytes to protein levels not observed in the data. Most proteins ( 80%) had no <LLOQ measurements, while about 5% had at least some <LLOQ measurements (between 0% and 50%). Finally the expression data for 121 analytes in AS, 122 in PSOR and 155 in PSA was used for the following analysis.

Calculation of additive combined expression-based score
Combined analyte scores were calculated in two ways: 1) unweighted penalized geometric mean, we used the combination of geometric mean for positively correlated analytes and penalized geometric mean for negatively correlated analytes: where k -number of positively correlated analytes and l -number of negatively correlated analytes.
2) weighted score, constructed and tested multivariate additive models using the significant analytes from univariate analyses. All results were reported by the second approach.

Correlation analyses with disease scores
Correlations for AS subjects were calculated based on 3 sets of disease scores: ASDAS, BASDAI, and BASFI. For PSOR patients, PASI total score was used. The clinical score used for analysis of the PSA trial was the DAS28 score. Univariate linear models predicting clinical score with covariates and additive models were applied using R package stats. P-values were considered to be significant when less than 0.05 after FDR correction. Analytes selected as significant biomarkers, correlated significantly with disease score and R 2 > 0.2 in each time point and with all time points together.

Differential expression analysis
Differential expression analysis was done using limma R package 54, 55 with adjusted p-value threshold 0.05 after FDR correction. Linear regression was applied using treatment arms, gender and, in the case of meta-analysis, time points as covariates.

Time-series analyses
Both CoGAPS and PatternMakers algorithms were applied using R package CoGAPS 37, 38, 55 . For each analyte, mean expression values across samples at each contrasting condition were appended by rows into a vector. Vectors of values from each contrasting condition were appended by columns into a matrix. Similarly, a matrix with corresponding standard errors was used to estimate goodness of fit by the CoGAPS module. We used 5000 iterations for each CoGAPS run, each enforcing a minimum of three patterns. Algorithm was run for 100 times with different random seeds for each comparison. On the next step, similar patterns in runs were merged. Pattern similarities were calculated using RMSD. The threshold for the distance function was set to allow 5% of deviation of relative strength in each condition. Analyte memberships belonging to each pattern were calculated and ranked by the PatternMarkers algorithm. Analytes were assigned to a pattern if the occurrence frequency exceeded 50% across all the runs.

Systematic literature survey of KLK-7 and diseases
Sentence-level co-occurrence of terms referring to KLK-7 and terms describing phenotypes from MEDLINE R abstracts from January, 1991 to June 2018 were extracted using SciBite c v6.2 TERMite Expressions (TExpress) in Batch mode. A TExpress pattern specifying entity types including GENE and HPO was employed to ensure that synonyms of KLK-7 and all phenotype concepts from the Human Phenotype Ontology (HPO) were considered. Following the co-occurrence extraction, a specificity score for each unique gene-phenotype co-occurrence was calculated based on mutual information between the gene and phenotype assessed by tf-idf scores, as described by Frijters et al. 56 A higher specificity score shows that the association between the phenotype and the gene is more relevant. Each unique phenotype was tagged with their highest HPO parent class(es) under phenotypic abnormality based on the HPO hierarchical relationships, in order to elucidate essential findings.

Systematic mRNA survey of KLK-7 patterns in public psoriasis gene expression datasets
Relevant datasets were compiled, downloaded and re-analyzed from Gene Expression Omnibus (GEO) 57 . For each dataset, batch effects were explored and removed by the SVA package 55,58 . limma package 54, 55 was used for transcriptome-wide differential analyses with FDR-corrected p-value.   Week 16 in placebo arm (PBO) and apremilast arm (APR). List of proteins is sorted in descending order of occurrence of the proteins in the pattern. Proteins in common between patterns are bolded and KLK-7 is marked in red. B) CoGAPS patters when comparing the protein expression levels across baseline, Week 4 and Week 16 in apremilast non-responders (ANR) and apremilast responders (AR). In bold -proteins with occurrence of association with pattern greater than 80%, in grey -with occurrence less than 60%.

12/13
. CC-BY-NC-ND 4.0 International license certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was not this version posted May 30, 2019. . https://doi.org/10.1101/652875 doi: bioRxiv preprint

13/13
. CC-BY-NC-ND 4.0 International license certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was not this version posted May 30, 2019. . https://doi.org/10.1101/652875 doi: bioRxiv preprint