Ai-based method for characterization of variants of uncertain significance using RNA-SEQ data

Abstract

A method for constructing a classifier to classify variants of uncertain significance (VUS) is proposed according to the present disclosure, wherein raw data collected from databases comprising public Mendelian and hereditary disease RNA-Seq and / or DNA-Seq data is used for training a random forest model. Training data is based on a feature table created using variant analysis based on difference to reference genome, variant annotation based on at least multiple measures, RNA-Seq expression analysis, as well as RNA-Seq splicing analysis. A system for classifying the pathogenicity of variants of uncertain significance based on the disclosed method is also disclosed according to the present disclosure.

Description

[0001] AI-BASED METHOD FOR CHARACTERIZATION OF VARIANTS OF UNCERTAIN SIGNIFICANCE USING RNA-SEQ DATA

[0002] Technical Field of the Present Invention

[0003] The present invention relates generally to machine learning-based systems that receive genome or multiomics and clinical data input and generate disease classifications or predictions, and more particularly to systems that process genome or multiomics sequences for detecting pathogenicity of genetic variants of uncertain significance.

[0004] Background of the Present Invention

[0005] The human genome is extensive, encompassing over 3 billion bases, and the wealth of experimental data mapping to its positions further compounds its complexity. Even if one were to represent longer k-mers (consecutive loci or substrings) as a single token, handling a deep learning model's complete attention matrix for these tokens becomes utterly unfeasible due to its sheer size. To put it into perspective, this would demand storage capacities in the terabytes or even petabytes to compute, given its intricate nature. Prior approaches aimed to alleviate this complexity by techniques such as clustering or concentrating the attention mechanism on only a subset of the input, thereby reducing the region under consideration and limiting the input sequence to specific loci of interest.

[0006] For instance, if the entire genome sequences of both a human cancer tissue and a normal tissue were served as input into a neural network, one would require a matrix with approximately 3 billion columns and 10 - 20 rows, assuming a dense concatenated (stacked) one-hot encoded format. This calculation considers each of the two or four sets of columns specifying the four nucleotides (A, C, G, T) and a placeholder for unknown values, both for cancer and normal tissues, across all positions and possibly for both pairs of chromosomes for each tissue type. Additionally, various other dimensions (rows) may be required to account for different genomic features, such as methylation sites, microsatellite variability, insertions, deletions, and more. Given that a network trained with stochastic gradient descent necessitates multiple such setups batched together, the memory and computational resources required increase exponentially during training.

[0007] Consequently, previous solutions have been suboptimal, focusing on relatively small network inputs and potentially relying on partial knowledge of the underlying biology to filter input loci. Furthermore, a more robust and veritable method is required for interpreting and classifying clinical variants.

[0008] Clinical class of variants obtained by variant analysis are estimated computationally according to ACMG-AMP criteria or retrieved from publicly accessible databases such as ClinVar, variants wherein have been characterized either experimentally or computationally. However, not all variants are classified as pathogenic or benign. Those non-classified variants are referred to as Variant of Uncertain Significance (VUS hereinafter) in the literature. Variants are prioritized by filtering based on certain criteria such as population frequency, exonic function clinical class of the variants and then interpreted by experts. VUSs are often not reported by clinicians due to prioritization and interpretation problems, and as such they are the one of the reasons for the low diagnostic rate, especially in rare hereditary diseases. In recent years, many studies have been published demonstrating that VUSs can be prioritized using expression and splicing information obtained from RNA-Seq data. Next Generation Sequencing (NGS) is a massively parallel and high-throughput sequencing technology. The NGS method has facilitated the sequencing of large amounts of genomic data in recent years. Whole Genome Sequencing (WGS), in which the entire DNA of an individual is sequenced, and Whole Exome Sequencing (WES), that targets only protein coding regions of the genome, are popular DNA-based NGS methods. With the decrease in sequencing cost, these two DNA-based sequencing methods, mostly WES, have started to be used in the routine clinical practice. The purpose of these two DNA-based methods is to detect variants that are disease related. Variants are defined as differences between a representative sequence called reference genome and a sample. The raw data obtained via sequencing methods is processed with bioinformatic algorithms and a list of variants is generated. Information whether these variants are disease-related (pathogenic) or benign, in other words the information about clinical class of the variants, are obtained in the annotation step according to ACMG-AMP criteria as in the 2015 study titled "Standards and guidelines for the interpretation of sequence variants: a joint consensus recommendation of the American College of Medical Genetics and Genomics and the Association for Molecular Pathology" by Richards et al. or using databases such as ClinVar as set forth in the 2014 study titled "ClinVar: public archive of relationships among sequence variation and human phenotype" by Landrum et al.

[0009] ACMG-AMP criteria, published by the American College of Medical Genetics and Genomics (ACMG) and the Association for Molecular Pathology (AMP), consists of 28 criteria whereby pieces of information such as the frequency of the variant in the population, the prediction scores published in the literature (SIFT, PolyPhen etc.), the effect on the protein are considered to evaluate variants' clinical class in a computational sense. Proposed criteria are specific to variants detected in inherited disorders, especially Mendelian diseases. The criteria are grouped as "very strong", "strong", "supporting" etc. for both pathogenic and benign cases, whereas a variant is characterized by one of the five main classes according to the number of criteria met for each group. The main classes are "pathogenic", "likely pathogenic", "benign", "likely benign" and "variant of uncertain significance". For example, a variant is characterized as pathogenic if it meets more than two criteria which are in pathogenic strong (PS) group. If the variant cannot be classified in one of the main classes or classified conflictingly as both being pathogenic and benign, it is characterized as VUS.

[0010] VUSs are variants whose impact on individual have not been determined and therefore cannot be classified as pathogenic or benign. The largest number of variants in the ClinVar database belongs to the variant of uncertain significance (VUS) class. There are a total of 1,539,888 variants in the ClinVar Miner database, 674,594 of which being classified as VUS.

[0011] Although WES and WGS are used for disease diagnosis in routine clinical practice, many variants found are VUSs. These variants are mostly not reported by the clinicians due to issues regarding prioritization and interpretation. Thus, many variants with the potential to cause disease are left unaddressed. As a result, the diagnostic rate of WGS and WES technologies in rare diseases, remains at 50%. One of the most important reasons why the diagnosis rate in hereditary diseases is low, especially in rare diseases, is the failure to characterize VUSs.

[0012] International Patent Application WO 2023154778 A2 concerns training and utilizing a neural network obtaining sequencing variant sample data and training a neural network using the sequencing variant sample data. The example method further includes training the neural network by altering values at one or more loci for one or more samples in the sequencing variant sample data to produce altered sequencing variant sample data and using the altered sequencing variant sample data as input during training of the neural network such that the neural network is trained to predict values from unaltered sequencing variant sample data and the one or more reference genome maps. Another example method includes obtaining subject sequencing variant sample data for a subject and generating predicted values of a genome using a neural network. The method further includes determining one or more predicted conditions for the subject.

