The LMSz method - an automatable scalable approach to constructing gene-specific growth charts in rare disorders

Statistical analysis

Growth centile charts are conventionally constructed using the LMS (lambda-mu-sigma) method (8), a special case of the family of Generalized Additive Models for Location Scale and Shape (GAMLSS) (9). Using anthropometry from a reference sample of children the LMS method estimates three smooth curves in age. M is the median or 50^th centile curve; the S curve is the coefficient of variation or fractional standard deviation (SD), and the L curve is a measure of skewness. Any required centile curve can be constructed from the LMS curves, and by reversing the process data for individuals can be converted to exact z-scores.

The LMS method is most effective with large datasets, ideally many thousands of individuals (10). The issue with charts for rare diseases is the inherent scarcity of data. In this case, the chart cannot be based solely on the available data, it must “borrow strength” from previous knowledge about how children grow, as contained in an existing or “baseline” growth reference. This baseline reference can then be adjusted using the rare disease data to obtain a reference for the rare disease. This is analogous to a Bayesian analysis where prior information (the baseline reference) is combined with rare disease data to give a posterior or updated reference relevant to the rare disease.

The process is as follows: the data are first converted to z-scores using the baseline reference, which here is the British 1990 reference (UK90) (11), the official UK chart for birth and age 4-20 and the reference used by the DECIPHER database. Switching to the z-score scale creates centile curves which are horizontal lines, as each centile curve connects points with the same z-score; they are also the same for both sexes. For representative children the z-transformed centiles are Normally distributed with mean (and median) 0 and SD 1 at all ages. Both mean and SD can then be viewed as straight lines plotted against age. Other centile curves are lines above or below the median, spaced according to the corresponding z-score (e.g. the 97.7^th centile is at z = 2, i.e. 2 SDs above the mean).

The next stage is to assume that the rare disease median curve is also a straight line but estimated from the rare disease data as the linear regression of UK90 z-score on age and sex. The SD curve is estimated in the same way from the linear regression of the z-score SD on age and sex. GAMLSS is used to estimate the mean and SD regression lines simultaneously, using the Normal or NO distribution family. The two lines are analogous to the M and S curves of the LMS method, but on the z-score scale. Note that there is no equivalent L line as the normal distribution is by definition not skew, and the conversion from measurement to z-score is assumed to have removed any skewness.

The GAMLSS regression equation for the M line is made up of four terms: intercept, age trend, mean sex difference and age by sex interaction (i.e. the age trend differing between the sexes), and the coefficient for each term is estimated: But if one or more of the fitted coefficients is small enough, the corresponding term(s) can be dropped from the regression, leading to a set of five progressively simpler equations: The SD or S line is modelled in the same way, except it is log transformed to ensure positivity. The interaction term is also omitted for simplicity. Note that the final equations in these lists are null models where the coefficient is 0. They correspond to M = 0 and S = 1, the mean and SD for the baseline UK90 reference, which are also the default values for the GAMLSS NO distribution. GAMLSS estimates the M line and S line simultaneously, and there are 6 M line equations and 5 S line equations, so there are 6x 5 = 30 possible model combinations. All 30 are fitted to the data in turn, and they are compared for goodness of fit using either the Akaike information criterion (AIC) or the Bayesian information criterion (BIC), the optimal model being the one with the lowest AIC or BIC. Note that the S line is linear on the log scale, but backtransformed it is slightly curved if the b_S coefficient is non-zero.

Lack of fit is measured by the deviance, or -2 log likelihood. The AIC penalises the deviance by adding two units for each coefficient in the model, while the BIC uses a larger penalty of log n units per coefficient, where n is the sample size. This penalises more complex models, which tend to fit better.

For a model where n = 7 the AIC and BIC are the same, but for larger n the BIC penalty is larger. For example, the most numerous measurement is weight in ANKRD11 where there are 488 points and log (488) ≈ 6, so here the BIC penalty is three times the AIC penalty. For this reason, the BIC models are either the same as or simpler than the AIC models.

For larger samples the set of models for the M line can be extended to include a cubic P-spline curve in age as an alternative to a straight line. This increases the number of possible fitted M lines from 6 to 9, and the total of combined M+S models from 30 to 45.

Once the optimal model has been identified, it is used to predict the required centile lines on the UK90 z-score scale. They are then backtransformed to the measurement scale using the UK90 reference, giving a set of centile curves appropriate for the rare disease.

Error-checking

To detect possible data outliers, the data in the form of UK90 z-scores were centred and scaled to gene-specific z-scores, and z-scores exceeding 3.5 in absolute value were excluded from the analysis. In this way several data errors were identified and either corrected or excluded.

