From Clinical Phenotype to Genotypic Modelling: Incidence and Prevalence of Recessive Dystrophic Epidermolysis Bullosa (RDEB)

Background

Recessive dystrophic epidermolysis bullosa (RDEB) is an inherited genetic disorder characterized by recurrent and chronic open wounds with significant morbidity, impaired quality of life, and early mortality. RDEB patients lack functional type VII collagen (C7) owing to mutations in the gene COL7A1, the main component of anchoring fibrils (AF) necessary to maintain epidermal-dermal cohesion. While over 700 alterations in COL7A1 have been reported to cause Dystrophic epidermolysis bullosa (DEB), which may be inherited in an autosomal dominant or autosomal recessive pattern,¹ the incidence and prevalence of RDEB are not well defined. To date, the widely estimated incidence (0.2–6.65 per million births) and prevalence (3.5–20.4 per million people) of RDEB has been primarily characterized by limited analyses of clinical databases or registries (Table 1).

Table 1 Reported Incidence and Prevalence of Dystrophic Epidermolysis Bullosa

Of note, patients face significant cost and delays in receiving an accurate and timely diagnosis due to the cost and availability of genetic testing, as well as a limited number of rare disease specialists. In the US, rare disease patients visit an average of 7.3 physicians over 7.6 years before receiving a diagnosis.^2,3 Also, a growing number of population genome sequencing efforts such as the 100,000 Genome Project in the UK are highlighting cases of missed diagnoses and developing more accurate genotype-phenotype correlations,⁴ including mutations associated with various degrees of pathogenicity.⁵ In this context, we explore whether genetic allele frequency estimation may enhance the under-diagnosis of rare genetic diseases generally, and RDEB specifically, towards improving the incidence and prevalence estimates of patients who may benefit from treatment.

Using a genetic modelling approach, we use whole exome and genome sequencing data to estimate the allele frequency of recessive pathogenic variants in the COL7A1 gene. Through ClinVar, the NCBI database of human genome variants and phenotypes, DEB Register, and analyzing premature COL7A1 termination variants, we identified 270 variants with documented clinical significance. From these, we built a model to predict the pathogenicity of previously unclassified variants. Applying the model to variant data in the Exome Aggregation Consortium (ExAC) database and the Genome Aggregation Database (gnomAD) in Appendix A, we identified the aggregate frequency of variants predicted to be pathogenic for RDEB and estimated the incidence of the disease.

Dystrophic Epidermolysis Bullosa (DEB)

Dystrophic epidermolysis bullosa (DEB) is one of four major types of epidermolysis bullosa (EB), a group of genetic disorders of the skin and mucous membranes which arise from the defects of basal keratinocyte attachment to the underlying dermis. The skin layer in which the separation occurs defines the EB type. In DEB, separation occurs in the sublamina densa, which attaches the epidermis to the papillary dermis. The separation results from reduction or alteration of C7, a major component of the anchoring fibrils that mediate dermal-epidermal cohesion.^1,6

Gene and Protein

Structure

The gene that encodes C7 is COL7A1. Each C7 molecule is composed of 3 procollagen α₁ chains; each procollagen chain contains 3 domains: an amino-terminal noncollagenous domain (NC1), a central collagenous, triple-helical domain, and a carboxyl-terminal noncollagenous domain (NC2).

Genetics

Over 700 alterations in COL7A1 have been reported to cause DEB, which may be inherited in an autosomal dominant or autosomal recessive pattern.⁷ Less common genetic mechanisms have also been observed. A patient carrying only one COL7A1 mutation not present in the parents’ peripheral blood leukocytes suggested de novo mutation or parental germline mosaicism. Such cases recurred in at least one family, strongly suggesting parental mosaicism.⁸ In two cases with apparently recessive inheritance, only one heterozygous mutation was detected after sequencing all 118 exons and flanking exon–intron borders, while the parent carrying the identified mutation was unaffected.⁹ Finally, one case of RDEB-GS was homozygous for a frameshift mutation, 345insG, carried by the mother due to maternal isodisomy of chromosome 3.⁸ Most mutations are family-specific, but a few recurrent mutations have been noted.^6,8