[0013] Several machine learning (ML) based tools are developed for the variant classification. LEAP variant classifier is specific to missense variants. The training set consists of variants in genes associated with hereditary cancer and cardiovascular diseases. Publicly available information such as functional impact predictors and splicing impact predictors are used to train Logistic Regression and Random Forest Classifier models. Logistic regression model as a result of cross validation test achieved an area under the receiver operating characteristic curve (AUROC) of 97.8% (cancer) and 98.8% (cardiovascular), while the random forest model achieved 98.3% (cancer) and 98.6% (cardiovascular). Hold-out data was used to test the Leap models, as well. Gene-holdout predictions achieved 96.8% AUROC. In another ML based variant classifier called ClinPred, non-synonymous variants in the ClinVar database are used to train random forest and gradient boosted decision tree (XGBoost) models. Functional scores such as SIFT, PolyPhen, frequency and conservation scores are used as features during model development. The AUC score of the ClinPred model on the ClinVartest set is 0.98. It is important to note that, these two variant classifiers are DNA-Seq based models. Therefore, in the art, there exists a lack of robust machine learning models for classifying VUS and those take advantage of RNA-Seq data. Objects of the Present Invention

[0014] Primary object of the present invention is to provide a robust machine learning method for predicting pathogenicity of VUS.

[0015] Another object of the present invention is to provide a machine learning model based on random forest that utilizes splicing and expression features, as well as features such as pathogenicity prediction scores and population frequency.

[0016] Another object of the present invention is to provide a machine learning model based on random forest that utilizes transcriptome sequencing, otherwise known as RNA-Seq data.

[0017] Another object of the present invention is to provide a machine learning model based on random forest that utilizes gene expression level and splicing information, in addition to various other features, to increase success rate of hereditary disease diagnosis.

[0018] Brief Description of the Present Invention

[0019] Disclosed invention mainly aims to achieve classification of VUSs and their pathogenic potential. In clinical practice, from most of the diagnostically utilized DNA sequencing methods, whole exome sequencing (WES) technique is generally preferred. Based on difference to a reference genome obtained by sequencing of genes from thousands of healthy individuals, base differences observed in the genome of an individual is determined to be a variant. Many such variants are obtained based on DNA-Seq methods clinically used. However, the pathogenicity of these variants is restricted to the ones that are found in public databases such as ClinVar, whereas many variants cannot be classified as either pathogenic or healthy, which are consequently dubbed as Variants of Uncertain Significance (VUS). Since VUSs do not immediately offer any insight as to the pathogenicity, they are mostly unreported by clinicians. Thus, many potentially pathogenic variants are left unaddressed, which is one of the reasons for diagnosis rate in rare diseases remaining low.

[0020] According to the present disclosure, features that can be obtained from RNA- Seq data such as gene expression level and splicing, in addition to DNA-Seq data, can strongly predict pathogenicity when used for a ML-based model when targeting VUS. Thus, the potential of many undiagnosed diseases to be correctly diagnosed based on such a ML-based model is significantly increased. Strength of the model proposed by the present disclosure has a threefold advantage based on the processing of data obtained from publicly available data found in the art, the type of machine learning approach used namely random forest, and finally the features that are used for training said random forest for classifying VUSs.

[0021] RNA-Seq data is very rich in terms of information it contains, as it can be processed with many different analysis types. These types of analysis mainly include (i) measurement of gene expression levels, (ii) identification of genomic variants, (iii) identification of splicing events, and (iv) detection of gene fusions. Fusion gene detection from RNA-Seq data mainly performed for non-hereditary (somatic) diseases. Other analyses from RNA-Seq data are performed for both hereditary (germline) and non-hereditary (somatic) diseases. Recent studies have revealed the power of RNA-Seq on variant interpretation and prioritization. In these studies, the possibility of VUSs to be pathogenic was evaluated according to the splicing and expression information obtained from RNA-Seq data. If a variant disrupts splicing or causes excessive increase or decrease in gene expression level, variant is more likely to be pathogenic. It was observed that the diagnosis rate by using only WGS and WES data is, increased by 10-35% with the contribution of RNA-Seq data.

[0022] Proposed invention makes use of raw data obtainable from publicly available databases such as GEO (Gene Expression Omnibus) and ENA (European Nucleotide Archive), which comprise Mendelian and hereditary disease RNA- Seq data. Said RNA-Seq data are then subjected to bioinformatics to create training data, said training data being obtained based on two distinct steps of creating a feature table, and labeling of variants based on information obtained from ClinVar and ACMG. Feature table is created with compiling the results of for different types of analyses, namely variant, annotation, RNA-Seq expression and RNA-Seq splicing.

[0023] With variant analysis, all data pertaining to the differences to the reference genome are obtained. Subsequently, said variants in the data belonging to the outcome from the variant analysis are annotated, wherein information pertaining to the incidence, effect on protein function etc. are generated. As a result of this annotation, variants known to be benign, and variants known to be pathogenic are filtered and a feature table is compiled. Following this, a set of RNA-Seq expression analysis and RNA-Seq splicing analysis are executed whereby expression levels of variants based on the genes they belong to are calculated and their rate of expression compared to healthy individuals are obtained, and added to the compiled feature table as another set of features.

[0024] Said annotated feature table is finally used for training a random forest classifier whereby pathogenetic potential of variants of uncertain significance are rendered classifiable with a mean accuracy of 97%, which marks a significant advantage versus the known state of the art when VUSs are concerned for their pathogenicity potential.

[0025] Brief Description of the Figures of the Present Invention

[0026] Accompanying drawings are given solely for the purpose of exemplifying a random forest based VUS classification model and a system using said model, whose advantages over prior art were outlined above and will be explained in brief hereinafter.

[0027] The drawings are not meant to delimit the scope of protection as identified in the claims nor should they be referred to alone in an effort to interpret the scope identified in said claims without recourse to the technical disclosure in the description of the present invention.

[0028] Figure 1 demonstrates an exemplary feature table according to an embodiment of the present invention.

[0029] Figure 2 demonstrates a classification flow diagram according to an embodiment of the present invention.

[0030] Detailed Description of the Present Invention

[0031] According to the present disclosure, a novel variant classifier model based on random forest machine learning using RNA-Seq data is proposed. The highlight of the present invention is the contribution of RNA-Seq VUS characterization in machine learning models, which is brought about by using a wide and varied set of features used for training a random forest. Utilization of expression and splicing features data distinguishes the model from its counterparts. The model developed using public RNA-Seq data displays a performance of 98% accuracy and %88 accuracy for predicting pathogenicity of VUSs (specificity) on a holdout test data.

[0032] Present invention, according to multiple embodiments, utilizes transcriptome sequencing. Ribonucleic acid (RNA) is a single-chain molecule found in the nucleus, cytoplasm, ribosome, and mitochondria in human cells. On the other hand, deoxyribonucleic acid (DNA) that is found in nucleus and stores genetic information of humans is made up of two strands, in which the base pairs complement each other as base Adenine (A) pairs with base Thymine (T) and base Cytosine (C) pairs with base Guanine (G). RNA has a Uracile (U) base instead of Thymine and a five-carbon ribose sugar instead of deoxyribose unlike DNA. RNAs are divided into two main classes, coding, and non-coding, according to their potential to be translated into proteins. Messenger RNAs (mRNAs) are the coding RNAs. Non-coding RNAs are functional but cannot translated into protein. They are either structural (carrier RNA, ribosomal RNA) or regulatory RNAs. Regulatory RNAs are further divided into two types as nonprotein coding RNAs shorter than 200 nucleotides (microRNA) and longer than 200 nucleotides (IncRNA).