For simplicity the analysis treats any repeated measures in the data as independent.

Specificity

To test the specificity of the method, z-scores of measurements with ANKRD11 (being the most numerous) were replaced by z-scores randomly sampled from a standard Normal distribution, with mean 0 and SD 1. In this way simulated z-score data were generated with the age and sex structure of the underlying data, where the optimal model ought to be M =0 and log (S) = 0, corresponding to the baseline UK90 reference. The model was then fitted repeatedly to newly sampled data, and the optimal model each time noted, 100 times each for height (n = 355 points), weight (n = 488), BMI (n = 343) and OFC (n = 231).

Z-score centile charts and growth charts

The results are presented as a rectangular array of centile charts by gene (ANKRD11, ARID1B, ASXL3, DDX3X, KMT2A and SATB2) for height, weight, BMI, and OFC (Figure 1). Three centiles are displayed, corresponding to z-scores -2, 0 and 2, i.e. the 2.3^rd, 50^th and 97.7^th centiles. For simplicity they are referred to here as the 2^nd, 50^th and 98^th centiles.

• Download figure
• Open in new tab

Figure 1:

mean ±2 SD lines representing the 2^nd, 50^th and 98^th gene-specific z-score “centiles” for height, weight, BMI and OFC across the genes ANKRD11, ARID1B, ASXL3, DDX3X, KMT2A and SATB2. The corresponding UK90 z-score centiles at -2, 0 and 2 are shown as dashed lines.

In the five plots where there is a sex difference the centile lines are coloured red for females and turquoise for males. The raw data are also shown, colour coded similarly. The ANKRD11 charts use the data from all four countries, while Supplement 2 discusses the appreciable inter-country differences in ANKRD11 height.

Figure 1 penalises the plots using the BIC, while Supplementary Figure 2 shows the same plots penalised with the AIC. As expected, the two sets of centiles are broadly similar, though 13 of the 24 plots in SF2 show a sex difference as against five in Figure 1. Similarly, Supplementary Figure 3 shows the BIC plots of Figure 1 with the option of a spline curve for the M line. Five of the 24 facets select a spline curve, and in all cases they show a dip in early life compared to later childhood.

The specificity of the method was assessed with the ANKRD11 data. Using the BIC models, respectively 98, 97, 97 and 97 times out of 100 the selected model was the null model, i.e. appropriate for the baseline growth reference. This corresponds to an average specificity of just over 97%, i.e. the percentage of times the procedure correctly chooses the baseline reference for data from non-syndromic children. For comparison, repeating the exercise using the AIC models gave a much lower specificity of 53%.

Overall, the differences between Figure 1 and Supplementary Figures 2 and 3 are small, and given their lower specificity there is no obvious reason to prefer the more complicated AIC centiles over the BIC-based centiles, so the focus from here on is the BIC centiles of Figure 1.

Supplementary Table 1 shows the regression coefficients and standard errors for the models in Figure 1. Except for weight and BMI for SATB2, all the models include one or more significant coefficients, some of them highly so, despite the small sample sizes involved. Conversely no model includes an age by sex interaction even though it was an option. Altogether there are 53 regression coefficients fitted, which with 24 distinct gene-measure models correspond to just over two coefficients per model. Supplementary Table 2 shows the 74 corresponding coefficients for the AIC models.

Figures 2 and 3 show the 2^nd, 50^th and 98^th measurement centiles obtained by back-transforming the z-score centiles of Figure 1, for females and males respectively, as a rectangular 4 x 6 plot by measurement and gene, with the corresponding UK90 centiles shown as dashed lines.

• Download figure
• Open in new tab

Figure 2:

measurement centiles for females by measurement and gene obtained by back- transforming the z-score centiles shown in Figure 1.

• Download figure
• Open in new tab

Figure 3:

measurement centiles for males by measurement and gene obtained by back-transforming the z-score centiles shown in Figure 1.

Gene-specific charts

In 19 of the 24 panels in Figure 1, the centiles for males and females are the same, including all four measurements for ANKRD11, ARID1B and KMT2A. The sex differences where they appear are generally small. None of the genes show an age trend in median OFC, and only KMT2A shows a (rising) age trend in median height. But for weight and BMI two and four genes respectively show rising trends, by up to +2 SDs from birth to age 18. This may reflect increasing adiposity with age, as the UK90 reference corresponds to child adiposity as seen in 1990, and obesity prevalence increases with age (12). In 20 of the panels the three centiles are parallel lines, indicating that the variability is constant across age. The most striking exception is for weight with DDX3X, where the variability is dramatically greater at age 18 than at birth.