RDEB Subtypes

Individuals with RDEB-GS produce little to no C7. They frequently have two COL7AI alleles with premature stop codons, resulting in mRNA decay, absence of mRNA expression, or formation of truncated C7 polypeptides that are structurally defective or rapidly degraded in cells. They often have generalized blistering from birth that results in extensive scarring, and alopecia. Patients with RDEB-GS also suffer from extracutaneous manifestations, resulting in an increased risk of corneal abrasions, mucous membrane blistering, oesophagal strictures, kidney problems and cardiomyopathy. Scarring can lead to difficulty eating, vision impairment, painful stools, and constipation. RDEB patients have a high risk of developing aggressive squamous cell carcinomas, which is associated with early mortality.⁶

Patients with RDEB-other (RDEB-O) have a similar, though less severe, phenotype compared to those with RDEB-GS. They produce some functional, albeit abnormal, C7. They also have a better prognosis, with some affected women capable of giving birth. RDEB-O patients are frequently compound heterozygotes: one allele is a missense, in-frame, or splice site mutation, while the second often contains a premature stop codon.^8,10

Genotype-Phenotype Correlation

RDEB generalized severe cases often result from null mutations. Most reported null mutations are nonsense or frameshift mutations that cause premature stop codons, but others destroy the methionine initiation code or completely disrupt splicing.^11,12 Approximately 12% of RDEB-GS cases are compound heterozygotes for one premature termination mutation and one missense mutations, or two missense mutations.¹²

Premature stop codons (PTCs) also play a significant role in RDEB-O. In one study, 34% of cases had PTC mutations in both alleles.¹² Over half of the remaining cases were compound heterozygotes for a PTC mutation and a missense mutation, primarily glycine substitutions. Glycine or arginine substitutions in the collagenous domain appear to cause RDEB inversa (RDEB-I). Some of these substitutions are specific to RDEB-I, but others have also been reported in non-inversa subtypes.¹³ Cases with one or two missense mutations usually have a milder phenotype than those with a PTC mutation. The presentation may be similar to DDEB in some cases.¹¹ It is not always clear why a case with two PTC alleles may present with the milder RDEB-O phenotype. Exon skipping, resulting in a shortened but functional protein, explains this phenomenon in some, but not all, cases.¹²

Epidemiology

Several authors have reported on the epidemiology of DEB (Table 1).^11,14–17 The incidence of DEB mutations reported in these studies ranged from 0.2 to 6.65 per million live births; the prevalence ranged from 3.0 to 20.4 per million people. All studies used multiple methods of ascertainment, including records from hospitals, dermatology clinics, paediatricians and general practitioners, newspaper advertisements, announcements in patient newsletters, and family networks. The two studies that examined ascertainment by source found that the majority of cases were not followed at referral centres or recognized as having EB by their primary care physician, highlighting the difficulty of complete ascertainment.^18,19

Methods

In this study, we estimate the incidence of recessive dystrophic epidermolysis bullosa cases by leveraging publicly available whole-exome sequencing (WES) and whole-genome sequencing (WGS) data to estimate the allele frequency of recessive pathogenic variants in the COL7A1 gene among healthy adult carriers.

Pathogenic Genetic Variants Observed in Cases

We developed an inventory of reported pathogenic COL7A1 variants from ClinVar²⁰ and the DEB Register,^12,21 an international registry for DEB patients and their COL7A1 mutations. The primary goal was to build a model to predict whether an allele could result in an RDEB phenotype when harboured in a homozygous state or compound heterozygous state with another pathogenic allele. To build our model, we constructed a training data set composed of established nonpathogenic variants in COL7A1 and variants that were clearly associated with some form of dystrophic EB. We classified alleles as pathogenic for RDEB if the majority of sources report the allele as causing or associated with RDEB; DDEB if 50% or more of the sources report the allele causing or in association with DDEB; and DEB, not otherwise specified (NOS) if no sources report inheritance pattern.