[0033] RNA-Seq data is very rich in terms of information it contains, as it can be processed with many different analysis types. These types of analysis mainly include (i) measurement of gene expression levels, (ii) identification of genomic variants, (iii) identification of splicing events, and (iv) detection of gene fusions. Fusion gene detection from RNA-Seq data mainly performed for non-hereditary (somatic) diseases. Other analyses from RNA-Seq data are performed for both hereditary (germline) and non-hereditary (somatic) diseases. Recent studies have revealed the power of RNA-Seq on variant interpretation and prioritization. In these studies, the possibility of VUSs to be pathogenic was evaluated according to the splicing and expression information obtained from RNA-Seq data. If a variant disrupts splicing or causes excessive increase or decrease in gene expression level, variant is more likely to be pathogenic.

[0034] Transforming genetic information in DNA into functional products (proteins) via RNA is called central dogma, which consists of two main stages of transcription and translation, with the former being the process of creating mRNA from DNA and the latter being the process of synthesizing protein from mRNA. Three main types of RNA are involved in protein synthesis, namely messenger RNA (mRNA), transfer RNA (tRNA) and ribosomal RNA (rRNA).

[0035] Transcriptome is the whole RNA transcripts in a sample, including proteincoding and other non-protein-coding transcripts. Transcriptome sequencing or RNA sequencing (RNA-Seq hereinafter) is a technique that provides comprehensive analysis of the transcriptome in a sample. Standard transcriptome sequencing steps are isolation of all RNA strands in the cells, selection of RNA molecule to be sequenced, conversion from RNA to cDNA, adapter ligation, amplification, and sequencing. The steps from conversion to cDNA to amplification is called library preparation. RNA isolation is the first step in transcriptome sequencing. The quality of the isolated RNA should be checked before proceeding to the sequencing step. The often-preferred criterion for assessing RNA quality is the RNA integrity number (RIN), which produces values between 1 and 10. The closer this value is to 10, the higher the quality of the RNA. Isolates are fragmented before being converted to cDNA. These fragments are then converted to cDNA, which is more stable than RNA. Adapter sequences that specific to the kit and platform are ligated to the fragments. Adaptor sequences mainly used for binding of fragments to the flow cells. After the fragments are amplified by methods such as polymerase chain reaction (PCR), the targeted RNA is sequenced. Next Generation Sequencing (NGS) is a sequencing technology known as the second-generation sequencing, whereas the first-generation "chain termination method" (or Sanger method) for determining the nucleotide sequence of DNA or RNA. In principle, the concepts behind Sanger and nextgeneration sequencing (NGS) technologies are similar. In both NGS and Sanger sequencing (also known as dideoxy or capillary electrophoresis sequencing), DNA polymerase adds own fluorescent nucleotides one at a time to a growing DNA sequence. Each included nucleotide is identified by its fluorescent label. The critical difference between Sanger sequencing and NGS is the sequencing size. While the Sanger method sequences only a single piece of DNA at a time, NGS is largely parallel. So, sequencing millions of pieces simultaneously at each run is possible.

[0036] NGS methods are divided into deoxyribonucleic acid (DNA) and ribonucleic acid (RNA) sequencing according to the genetic material of the organism to be sequenced. DNA sequencing methods are divided into whole genome sequencing (WGS), whole exome sequencing (WES) and panel sequencing according to the targeted region. WGS is the method by which all the nucleotides that make up the genome are sequenced, which expresses the entire genetic information of an organism. WES is a method in which proteincoding regions corresponding to 2% of the genome are sequenced. Panel sequencing is the method that is mostly used to search for a specific mutation associated with the disease and to target certain genes. This method is cost effective compared to other DNA sequencing methods. If it is aimed to discover a new gene associated with a disease, other DNA sequencing methods with higher cost and lower sequencing depth should be preferred instead of panel sequencing. General process workflow for all NGS methods are as follows: Sample pre-processing, library preparation, sequencing and bioinformatics analysis. NGS methods commence with sample isolation. Longer reads are obtained with NGS than with Sanger sequencing, but NGS devices are not capable of sequencing the entire sequence. Therefore, isolated DNA or RNA molecules are fragmented. The library preparation step consists of amplification and adapter ligation. Adaptors are short sequences of 25-30 nucleotides and are designed to complement sequences embedded on the flow cell in sequencing devices. In the adapter ligation step, the fragments can be attached to the flow cell by means of adapters that ligate the fragments. The purpose of the amplification process is to increase the number of fragments, this process is performed by PCR. The library preparation step is followed by sequencing. The raw data used in the analysis and obtained by the sequencing is FASTQ file. The final step in the NGS workflow is bioinformatic analysis where the sequence is interpreted. For example, if WES is performed, the aim of the bioinformatic analysis would be discovery of disease-related variants.

[0037] Bioinformatic analysis is the process of extracting useful information from biological data. The analysis workflow or pipeline refers to all the steps performed in a certain order to process said biological data. First few basic steps of the bioinformatic analysis pipeline are common for both DNA and RNA- based processing, which are data quality control, trimming, and alignment to the reference genome respectively. However, RNA-based data provide more information over DNA data via downstream analysis. While variant detection is the main bioinformatic analysis type with genomic data, various analyses such as measurement of expression levels, detection of variants, and splicing are performed using RNA-Seq data.

[0038] The bioinformatics pipeline of NGS data begins with quality control of reads. The quality control (QC) process typically involves assessing the actual quality of raw reads using metrics generated by the sequencing platform (e.g. quality scores) or calculated directly from raw reads. The purpose of the quality control step is to detect many sequencing artifacts such as base calling errors and primer / adaptor contamination. The quality of data is crucially important for various downstream analyses, such as variant detection or gene expression studies. For example, if adapter sequences are not removed before variant analysis, they will be mapped to regions that are like these sequences in the reference genome and differences will be reported as fictional variants.

[0039] There exist many tools in the art for assessing the quality of NGS reads such as NGS QC Toolkit, and NGSQC. The most popular one is FastQC that published by Babraham Institute by Trivedi et al. FastQC tool reports basic statistics and metrics such as base quality, overrepresented sequences, adapters, GC content. Duplication, k-mer or GC content levels in the FastQC output depend on experiment and organism.

[0040] Poor quality reads need to be removed besides adapters and barcodes. Common tools for this trimming operation are known in the art by names of Trimmomatic and Cutadapt. When a trimming produces a base of less than 20, those bases are excluded from the reads.

[0041] Clean reads that are obtained after trimming are aligned to the reference genome. The reference genome is the representative, ideal DNA sequence of all individuals of the species. Between the reads and the reference genome, alignment or mapping is the process of obtaining a complete genome of an organism by putting together the matching parts of the reference genome of the clean reads. The actual locations of the reads from the sequencing machine are tried to be located in the reference genome.