Supplementary Table 3 shows the UK90 centiles at birth corresponding to the 2^nd, 50^th and 98^th centiles for each measurement, averaged across sex, by gene.

50 centile

The median for SATB2 (height, weight and BMI) and for DDX3X (height and BMI) correspond to the UK90 median. However for all other measurements and genes the median is much lower than reference, ranging from 0.08th (height KMT2A) to 28^th (BMI ANKRD11) UK90 centile.

98th centile

The 98^th centiles are nearly all above the UK90 90^th centile. However 98^th centile for height in KMT2A is on the 44^th UK90 centile. Very few KMT2A cases achieve even median height on the UK90 chart. Other relatively low 98^th centiles are for height in ASXL3 (79^th) and weight in DDX3X (76^th).

2nd centile as z-score

The 2^nd centiles at birth are shown as UK90 z-scores because many are far below the 1^st UK90 centile. The most extreme cases are KMT2A (-6 SD) and ARID1B (-5 SD), while all the 2^nd centiles are on or below the UK90 2^nd centile.

Mowat Wilson Syndrome (MWS) – ZEB2

Given the large sample size, we estimated MWS centiles in two ways: the LMSz method as applied to the z-scores (Figure 4A), and the LMS method as applied to the raw data (Figure 4B). On the z-score scale (Figure 4A) all four measurement medians are below the UK90 median at birth, with the OFC median close to the UK90 2^nd centile. The medians all fall further with increasing age, and at faster rates for weight and OFC among the females. Figure 4B compares the 2^nd, 50^th and 98^th centiles by the two methods, for the sexes separately, where the solid curves correspond to the LMS method while the dashed curves are for the LMSz method. For height, weight and BMI there is reasonable agreement, though the BMI 98^th centiles agree less well. The height and weight centiles are consistent by sex, being similar until puberty when the males become taller and heavier. However, unlike the LMSz centiles, the LMS centiles continue rising after puberty despite growth ending, because there are too few data at older ages to pull them down. OFC performs poorly, with all three LMS-based centiles rising and falling at different ages, while the LMSz centiles rise linearly for males but fall in later childhood for females.

• Download figure
• Open in new tab

Figure 4:

A) MWS data on the z-score scale, in the same form as Figure 1. B) Comparison of the 2^nd, 50^th and 98^th centiles by the two methods, for the sexes separately - solid curves correspond to the LMS method and dashed curves to the LMSz method.

Sample size

Growth charts typically require many thousands of data points to ensure adequate precision for the outer centiles (10). This does not work for the small samples available in rare disease cohorts. The LMSz method instead considers the universe of all possible simple models, such as all combinations of sex and linear age for the mean and SD, and selects the best model, penalised for complexity. This works because the universe of all possible linear models is small, between 30 and 45 alternative models, and the model selection does not depend on statistical significance.

Fitting of centiles

Due to the larger penalty per extra df, the BIC models ought to be simpler than the AIC models, and indeed they are—the BIC models underlying the 24 facets in Figure 1 involve 53 coefficients, i.e. ∼2 coefficients per facet. In contrast the corresponding AIC models in Supplementary Figure 2 involve 74 coefficients, averaging over 3 per facet. Despite this the centiles in Figure 1 and Supplementary Figure 2 are not materially different, probably because the coefficients are all small. The specificity calculation showed that the AIC models had only 53% specificity, compared to 97% for the BIC models, for selecting the correct null model when the data were random Normal. This shows that AIC tends to choose over-complex models, which in turn supports the use of the BIC models. Note that the sensitivity of the method cannot be tested in the same way as the specificity. It involves knowing in advance the centiles appropriate for the particular genotype, which are by definition unknown.

Z-scores versus growth centiles

The z-score centile plots show clearly the differences in growth pattern for a particular gene compared to UK90. The z-score plots compactly show growth trends over time, and sex differences where present. However, in clinical practice health professionals measure and plot their patients’ growth on centile charts as part of standard clinical care. For this reason, it is important to provide the back-transformed centile charts for each gene, as in Figure 2 and 3, as clinicians will find them more familiar in measurement units.

Comparison with the literature and published charts