We considered variants as causal if they were consistently classified as “Pathogenic” or “Likely Pathogenic” in ClinVar or the majority of DEB Register participants that carried the allele had a dystrophic EB phenotype and there was no conflicting evidence in ClinVar. We identified a total of 86 variants that cause DDEB, 155 variants causing RDEB, and 15 DEB-NOS variants. All variants that were removed from the training set due to conflicting evidence annotations were retained in the test set and subsequently evaluated for their potential pathogenic influence with respect to RDEB.

Population Genetics Database

In order to estimate the allele frequency of pathogenic variants, we leveraged the Exome Aggregation Consortium (ExAC) database and the Genome Aggregation Database (gnomAD).²² ExAC contains over 60,000 genotypes from 14 datasets and is ethnically diverse (Appendix A). The only individuals excluded from the final database are those with severe pediatric diseases, such as DEB and other genetic diseases, and those without adequate informed consent. The inclusion criteria allow for the detection of heterozygous carriers of COL7A1 variants, but should exclude affected individuals with homozygous alleles for pathogenic variants. Therefore, COL7A1 variants that are homozygous within ExAC samples are unlikely to be causal for RDEB.

The ExAC database is a subset of gnomAD. GnomAD contains 123,136 exomes and 15,496 genomes from healthy human donors. As the first degree relatives of individuals with severe congenital disorders are excluded in gnomAD, using it risks a downward bias in the estimated allele frequency of pathogenic COL7A1 variants. Nonetheless, as RDEB is a rare disorder, the likelihood that any carrier has a first degree relative with the disease is still quite low. Therefore, while we preferentially report estimates derived from ExAC when a pathogenic variant is listed in both databases, pathogenic and predicted pathogenic variants found only in gnomAD are still reported.

Predicted Pathogenic Alleles

We predicted the pathogenicity of previously unclassified missense variants observed in whole-exome sequences using two automated classifiers, one for missense mutations and one for frameshift mutations. The missense classifier predicts the clinical impact of amino acid substitutions caused by missense mutations by synthesizing the predictions from six variant functional prediction tools: FATHMM,²³ MutationAssessor,²⁴ PolyPhen-2,²⁵ PROVEAN,¹⁸ SIFT,²⁶ and CONDEL.²⁷ Variants for training the model were selected from variants with known clinical relationships to dystrophic EB as described in section 2.1.

We first compiled a training set of 256 known pathogenic and 14 known benign missense mutations based on well-established pathogenic status from ClinVar²⁸ or the DEB Register.²⁹ We applied six prediction tools to these variants, and from their scores fitted a multivariable logistic regression model to predict the disease state of a hypothetical homozygous carrier for a given variant (Figure 1).

Figure 1 Model Classification of Training Data. The 3D scatter plot shown here demonstrates the ability of the classifier and three components to stratify the variants into similar functional clusters.

Based on the clustering of pathogenic variants and distribution of prediction scores (Figure 2) in the training database, we classified variants as pathogenic if the prediction score was 0.95 or higher. With this cut-point, the model correctly classified 14 of 14 benign variants and 225 of 256 pathogenic variants, resulting in a positive predictive value of 100%, sensitivity of 88% and specificity of 100%. To guard against the final estimate of net pathogenic allele frequency being biased upward by likely false positives we filtered our list of previously reported and newly predicted pathogenic variants to remove variants for which homozygous carriers are observed in either ExAC or gnomAD (both of which should only contain unaffected individuals), variants with low confidence loss of function flags in either sequencing database, and any variant with an allele frequency greater than 0.0001 as it would suggest a disease prevalence that is unrealistically high for this rare disorder.

Figure 2 Distribution of Pathogenicity Scores in Training Data.

The splice site classifier predicts the impact of non-coding variants on intron-exon splicing. The classifier is based on a logistic regression model, which integrates the functional annotation scores from the CADD¹⁹ and EIGEN²⁸ algorithms into a composite score. As with the missense classifier, the predictions from the splice site model are refined using a maximum cutoff 0.0001 for allele frequency. We applied the model to a training dataset of 16 known pathogenic variants and 17 known benign variants. We set the threshold of pathogenicity conservatively at a composite score of 0.8, resulting in a sensitivity of 0.44 and a specificity of 1.