[0042] There also exist many aligner tools in the art for NGS data analysis. Some of the criteria that are to be considered when selecting the appropriate aligner that affect the performance of the aligner are as follows, i) length of reads obtained after sequencing, ii) sequencing quality, iii) errors due to sequencing, iv) memory requirement, v) sequencing depth. In the light of these criteria, the appropriate algorithm for analysis and the tool that adopts this algorithm can be selected as appreciated by the person skilled in the art. The alignment algorithms can be divided into two main groups, called as global alignment and local alignment. The purpose of the global alignment algorithm is to align over the entire sequence, and the purpose of the local alignment is to find the regions with most similarity and align to those regions. The most well-known examples of these algorithms are Needleman-Wunsch and Smith-Waterman respectively. Many aligners including Bowtie and BWA are also known in the art that are developed based on the theory of Burrows-Wheeler Transform (BWT). BWA and Bowtie tools are designed for alignment of DNA sequences. Splice aware aligner is needed for alignment of transcripts since tools like BWA do not handle intron-sized gaps. Two of the well-known and mostly used transcriptomic aligners are HISAT and STAR. According to a preferred embodiment, STAR aligner can be used for RNA-Seq analysis.

[0043] The reference genome must be indexed before the alignment with the STAR tool, as with other alignment tools. STAR alignment algorithm consists of seed searching and clustering / stitching / scoring steps. In the seed searching step, STAR searches for the longest sequence that exactly matches one or more locations. The longest matching sequences are called Maximal Mappable Prefixes. The part of the read that aligns with the reference genome is called the seed. The seeds are clustered and scored based on their proximity to the seed in the original read, and scored based on mismatches, indels, gaps after seed searching step. Seeds showing the best alignment for reading are stitched. At the end of the alignment step, the sequences mapped to the reference genome are stored in the Sequence Alignment Map (SAM), a textbased file. The file in SAM format consists of header and alignment parts. The parts that start with in the SAM file represent the headers. Each line in the alignment section represents the alignment of a segment and consists of eleven or more fields.

[0044] To increase the performance in the analysis and to reduce memory usage, the analysis continues with the binary alignment map (BAM) extension, which is the compressed binary form of the SAM file.

[0045] Reads can be aligned to the reference genome or to the reference transcriptome during the alignment of RNA data. While it is more advantageous to align to the transcriptome to speed up the analysis, the disadvantage is that de novo transcripts cannot be detected.

[0046] Variants are the difference from the reference genome in at least one base. It can have different names according to the rate of occurrence in the population. Variants with a population frequency below 1% are called mutations, and variants with a population frequency above 1% are called polymorphisms. Variants are divided into two main classes, single nucleotide variation (SNV) and insertion / deletion (INDEL) or structural variants such as copy number variations and fusion genes. A single base change is called a single nucleotide variant, while the addition of more than one base is called insertion, and the deletion of more than one base is called deletion. In the literature, insertion and deletion events are mostly expressed as INDEL.

[0047] The process of detecting variants from NGS data is dubbed variant calling. Different algorithms have been developed for somatic and germline variant calling. Germline variants occur in the germ cells and that are passed onto offspring. The germline variant affects all body cells, so it can be detected from any tissue. Somatic variants occur in somatic cells and are non-inherited. They are detected from the target tissue. The germline variants are detected easier than somatic variants. One of the major problems of somatic variant detection is caused by tumor heterogeneity, which is a result of low variant frequency. There are many tools developed to detect germline and somatic variants known in the art, including Genome Analysis Toolkit (GATK), VarScan, and SAMtools. The output of the tools is called variant call format (VCF) file, which consist of a header section and variant call records. The definitions of the columns are in the header section. Filter column is empty before the variant filtration step. The filtration step is carried out to filter out false positives. The filtering criteria depends on the data and sequencing technology, and filtering is done according to the quality metrics such as fisher strand, sequencing depth in the INFO section.

[0048] The variant analysis is generally performed for DNA-Seq data previously stated as WGS and WES. Whole exome sequencing (WES) variant analysis is mostly used in clinical implementation. WES increased the diagnostic efficiency of Mendelian diseases and enabled the discovery of new genes. Nevertheless, the diagnostic yield of WES analysis for rare diseases remains at 50%, mainly due to prioritization and interpretation problems. Although whole genome sequencing (WGS) allows detection of all variants, the diagnosis rate in clinical practice is close to WES. Furthermore, prioritizing and interpreting the total 3 million of variants obtained from human WGS data is a challenge.

[0049] Identification of genomic variants from RNA-Seq data is intrinsically complex. However, it has been shown that variant analysis can be performed from RNA- Seq data due to its cost-effectiveness compared to WGS and the detection of variants such as splice or intronic variants which cannot be detected by WES. Variant detection from RNA-Seq data according to the present invention is performed according the following workflow: After the alignment step, duplicates are marked. After the this step an additional step, compared to WES variant analysis, is applied which is performed with the SplitNCigarReads tool. SplitNCigarReads tool is used when RNA reads are aligned against the reference genome. This tool splits reads that contain Ns in their cigar string (e.g., spanning splicing events in RNA-Seq data). The aim of this process is hard clipping any sequences overhanging to the intronic regions.

[0050] The workflow continues with variant calling and variant filtration after SplitNCigarReads. Variant filtration step is applied to filter false positive variants. Variant calling tools report all differences as variants. However, variants that pass certain criteria such as alignment quality, and strand bias are considered true positive variants. If the variant is found on only one thread in paired reads (that is, only forward or reverse), this variant is considered false positive. The Fisher strand score represents this criterion, and variants with FS > 30 for RNA-Seq data were considered false positives. For the quality by depth (QD) criterion, which expresses variant reliability, variants with QD < 2 are considered false positives.

[0051] Annotation is a process where information stored in databases is assigned to variants, subsequent to which the variants are interpreted and prioritized with information such as the population frequency of the variant, its effect on protein function, etc.

[0052] There exist many annotation tools known in the art including Annovar, GATK Variant Annotator, SnpEFF etc. Annotation of the variants detected according to the present invention can be, in at least one embodiment, carried out with the Variant Effect Predictor (VEP v.104) tool ENSEMBL. VEP is an annotation tool that provides information at the transcript level and includes various plugins such as non-mediated decay and loFtool. This tool also allows for adding custom databases that can be used for annotation. After annotation, details such as genotype and gene-level consequences obtained from VEP are sourced from RefSeq and Ensembl.

[0053] According to at least one embodiment of the present invention, Transcript support level (TSL) is used: GENCODE TSL offers a standardized approach to evaluate the degree of support indicating the actual expression of a GENCODE transcript annotation in humans.

[0054] According to at least one embodiment of the present invention, RNA Editing is used: RNA editing is the posttransciptional modification, RNA editing can be observed in two general types, substitution editing and insertion / deletion editing respectively.

[0055] According to at least one embodiment of the present invention, SIFT (Sorting Intolerant from Tolerant) tool is used: It is the score showing the effect of amino acid substitution on protein function. SIFT algorithm uses sequence homology considering the position at which the change occurred and the type of amino acid change.

