Transcriptome-based variant calling and aberrant mRNA discovery enhance diagnostic efficiency for neuromuscular diseases

Background Whole-exome sequencing-based diagnosis of rare diseases typically yields 40%–50% of success rate. Precise diagnosis of the patients with neuromuscular disorders (NMDs) has been hampered by locus heterogeneity or phenotypic heterogeneity. We evaluated the utility of transcriptome sequencing as an independent approach in diagnosing NMDs. Methods The RNA sequencing (RNA-Seq) of muscle tissues from 117 Korean patients with suspected Mendelian NMD was performed to evaluate the ability to detect pathogenic variants. Aberrant splicing and CNVs were inspected to identify additional causal genetic factors for NMD. Aberrant splicing events in Dystrophin (DMD) were investigated by using antisense oligonucleotides (ASOs). A non-negative matrix factorisation analysis of the transcriptome data followed by cell type deconvolution was performed to cluster samples by expression-based signatures and identify cluster-specific gene ontologies. Results Our pipeline called 38.1% of pathogenic variants exclusively from the muscle transcriptomes, demonstrating a higher diagnostic rate than that achieved via exome analysis (34.9%). The discovery of variants causing aberrant splicing allowed the application of ASOs to the patient-derived cells, providing a therapeutic approach tailored to individual patients. RNA-Seq data further enabled sample clustering by distinct gene expression profiles that corresponded to clinical parameters, conferring additional advantages over exome sequencing. Conclusion The RNA-Seq-based diagnosis of NMDs achieves an increased diagnostic rate and provided pathogenic status information, which is not easily accessible through exome analysis.


RNA-seq variant calling and prioritization
RNA-seq variant calling was performed after GATK's best practices with modification. STAR aligner (v.2.7.1a) with two-pass mode was used to map raw reads into the hg19 reference genome.
Picard tools marked duplicate reads in the resulting BAM files. GATK's (v.3.7.0) SplitNCigarReads was utilized to extract exon segments and hard clip intronic regions along with reassigning read mapping quality. Base quality recalibration was performed with GATK's BaseRecalibrator.
Variant prioritization was performed with Exomiser (v.12.1.0), after its performance was compared with those of Divine, and DeepPVP (Table S3). First, using GATK Select Variants and bcftools, indels exceeding 70bp and variants with RNA-seq coverage < 5 were excluded. Phenotype information including clinical symptoms, laboratory test results, and pathology reports was converted into the HPO terms. Exomiser was used to further filter variants and prioritize genes.
Resulting variants with allele frequency > 1% in 1000Genomes, TOPMED, UK10K, EXAC, or GNOMAD exomes were removed. Variants annotated in 5'/3'UTR, noncoding transcript exon, upstream/downstream gene, intergenic, or intron were removed. Inheritance filter was set as default (UNDEFINED). Both omimPrioritiser and hiPhivePrioritiser in Exomiser were used to score each gene. We evaluated whether pathogenic gene identified by WES was prioritized within the top 10 candidates ranked by EXOMISER_GENE_COMBINED_SCORE (Table S3).

Aberrant splicing analysis
The LeafcutterMD software was employed to detect aberrant splicing events. RNA-seq BAM files were converted into junction files by quantifying the usage of each intron. Introns in junction files were clustered into intron clusters with default parameters. Outlier intron excision events were detected by modeling intron usage into mathematical model (Dirichlet-multinomial model) and outlier detection process. Cluster P-values were multiple-corrected by Benjamini and Hochberg method and adjusted P-values < 0.05 were selected as a candidate of aberrant splicing events.
Finally, calls residing in muscle disease-related genes were manually inspected and ggsashimi was used for visualization.

Allele specific expression analysis
In order to perform allele specific expression analysis, we adapted the python script presented by Castel et al. In short, heterozygous variants with a coverage depth >20 and a minor allele frequency between 0.3 and 0.7 were extracted from the normalized unfiltered WES VCF file. The positions of these variants were converted into a BED file, and Samtools mpileup was used on RNA-seq BAM files. Information from the WES VCF file and mpileup output file carrying the expression information for all heterozygous exome variants covered in RNA-seq were combined into a final file. Minor allele frequency of variants in RNA seq VCF (MAF < 0.3 or MAF > 0.7) were used to define allelic imbalance. Variants showing allele specific expression patterns were manually examined for clinical relevance.

Gene expression outlier analysis
Gene expression outliers were detected by OUTRIDER package in R. RNA-seq raw counts matrix were combined into one matrix and used as an input file. After removing low expressed genes BMJ Publishing Group Limited (BMJ) disclaims all liability and responsibility arising from any reliance Supplemental material placed on this supplemental material which has been supplied by the author(s)

ASO design, primary myoblast culture, and transfection
Twenty-nine antisense oligonucleotides (ASOs), 20-24 nucleotides (nts) in length, were designed (13 ASOs for CDC_NM54.1 and 16 ASOs for CDC_NM55.1) to target splice donors, splice acceptors, or exonic splicing enhancers of each cryptic exon using the published splicing prediction algorithms (Table S8) hours after transfection, total RNA was isolated and reverse transcribed with oligo dT and random hexamers for subsequent PCR and Sanger sequencing confirmation.

Unsupervised clustering and cell type deconvolution
Non-negative matrix factorization (NMF) was performed following a previous study utilizing R package NMF. The trimmed mean of M-values was used for normalization of the raw counts matrix. The top 25% most-variable genes were selected based on absolute deviation from the median. To optimize cluster number k, 50 random initialization and NMF runs were performed with k ranging from 2 to 20. Stability of clustering was evaluated following the cophenetic correlation coefficient and dispersion coefficient. The chosen cluster number was k = 11 with which 500 random initialization and NMF runs were performed and the best factorization was used for further analysis. Marker genes for each cluster were extracted using a gene scoring algorithm implemented in the NMF package. GO analysis was performed using Goseq. Muscle cell types were deconvoluted using Bisque with input marker genes from a published muscle single-cell RNA-seq study (Rubenstein AB, Smith GR, Raue U et al. Single-cell transcriptional profiles in human skeletal muscle. Sci Rep 2020; 10: 229). R 2 and P-value were calculated from linear regression model for fibrosis/adipose score by fibroadipogenic progenitor (FAP) cell abundance, and reciprocal of (age + 1) by satellite cell abundance. Custom Python and R scripts were used to parse and visualize results.