We estimated the net carrier allele frequency as the sum of the individual pathogenic allele frequencies from the ExAC database, or when a variant was not found in ExAC, from the larger gnomAD database. We calculated the frequency estimates for pathogenic alleles using all available genotypes in each respective database, as we found no differences in the net carrier allele frequency by ancestry.

The assumption was made that individuals with two null alleles would have RDEB-GS and those with any other two pathogenic alleles would have RDEB-O. We estimated genotype frequency from observed allele frequencies assuming homozygote/heterozygote proportions expected under Hardy–Weinberg equilibrium.

The incidence of RDEB was estimated as the number of births in the US per year times the total expected frequency of pathogenic genotypes. Next, the number of cases per group were summed to get the total number of cases for a birth cohort.²⁹ We estimated prevalence using the estimated number of cases born in a birth cohort and the cumulative probability of death for RDEB subtype (Table 2).³⁰ The calculations were done using cases born in 1960 or later.

Table 2 Cumulative Probability of Mortality by Age and RDEB Subtype

Results

Frequency of RDEB Variants

The ExAC and gnomAD databases contained 1620 COL7A1 exonic variants that cause amino acid substitutions absent in our training set. In total, 523 variants were classified as pathogenic for RDEB, encompassing 193 previously reported pathogenic variants, including 5 dominant mutations. The variants included 128 premature termination codons, 323 missense and 67 splice site mutations. The vast majority of both known and predicted pathogenic alleles had frequencies of less than the 0.0001 inclusion criteria (Figure 3). In all, we excluded 12 variants from our final inventory of pathogenic variants (2 previously reported and 10 newly predicted pathogenic alleles) whose observed allele frequencies exceeded this cutoff.

Figure 3 Distribution of Allele Frequencies.

After adjusting the model to account for the overrepresentation of pathogenic variants in the training dataset, variant classification was more conservative and less sensitive to the threshold used to classify a variant as pathogenic (Table 3). To adjust for the pathogenic over-representation in the training set, we first estimated the true proportion of pathogenic variants in the population. We estimate this proportion as a function of the mean pathogenicity scores as follows:

Table 3 Allele Frequency of DEB Pathogenic Variants in Col7A1 Gene

Where

µ_b = mean composite score of the benign variants;

µ_p = mean composite score of the pathogenic variants;

µ_u = mean composite score of the unclassified variants; and

p = proportion of pathogenic variants in the unclassified set

We then performed weighted logistic regression similar to the methods previously described by Rose et al³² Since we are controlling for an overrepresentation of variants as opposed to patients, we modified Rose’s approach by treating pathogenic variants as cases, benign variants as controls, and the weight for each variant as the ratio of the proportion of the variant class in the unclassified set to the proportion of the variant class in the training set. After adjusting the model to account for the overrepresentation of pathogenic variants in the training dataset, variant classification was more conservative and less sensitive to the threshold used to classify a variant as pathogenic (Table 3). The adjustment also resulted in less dispersal in the distribution of the pathogenicity scores.

We further explored whether the performance of the bias-corrected model varied due to the choice of benign and pathogenic variants in the training datasets using 3-fold cross-validation. We divided the training datasets of known benign and pathogenic variants into three distinct subsets. In each of the three cross-validations, two subsets formed the training dataset and one subset formed the testing dataset. Because the number of benign variants (n=14) and pathogenic variants (n=256) are not divisible by three, 2 of 3 cross-validation testing datasets will contain five benign variants and one testing dataset will contain only four benign variants, with an additional pathogenic variant. We conducted 100 iterations of the 3-fold cross-validation, maintaining the same weights and pathogenicity score cutoff (0.95) for each cross-validation. For each iteration of the 3-fold cross-validation, we calculated the sensitivity, specificity, positive predictive value and negative predictive value of the model. We also calculated the average sensitivity, specificity, positive predictive value and negative predictive value across all 100 iterations.