[0056] According to at least one embodiment of the present invention, LRT (Likelihood Ratio Test) is used: LRT shows the effect of missense changes on the function of human protein. LRT identifies a subset of deleterious mutations using comparative genomics dataset of 32 vertebrate species.

[0057] According to at least one embodiment of the present invention, MutationTaster tool is used: MutationTaster scores show the disease-causing potential of DNA sequence alterations.

[0058] According to at least one embodiment of the present invention, MetaLR and MetaSVM tools are used: These are used to evaluate other prediction scores using two ensemble-based approaches, support vector machine (SVM) and logistic regression (LR).

[0059] According to at least one embodiment of the present invention, PROVEAN (Protein Variation Effect Analyzer) is used: PROVEAN predicts the functional effects of protein sequence variations including single or multiple amino acid substitutions, and in-frame insertions and deletions. PROVEAN consists of two main steps: The first step is the collection of homologous and distantly related sequences (with BLASTP) and the second step is clustering.

[0060] According to at least one embodiment of the present invention, FATHMM (Functional Analysis Through Hidden Markov Models) is used: It is used to predict the functional effects of protein missense variants.

[0061] According to at least one embodiment of the present invention, FATHMM- coding is used: It is based on a machine learning approach (FATHMM-MKL (Multiple Kernel Learning)) that integrates functional annotations from ENCODE with nucleotide-based sequence conservation measures and predicts the functional consequences of coding variants. Pathogenic dataset was constructed using heritable germ-line mutations from the Human Gene Mutation Database.

[0062] According to at least one embodiment of the present invention, FATHMM-XF (FATHMM-Extended Feature) is used. FATHMM-XF is a model based on machine learning. The training set of the model consists of 156,775 coding examples and 25,720 non-coding examples. SNVs are characterized using features from datasets.

[0063] According to at least one embodiment of the present invention, DEOGEN2 is used: DEOGEN2 integrates information such as the molecular effect of the variant, the relevance of the mutated gene, its known interactions, and the pathways in which it is involved. The Random Forest Classifier from sckit-learn is used in DEOGEN2 to predict the effect of single amino acid substitution on protein.

[0064] According to at least one embodiment of the present invention, BayesDel is used. BayesDel is meta score included PolyPhen, SIFT, FATHMM, LRT, Mutation Taster, Mutation Assessor, PhyloP, GERP++, and SiPhy. This score can be used coding / non coding and SNP / INDEL variants. There are two different type scores as included population frequency information and not population frequency information. As a more prominent example, BayesDel score, which does not include frequency information, is used to prevent double counting.

[0065] According to at least one embodiment of the present invention, primateAI is used: Primate Al is a prediction score based on deep learning. The model of the scores was trained with ~380,000 common missense variants from humans and six non-human primate species.

[0066] According to at least one embodiment of the present invention, PhyloP is used: PhloP means phylogenetic p value. A positive value indicates evolutionary conservation, and a negative value indicates fast evolving.

[0067] Various plugins can also be used according to the present invention when running the VEP tool, and can be listed as follows: EVE (Evolutionary model of Variant Effect) score: EVE is a pathogenicity prediction score. Amino acid sequences of over 140K species were used while model training.

[0068] NMD (Nonsense Mediated mRNA Decay): Nonsense-mediated mRNA decay (NMD) is a biological mechanism that eliminates mRNAs containing premature stop codons and prevents production of truncated protein. The plugin indicates whether the variant can escape this mechanism or not. loFtool: Loss of function variants disrupt the function of protein coding genes.

[0069] LoFtool ranks genes by calculating the ratio of LOF variants to synonymous variants. The gene's low LoFtool percentage value indicates that it is the gene most intolerant to functional variation.

[0070] LOEUF (Loss-of-function Observed / Expected Upper-bound Fraction): LOEUF is another loss of function score, wherein low scores are associated with disease genes. gnomAD (Genome Aggregation Database): GnomAD is a database containing the population frequency of the variant. The variants in gnomAD_exome (v.2.1.1) are derived from 125,748 exome sequencing data. The variants in gnomAD_genome (v.3.1.2) are derived from 76,156 genome sequencing data.

[0071] ClinVar: ClinVar is a database that holds the clinical classification information for the variants.

[0072] RNA-Seq data is mostly used to measure gene or transcript level in sample and to detect differentially expressed genes between two or more conditions. The expressed gene refers to the active gene. In other words, a functional product is created from DNA. RNA-Seq expression analysis commence with quantification. Two different approaches can be implemented to quantify at the gene level. In the first approach, the number of reads belonging to a gene is determined by counting the reads that overlap with the gene location. In the second approach, the number of transcripts corresponding to the gene is calculated and then the number of reads per gene is determined.

[0073] There are two major quantification methods. Alignment-free quantification methods quantify transcripts with the pseudo alignment in k-mer space. Salmon and Kallisto are the popular alignment free and transcript level quantification tools known in the art. Alignment based quantification is a method to quantify only reads aligned to reference genome or transcript, wherein reads that are aligned to two or more gene locations are ignored. Quantification tools such as Featurecounts or htseq-count are configured to ignore reads overlapping multiple genes in their default operating mode and do not count them. These two tools are the most common used tools in the quantification step. These tools are mostly preferred because of they are compatible with the tools used in differential expression analysis without the need for any processing. In exemplary read count matrices, there are raw read counts of the genes obtained from the quantification step. In a read count matrix, each row corresponds to a gene while each column corresponds to a sample.

[0074] Gene read counts in the read count matrix are not normalized. Raw read counts cannot be used to compare expression levels between samples. These counts need to be normalized regarding factors such as gene length. Genes with longer lengths are expected to have more read counts. In addition, reading depth is another factor that affects read counts. As the reading depth increases, the number of aligned reads also tend to increase. Therefore, comparison of raw counts of samples with different read depths will not be reliable. Different methods have been developed for the normalization of read counts. TPM (transcript per million), RPKM (reads per kilobase of transcript per million reads mapped) and FPKM (fragments per kilobase of transcript per million reads mapped) are the normalization methods based on gene length and depth. FPKM is very similar to RPKM. RPKM is used for single-end RNA-seq while FPKM is used for paired-end RNA-seq. The formulas of the normalization methods are as given below, per Formula 1 and Formula 2.

[0075] Number of reads mapped to gene x 103x 106

[0076] FPKM or RPKM= (1)

[0077] Total number of mapped reads x Gene length (bp)

[0078] Total reads mapped to gene x 103

[0079] TPM = E(A) Gene length (bp)

[0080] Unlike FPKM, TPM is first normalized to the gene length and then to the total number of mapped read(depth). As a result, the sum of all TPM values is the same in all samples. These three mentioned methods are gene length normalization techniques. When considering the differences in library sizes between samples, it is crucial to take into account the library size when comparing the samples. TMM (trimmed mean of m-values) is a library size normalization method known in the art. In this method, calculations are performed by using both case and control samples together.

[0081] Differentially expressed genes (DEG) can be detected using packages such as DESeq. There are many tools that make bioinformatic analysis practical. DeSeq2, EdgeR and limma tools are among the packages coded in R and widely used in DEG analysis as known in the art. These tools take the raw read counts as input and output the normalized expression values, statistical p values and fold change for each gene.