There is a published centile chart for Coffin-Siris Syndrome (CSS; OMIM #135900), caused by variants in BAF complex genes (4). Direct comparison with ARID1B is not possible given the greater genetic heterogeneity and differing endpoints (height to age 10; OFC to 36 months). In fact, we did not see a sex difference for weight, but the 50^th centile for height at age 10, and weight and OFC were consistent.

We were able to analyse the MWS dataset in two ways: LMSz as described here and the conventional LMS method applied to the raw data. The two sets of centiles are reasonably similar though the z-score-based centiles look more convincing for height and weight, particularly at older ages. For BMI, the 2^nd and 50^th centiles agree well, but the 98^th centiles are more discordant. This is probably because the z-score conversion applies the UK90 skewness adjustment to the MWS data, but the degree of skewness in the MWS data is different. With right skewness this would tend to affect the upper rather than the lower centiles of the distribution, as seen here.

MWS centiles have already been published by the MWS Consortium(3), but superficially they look different from the centiles in Figure 4. This may be because we excluded different outliers and fitted a different model–see Supplement 1 for our gamlss code.

There are no published growth charts for ANKRD11, KMT2A, DDX3X or ASXL3. Variants in ANKRD11 cause KBG Syndrome (MIM # 148050) (13). Our growth charts are consistent with the published phenotypes. Growth hormone has been trialled as a therapeutic option for a small number of individuals with short stature, a review of which is underway by Dr Ockeloen’s group (Aukema S. et al, submitted). These growth charts will be an essential tool in monitoring children both in determining need for potential treatment and judging its success.

We selected DDX3X syndrome (MIM #300160), a neurodevelopmental disorder predominantly in females (14) to investigate how our charts would work for a very rare gene, where n= 5 for the males. Our charts agree with the literature of fairly normal heights and weights in females, with a

proportion having borderline microcephaly. The phenotype in boys is less well understood and the data are limited. Nonetheless our method enables some form of chart to be produced. Considered with caution this could still provide a useful adjunct in a clinical setting, particularly if the data points are viewable in the chart to alert the clinician to the small sample size.

Alterations in KMT2A cause Wiedemann-Steiner Syndrome (WSS) (MIM # 605130 (15). WSS is associated with short stature in about 60% of individuals, microcephaly in a third, and weight below the 5^th centile in a third. Our charts are concordant with published descriptions depicting median height, weight, and OFC for KMT2A at the reference 2^nd, 2^nd and 25^th centiles respectively.

ASXL3-related disorder (Bainbridge-Ropers Syndrome, MIM # 615485) is associated with normal birth weight but poor postnatal growth due to feeding issues in infancy, which stabilise following feeding intervention (16–18). BMI rises with age which may be explained by sustained feeding intervention or because of dysregulated or impulsive eating behaviours that can develop later in childhood. Feeding interventions should not target reference height (ASXL3 median height = 12^th centile) and our charts will be important for growth monitoring. Median OFC is on the 3^rd centile, consistent with published literature of postnatal microcephaly.

Pre- and postnatal growth restriction, sometimes with microcephaly, can be found in up to 34% of individuals with variants in SATB2 (19, 20) (MIM # 608148). Growth charts were produced recently (1508 data points) but they are hard to compare with ours as they included single variants together with chromosomal microdeletions; and they stopped at age 10 (5). They presented z-score heatmaps which indicate that the microdeletion group was most affected by growth restriction, and this probably explains the difference with our charts.

Limitations

Most individuals were likely Caucasian although ethnicity was not formally recorded. Data were collated from various sources and some outliers may be due to measurement error, though we sought out extreme measurements and either corrected or excluded them. Fewer OFC measurements were available and birth OFC measurements were frequently missing. UK90 was used as the baseline reference for three reasons: most of the data were from the UK (except for the European ANKRD11 and MOWS and International ASXL3 cohorts), UK90 is the official UK growth reference for data at birth and over 4 years of age (21), and the DECIPHER database uses UK90 to display the growth data. However, there are inter-country differences in height for any child – see Supplement 2 – and this is true for neurodevelopmental disorders as shown with ANKRD11 - clinical interpretation needs to take this into account.

The LMSz method converts raw data to z-scores using the baseline reference, and then models the z- scores. It assumes that any skewness in the raw data matches the skewness in the reference, and this may not be the case. If the two skewness patterns differ it will introduce bias, as seen in Figure 4B for BMI, where the 98^th centiles based on the LMS and LMSz methods are very different, particularly in males. BMI is markedly right skew at older ages (i.e. the L value is large and negative), which affects the 98^th centile more than the 2^nd and 50^th centiles.

Implementation

We plan to automate LMSz within the DECIPHER platform (22, 23) using its open access datasets alongside further growth measurements being collated through the GenROC study (7). Gene-specific growth charts will be viewable within DECIPHER. A future development would be to variant subtype specific charts within our proposed pipeline.