Approximately 4 million infants are born in the United States annually.³³ Based on this number of births and the frequency of pathogenic allele variants, we estimate that 251 infants with RDEB are born each year in the US (Table 4). Of these, 10 are predicted to have the severe form of RDEB. During the last decade, the European Union averaged a total of 5.23 million births per year.³³ Based on the above frequency of RDEB pathogenic alleles, an estimated 326 cases of RDEB are born in the European Union each year, of which 13 are predicted to be RDEB-GS (Table 4).

Table 4 Estimated Incidence of Recessive Dystrophic Epidermolysis Bullosa

Applying life table analysis to the number of cases born per year, an estimated 12,562 individuals affected with RDEB who were born since 1960 are living in the US, and 16,290 are living in the European Union (Table 5). Fewer than 2% of those individuals are predicted to have the RDEB-GS subtype.

Table 5 Estimated Prevalence of Recessive Dystrophic Epidermolysis Bullosa

Simulation

Given the challenges noted in clinical RDEB ascertainment,^18,19 we further performed a Monte Carlo simulation to estimate the expected number of RDEB cases who may benefit from COL7A1-mediated treatments (e.g., gene therapy, gene editing, mRNA) using MATLAB^® and Simulink^®. The Monte Carlo method uses computational algorithms based on iterations of random sampling from a positively skewed probability distribution (i.e., the probability is inversely related to larger population-based inferences versus known validated patient registries) in order to acquire numerical results for the goals of sampling, optimization, and estimations.^34–36 In other words, our population confidence decreased as we increased our reliance on sequencing point estimates, resulting in positive skewness for the probability distribution. Using a minimum amount defined as the mean of RDEB-GS of 343 and a maximum population estimate of 12,562 of any RDEB, we conducted 10,000 repetitions resulting in a mean estimate of approximately 3850 patients in the US who may benefit from COL7A1-mediated treatments in the US.

Discussion

Genetic modelling applied to the whole exome and genome sequencing data resulted in the identification of predicted RDEB pathogenic alleles, from which our estimate of the incidence of RDEB is 95 per million live births, 30 times the 3.05 per million live birth incidence estimated by the National Epidermolysis Bullosa Registry (NEBR). The NEBR may have under-ascertained less severely affected cases, because at least two studies have reported that a majority of the cases in their studies were not evaluated and treated at EB referral centres.^18,19 Alternatively, we may have overestimated the incidence of RDEB. Our estimates are based on allele frequencies, not clinical symptoms. Also, while a robust database of pathogenic mutations is available, there are a limited number of benign missense mutations available²³ to utilize as a training set. The pathogenicity of splice site mutations is also difficult to predict, making them particularly susceptible to misclassification; even though a conservative threshold was used, some benign variants may be misclassified as pathogenic. Our methodology also overestimates the incidence of RDEB for those genotypes that are incompletely penetrant³⁷ and if fetuses with RDEB are at an increased risk of fetal death. DEB is a complex disorder; clinical symptoms manifest along a continuum of mild to extremely severe. Genotype does not always predict either protein function or disease severity well.¹³ Our estimates of incidence and prevalence apply to a range of RDEB phenotypes and may be limited in the ability to estimate the prevalence of specific symptoms.

A further complication to estimating the incidence of RDEB is the reproductive potential of moderately affected patients. Shinkuma⁶ noted that some RDEB-GO patients are capable of giving birth, but we did not find any information on male fertility in RDEB patients. The methodology we used assumes alleles are in Hardy–Weinberg equilibrium. It remains valid if cases of either sex reproduce, as long as this assumption is met. Notwithstanding, consensus panels recognize that RDEB-GO is a diverse group of RDEB subtypes ranging from RDEB generalized intermediate, RDEB pruiginosa, RDEB bullous dermolysis of the newborn, RDEB pretibial, RDEB centripetalis, RDEB inversa and RDEB localized which may be explored further.³⁸

Conclusion

In sum, this study evaluates the incidence and prevalence of RDEB using publicly available whole-exome sequencing and whole-genome sequencing databases. We established a range estimate, as well as a simulation model, of RDEB patients who may benefit from COL7A1-directed treatments. We conclude that genetic allele frequency estimation may enhance the underdiagnosis of rare genetic diseases generally, and RDEB specifically, which may improve incidence and prevalence estimates of patients who may benefit from treatment.