[0082] Splicing is a process in which introns are removed from the pre-mRNA and exons are ligated, resulting in the transcript being converted into mature mRNA. Alternative splicing is the process of ligating different exon combinations to transform the final mature mRNA. It is possible to produce more than one protein from one gene with alternative splicing processes. In constitutive splicing, introns are removed, and exons are ligated. In exon skipping, cutting occurs from the receiving region of the next exon instead of successive exons, and introns are removed along with the exon in between. The remaining exons are ligated. Alternative 5'splice site is upstream donor site and intron, with a part of exon are removed. The alternative 3' splice site is downstream of the acceptor site and intron, with a part of exon are removed. Intron retention is an alternative splicing event in which mature mRNA containing intron is occurred. Intron often contain premature termination codons, intron retaining isoforms are often degraded. If intron retaining isoforms are not degraded and converted to truncated protein isoforms. Splicing process can be repeated if intron retaining isoforms are retained in the nucleus or cytoplasm.

[0083] The removal of the intron is carried out by the ribonucleoprotein complex called the spliceosome. Splicing is regulated trans-acting proteins (repressors and activators) and cis-acting regulatory sites (enhancer and silencer). When the repressor proteins bind to the silencer site, the splicing process is inactivated. When the activator proteins bind to the enhancer site, the splicing process is activated. The boundaries between introns and exons are called splice sites. Upstream part of the intron is called the donor splice site and downstream part of the introns is called the acceptor splice site. The first two nucleotides GT (donor splice site) at the beginning of the intron and the last two nucleotides AG (acceptor splice site) represent the canonical splice site. 98.7% of splice site pairs are canonical. Splice site pairs, AT-AC and GC and AC, respectively, are called non-canonical. The canonical and non-canonical splice sites are known for their different snRNP combinations.

[0084] A variant that occurs at the boundary of an exon and an intron (splice site) is called splice site mutation. 15-30% of variants causing hereditary diseases are observed to be splice site mutations. Furthermore, about 62% of all pathogenic single nucleotide variants are thought to affect RNA splicing. It is also known that variants occurring in non-coding regions may cause aberrant splicing. Studies in the literature have also shown that RNA-Seq splicing analysis contributes to the interpretation and prioritization of variants and increases diagnosis rate. There is no best practice for splicing analysis, however, there are several specialized tools for predicting splice changes. Early published splicing tools such as MaxEntScan focus on the close neighbourhood of splicing junctions. In the last few years, the number of published artificial intelligence based splicing tools has increased with the use of machine learning for biological data. Recently, deep neural networks (DNNs) have also achieved good results in predicting genome-wide splice variants. SPANR / SPIDEX is the first deep learning (DL) based tool. SPANR / SPIDEX model trained with the experimentally observed exon skipping events. According to at least one aspect of the present disclosure, splicing analysis can be performed with the SpliceAI tool released by Illumina, which is another deep learning-based tool.

[0085] SpliceAI analyses each position in a pre-mRNA transcript and evaluates whether it is likely to be a splice donor, acceptor, or neither. The model is only trained on reference transcript sequences and splice junction annotations, and never saw variant data during training. Training set can be constructed with protein-coding transcripts from the GENCODE v24 annotation. The effect of the variant on splicing is evaluated according to the delta score produced by the tool. The tool is used to predict firstly exon-intron boundaries for both the reference pre-mRNA transcript sequence and the alternative transcript sequence containing the variant and then takes the difference between the scores.

[0086] The SpliceAI tool uses a .vcf file, a reference genome, and an annotation file (.txt) as input. Annotation file contains chromosome, strand, transcript start (TX_START) and end (TX_END) position and exon start end positions. The tool cannot calculate for genes that are not included in the annotation file in the tool. The file should be updated by adding the relevant information of the genes for which calculation is requested to this file. The tool generates a vcf file that includes the SpliceAI scores as output.

[0087] The output file contains eight values for each variant, of which four are delta scores and four are positions, the scores show the effect of the variant on splicing, while the positions indicate the position where the splicing changes relative to the variant position.

[0088] The American College of Medical Genetics and Genomics (ACMG) collaboration with the Association of Molecular Pathology (AMP) have published a guideline consisting of 28 criteria to be used for interpretation of germline variants in 2015 known in the art. The criteria are grouped according to biological impact, in silico predictions, presence in the control cohort, familial information, and inheritance mode. Variants are evaluated hierarchically and divided into 5 main classes as pathogenic, likely pathogenic, benign, likely benign and VUS according to the proposed criteria. Manual evaluation of these criteria on a one- by-one basis is challenging as it is time consuming and requires expertise. Several tools have been developed that automate the ACMG criteria, while neither of them being machine learning based, these tools are limited due to their gene specific nature. Machine learning based tools developed in line with the ACMG criteria and published in the literature to automate the variant classification with additional features, but since some of these criteria evaluate variants according to familial information, not all criteria are suitable for general use. Generalizable criteria such as variant population frequency from ACMG-AMP criteria is used as features in ML based tools. The success of ML-based methods depends on the selected training dataset and features. The features used in the models developed so far can be obtained from databases or can be calculated by DNA- Seq-based bioinformatic analysis methods. These methods are still insufficient for characterizing remaining variants in the deep intronic region. Therefore, there is a need for a tool that uses features obtained from RNA-Seq data.

[0089] In unsupervised machine learning algorithms, it is unclear which class the input data belongs to. These algorithms utilize a function to extract an unknown structure over unlabeled data. Clustering and association are examples of tasks that can be accomplished by unsupervised machine learning algorithms. Clustering is the process of assigning similar observations to the same clusters. An association algorithm analyses the co-occurrence of events and establishes relationships between data.

[0090] Semi-supervised learning is a method that falls between supervised and unsupervised learning. These types of algorithms are used for improving the accuracy of learning generally when there is small amount of labelled data for training and a large amount of unlabeled data. Reinforcement learning is a method that enables an agent to learn through trial and error in an interactive environment using feedback from their own actions and experiences. Reinforcement learning is mostly used to train robot motions.

[0091] In machine learning, the input is known as the instance (data) and the output is known as the target (label). The input used for training is a series of data called set of features. The features can either be continuous, binary, or categorical. Categorical features must be converted to numeric due to machine learning models require all input and output variables to be numeric. The terms, validation data and test data are frequently handled along with training data. The validation data is used to set the model parameters correctly during model training. Validation data is not used in unsupervised methods because it is not strictly correct. The test data is the input that the model does not encounter during the training. Model performance is estimated using the test data.

[0092] Data and model fit are extremely critical for the success of the model. Compared to data complexity, overly complex or simplistic models would reveal overfitting or underfitting problems, which is reflected in the model selection process which requires a trade-off between bias and variance. The bias the value that reflects the distance between the estimated data and the actual data because of the modelling. The variance is the variability of the model estimate for a given data point, or the value that indicates how the data is spread out. High bias and low variance are interpreted as underfitting. In other words, the model is simplistic when it is compared to the data complexity. In this case, the accuracy is low since the model cannot learn the data sufficiently. Low bias and high variance are interpreted as overfitting, meaning that the model is more complex than necessary. The model picks up important relationships between variables, as well as relationships that turn out to be just the result of a noise. Hence the model cannot perform well on the test data because it memorizes the training data.

