Advertisement

Research Techniques Made Simple: Whole-Transcriptome Sequencing by RNA-Seq for Diagnosis of Monogenic Disorders

  • Author Footnotes
    3,4 These authors contributed equally to this work.
    Amir Hossein Saeidian
    Footnotes
    3,4 These authors contributed equally to this work.
    Affiliations
    Department of Dermatology and Cutaneous Biology, Sidney Kimmel Medical College, and Jefferson Institute of Molecular Medicine, Thomas Jefferson University, Philadelphia, Pennsylvania, USA

    Genetics, Genomics and Cancer Biology PhD Program, Thomas Jefferson University, Philadelphia, Pennsylvania, USA
    Search for articles by this author
  • Author Footnotes
    3,4 These authors contributed equally to this work.
    Leila Youssefian
    Footnotes
    3,4 These authors contributed equally to this work.
    Affiliations
    Department of Dermatology and Cutaneous Biology, Sidney Kimmel Medical College, and Jefferson Institute of Molecular Medicine, Thomas Jefferson University, Philadelphia, Pennsylvania, USA

    Genetics, Genomics and Cancer Biology PhD Program, Thomas Jefferson University, Philadelphia, Pennsylvania, USA
    Search for articles by this author
  • Author Footnotes
    3,4 These authors contributed equally to this work.
    Hassan Vahidnezhad
    Footnotes
    3,4 These authors contributed equally to this work.
    Affiliations
    Department of Dermatology and Cutaneous Biology, Sidney Kimmel Medical College, and Jefferson Institute of Molecular Medicine, Thomas Jefferson University, Philadelphia, Pennsylvania, USA
    Search for articles by this author
  • Author Footnotes
    3,4 These authors contributed equally to this work.
    Jouni Uitto
    Correspondence
    Correspondence: Jouni Uitto, Department of Dermatology and Cutaneous Biology, Sidney Kimmel Medical College at Thomas Jefferson University, 233 S. 10th Street, Suite 450 BLSB, Philadelphia, Pennsylvania 19107, USA.
    Footnotes
    3,4 These authors contributed equally to this work.
    Affiliations
    Department of Dermatology and Cutaneous Biology, Sidney Kimmel Medical College, and Jefferson Institute of Molecular Medicine, Thomas Jefferson University, Philadelphia, Pennsylvania, USA
    Search for articles by this author
  • Author Footnotes
    3,4 These authors contributed equally to this work.
      Mendelian disorders with cutaneous manifestations comprise a genotypically heterogeneous group of over 1,000 diseases, and in most of them mutant genes have been identified. Mutation detection approaches in these diseases have largely focused on DNA analysis by next-generation sequencing techniques, including gene-targeted sequencing panels as well as whole-exome and whole-genome sequencing. Genome-wide homozygosity mapping (HM), based on DNA polymorphism, has also assisted in the identification of candidate genes in families with consanguinity. However, specific pathogenic variants have not been disclosed in many individual patients when analyzed by next-generation sequencing, and in particular, DNA-based analysis failed to identify many of the mutations impacting on splicing or gene expression. Whole-transcriptome sequencing by RNA sequencing (RNA-Seq), with appropriate bioinformatics, provides a robust tool to identify additional mutations to facilitate genetic diagnosis in genodermatoses. RNA-Seq can be used for variant calling and HM similar to DNA-based approaches, but it also allows for the identification of mutations that result in aberrant transcriptome expression, as displayed by heatmap analysis, and altered splicing patterns of RNA, as visualized by Sashimi plots. Thus, clinical RNA-Seq extends molecular diagnostics of rare genodermatoses, and it could provide a reliable first-tier diagnostic approach to extend mutation databases in patients with heritable skin diseases.

      Abbreviations:

      DEG (differentially expressed genes), EB (epidermolysis bullosa), HM (homozygosity mapping), MAF (minor allele frequency), RNA-Seq (RNA sequencing), ROH (run of homozygosity), WES (whole-exome sequencing), WGS (whole-genome sequencing)

      Introduction

      Monogenic heritable diseases comprise a highly heterogeneous group of as many as 10,000 disorders, and some 1,000 of them have cutaneous manifestations (OMIM: http://www.OMIM.org; World Health Organization: https://www.who.int/genomics/public/geneticdiseases/en/index2.html). In some diseases, the manifestations are present only in the skin, thus being nonsyndromic, whereas in some cases the cutaneous manifestations are associated with a number of extracutaneous manifestations being syndromic (
      • Vahidnezhad H.
      • Youssefian L.
      • Saeidian A.H.
      • Uitto J.
      Phenotypic spectrum of epidermolysis bullosa: the paradigm of syndromic versus non-syndromic skin fragility disorders.
      ). The genetic diagnosis has been established in some of these diseases, many being caused by mutations in several distinct genes. Identification of mutated genes and specific mutations has greatly advanced our understanding of the pathomechanisms of these often complex disorders, and such information can be used for confirmation of the diagnosis with subclassification and prognostication. Such rare monogenic disorders can also serve as simplified models for better understanding of common multisystem disorders (
      • Youssefian L.
      • Vahidnezhad H.
      • Saeidian A.H.
      • Pajouhanfar S.
      • Sotoudeh S.
      • Mansouri P.
      • et al.
      Inherited non-alcoholic fatty liver disease and dyslipidemia due to monoallelic ABHD5 mutations.
      ). Knowledge of the genetic defect also forms the basis for prenatal testing and preimplantation genetic diagnosis. Furthermore, information about the specific mutations is required for potential application of allele-specific treatment approaches that are currently in the developmental pipeline, some of which are already in early clinical trials, for heritable skin diseases (
      • Has C.
      • South A.P.
      • Uitto J.
      Molecular therapeutics in development for epidermolysis bullosa: update 2020.
      ).
      Conventional mutation detection strategies have focused on DNA-based analyses, including next-generation sequencing in the form of gene-targeted arrays or whole-exome sequencing (WES) and whole-genome sequencing (WGS) (
      • Adams D.R.
      • Eng C.M.
      Next-generation sequencing to diagnose suspected genetic disorders.
      ,
      • Bamshad M.J.
      • Ng S.B.
      • Bigham A.W.
      • Tabor H.K.
      • Emond M.J.
      • Nickerson D.A.
      • et al.
      Exome sequencing as a tool for Mendelian disease gene discovery.
      ,
      • Yang Y.
      • Muzny D.M.
      • Reid J.G.
      • Bainbridge M.N.
      • Willis A.
      • Ward P.A.
      • et al.
      Clinical whole-exome sequencing for the diagnosis of Mendelian disorders.
      ). Information derived from these genomic approaches has also been used for homozygosity mapping (HM) to identify putative candidate genes in consanguineous families (
      • Vahidnezhad H.
      • Youssefian L.
      • Jazayeri A.
      • Uitto J.
      Research techniques made simple: genome-wide homozygosity/autozygosity mapping is a powerful tool to identify candidate genes in autosomal recessive genetic diseases.
      ,
      • Vahidnezhad H.
      • Youssefian L.
      • Saeidian A.H.
      • Zeinali S.
      • Touati A.
      • Abiri M.
      • et al.
      Genome-wide single nucleotide polymorphism-based autozygosity mapping facilitates identification of mutations in consanguineous families with epidermolysis bullosa.
      ). However, less than 50% of all cases received genetic diagnosis by DNA-based sequencing techniques, as the DNA analysis does not capture many of the disease-causing variants located in the noncoding regions of the genes or may overlook the consequences of certain types of mutations (
      • Wright C.F.
      • McRae J.F.
      • Clayton S.
      • Gallone G.
      • Aitken S.
      • FitzGerald T.W.
      • et al.
      Making new genetic diagnoses with old data: iterative reanalysis and reporting from genome-wide data in 1,133 families with developmental disorders.
      ). The latter situation is exemplified by synonymous or silent nucleotide substitutions in exons or in noncanonical splicing sequences both impacting on splicing. In such cases, whole-transcriptome sequencing by RNA sequencing (RNA-Seq), with appropriate bioinformatics analysis steps, provides a complementary tool to identify additional mutations to facilitate genetic diagnosis of genodermatoses.
      RNA-Seq is a multifaceted technique that can be used for variant calling and HM, similar to DNA-based approaches. In addition, RNA-Seq is a powerful tool to identify mutant genes with aberrant expression and perturbed splicing patterns (Figure 1). The importance of this ability of RNA-Seq to identify pathogenic sequence variants is emphasized by the fact that (i) splicing defects are among the major causes of Mendelian disorders and they can be located either deeply in intronic sequences, not captured by WES, or as silent variants inside the exons (
      • Chmel N.
      • Danescu S.
      • Gruler A.
      • Kiritsi D.
      • Bruckner-Tuderman L.
      • Kreuter A.
      • et al.
      A deep-intronic FERMT1 mutation causes Kindler syndrome: an explanation for genetically unsolved cases.
      ); and (ii) pathogenic mutations leading to a premature termination codon of translation can cause dramatically reduced gene expression at the mRNA level (
      • Frésard L.
      • Smail C.
      • Ferraro N.M.
      • Teran N.A.
      • Li X.
      • Smith K.S.
      • et al.
      Identification of rare-disease genes using blood transcriptome sequencing and large control cohorts.
      ,
      • Kremer L.S.
      • Bader D.M.
      • Mertes C.
      • Kopajtich R.
      • Pichler G.
      • Iuso A.
      • et al.
      Genetic diagnosis of Mendelian disorders via RNA sequencing.
      ). Thus, RNA-Seq can be utilized for interrogation of culprit genes. Furthermore, filtering of the WES results by frequencies is highly efficient for coding sequences of the gene but does not capture intronic or intergenic variants. Some regions in the genome are difficult to sequence, and RNA-Seq can be helpful to find the causative variants in those regions that are not well characterized by DNA-based genome sequencing methods. Collectively, these observations argue that clinical RNA-Seq is a powerful tool to facilitate and extend diagnostics of rare sequence variants in genodermatoses.
      Figure thumbnail gr1
      Figure 1RNA-Seq technique workflow, the bioinformatics analysis steps, and potential applications of the information. (a) RNA is extracted from whole skin biopsies, dermal fibroblasts, or epidermal keratinocytes and is then sequenced according to protocols and platforms, such as Illumina TruSeq or Takara SMARTer library preparation. The FASTQ file will be analyzed by bioinformatics pipelines. (b) The workflow of bioinformatics analysis includes the steps in the diagram. For details, see where the corresponding steps are color coded for clarity. For details of the software packages and databases, see . (c) Six potential applications of the RNA-Seq technique are shown; they include variant calling, variant prioritization, homozygosity mapping, validation of variants of unknown significance in the genome, quantitative gene expression analysis, and differential gene expression.
      We recognize that the whole-transcriptome sequencing by RNA-Seq with the bioinformatics pipeline is a complex process requiring special expertise and equipment. Therefore, the primary purpose of this summary is to familiarize noncognoscenti with the principles of this technique and alert the readers of the availability of these contemporary approaches for mutation detection in heritable skin diseases. Also, to enhance readability of this review, a glossary is enclosed as Supplementary Material.

      RNA-Seq: Workflow and bioinformatics

      RNA-Seq is initiated by isolation of RNA from tissues or cells that actively express their genome. In skin, RNA can be isolated from a relatively small (3 mm) whole skin biopsy or from cultured cells, such as dermal fibroblasts and epidermal keratinocytes (Figure 1). In this context, the selection of correct cell types is important depending on their gene expression profile. For example, many of the genes involved in cutaneous blistering or cornification disorders are expressed in keratinocytes but not in fibroblasts. Also, it is important to isolate RNA by procedures that preserve the quality of RNA to ensure high quality sequence reads (
      • Vahidnezhad H.
      • Youssefian L.
      • Sotoudeh S.
      • Liu L.
      • Guy A.
      • Lovell P.A.
      • et al.
      Genomics-based treatment in a patient with two overlapping heritable skin disorders: epidermolysis bullosa and acrodermatitis enteropathica.
      ). The biopsies can be transferred to the laboratory in RNA stabilization solution, such as RNAlater transport medium, in which RNA is stable for at least one week at room temperature, one month at 4 °C, and several months at −20 °C. Conveniently, there is no need to freeze samples in liquid nitrogen or rush the samples to the laboratory freezer. RNA is extracted and then subjected to sequencing, initially synthesizing a cDNA library dedicated to RNA-Seq; this approach is different from the conventional cDNA synthesis for Sanger sequencing. RNA-Seq then provides data for bioinformatics analysis. The different steps of data analysis, as shown in Figure 1, consist of mapping the raw data reads to genome and transcriptome references, followed by variant calling and callset refinement, annotation of the variants for HM, and counting and normalization of the reads (Figure 1b). Details of this stepwise process for filtering, including information on bioinformatics software and packages and the output files, are shown in Figure 2 (For details about some examples of the software packages and databases used for data analysis, variant calling and differentially expressed gene [DEG] analysis, see Table 1). The information on these endpoints, specifically variant detection and prioritization as well as HM, is complementary to the information provided by DNA-based techniques. Filtering of the annotated variants for prioritization assists in identification of candidate genes by first focusing on exonic sequence variants and removal of benign synonymous variants with a Combined Annotation Dependent Depletion score, CADD < 20. In the case of rare heritable diseases, filtering to include variants with a minor allele frequency (MAF) < 0.001 and removal of benign variants by bioinformatics prediction programs, followed by alignment of the candidate genes harboring homozygous sequence variants with runs of homozygosity (ROHs), can significantly reduce the number of variants to be considered as pathogenic. As an example, shown in Figure 1c (Variant Prioritization), this approach allowed for the reduction of the number of variants under consideration from 53,035 to 50. The remaining variants, when matched with the clinical phenotypes, identified a single gene (Gene X) as a candidate gene, which was further confirmed by segregation analysis in the family and by modeling of the consequences of the mutation at the protein levels, and further corroborated by datasets at the mRNA level derived from RNA-Seq.
      Figure thumbnail gr2
      Figure 2The flowchart of bioinformatics analysis. Preprocessing step of the raw reads (FASTQ files) includes quality control, trimming the low-quality sequences, align and map high quality reads to the genome, and/or transcriptome reference sequences by programs such as STAR, HISAT2, TopHat, Salmon and kallisto. Two different analyses, including variant calling and extracting DEGs, can be done on the output file. For variant calling by GATK pipelines, the duplicate reads are marked, and the variants are called and filtered. The.vcf output file is annotated with ANNOVAR software for mutation detection, and for homozygosity mapping with the Plink software. The second analysis for identification of DEGs involves a software called SAMtools (Sequence Alignment/Map) for converting SAM to BAM files, followed by assembly of aligned reads by StringTie. The output files are used for visualization of DEGs, or other analysis such as gene set enrichment analysis, and ingenuity pathways analysis (IPA). DEG, differentially expressed gene.
      Table 1Description of Examples of Computer Software Programs and Online Tools for Bioinformatics Analyses of NGS Data and Pathway Analysis
      SoftwareDescription and purposeURLReferences
      BWAA package used to map read sequences to a reference genome.http://bio-bwa.sourceforge.net/(
      • Li H.
      • Durbin R.
      Fast and accurate short read alignment with Burrows-Wheeler transform.
      )
      GATKThe Genome Analysis Toolkit is a multi-purpose variant discovery and genotyping tool.https://software.broadinstitute.org/gatk/(
      • McKenna A.
      • Hanna M.
      • Banks E.
      • Sivachenko A.
      • Cibulskis K.
      • Kernytsky A.
      • et al.
      The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data.
      )
      PicardTools used for manipulating sequencing data in different formats such as BAM and VCF files.http://broadinstitute.github.io/picardRefer to corresponding URL
      PLINKTools used for whole-genome association studies.

      It can be used for homozygosity mapping and IBD estimation.
      http://zzz.bwh.harvard.edu/plink/(
      • Purcell S.
      • Neale B.
      • Todd-Brown K.
      • Thomas L.
      • Ferreira M.A.
      • Bender D.
      • et al.
      PLINK: a tool set for whole-genome association and population-based linkage analyses.
      )
      RA general-purpose software environment for statistical data analysis and visualization.https://www.r-project.org/Refer to corresponding URL
      SAMtoolsTools can be used for manipulating files with SAM/BAM formats.http://samtools.sourceforge.net/(
      • Li H.
      • Handsaker B.
      • Wysoker A.
      • Fennell T.
      • Ruan J.
      • Homer N.
      • et al.
      The sequence alignment/map format and SAMtools.
      )
      FastQCA quality control tool for high throughput sequence data.http://www.bioinformatics.babraham.ac.uk/projects/fastqc/(
      • Andrews S.
      FastQC: A quality control tool for high throughput sequence data
      )
      TrimmomaticA flexible read trimming tool for Illumina NGS data.http://www.usadellab.org/cms/?page=trimmomatic(
      • Bolger A.M.
      • Lohse M.
      • Usadel B.
      Trimmomatic: a flexible trimmer for Illumina sequence data.
      )
      STARAn aligner designed to specifically address many of the challenges of RNA-Seq data mapping using a strategy to account for spliced alignments.https://github.com/alexdobin/STAR/releases(
      • Dobin A.
      • Davis C.A.
      • Schlesinger F.
      • Drenkow J.
      • Zaleski C.
      • Jha S.
      • et al.
      STAR: ultrafast universal RNA-seq aligner.
      )
      HISAT2A fast and sensitive alignment program for mapping next-generation sequencing reads (both DNA and RNA) to a population of human genomes.http://ccb.jhu.edu/software/hisat2/index.shtml(
      • Kim D.
      • Paggi J.M.
      • Park C.
      • Bennett C.
      • Salzberg S.L.
      Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype.
      )
      TopHat2A fast splice junction mapper for RNA-Seq reads.https://ccb.jhu.edu/software/tophat/index.shtml(
      • Kim D.
      • Pertea G.
      • Trapnell C.
      • Pimentel H.
      • Kelley R.
      • Salzberg S.L.
      TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions.
      )
      StringTieA fast and highly efficient assembler of RNA-Seq alignments into potential transcripts.https://ccb.jhu.edu/software/stringtie/(
      • Kovaka S.
      • Zimin A.V.
      • Pertea G.M.
      • Razaghi R.
      • Salzberg S.L.
      • Pertea M.
      Transcriptome assembly from long-read RNA-seq alignments with StringTie2.
      )
      Pheatmap packageA useful R package that implement a heatmaps that offers more control over dimensions and appearance.https://CRAN.R-project.org/package=pheatmapRefer to corresponding URL
      EdgeRA Bioconductor package for differential expression analysis of digital gene expression data.https://bioconductor.org/packages/edgeR/(
      • Robinson M.D.
      • McCarthy D.J.
      • Smyth G.K.
      edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.
      )
      DESeq2A Bioconductor package that provides methods to test for differential expression by use of negative binomial generalized linear models.http://bioconductor.org/packages/release/bioc/html/DESeq2.html(
      • Love M.I.
      • Huber W.
      • Anders S.
      Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.
      )
      EnrichrAn easy to use intuitive enrichment analysis web-based tool providing various types of visualization summaries of collective functions of gene lists/https://amp.pharm.mssm.edu/Enrichr/(
      • Kuleshov M.V.
      • Jones M.R.
      • Rouillard A.D.
      • Fernandez N.F.
      • Duan Q.
      • Wang Z.
      • et al.
      Enrichr: a comprehensive gene set enrichment analysis web server 2016 update.
      )
      GOAn up-to-date and useful database for enrichment and pathway analysis using RNA-Seq data.http://geneontology.org/(
      • Ashburner M.
      • Ball C.A.
      • Blake J.A.
      • Botstein D.
      • Butler H.
      • Cherry J.M.
      • et al.
      Gene ontology: tool for the unification of biology. The Gene Ontology Consortium.
      )
      CytoscapeAn open source software platform for visualizing complex networks and integrating these with any type of attribute data.https://cytoscape.org/(
      • Shannon P.
      • Markiel A.
      • Ozier O.
      • Baliga N.S.
      • Wang J.T.
      • Ramage D.
      • et al.
      Cytoscape: a software environment for integrated models of biomolecular interaction networks.
      )
      Abbreviations: NGS, next-generation sequencing; RNA-Seq, RNA sequencing.
      As noted in the variant prioritization flowchart (Figure 1c), this variant calling would remove synonymous variants that could result in aberrant splicing, including those residing at the very end of exons at the exon-intron border within a canonical splice-site sequence. The confirmation of the pathogenicity resulting in aberrant splicing can be facilitated by RNA-Seq visualized by Sashimi plots demonstrating the pattern of splicing in qualitative and semi-quantitative terms (Figure 1c). As a consequence of such splice junction mutations, several studies have demonstrated exon skipping, partial or complete intron retention, and utilization of alternative splice sequences, which are predicted to lead to alterations in translation, with synthesis of defective protein or absence of the protein expression (
      • Kremer L.S.
      • Bader D.M.
      • Mertes C.
      • Kopajtich R.
      • Pichler G.
      • Iuso A.
      • et al.
      Genetic diagnosis of Mendelian disorders via RNA sequencing.
      ). RNA-Seq utilizing Sashimi plots can also demonstrate a lack of gene expression at the mRNA level in the case of promoter mutations or large gene deletions. DEGs can be quantitated from transcriptome data by heatmap visualization of the expression profiles, which can then direct the attention to the most likely pathogenic gene variants among those under consideration (Figure 1c). It should be noted that gene expression levels in cultured cells can be affected by a number of factors, and analysis of skin biopsies may more accurately reflect the expression profiles in situ.

      Examples of the utility of RNA-Seq in facilitating the identification of mutated genes

      Over the past decade, our laboratory has focused on mutation detection in large cohorts of heritable skin disorders, including epidermolysis bullosa (EB), Mendelian disorders of cornification (ichthyosis and keratodermas), pseudoxanthoma elasticum, and more recently epidermodysplasia verruciformis. In each case, DNA and RNA were isolated after obtaining a written, informed consent by the patient or his/her parents or guardians, who also consented to the publication of the patient’s images (These studies were approved by the institutional review boards of the Pasteur Institute of Tehran, Iran and Thomas Jefferson University, Philadelphia, PA).

      Case 1

      A 1.5-year-old patient with scaly skin, with clinical features consistent with lamellar ichthyosis (Figure 3a), was analyzed for the underlying genetic mutations by WES, which was inconclusive (
      • Youssefian L.
      • Vahidnezhad H.
      • Saeidian A.H.
      • Touati A.
      • Sotoudeh S.
      • Mahmoudi H.
      • et al.
      Autosomal recessive congenital ichthyosis: genomic landscape and phenotypic spectrum in a cohort of 125 consanguineous families.
      ). Subsequently, RNA-Seq of the patient’s skin biopsy revealed a homozygous G>A substitution in the last nucleotide of exon 10 of the TGM1 gene within the canonical splice-site sequence, which changed the wild-type AG-gt to AA-gt (Figure 3b). Because this nucleotide substitution (GAG>GAA) did not change the corresponding amino acid, both codons encoding glutamic acid in transglutaminase 1 protein, conventional bioinformatics filtering of DNA data overlooked this variant and did not prioritize it as a pathogenic variant. However, a Sashimi plot of the RNA-Seq data revealed that this nucleotide change abolished the canonical splice site at the exon 10/intron 10 border, and instead, the patient’s mRNA was processed from two alternate splice sites, the major one residing within intron 10 and the other one within exon 10, resulting in partial retention of intron 10 sequences (Figure 3b). These altered transcripts were predicted to result in frameshift and abolish the synthesis of functional, full-length transglutaminase 1, explaining the patient’s phenotype. Heatmap analysis of the patient’s RNA transcripts, as compared with average transcript levels in three controls, revealed that among the 57 tested genes known to be associated with ichthyosis phenotypes, the level of expression of TGM1 was the lowest (Figure 3c). This is apparently a reflection of the nonsense-mediated mRNA decay because of the splicing mutation disclosed by RNA-Seq.
      Figure thumbnail gr3
      Figure 3Examples of the utility of RNA-Seq in mutation detection in challenging cases of heritable skin diseases. Left panel: Identification and confirmation of a homozygous synonymous/splice-site mutation in TGM1 in a 1.5-year-old patient with lamellar ichthyosis. (a) The patient manifested with scaly skin, sparse hair, and mild erythroderma consistent with diagnosis of lamellar ichthyosis. (b) Sashimi plot of RNA-Seq revealed that a homozygous synonymous mutation results in aberrant splicing and partial intron retention (red), as compared with splicing in control RNA (blue) (upper panel). Screen shot of the genomic sequence visualized by IGV demonstrating a homozygous mutations in TGM1: c.1491G>A at the border of exon 10 and intron 10 (lower panel). (c) Heatmap analysis of the patient’s RNA in comparison with the average quantity of three healthy controls showed that TGM1 had the lowest level of gene expression (red arrow). Right panel: Identification and confirmation of a canonical splice-site mutation in KRT5 in a neonate with EBS. (d) The neonate with severe generalized EBS presented with fragile skin and mucosa. Milia was present at one week of life. The patient died shortly after birth. (e) Sashimi plot of the transcriptome profile of the mutant KRT5 RNA revealed complex aberrant splicing because of the canonical splice-site mutation of c.1440-1G>A that results in retention of intron 7 and intron 8 sequences in patient (red, asterisks in yellow highlighted area) as compared with splicing in control RNA (blue). Screenshot of the genomic sequence visualized by IGV demonstrating canonical splice-site mutation of KRT5: c.1440-1G>A at the border of exon 8 and intron 7 (lower panel). (f) Differential gene expression of 21 genes associated with a blistering phenotype by heatmap analysis revealed that KRT5 was the most downregulated gene among those associated with EB. Permission to publish the patients’ images was provided by the patient and/or his/her parents (guardians). EB, epidermolysis bullosa; EBS, epidermolysis bullosa simplex; IGV, Integrative Genomics Viewer; RNA-Seq, RNA sequencing.

      Case 2

      The utility of heatmap analysis as a tool to guide in the calling of pathogenic variants by RNA-Seq was further illustrated in a neonate with clinical diagnosis consistent with a lethal form of EB (Figure 3d) (
      • Vahidnezhad H.
      • Youssefian L.
      • Daneshpazhooh M.
      • Mahmoudi H.
      • Kariminejad A.
      • Fischer J.
      • et al.
      Biallelic KRT5 mutations in autosomal recessive epidermolysis bullosa simplex, including a complete human keratin 5 "knock-out".
      ). Analysis of the transcriptome data from RNA-Seq revealed a G>A transition in position −1 of intron 7 just preceding exon 8 in the KRT5 gene (Figure 3e). A Sashimi plot revealed that this homozygous mutation rendered the canonical splice site at the intron 7/exon 8 border nonfunctional, and instead an AG sequence six nucleotides upstream within exon 8 was used as an alternate acceptor splice site and led to retention of intron 7 and 8 sequences in the patient’s cells. Heatmap analysis of 21 genes associated with the blistering phenotype in EB revealed that KRT5 expression was the most downregulated among all these genes (Figure 3f). Thus, prioritization of the candidate genes was clearly facilitated by quantitative assessment of the transcript levels by heatmap analysis, further supporting the notion of pathogenicity of the sequence variant in KRT5.

      Case 3

      The utility of RNA-Seq in confirming the pathogenic consequences of nucleotide changes identified by NGS is further illustrated in a 3-year-old patient with EB blistering phenotype (Figure 4a–c) (
      • Vahidnezhad H.
      • Youssefian L.
      • Saeidian A.H.
      • Touati A.
      • Sotoudeh S.
      • Jazayeri A.
      • et al.
      Next generation sequencing identifies double homozygous mutations in two distinct genes (EXPH5 and COL17A1) in a patient with concomitant simplex and junctional epidermolysis bullosa.
      ). Initial DNA sequencing with an EB-associated gene panel of 21 genes revealed two homozygous mutations in this proband from a consanguineous family. One of the mutations, c.5422C>T, in the EXPH5 gene resulted in a nonsense codon p.Arg1808Ter and was predicted to result in the synthesis of a truncated exophilin 5 protein and nonsense-mediated RNA decay (Figure 4a). The second mutation in the COL17A1 gene was a homozygous delA that was predicted to result in a frameshift and termination codon 106 nucleotides downstream of the site of deletion. The mutation in COL17A1 was indeed shown to result in aberrant splicing, and a Sashimi plot revealed partial skipping of exon 4 as well as exon 6, with retention of intron 5 sequences (asterisks in Figure 4b). Heatmap analysis of these two genes placed them on the top of the three most downregulated genes among the 21 genes known to be associated with skin fragility in the spectrum of EB (Figure 4c).
      Figure thumbnail gr4
      Figure 4Benefits and limitations of RNA-Seq in molecular diagnostic settings in confirmation of the pathogenic consequences of mutations. Left panel: Identification of double homozygous mutations in two distinct genes, EXPH5 and COL17A1, in a patient with EB. (a) Next-generation sequencing panel of 21 genes identified homozygous mutations in COL17A1 and EXPH5, which were confirmed by Sanger sequencing. Clinical findings in the proband at 11 months of age consisted of blistering and erosions in the fingers. Immunofluorescence staining for exophilin 5 and type XVII collagen demonstrated complete absence (exophilin 5) and/or a markedly attenuated and discontinuous pattern (type XVII collagen) in the patient's skin, as compared with the control skin. (b) Screenshot of the genomic sequence visualized by IGV demonstrating deletion of a nucleotide within exon 4 of COL17A1: c.202_202delA; this frameshift mutation was predicted to result in synthesis of a truncated polypeptide (p.Thr68LeufsTer106). Sashimi plot of RNA-Seq revealed complex aberrant splicing because of this frameshift mutation. The donor splice site at the 3′ end of exon 3 of COL17A1 utilized the acceptor splice site at the intron 3 and exon 4 border, resulting in partial skipping of exon 4 and 6 (asterisks). (c) Heatmap analysis of the patient’s RNA in comparison with the average quantity of three healthy controls showed that COL17A1 and EXPH5 genes were among the top three most downregulated among the 21 skin fragility–associated genes (Adopted from
      • Vahidnezhad H.
      • Youssefian L.
      • Saeidian A.H.
      • Zeinali S.
      • Touati A.
      • Abiri M.
      • et al.
      Genome-wide single nucleotide polymorphism-based autozygosity mapping facilitates identification of mutations in consanguineous families with epidermolysis bullosa.
      , with permission). Right panel: A missense variant in a novel candidate gene PLOD3 as the cause of syndromic EB identified by WES did not alter transcriptome profiling. (d) The proband manifested with severe scoliosis; joint contractures; and the presence of a tense, hemorrhagic blister on the fifth toe and small erosions on the arm (white arrow) with atrophic scarring at the elbow. (e) WES identified a total of approximately 134,000 annotated sequence variants, which were filtered by the steps indicated to yield eight variants with MAF < 1:1,000 and residing within regions of homozygosity. Matching of the patient’s phenotype identified a mutation in PLOD3, encoding LH3. (f) Immunofluorescence reveals the complete absence of LH3 in the patient's skin with a blister (asterisk), whereas in the control skin, a punctate pattern at the dermal-epidermal junction is noted. Type VII collagen expression is significantly reduced, the remaining protein being primarily in the roof of the blister, as compared with the control skin. (g) Western blotting revealed markedly reduced LH3 protein levels in fibroblasts cultured from the skin of the patient (Pt) as compared with controls (C1–3). Reprobing the filter with an anti–β-tubulin antibody revealed equal protein loading. (h) The LH3 protein levels were quantitated by scanning the bands in western blots in three separate experiments (mean ± SEM; ∗∗P < 0.01, ∗∗∗P < 0.001). (i) Histopathology revealed separation at the dermal-epidermal junction (Richardson's stain). (j) Differential gene expression of 21 genes associated with blistering phenotype by heatmap analysis revealed that PLOD3 gene expression was among the topmost expressed genes. Bar = 50 μm. Permission to publish the patients’ images was provided by the patient and/or his/her parents (guardians) (Adopted from
      • Vahidnezhad H.
      • Youssefian L.
      • Saeidian A.H.
      • Uitto J.
      Phenotypic spectrum of epidermolysis bullosa: the paradigm of syndromic versus non-syndromic skin fragility disorders.
      , with permission). C, control; EB, epidermolysis bullosa; IGV, Integrative Genomics Viewer; MAF, minor allele frequency; Pt, patient; RNA-Seq, RNA sequencing; WES, whole-exome sequencing.
      It is of interest to note that EXPH5 and COL17A1 are associated with two different subtypes of EB, simplex and junctional, respectively. In support of pathogenicity of both mutations, the corresponding genes were both located inside of ROHs derived from RNA-Seq data. The combination of these mutations in this consanguineous family explains the blended phenotype, which initially made subclassification on the clinical basis difficult. This case also illustrates the importance of genome-wide approaches for mutation detection. Specifically, if methodologies that sequence candidate genes one at a time would have been applied, one of the mutations would probably have been missed. This clearly would have impacted the ability to provide accurate genetic counseling to the family in terms of risk of inheritance and its use for prenatal testing and preimplantation genetic diagnosis.
      Collectively, as illustrated by Cases 1–3, RNA-Seq provides a robust tool to explore pathogenicity of heritable skin diseases by identifying variants at both noncoding and coding regions of the gene. RNA-Seq allows prioritization and interpretation of identified sequence variants and provides information complementary to DNA sequencing by transcriptome analysis at the level of RNA splicing and expression. In this regard, uncovering transcriptional consequences of genetic variants allows prioritization or identifies variants that were missed by the applied filters in the bioinformatics pipeline when analyzing DNA-derived data. One of the considerations for RNA-Seq is also its ability to detect genes with monoallelic expression in search of causal variants, particularly for diseases with a recessive mode of inheritance that would not be captured by a single heterozygous variant identified by WES or WGS. Thus, RNA-Seq provides substantial potential to reliably identify pathogenic variants in both known and new disease genes.
      This presentation highlights the utility and limitations of quantitative transcriptome sequencing as visualized by heatmap analysis in identifying candidate genes with particular emphasis on frameshift, loss-of-function, and splicing mutations. It is prudent to point out that these approaches examining the transcriptome levels may not be applicable to missense mutations. Emphasizing this point is Case 4, a 4.5-year-old patient with complex connective tissue disorder, including a syndromic form of dystrophic EB, was shown to harbor a missense mutation, c.1880T>C; p.Leu627Pro, in the PLOD3 gene encoding LH3, a critical enzyme for post-translational modification of collagens (Figure 4d–j) (
      • Vahidnezhad H.
      • Youssefian L.
      • Saeidian A.H.
      • Touati A.
      • Pajouhanfar S.
      • Baghdadi T.
      • et al.
      Mutations in PLOD3, encoding lysyl hydroxylase 3, cause a complex connective tissue disorder including recessive dystrophic epidermolyis bullosa-like blistering phenotype with abnormal anchoring fibrils and type VII collagen deficiency.
      ). The study demonstrated that this missense mutation caused the absence of LH3 at the protein level but also resulted in reduced collagen VII expression and reduction in anchoring fibrils, apparently because of deficient post-translational modification of collagen VII. Although the primary genetic defect was clearly shown to be in PLOD3, the heatmap profile showed that this gene was expressed among the top three highest levels among the genes associated with skin fragility in the spectrum of EB. Although the reason for upregulation of PLOD3 gene expression is not entirely clear, it may reflect compensatory changes in mRNA levels in response to absent protein (Figure 4j). Thus, transcript quantitation of DEGs by heatmap in this case of homozygous missense mutations was not helpful in prioritizing or identifying mutant genes.

      RNA-Seq reveals an unexpected splicing event in COL7A1

      Transcriptome profiling is also able to evaluate tissue- and cell type–dependent expression and splicing profiles using the corresponding material as a source of RNA selected at the beginning of the analysis. This consideration highlights one of the limitations of RNA-Seq, which requires tissues or cells that express the corresponding genes as starting material. In this context, it is important to note that whole skin biopsy consists not only of epidermal keratinocytes and dermal fibroblasts but also has a large number of other cell types with unique expression profiles. This issue is highlighted by the demonstration that COL7A1, which encodes type VII collagen, was shown to have a different splicing pattern in RNA isolated from whole skin as compared with cultured dermal fibroblasts (Figure 5a). Specifically, a Sashimi plot revealed alternative splicing in control skin, which was generated by a different exon 18 acceptor site 27 base pairs upstream from the canonical acceptor site as compared with some, but not all, fibroblast cultures (Figure 5a). This resulted in insertion of nine amino acid residues into the fibronectin type III linker domain of the noncollagenous NC-1 region of type VII collagen (Figure 5b). This splice variant has been shown previously to be differentially expressed in wound edges of patients with epithelizing skin ulcers in patients with EB as compared with normal keratinocytes from steady-state body sites (Figure 5c) (
      • Sawamura D.
      • Goto M.
      • Yasukawa K.
      • Kon A.
      • Akiyama M.
      • Shimizu H.
      Identification of COL7A1 alternative splicing inserting 9 amino acid residues into the fibronectin type III linker domain.
      ). Our data demonstrating differential expression of this alternative spliced insertion in keratinocytes and fibroblasts provide further evidence of the importance of splice variants in disease processes, and such variants can be visualized by RNA-Seq. This demonstration also emphasized previous observations that cell culture conditions can profoundly alter the gene expression profiles and, for example, the expression levels and transcriptome profiles of fibroblasts in culture may not necessarily reflect those in intact skin (
      • Mahmoudi S.
      • Mancini E.
      • Xu L.
      • Moore A.
      • Jahanbani F.
      • Hebestreit K.
      • et al.
      Heterogeneity in old fibroblasts is linked to variability in reprogramming and wound healing.
      ).
      Figure thumbnail gr5
      Figure 5Confirmation of alternative splicing in COL7A1. (a) Comparison of COL7A1 gene expression in different sample types, including normal keratinocytes, fibroblasts, and whole skin biopsies. The Sashimi plot of RNA-Seq revealed in control skin alternative splicing, which was generated by a different exon 18 acceptor site 27 bp upstream from the canonical acceptor site as compared with some of the fibroblast cultures. (b) Screenshot of the genomic sequence visualized by IGV demonstrating the 27 bp retention from intron 17 that would result in insertion of nine amino acid residues into the fibronectin type III linker domain of the noncollagenous NC-1 region of type VII collagen. (c) This alternative splicing has been shown previously to be differentially expressed in wounds of EB patients as compared with normal keratinocytes. The cDNA containing exon 18 and with the 27 bp intron 17 retention is shown by arrow (lower band, exon 18; upper band, exon 18 plus 9 amino acids) (Adopted from
      • Sawamura D.
      • Goto M.
      • Yasukawa K.
      • Kon A.
      • Akiyama M.
      • Shimizu H.
      Identification of COL7A1 alternative splicing inserting 9 amino acid residues into the fibronectin type III linker domain.
      , with permission). bp, base pair; EB, epidermolysis bullosa; IGV, Integrative Genomics Viewer; RNA-Seq, RNA sequencing.

      Clinical utility of transcriptome profiling by RNA-Seq

      RNA-Seq is increasingly becoming a complementary way of identifying genes underlying rare heritable diseases, bypassing hurdles in interpreting intronic or splice-altering variants (
      • Li D.
      • Tian L.
      • Hakonarson H.
      Increasing diagnostic yield by RNA-sequencing in rare disease-bypass hurdles of interpreting intronic or splice-altering variants.
      ). It has been indicated that the genetic diagnostic rate for Mendelian diseases by WES is typically only in the range of 20–40%, in part because of the fact that WES misses deep intronic, silent, or synonymous exonic variants leading to aberrant splicing. The power of RNA-Seq in solving unrecognized pathogenic variants (
      • Hamanaka K.
      • Miyatake S.
      • Koshimizu E.
      • Tsurusaki Y.
      • Mitsuhashi S.
      • Iwama K.
      • et al.
      RNA sequencing solved the most common but unrecognized NEB pathogenic variant in Japanese nemaline myopathy.
      ) has been demonstrated in studies on a number of heritable skin diseases; some of them were highlighted in this review with the focus on EB and ichthyosis, two heterogeneous disorders that have been associated with 21 and over 67 mutant genes, respectively. Finally, it should be noted that although immunohistochemical and fluorescent staining of the skin samples is often informative and indicates the candidate gene through lack of immunoreactivity for a specific protein, this approach does not provide information about the specific mutations (
      • Chmel N.
      • Danescu S.
      • Gruler A.
      • Kiritsi D.
      • Bruckner-Tuderman L.
      • Kreuter A.
      • et al.
      A deep-intronic FERMT1 mutation causes Kindler syndrome: an explanation for genetically unsolved cases.
      ,
      • He Y.
      • Balasubramanian M.
      • Humphreys N.
      • Waruiru C.
      • Brauner M.
      • Kohlhase J.
      • et al.
      Intronic ITGA3 mutation impacts splicing regulation and causes interstitial lung disease, nephrotic syndrome, and epidermolysis bullosa.
      ).

      Conclusions

      The feasibility of RNA-Seq applications for heritable skin diseases is emphasized by the ready availability of skin biopsies and the capability to culture epidermal keratinocytes or dermal fibroblasts for RNA isolation. This situation contrasts with many heritable disorders, such as those manifesting with neurologic or internal organ involvement. Although transcriptome sequencing has been recognized as complementary to DNA-based next-generation sequencing approaches, one could argue that RNA-Seq, capable of variant calling and HM analogous to DNA-based analysis, with additional capability to reveal aberrant gene splicing and determine the gene expression levels visualized by heatmapping, could be considered as an expedient, reliable, and affordable first-tier diagnostic approach for finding mutations in patients with heritable skin diseases.

      Conflict of Interest

      The authors state no conflict of interest.

      Multiple Choice Questions

      • 1.
        What type of sequence variants may not be detected by quantitative whole-transcriptome sequencing by heatmap analysis?
        • A.
          Missense mutations
        • B.
          Nonsense mutations
        • C.
          Canonical splicing site mutations
        • D.
          Frameshift mutations
      • 2.
        Which of the following statements is true when comparing RNA sequencing (RNA-Seq) with whole-exome sequencing?
        • A.
          Germline mutations are detectable both in RNA and DNA samples.
        • B.
          Whole-transcriptome sequencing provides tissue-specific gene expression information and can help in variant prioritization.
        • C.
          For splicing variants, RNA-Seq provides information that can be used for pathogenicity assessment.
        • D.
          All of the above.
      • 3.
        Which statement is NOT true regarding RNA-Seq technique?
        • A.
          Sampling is invasive and, in some cases, impossible.
        • B.
          Sometimes tissue availability is an issue.
        • C.
          RNA, if not preserved in RNAlater, is unstable and readily subject to degradation.
        • D.
          Because of low coverage of exonic sequences, it cannot detect substitutions.
      • 4.
        For rare heritable skin disorders, which type of tissue(s) or cells could be used for mutation detection by whole-transcriptome sequencing?
        • A.
          Fibroblasts
        • B.
          Keratinocytes
        • C.
          Whole skin
        • D.
          All of the above
      • 5.
        Which statement is NOT correct regarding the combination of RNA-Seq and homozygosity mapping techniques?
        • A.
          It is useful to detect germline mutations.
        • B.
          It is possible to find loss-of-function pathogenic mutations in novel candidate genes.
        • C.
          Combination of these methods is a robust approach to detect exonic synonymous mutations affecting splicing.
        • D.
          It is applicable only to the families affected by diseases with autosomal dominant inheritance.See online version of this article for a detailed explanation of correct answers.

      Acknowledgments

      Correspondence regarding the technical aspects of this study should be addressed to HV, whereas those addressing the clinical application should be directed to JU. The authors’ original studies were supported by the United States National Institutes of Health (Grant no. R01AI143810 ) and Debra International , Austria. Carol Kelly assisted in manuscript preparation. This study is in partial fulfillment of the PhD Thesis of Amir Hossein Saeidian.

      Author contributions

      Conceptualization: JU, HV; Data Curation: AHS, LY; Formal Analysis: AHS, LY, HV; Funding Acquisition: JU; Investigation: AHS, LY, HV; Project Administration: JU; Supervision: JU, HV; Visualization: AHS, LY; Writing - Original Draft Preparation: AHS, JU; Writing - Review and Editing: AHS, LY, JU, HV

      Detailed Answers

      • 1.
        What type of sequence variants may not be detected by quantitative whole-transcriptome sequencing by heatmap analysis?
      • Answer: A. All mutations except missense mutations definitively cause quantitative changes in gene expression pattern. Depending on the effect of the missense mutation on protein function, it may or may not alter the gene expression.
      • 2.
        Which of the following statements is true when comparing RNA sequencing (RNA-Seq) with whole-exome sequencing?
      • Answer: D. In the comparison of RNA-Seq and whole-exome sequencing, all mentioned statements are correct.
      • 3.
        Which statement is NOT true regarding RNA-Seq technique?
      • Answer: D. Regarding RNA-Seq technique, all choices are correct except choice D. In the RNA-Seq method, exons that are coding sequences of the gene are covered and all substitutions should be detected (if there is no technical pitfall or low RNA sample quality).
      • 4.
        For rare heritable skin disorders, which type of tissue(s) or cells could be used for mutation detection by whole-transcriptome sequencing?
      • Answer: D. For heritable skin disorders, all mentioned samples could be used for transcriptomic studies.
      • 5.
        Which statement is NOT correct regarding the combination of RNA-Seq and homozygosity mapping techniques?
      • Answer: D. The combination of RNA-Seq and homozygosity mapping can detect germline mutations, loss-of-function mutations, and exonic synonymous mutations affecting splicing. However, it is not useful for diseases with autosomal dominant inheritance.

      Supplementary Material

      References

        • Adams D.R.
        • Eng C.M.
        Next-generation sequencing to diagnose suspected genetic disorders.
        N Engl J Med. 2018; 379: 1353-1362
        • Andrews S.
        • FastQC: A quality control tool for high throughput sequence data
        (2010 (accessed 1 December 2019))
        • Ashburner M.
        • Ball C.A.
        • Blake J.A.
        • Botstein D.
        • Butler H.
        • Cherry J.M.
        • et al.
        Gene ontology: tool for the unification of biology. The Gene Ontology Consortium.
        Nat Genet. 2000; 25: 25-29
        • Bamshad M.J.
        • Ng S.B.
        • Bigham A.W.
        • Tabor H.K.
        • Emond M.J.
        • Nickerson D.A.
        • et al.
        Exome sequencing as a tool for Mendelian disease gene discovery.
        Nat Rev Genet. 2011; 12: 745-755
        • Bolger A.M.
        • Lohse M.
        • Usadel B.
        Trimmomatic: a flexible trimmer for Illumina sequence data.
        Bioinformatics. 2014; 30: 2114-2120
        • Chmel N.
        • Danescu S.
        • Gruler A.
        • Kiritsi D.
        • Bruckner-Tuderman L.
        • Kreuter A.
        • et al.
        A deep-intronic FERMT1 mutation causes Kindler syndrome: an explanation for genetically unsolved cases.
        J Invest Dermatol. 2015; 135: 2876-2879
        • Dobin A.
        • Davis C.A.
        • Schlesinger F.
        • Drenkow J.
        • Zaleski C.
        • Jha S.
        • et al.
        STAR: ultrafast universal RNA-seq aligner.
        Bioinformatics. 2013; 29: 15-21
        • Frésard L.
        • Smail C.
        • Ferraro N.M.
        • Teran N.A.
        • Li X.
        • Smith K.S.
        • et al.
        Identification of rare-disease genes using blood transcriptome sequencing and large control cohorts.
        Nat Med. 2019; 25: 911-919
        • Hamanaka K.
        • Miyatake S.
        • Koshimizu E.
        • Tsurusaki Y.
        • Mitsuhashi S.
        • Iwama K.
        • et al.
        RNA sequencing solved the most common but unrecognized NEB pathogenic variant in Japanese nemaline myopathy.
        Genet Med. 2019; 21: 1629-1638
        • Has C.
        • South A.P.
        • Uitto J.
        Molecular therapeutics in development for epidermolysis bullosa: update 2020.
        Mol Diagn Ther. 2020;
        • He Y.
        • Balasubramanian M.
        • Humphreys N.
        • Waruiru C.
        • Brauner M.
        • Kohlhase J.
        • et al.
        Intronic ITGA3 mutation impacts splicing regulation and causes interstitial lung disease, nephrotic syndrome, and epidermolysis bullosa.
        J Invest Dermatol. 2016; 136: 1056-1059
        • Kim D.
        • Paggi J.M.
        • Park C.
        • Bennett C.
        • Salzberg S.L.
        Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype.
        Nat Biotechnol. 2019; 37: 907-915
        • Kim D.
        • Pertea G.
        • Trapnell C.
        • Pimentel H.
        • Kelley R.
        • Salzberg S.L.
        TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions.
        Genome Biol. 2013; 14: R36
        • Kovaka S.
        • Zimin A.V.
        • Pertea G.M.
        • Razaghi R.
        • Salzberg S.L.
        • Pertea M.
        Transcriptome assembly from long-read RNA-seq alignments with StringTie2.
        Genome Biol. 2019; 20: 278
        • Kremer L.S.
        • Bader D.M.
        • Mertes C.
        • Kopajtich R.
        • Pichler G.
        • Iuso A.
        • et al.
        Genetic diagnosis of Mendelian disorders via RNA sequencing.
        Nat Commun. 2017; 8: 15824
        • Kuleshov M.V.
        • Jones M.R.
        • Rouillard A.D.
        • Fernandez N.F.
        • Duan Q.
        • Wang Z.
        • et al.
        Enrichr: a comprehensive gene set enrichment analysis web server 2016 update.
        Nucleic Acids Res. 2016; 44: W90-W97
        • Li D.
        • Tian L.
        • Hakonarson H.
        Increasing diagnostic yield by RNA-sequencing in rare disease-bypass hurdles of interpreting intronic or splice-altering variants.
        Ann Transl Med. 2018; 6: 126
        • Li H.
        • Durbin R.
        Fast and accurate short read alignment with Burrows-Wheeler transform.
        Bioinformatics. 2009; 25: 1754-1760
        • Li H.
        • Handsaker B.
        • Wysoker A.
        • Fennell T.
        • Ruan J.
        • Homer N.
        • et al.
        The sequence alignment/map format and SAMtools.
        Bioinformatics. 2009; 25: 2078-2079
        • Love M.I.
        • Huber W.
        • Anders S.
        Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.
        Genome Biol. 2014; 15: 550
        • Mahmoudi S.
        • Mancini E.
        • Xu L.
        • Moore A.
        • Jahanbani F.
        • Hebestreit K.
        • et al.
        Heterogeneity in old fibroblasts is linked to variability in reprogramming and wound healing.
        Nature. 2019; 574: 553-558
        • McKenna A.
        • Hanna M.
        • Banks E.
        • Sivachenko A.
        • Cibulskis K.
        • Kernytsky A.
        • et al.
        The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data.
        Genome Res. 2010; 20: 1297-1303
        • Purcell S.
        • Neale B.
        • Todd-Brown K.
        • Thomas L.
        • Ferreira M.A.
        • Bender D.
        • et al.
        PLINK: a tool set for whole-genome association and population-based linkage analyses.
        Am J Hum Genet. 2007; 81: 559-575
        • Robinson M.D.
        • McCarthy D.J.
        • Smyth G.K.
        edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.
        Bioinformatics. 2010; 26: 139-140
        • Sawamura D.
        • Goto M.
        • Yasukawa K.
        • Kon A.
        • Akiyama M.
        • Shimizu H.
        Identification of COL7A1 alternative splicing inserting 9 amino acid residues into the fibronectin type III linker domain.
        J Invest Dermatol. 2003; 120: 942-948
        • Shannon P.
        • Markiel A.
        • Ozier O.
        • Baliga N.S.
        • Wang J.T.
        • Ramage D.
        • et al.
        Cytoscape: a software environment for integrated models of biomolecular interaction networks.
        Genome Res. 2003; 13: 2498-2504
        • Vahidnezhad H.
        • Youssefian L.
        • Daneshpazhooh M.
        • Mahmoudi H.
        • Kariminejad A.
        • Fischer J.
        • et al.
        Biallelic KRT5 mutations in autosomal recessive epidermolysis bullosa simplex, including a complete human keratin 5 "knock-out".
        Matrix Biol. 2019; 83: 48-59
        • Vahidnezhad H.
        • Youssefian L.
        • Jazayeri A.
        • Uitto J.
        Research techniques made simple: genome-wide homozygosity/autozygosity mapping is a powerful tool to identify candidate genes in autosomal recessive genetic diseases.
        J Invest Dermatol. 2018; 138: 1893-1900
        • Vahidnezhad H.
        • Youssefian L.
        • Saeidian A.H.
        • Touati A.
        • Pajouhanfar S.
        • Baghdadi T.
        • et al.
        Mutations in PLOD3, encoding lysyl hydroxylase 3, cause a complex connective tissue disorder including recessive dystrophic epidermolyis bullosa-like blistering phenotype with abnormal anchoring fibrils and type VII collagen deficiency.
        Matrix Biol. 2019; 81: 91-106
        • Vahidnezhad H.
        • Youssefian L.
        • Saeidian A.H.
        • Touati A.
        • Sotoudeh S.
        • Jazayeri A.
        • et al.
        Next generation sequencing identifies double homozygous mutations in two distinct genes (EXPH5 and COL17A1) in a patient with concomitant simplex and junctional epidermolysis bullosa.
        Hum Mutat. 2018; 39: 1349-1354
        • Vahidnezhad H.
        • Youssefian L.
        • Saeidian A.H.
        • Uitto J.
        Phenotypic spectrum of epidermolysis bullosa: the paradigm of syndromic versus non-syndromic skin fragility disorders.
        J Invest Dermatol. 2019; 139: 522-527
        • Vahidnezhad H.
        • Youssefian L.
        • Saeidian A.H.
        • Zeinali S.
        • Touati A.
        • Abiri M.
        • et al.
        Genome-wide single nucleotide polymorphism-based autozygosity mapping facilitates identification of mutations in consanguineous families with epidermolysis bullosa.
        Exp Dermatol. 2019; 28: 1118-1121
        • Vahidnezhad H.
        • Youssefian L.
        • Sotoudeh S.
        • Liu L.
        • Guy A.
        • Lovell P.A.
        • et al.
        Genomics-based treatment in a patient with two overlapping heritable skin disorders: epidermolysis bullosa and acrodermatitis enteropathica.
        Hum Mutat. 2020; 41: 906-912
        • Wright C.F.
        • McRae J.F.
        • Clayton S.
        • Gallone G.
        • Aitken S.
        • FitzGerald T.W.
        • et al.
        Making new genetic diagnoses with old data: iterative reanalysis and reporting from genome-wide data in 1,133 families with developmental disorders.
        Genet Med. 2018; 20: 1216-1223
        • Yang Y.
        • Muzny D.M.
        • Reid J.G.
        • Bainbridge M.N.
        • Willis A.
        • Ward P.A.
        • et al.
        Clinical whole-exome sequencing for the diagnosis of Mendelian disorders.
        N Engl J Med. 2013; 369: 1502-1511
        • Youssefian L.
        • Vahidnezhad H.
        • Saeidian A.H.
        • Pajouhanfar S.
        • Sotoudeh S.
        • Mansouri P.
        • et al.
        Inherited non-alcoholic fatty liver disease and dyslipidemia due to monoallelic ABHD5 mutations.
        J Hepatol. 2019; 71: 366-370
        • Youssefian L.
        • Vahidnezhad H.
        • Saeidian A.H.
        • Touati A.
        • Sotoudeh S.
        • Mahmoudi H.
        • et al.
        Autosomal recessive congenital ichthyosis: genomic landscape and phenotypic spectrum in a cohort of 125 consanguineous families.
        Hum Mutat. 2019; 40: 288-298