Gene expression signatures predict response to therapy with growth hormone

Abstract

Recombinant human growth hormone (r-hGH) is used as a therapeutic agent for disorders of growth including growth hormone deficiency (GHD) and Turner syndrome (TS). Treatment is costly and current methods to model response are inexact. GHD (n = 71) and TS patients (n = 43) were recruited to study response to r-hGH over 5 years. Analysis was performed using 1219 genetic markers and baseline (pre-treatment) blood transcriptome. Random forest was used to determine predictive value of transcriptomic data associated with growth response. No genetic marker passed the stringency criteria for prediction. However, we identified an identical set of genes in both GHD and TS whose expression could be used to classify therapeutic response to r-hGH with a high accuracy (AUC > 0.9). Combining transcriptomic markers with clinical phenotype was shown to significantly reduce predictive error. This work could be translated into a single genomic test linked to a prediction algorithm to improve clinical management. Trial registration numbers: NCT00256126 and NCT00699855.

Conflict of interest statement

AS and PM have received speaker honoraria from Merck KGaA, Darmstadt, Germany. P Ch has received investigator research support, consultant and speaker honoraria from Merck KGaA, Darmstadt, Germany. P Cl had received research investigator support and speaker honoraria from Merck KGaA, Darmstadt, Germany. EK is an employee of Merck KGaA, Darmstadt, Germany.

© 2021. The Author(s).

Figures

Fig. 1. The association of whole blood gene expression at baseline with response to recombinant human growth hormone (r-hGH) over all 5 years of therapy in patients with growth hormone deficiency (GHD) and Turner syndrome (TS).

Comparison of patient response to r-hGH using Discriminant Analysis of Principal Components (DAPC). Low quartile (green, LoQ) and high quartile (red, HiQ) of growth response over 5 years of therapy (cms grown) compared to the remaining patients (orange) in GHD (N = 50) and TS (N = 22). Unsupervised transcriptomic data with no normalisation for phenotype are shown, GHD = 8875 and TS = 8455 gene probesets. DAPC generates a discriminant function, a synthetic variable that optimises the variation between the groups whilst minimising the variation within a group. The frequency of the discriminant function of DAPC is plotted (colour figure online).

Fig. 2. Whole blood gene expression is associated with response to recombinant human growth hormone (r-hGH) over 5 years of therapy in patients with growth hormone deficiency (GHD) and Turner syndrome (TS).

Partial least squares discriminant analysis (PLS-DA) of unsupervised transcriptome using three components. The low and high quartiles of growth response are shown for response to r-hGH (cm) over 5 years in A GHD and B TS. Star plot shows sample distance from the centroid, the arithmetic mean position of all the points in each group.

Fig. 3. Predictive value of whole blood gene expression associated with response to recombinant human growth hormone (r-hGH) in patients with growth hormone deficiency (GHD).

Classification of low quartile (LoQ) and high quartile (HiQ) of growth response (height velocity, cm/year) over each of 5 years of therapy with r-hGH (Y1-Y5) was performed in GHD patients and TS patients. Gene expression associated with growth response was determined using rank regression (p < 0.01) and Partial Least Squares Discriminant Analysis (PLS-DA) with two components (X-variate 1 and 2) was used to visualise response groups; PLS-DA is an analytical approach that determines the similarity between individual patients whilst maximising the difference between patient groups. Low quartile (green) and high quartile (red) compared to the rest of the data (orange) are shown for first year growth response to r-hGH in GHD (N = 71, 330 gene probesets with rank regression p < 0.01). Similarity between samples is represented by their proximity. The star plot shows sample distance from the centroid, the arithmetic mean position of all the points in each group (colour figure online).

Fig. 4. Overlap of the core interactome models of height velocity-related gene expression in GHD and TS.

Interactome models were generated from the gene expression associated (p < 0.01) with the height velocity at each year of the study. The functional hierarchy of gene interaction modules within the interactome models was determined using the Moduland algorithm, and the core of the interactome model was defined as the unique sum of the top ten elements of the modules as ranked by network centrality. The overlap of the core of the interactome models between GHD and TS was then determined and visualised as a Venn diagram.

Fig. 5. Network structure of the common core network module shared in patients with growth hormone deficiency (GHD) and Turner syndrome (TS) related to response to recombinant human growth hormone (r-hGH).

A Similarities in the interactome models of the response of GHD and TS to r-hGH were identified by overlap at each year of therapy. Genes were selected that were significantly related to growth response in either or both GHD and TS. The genes related to each year of therapy were combined into a set of 58 uniquely identified genes and this set was used to generate an interactome module (reactome plugin for Cytoscape 3.6.0). Genes with a dark border also have a genetic association with growth response in either GHD or TS. Connecting lines represent known protein:protein interactions, size of the node is proportional to the number of connections made. B The clustering coefficient of the group of genes in the network module associated with each year of therapy was determined and presented as a histogram (average ± standard error of the mean). The clustering coefficient measures the tendency of nodes to cluster together within a network. C The correlation coefficient linking gene expression with growth response at each year of therapy was mapped to the network model, red = positive correlation, green = negative correlation. Genes with a thick border also have a genetic association with growth response in either GHD or TS (colour figure online).

Fig. 6. Predictive value of an identical set of blood gene expression markers identified by network analysis in the classification of response to recombinant human growth hormone (r-hGH) in patients with growth hormone deficiency (GHD) and Turner syndrome (TS).

First year growth response is used as an example. Similarities in the interactome models of the response of GHD and TS to r-hGH were identified by overlap at each year of therapy. Genes were selected that were significantly related to growth response in either or both GHD and TS, generating an identical set of gene probesets used for prediction of both high and low response in both GHD and TS. BORUTA, an all relevant feature selection wrapper random forest-based algorithm, was used to confirm the importance of gene expression probe-sets used for classification of response to r-hGH. The BORUTA algorithm uses a 100-fold permutation to define the noise present in the data; the noise is modelled as shadow variables and used as a basis to assess confidence in the data. Green = confirmed gene probeset, yellow = tentative gene probeset, red = rejected gene probeset, blue = shadow variables (high, medium and low shadow variables are derived to define the noise within the dataset). Low quartile (left column—LoQ) and high quartile (right column—HiQ) are shown for first year growth response to r-hGH in GHD and TS. The same group of gene probesets are used in each case (colour figure online).

Fig. 7. Gene level summary of DNA methylation in GHD patients is related to growth response as measured by Knemometry.

Whole epigenome measurements of six GHD patients with growth response after 4 days of r-hGH therapy measured by knemometry were available from previously published data (GSE57107). A gene level summary of DNA methylation was conducted using median values in Qlucore Omics Explorer (version 3.3) (n = 20,618). A Rank regression of whole-genome DNA methylation against growth response after 4 days of r-hGH as measured by knemometry (p < 0.01) found 497 genes with differential methylation the majority of which showed increased methylation at low rates of growth (negative correlation). B Whole-genome methylation in the six GHD patients ordered by growth response in the sets of genes identified as predicting response to r-hGH in the first year of therapy in both GHD and TS.