[0093] Random Forest is a supervised machine learning model developed by Breiman, whose algorithm is built on decision tree models. The nodes with the first split and representing the entire sample is called the root node. If a sub-node is divided into further sub-nodes, it is called a decision node. Nodes that represent the class or label that are not divided into further sub-nodes are called leaves or terminal nodes. Each node in the tree represents some features edge represent possible answer. The purpose of dividing into sub nodes is to increase homogeneity. The purest node or most homogeneous node has samples from only one class. Node splitting is one of the factors affecting the performance of decision trees. There are different node splitting strategies depending on whether the target variable is categorical or continuous. If the target variable is categorical (such as classification problem), information gain or entropy is used, otherwise (such as regression problem) reduction in variance is used. The variance is calculated to decide the homogeneity of the node in reduction in variance. A completely homogeneous node has zero variance.

[0094] Random forest consists of many decision trees, wherein the predictions from individual trees are collected, and the final class is estimated by combining the individual predictions using a method such as the majority voting for classification and averaging for regression. In fact, each prediction is the output of a model. This process of combining different model results is called ensemble learning. The ensemble learning method used by the random forest algorithm for training is bagging or bootstrap aggregation. The concepts of using multiple trees in parallel and bootstrapping results in overcoming of overfitting issues by random forests, whereby random samples are selected from the dataset with replacement.

[0095] Model performance is calculated according to the differences between the correct prediction by the model and the actual class in supervised machine learning algorithms. The closer the model predicts the actual class, the higher the performance of the model. There are different performance metrics that are used for model evaluation purposes. To calculate these metrics, the numbers of true positive (TP), false positive (FP), true negative (TN), and false negative (FN) predictions are employed. In a confusion matrix produced for assessing the classes, the actual classes are assumed to be positive and negative. It is called TP, when the model predicts a member of positive class as positive, and FN when it predicts the member as negative. It is called TN when the model predicts a member of negative class as negative, and FP when it predicts the member as positive.

[0096] The performance of classification problems can be evaluated using accuracy, Fl score, sensitivity, or specificity. Accuracy represents the number of correctly classified samples over the total number of data samples. Accuracy considers the number of all correctly classified samples and do not concentrate on a specific class. Accuracy calculation is demonstrated according to the formula given below:

[0097] Accuracy = TP+TN / (TP+TN+FP+FN) (3)

[0098] Since there are different numbers of negative and positive samples in the imbalanced data, the Fl score is a better measure to evaluate model performance instead of accuracy. Fl score calculation is demonstrated in the formula below. If the positive class mattered the most as compared to negative, both precision and recall are crucial for the evaluation of model performance. Fl score is balancing precision and recall on the positive class while accuracy looks at correctly classified observations of both positive and negative set. Fl score = 2x (Precision xRecall) / (Precision + Recall) (4)

[0099] Sensitivity is the positive class detection performance of the model, that is, it shows how many of the members in the positive class are correctly classified. Specificity, the opposite of sensitivity, is the negative class detection performance of the model, that is, it shows how much of it is classified as negative within the negative class.

[0100] Random forest is a popular supervised machine learning algorithm that can be used for both regression and classification problems. It is an easy-to-use algorithm that gives fast results, in addition to giving good results even without an extensive parameter setting process. Therefore, the random forest classifier is used while developing the clinical characterization model of germline variants according to the present disclosure.

[0101] The feature table which is input to the random forest classifier model is, according to an exemplary embodiment, prepared as follows: First, public RNA- Seq data is collected to be used in bioinformatic analyses. All studies on Mendelian or hereditary cancer with published raw RNA-Seq data are considered, during data collection efforts. The publicly available datasets of the studies are, in an embodiment, obtained from GEO database. The datasets consist of patients (with Mendelian and hereditary cancers) and control individuals.

[0102] Some of the diseases that can be targeted with the aid of the present invention can be exemplified as follows: Myotonic dystrophy type 1 (DM1) is an inherited muscle disease that causes progressive muscle weakness and loss. When one copy of the gene is altered, the disease appears in the phenotype. This indicates that the inheritance pattern is autosomal dominant. Cystic fibrosis (CF) is an important genetic disease that affects various biological systems, especially the lungs and digestive system. Inheritance pattern of CF is autosomal recessive. Mutations occur in both copies of the gene in autosomal recessive diseases. CF affects cells that produce mucus, sweat, and digestive juices. Lynch Syndrome (LS) and Familial Adenomatous Polyposis (FAP) are autosomal dominant (dominant) diseases. The most common hereditary colon cancer is non-polyposis colorectal cancer. This disease is also known as Hereditary Non-polyposis Colorectal cancer (HNPCC) or Lynch Syndrome. Unlike familial adenomatous polyposis (FAP) syndrome, hundreds of polyps do not occur in these patients. Breast cancer is a type of cancer that is very common in women. 70% of breast cancer is sporadic, 10% is hereditary and 20% is familial. Hereditary breast cancer exhibits dominant autosomal inheritance.

[0103] Paired RNA-Seq data collected from samples with disease (case samples) in the datasets are used for variant calling, variant annotation, expression, and splicing analyses, respectively. The goal of variant calling analysis is to detect germline variants related to the hereditary diseases. For this purpose, HaplotypeCaller tool is used in the variant calling step and the variant list was obtained in vcf format.

[0104] According to one particular embodiment, the annotation tool used is VEP. VEP tool (v. 104) is used in the invention for the purpose of annotating variants, plugins of which are explained previously. ClinVar (v. 10.07.23) and ACMG 2015 criteria are used for clinical characterization of variants. Classification based on ACMG criteria is used together with the ClinVar database for characterizing variants. ACMG criteria that variants meet are also added to the feature table as a new feature. ClinVar database does not divide variants into 5 main classes as pathogenic, benign, likely pathogenic, likely benign and VUS, as is the case with ACMG criteria. Rather, ClinVar characterizes variants in different terms, including the main classes found in the link under Disagreement.

[0105] Among these terms is the 'Conflicting Interpretation of Pathogenicity'.

[0106] Conflicting Interpretation of Pathogenicity variant label represents the variants reported differently by different submitters of ClinVar database. For example, one such variant may be reported as possibly pathogenic by one submitter, another submitter may report the same variant as likely benign. Since there is no common decision on such variants, they are to be ignored under the scope of the present disclosure, and treated as VUS.

[0107] In expression analysis, raw counts at the gene level are initially obtained from RNA-Seq data using the FeatureCount tool. For data sets without control samples, the number of control samples is obtained from the GTEx database. Raw counts are normalized using the EdgeR R package, with TMM (trimmed mean of M values) employed as the normalization method. PCA analysis is performed for each data set to determine whether there are outlier samples before differential expression analysis (DEG). Outlier samples are then removed from the analysis. DEG is performed using the EdgeR R package. The decision on whether gene expression significantly changed is based on the p- value obtained from the analysis. Upregulated and downregulated genes are determined based on the fold change (FC). The p value and log FC values obtained as a result of the DEG analysis are also added to the feature table.

[0108] Splicing effects of the variants obtained in the variant calling step are determined with the precomputed SpliceAI scores available in VEP and released by illumina. SpliceAI tool produces four delta scores for acceptor and donor site as output. These scores are AG (acceptor gain), AL (acceptor loss), DG (donor gain) and DL (donor loss). Acceptor and donor site are the splice site region in the genome. Delta scores in the relevant regions are the difference of the maximum scores calculated according to the reference and alternative variant positions. As a result of the splicing analysis, it is determined whether the splicing process influenced by the variant or not according to values of the delta scores. A value of 0.2 is used as threshold. The variants with delta scores greater than 0.2 are consequently determined as splice-affecting variants. Such variants are labelled with 'yes' while variants not affecting splicing processes are labelled as 'no'. Two example variants with and without splicing effect are determined based on the 0.2 cut-off. The effect of the variant on the splicing process is also added to the feature table as a new feature.

[0109] Expression information is tissue-specific information, and the RNA-Seq data of the samples in the datasets used for the model are isolated from different tissues. As a solution to this problem, the distance of the genes to the housekeeping genes necessary for the maintenance of basic cellular functions is calculated. Housekeeping genes are founder genes expressed in all cells of an organism. Tissue specific housekeeping gene list was obtained from the database called Housekeeping and Reference Transcript (HRT) Atlas. This distance calculated was added to feature table.

[0110] In various embodiments, a feature table is prepared with the information obtained by bioinformatic analyses. The class labels are given prior to the training step due to the random forest model used according to this invention being a supervised machine learning algorithm. Variants are labelled as pathogenic and benign according to ClinVar and ACMG criteria. The pathogenic and benign variants are selected as follows. If a variant is pathogenic / likely pathogenic according to both sources (ClinVar and ACMG) or pathogenic / likely pathogenic according to one database and not benign / likely benign according to the other database, it is included in the training set. The opposite is done for benign ones. Selected variants are in canonical transcript or pathogenic variants in other transcripts. Accordingly, pathogenic and likely pathogenic variants are considered to be pathogenic and benign and likely benign variants are considered to be benign. Benign variants labelled with "0" in datasets and pathogenic variants labelled with "1". After identifying the pathogenic and benign variants according to the ClinVar and ACMG, filters for population frequency and consequence features are applied. For the first filter on population frequency, the threshold value of 0.05 was used due to the most applied threshold value being 0.05. As a result, a total of 7060 pathogenic and 420,373 benign variants are obtained.

[0111] SNP variants are used for training since the feature information used being mostly available for SNP variants. The number of benign variants is always greater than the number of pathogenic variants. To ensure the balance between pathogenic and benign variants in the training set and to prevent bias in terms of a feature, as many benign variants as the number of pathogenic variants are to be selected. Clustering is performed for benign variants selection. The star rating in ClinVar review status reflects the reliability of the variant's classification. Benign variants used in clustering analysis are determined to have at least one star according to ClinVar review status. Following the selection process, a total of 3940 variants (benign / pathogenic) are utilized for training.

[0112] As explained above, using the public RNA-Seq data, the features required for variant characterization are obtained via bioinformatic analyses and the datasets are formed. In the datasets, the feature table contains both numeric (e.g., SIFT score and population frequency) and categorical data (e.g., Gene Region and consequence). Also, some variants do not have value (missing value or NA) for several features. For example, there are variants without population frequency information. So, the features in the input data needs to be converted into a format suitable for the model. Missing data can be handled by deleting or imputing. Deleting missing values is expected to reduce the size of the dataset. However, since the purpose of the present disclosure is to accurately detect pathogenic variants, variants with missing information are not deleted so as not to reduce the number of pathogenic variants. Also, since features will be prioritized by the algorithm, the missing features for most variants are not deleted. Instead, missing values are sent as inputs to the present model.

[0113] There are various techniques in statistics such as mean imputation and median imputation to impute missing data. In mean or median imputation, missing values are replaced with the mean or median of the values for the corresponding feature. However, this kind of imputation is not reliable for biological features. Therefore, in the present invention data imputation is avoided as much as possible. Values that can be converted to categorical in the feature table are used as categorical, and "unKnown" is used for missing values. For example, If the SIFT score is less than 0.05, it is denoted by the letter "D" representing "Deleterious", otherwise it is denoted by a "T" representing "Tolerated". The thresholds from the literature are used to categorize other in silica prediction scores. In numeric features, "-1" is used for missing values in features that do not take negative values and cannot be converted to categorical. Missing values in features that contain negative values and cannot be converted to categorical are imputed according to the median. For features that are categorical and converted to categories are converted to numeric by the "one hot encoding" method, which is a technique that is widely utilized in the art. In one hot encoding, binary representation of categorical variables is used. The binary vector representing a categorical value is obtained by setting values in the vector to zero except the value at the integer index which is marked with one. In this representation, missing data is another categorical value, and it is converted to numeric value similarly, so the need for an additional data imputation is ameliorated.

[0114] After the model training set was prepared, the determination of the model parameters are carried out using Stratified K-fold (K:5) cross-validation. The n_estimator parameter was determined as 10. As a result of cross validation, the model mean accuracy is 97%.

Claims

CLAIMS1) A method for constructing a classifier to classify variants of uncertain significance (VUS), comprising steps of, obtaining raw data pertaining to information about clinical classes of variants, wherein data are collectible from databases comprising said info; said data comprising public Mendelian and hereditary disease RNA-Seq and / or DNA-Seq data and, training a random forest model using obtained variant clinical class data, wherein training said random forest model includes: generating training data wherein obtained data are filtered to highlight features and creating a feature table, said filtering being based on variant analysis based on difference to reference genome whereby pathogenic and benign variants according to ClinVar and ACMG data are included , variant annotation, based on at least multiple measures including distance to housekeeping genes, a loss-of-function score, nonsense mediated decay score, ACMG-AMP criteria,RNA-Seq expression analysis based on at least multiple measures including normalized gene expression levels, log-transformed expression values, and a p-value calculation based on normalization using trimmed mean of M values (TPM) and a coefficient change calculated based on change in gene expression of tissue-specific patient samples; and,RNA-Seq splicing analysis based on splicing scores predicting the variant's effect on splicing, using the created feature table as input during training of a random forest model such that the random forest model is trained to predict pathogenesis from unfiltered Mendelian and hereditary disease RNA-Seq data.2) A method for constructing a classifier to classify variants of uncertain significance (VUS) as set forth according to Claim 1, characterized in that variant annotation step is done based on a measure of post transcriptional RNA editing obtainable from a variant effect predictor tool, such as VEP.3) A method for constructing a classifier to classify variants of uncertain significance (VUS) as set forth according to Claim 1, characterized in that said raw data pertaining to information about clinical classes of variants are obtained from databases selectable from a group including ClinVar, ACMG, GEO and ENA.4) A system for classifying pathogenicity of variants of uncertain significance (VUS) comprising, an input means, whereby raw data comprising variants of uncertain significance are accepted, a processor configured to accept raw data via said an input means, wherein said processor being further configured to execute a random forest classification model as set forth according to Claims 1 to 4.

Images