Validation of predicted mRNA splicing mutations using high-throughput transcriptome data

Interpretation of variants present in complete genomes or exomes reveals numerous sequence changes, only a fraction of which are likely to be pathogenic. Mutations have been traditionally inferred from allele frequencies and inheritance patterns in such data. Variants predicted to alter mRNA splicing can be validated by manual inspection of transcriptome sequencing data, however this approach is intractable for large datasets. These abnormal mRNA splicing patterns are characterized by reads demonstrating either exon skipping, cryptic splice site use, and high levels of intron inclusion, or combinations of these properties. We present, Veridical, an method for in silico the automatic validation of DNA sequencing variants that alter mRNA splicing. Veridical performs statistically valid comparisons of the normalized read counts of abnormal RNA species in mutant versus non-mutant tissues. This leverages large numbers of control samples to corroborate the consequences of predicted splicing variants in complete genomes and exomes. 1 2 3


Introduction
DNA variant analysis of complete genome or exome data has typically relied on filtering of alleles according to population frequency and alterations in coding of amino acids. Numerous variants of unknown significance (VUS) in both coding and non-coding gene regions cannot be categorized with these approaches. To address these limitations, in silico methods that predict biological impact of individual sequence variants on protein coding and gene expression have been developed, which exhibit varying degrees of sensitivity and specificity 1 . These approaches have generally not been capable of objective, efficient variant analysis on a genome-scale.
Splicing variants, in particular, are known to be a significant cause of human disease 2-5 and indeed have even been hypothesized to be the most frequent cause of hereditary disease 6 . Computational identification of mRNA splicing mutations within DNA sequencing (DNA-Seq) data has been implemented to varying degrees of sensitivity, with most software only evaluating conservation solely at the intronic dinucleotides adjacent to the junction (i.e.7). Other approaches are capable of detecting significant mutations at other positions with constitutive, and in certain instances, cryptic, splice sites 5,8,9 which can result in aberrations in mRNA splicing. Presently, only information theory-based mRNA splicing mutation analysis has been implemented on a genome scale 10 . Splicing mutations can abrogate recognition of natural, constitutive splice sites (inactivating mutation), weaken their binding affinity (leaky mutation), or alter splicing regulatory protein binding sites that participate in exon definition. The abnormal molecular phenotypes of these mutations comprise: (a) complete exon skipping, (b) reduced efficiency of splicing, (c) failure to remove introns (also termed intron retention or intron inclusion), or (d) cryptic splice site activation, which may define abnormal exon boundaries in transcripts using nonconstitutive, proximate sequences, extending or truncating the exon. Some mutations may result in combinations of these molecular phenotypes. Nevertheless, novel or strengthened cryptic sites can be activated independently of any direct effect on the corresponding natural splice site. The prevalence of these splicing events has been determined by ourselves and others 5,11-13 . The diversity of possible molecular phenotypes makes such aberrant splicing challenging to corroborate at the scale required for complete genome (or exome) analyses. This has motivated the development of statistically robust algorithms and software to comprehensively validate the predicted outcomes of splicing mutation analysis.
Putative splicing variants require empirical confirmation based on expression studies from appropriate tissues carrying the mutation, compared with control samples lacking the mutation. In mutations identified from complete genome or exome sequences, corresponding transcriptome analysis based on RNA sequencing (RNA-Seq) is performed to corroborate variants predicted to alter splicing. Manually inspecting a large set of splicing variants of interest with reference to the experimental samples' RNA-Seq data in a program like the Integrative Genomics Viewer (IGV) 14 , or simply performing database searches to find existing evidence would be time-consuming for large-scale analyses. Checking control samples would be required to ensure that the variant is not a result of alternative splicing, but is actually causally linked to the variant of interest. Manual inspection of the number of control samples required for statistical power to verify that each displays normal splicing would be laborious and does not easily lend itself to statistical analyses. This may lead to either missing contradictory evidence or to discarding a variant due to the perceived observation of statistically insignificant altered splicing within control samples. In addition, a list of putative splicing variants returned by variant prediction software can often be extremely large. The validation of such a significant quantity of variants may not be feasible, for example, in certain types of cancer, in instances where the genomic mutational load is high and only manual annotation is performed. We have therefore developed Veridical, a software program that automatically searches all given experimental and control RNA-Seq data to validate DNA-derived splicing variants. When adequate expression data are available at the locus carrying the mutation, this approach reveals a comprehensive set of genes exhibiting mRNA splicing defects in complete genomes and exomes. Veridical and its associated software programs are available at: www.veridical.org.

Methods
The program Veridical was developed to allow high-throughput validation of predicted splicing mutations using RNA sequencing data. Veridical requires at least three files to operate: a DNA variant file containing putative mRNA splicing mutations, a file listing of corresponding transcriptome (RNA-Seq) BAM files, and a file annotating exome structure. A separate file listing RNA-Seq BAM files for control samples (i.e. normal tissue) can also be provided. Here, we predict mutations in a set of breast tumours which are validated with RNA-Seq data from the same individuals, with RNA-Seq data from normal breast tissues as controls. However, in principle, potential splicing mutations for any disease state with available RNASeq data can be investigated. In each tumour, every variant is analyzed by checking the informative sequencing reads from the corresponding RNA-Seq experiment for non-constitutive splice isoforms, and comparing these results with the same type of data from all other tumour and normal samples that do not carry the variant in their exomes.
Veridical concomitantly evaluates control samples, providing for an unbiased assessment of splicing variants of potentially diverse phenotypic consequences. Note that control samples include all non-variant containing files, as well any normal samples provided. Maximizing the set of control samples, while computationally more expensive, increases the statistical robustness of the results obtained.
For each variant, Veridical directly analyzes sequence reads aligned to the exons and introns that are predicted to be affected by the genomic variant. We elected to avoid indirect measures of exon skipping, such as loss of heterozygosity in the transcript, because of the possibility of confusion with other molecular etiologies (i.e. deletion or gene conversion), unrelated to the splicing mutations. The nearest natural site is found using the exome annotation file provided, based upon the directionality of the variant, as defined within Table 1. The genomic coordinates of the neighboring exon boundaries are then found and the program proceeds, iterating over all known transcript variants for the given gene. A diagram of this procedure is provided in Figure 1. The variant location, C, is specifically referring to the, variant-induced, location of the predicted mRNA splice site, which is often proximate to, but distinct from the coordinate of the actual genomic mutation itself.
The program uses the BamTools API 15 to iterate over all of the reads within a given genomic region across experimental and control samples. Individual reads are then assessed for their corroborating value towards the analysis of the variant being processed, as outlined in the flowchart in Figure 2. Validating reads are based on whether they alter either the location of the splice junction (i.e. junction-spanning) or the abundance of the transcript, particularly in intronic regions (i.e. read-abundance). Junction-spanning reads contain DNA sequences from two adjacent exons or are reads that extend into the intron. These reads directly show whether the intronic sequence is removed or retained by the spliceosome, respectively. Read-abundance validated reads are based upon sequences predicted to be found in the mutated transcript in comparison with sequences that are expected to be excised from the mature transcript in the absence of a mutation. Both types of reads can be used to validate cryptic splicing, exon skipping, or intron inclusion. A read is only said to corroborate cryptic splicing if and only if the variant under consideration is expected to activate cryptic splicing. Junction-spanning, cryptic splicing reads are those in which a read is exactly split from the cryptic splice site to the adjacent exon junction. For read-abundance cryptic splicing, we define the concept of a read fraction, which is the ratio of the number of reads corroborating the cryptically spliced isoform and the number of reads that do not support the use of the cryptic splice site (i.e. non-cryptic corroborating) in the same genomic region of a sample. Cryptic corroborating reads are those which occur within the expected region where cryptic splicing occurs (i.e. spliced-in regions). This region is bounded by the variant splice site location and the adjacent (direction dependent) splice junction. Non-cryptic corroborating reads, which we also term "anti-cryptic" reads, are those that do not lie within this region, but would still be retained within the portion that would be excised, had cryptic splicing occurred. To identify instances of exon skipping, Veridical only employs junction-spanning reads. A read is considered to corroborate exon skipping if the connecting read segments are split such that it connects two exon boundaries, skipping an exon in between. A read is considered to corroborate intron inclusion when the read is continuous and either overlaps with the intron-exon boundary (and is then said to be junction-spanning) or if the read is within an intron (and is then said to be based upon read-abundance).
We proceed to formalize the above descriptions as follows. A given read is denoted by r, with start and end coordinates (r s , r e ), if the read is continuous, as diagrammed within Figure 1(c), or otherwise, with start and end coordinate pairs, (r s 1 , r e 1 ) and (r s 2 , r e 2 ), as diagrammed within Figure 1(d). Let l be the length of the read. The set ζ denotes the totality of validating reads. The criterion for r ∈ ζ is detailed below. It is important to note that validating reads are necessary but not sufficient to validate a variant. Sufficiency is achieved only if the number of validating reads is statistically significant relative to those present in control samples. ζ itself is partitioned into three sets: ζ c , ζ e , and ζ i for evidence of cryptic splicing, exon skipping, and intron inclusion, respectively. We allow partitions to be empty. Let S denote the relevant splice junctions as defined by Figure 1 and Table 1. Without loss of generality, we consider only the red (i.e. direction is right) set of labels within Figure 1(b). Then the (splice consequence) partitions of ζ are given by:  Figure 1(b).

A B Strand Direction
Exonic    Table 1 for definitions concerning direction and Figure 1 for variable definitions.
We separately partition ζ by its evidence type, the set of junctionspanning reads, δ and read-abundance reads, α, as follows: Once all validating reads are tallied for both the experimental and control samples, a p-value is computed. This is determined by computing a z-score upon Yeo-Johnson (YJ) 16 transformed data. This transformation, shown in Equation 1, ensures that the data is sufficiently normally distributed to be amenable to parametric testing.
The transform is similar to the Box-Cox power transformation, but obviates the requirement of inputting strictly positive values and has more desirable statistical properties. Furthermore, this transformation allowed us to avoid the use of non-parametric testing, which has its own pitfalls regarding assumptions of the underlying data distribution 17 . We selected λ = 1 -2, because Veridical's untransformed output is skewed left, due to their being, in general, less validating reads in control samples and the fact that there are, by design, vastly more control samples than experimental samples. We found that this value for λ generally made the distribution much more normal.
A comparison of the distributions of untransformed and transformed data is provided in Figure S1. We were not concerned about small departures from normality as a z-test with a large number of samples is robust to such deviations 18 . It is important to realize, therefore, that the p-values given by Veridical are much more robust when the program is provided with a large number of samples.
Thus, we can compute the p-value of the pairwise unions of the two sets of partitions of ζ, except the irrelevant ζ e ∪ α = ∅. We only provide p-values for these pairwise unions and do not attempt to provide p-values for the partitions for the different consequences of the mutations on splicing. While such values would be useful, we do not currently have a robust means to compute them. Our previous work provides guidance on interpretation of splicing mutation outcomes 3-5,10 . Thus for ζ x ∈ {ζ c , ζ e , ζ i }, let Φ Z (z) represent the cumulative distribution function of the one-sided (right-tailedi.e. P[X > x]) standard normal distribution. Let N represent the total number of samples and let V represent the set of all ζ x validations, across all samples. Then: The program outputs two tables, along with summaries thereof. The first table lists all validated read counts across all categories for experimental samples, while the second table does the same for the control samples. P-values are shown in parentheses within the experimental table, which refer to the column-dependent (i.e. the category is given in the column header) p-value for that category with respect to that same category in control samples. The program produces three files: a log file containing all details regarding validated variants, an output file with the programs progress reports and summaries, and a filtered validated variant file. The filtered file contains all validated variants of statistical significance (set as p < 0.05, by default), defined as variants with one or more validating read categories which achieve statistical significance in a relevant category (i.e. a cryptic variant for which p = 0.04 in the junction-spanning cryptic column would meet this criteria).
We elected to use RefSeq 19 genes for the exome annotation, as opposed to, the more permissive exome annotation sets, UCSC Known Genes 20 or Ensembl 21 . The large number of transcript variants within Ensembl, in particular, caused many spurious intron inclusion validation events. This occurred because reads were found to be intronic in many cases, when in actuality they were exonic with respect to the more common transcript variant. In addition, the inclusion of the large number of rare transcripts in Ensembl significantly increased program runtime and made validation events much more challenging to interpret unequivocally. The use of RefSeq, which is a conservative annotation of the human exome, resolves these issues.
We also provide an R program 22 which produces publication quality histograms displaying embedded Q-Q plots and p-values, to evaluate for normality of the read distribution and statistical significance, respectively. The R program performs the YJ transformation as implemented in the car package 23 . The histograms generated by the program use the Freedman-Draconis 24 rule for break determination, and the Q-Q plots use algorithm Type 8 for their quantile function, as recommended by Hyndman and Fan 25 . Lastly, a Perl program was implemented to automatically retrieve and correctly format an exome annotation file from the UCSC database 20 for use in Veridical. All data uses hg19/GRCh37, however when new versions of the genome become available, this program can be used to update the annotation file.

Results
Veridical validates predicted mRNA splicing mutations using highthroughput RNA sequencing data. The performance of the software is affected by the number of predicted splicing mutations, the number of abnormal samples containing mutations and control samples and the corresponding RNA-Seq data for each type of sample.
Veridical has the ability to analyze approximately 3000 variants in approximately 4 hours assuming an input of 100 BAM files of RNA-Seq data. The relationship between time and numbers of BAM files and variants are plotted in Figure 3 for a 2.27 GHz processor. Veridical uses memory in linear proportion to the number and size of the input BAM files. In our tests, using RNA-Seq BAM files with an average size of approximately 6 GB, Veridical used approximately 0.7 GB for ten files to 1 GB for 100 files. Currently, splicing consequences that are reported include intron inclusion, exon skipping, and cryptic splicing, which are validated through junction-spanning reads, or based on read-abundance in the region circumscribing the variant (see Methods for details). For example, a cryptic splicing junction-spanning read will show that the mRNA contains a truncated or extended exon at the predicted location, which is directly attached to the sequence of the corresponding adjacent exon. For mutations that alter read-abundance, each read within the genomic location assessed (i.e. intron for intron inclusion) is counted for the variant-containing samples and then compared with the number of reads in the control files. For each input variant, Veridical outputs the number of validating reads (i.e. RNA-Seq reads which corroborate the predicted splicing consequence) for a given splice consequence within the variant-containing tumour samples and within control samples (i.e. non-variant containing tumour samples and normal samples). The program provides read counts for the different categories for all experimental and control samples as tab-delimited tables, along with the relevant p-values, indicating the statistical probability that the predicted mutation exhibits a normal expression pattern.
We demonstrate how Veridical and its associated R program are used to validate predicted splicing mutations in somatic breast cancer. Each example depicts a particular variant-induced splicing consequence, analyzed by Veridical, with its corresponding significance level. The relevant primary RNA-Seq data are displayed in IGV, along with histograms and Q-Q plots showing the read distributions for each example. The source data are obtained from controlled-access breast carcinoma data from The Cancer Genome Atlas (TCGA) 28 . Tumour-normal matched DNA sequencing data from the TCGA consortium was used to predict a set of splicing mutations, and a subset of corresponding RNA sequencing data was analyzed to confirm these predictions with Veridical. The following examples demonstrate the utility of Veridical to identify potentially pathogenic mutations from a much larger subset of predicted variants. Leaky mutations are those variants that reduce, but not abolish, the spliceosome's ability to recognize the intron/exon boundary 3 . This can lead to the mis-splicing (intron inclusion and/or exon skipping) of many but not all transcripts. An example of a predicted leaky mutation (chr5:162905690 G>T) in the HMMR gene in which both junction-spanning exon skipping (p < 0.01) and read-abundancebased intron inclusion (p = 0.04) are observed is provided within Figure 4. We predict this mutation to be leaky because its final R i exceeds 1.6 bits -the minimal individual information required to recognize a splice site and produce correctly spliced mRNA 4 . Indeed, the natural site, while weakened by 2.16 bits, remains strong -10.67 bits. This prediction is validated by the variantcontaining sample's RNA-Seq data (Figure 4), in which both exon skipping (5 reads) and intron inclusion (14 reads) are observed, along with 70 reads portraying wild-type splicing. Only a single normally spliced read contains the G→T mutation. These results are consistent with an imbalance of expression of the two alleles,  as expected for a leaky variant. Figure 5 shows that for the distribution of read-abundance-based intron inclusion is statistically significant (p = 0.04).

Input, output, and explanatory files for Veridical
Variants that inactivate splice sites have negative final R i values 3 with only rare exceptions 4 , indicating that splice site recognition is essentially abolished in these cases. We present the analysis of two inactivating mutations within the PTEN and TMTC2 genes from different tumour exomes, namely: chr10:89711873 A>G and chr12:83359523 G>A, respectively. The PTEN variant displays junction-spanning exon skipping events (p < 0.01), while the TMTC2 gene portrays both junction-spanning and read-abundancebased intron inclusion (both splicing consequences with p < 0.01). In addition, all intron inclusion reads in the experimental sample contain the mutation itself, while only one such read exists across all control samples analyzed (p < 0.01). The PTEN variant contains numerous exon skipping reads (32 versus an average of 2.466 such reads per control sample). The TMTC2 variant contains many junction-spanning intron inclusion reads with the G→A mutation (all of its junction-spanning intron inclusion reads: 22 versus an average of 0.002 such reads per control sample). IGV screenshots for these variants are provided within Figure 6. This figure also shows an example of junction-spanning cryptic splice site activated by the mutation (chr1:985377 C>T) within the AGRN gene. The concordance between the splicing outcomes generated by these mutations and the Veridical results indicates that the proposed method detects both mutations that inactivate splice sites and cryptic splice site activation.  tumours, two had significant levels of junction-spanning intron inclusion, and one showed statistically significant read-abundancebased intron inclusion. Details for all of the aforementioned variants, including a summary of read counts pertaining to each relevant splicing consequence, for experimental versus control samples, are provided in Table 2.

Discussion
We have implemented Veridical, a software program that automates confirmation of mRNA splicing mutations by comparing sequence read-mapped expression data from samples containing variants that are predicted to cause defective splicing with control samples lacking these mutations. The program objectively evaluates each mutation with statistical tests that determine the likelihood of and exclude normal splicing. To our knowledge, no other software currently validates splicing mutations with RNA-Seq data on a genome-wide scale, although many applications can accurately detect conventional alternative splice isoforms (i.e.29). Veridical is intended for use with large data sets derived from many samples, each containing several hundred variants that have been previously prioritized as likely splicing mutations, regardless of how the candidate mutations are selected. It is not practical to analyze all variants present in an exome or genome, rather only a filtered subset, due to the extensive computations required for statistical validation. As such, Veridical is a key component of an end-to-end, hypothesis-based, splicing mutation analysis framework that also includes the Shannon splicing mutation pipeline 10 and the Automated Splice Site Analysis and Exon Definition server 5 .
There is a trade-off between lengthy run-times and statistical robustness of Veridical, especially when there are either a large number of variants or a large number of RNA-Seq files. As with most statistical methods, those employed here are not amenable to small sample sets, but become quite powerful when a large number of controls are employed. In order to ensure that mutations can be validated, we recommend an excess of control transcriptome data relative to those from samples containing mutations (> 5 : 1), regardless of the computational expense. We do not recommend the use of a single control to corroborate a sample containing a putative mutation. Not surprisingly, we have found that junction-spanning reads have the greatest value for corroborating cryptic splicing and exon skipping. Even a single such read is almost always sufficient to merit the validation a variant, provided that sufficient control samples are used. For intron inclusion, both junction-spanning and read-abundancebased reads are useful and a variant can readily be validated with either, provided that the variant-containing experimental sample(s) show a statistically significant increase in the presence of either form of intron inclusion corroborating reads.
Veridical is able to automatically process variants from multiple different experimental samples, and can group the variant information if any given mutation is present in more than one sample. The use of a large sample size allows for robust statistical analyses to be performed, which aid significantly in the interpretation of results. The main utility of Veridical is to filter through large data sets of predicted splicing mutations to prioritize the variants. This helps to predict which variants will have a deleterious effect upon the protein product. Veridical is able to avoid reporting splicing changes that are naturally occurring through checking all variant-containing and non-containing control samples for the predicted splicing consequence. In addition, running multiple tumour samples at once allows for manual inspection to discover samples that contained the alternative splicing pattern, and consequently, permits the identification of DNA mutations in the same location which went undetected during genome sequencing.
The statistical power of Veridical is dependent upon the quality of the RNA-Seq data used to validate putative variants. In particular, a lack of sufficient coverage at a particular locus will cause Veridical to be unable to report any significant results. A coverage of at least 20 reads should be sufficient. This estimate is based upon alternative splicing analyses in which this threshold was found to imply concordance with microarray and RT-PCR measurements 30-33 . There are many potential legitimate reasons why a mutation may not be validated: (a) nonsense-mediated decay may result in a loss of expression of the entire transcript, (b) the gene itself may have multiple paralogs and reads may not be unambiguously mapped, (c) other nonsplicing mutations could account for a loss of expression, and (d) confounding natural alternative splicing isoforms may result in a loss of statistical significance during read mapping of the control samples. The prevalence of loci with insufficient data is dependent upon the coverage of the sequencing technology used. As sequencing technologies improve, the proportion of validated mutations is expected to increase. Such an increase would mirror that observed for the prevalence of alternative splicing events 34 . It is important to note that acceptance of the null hypothesis, due to an absence of evidence required to disprove it, does not imply that the underlying prediction of a mutation at a particular locus is incorrect, but merely that the current empirical methods employed were insufficient to corroborate it.
While there is considerable prior evidence for splicing mutations that alter natural and cryptic splice site recognition, we were somewhat surprised at the apparent high frequency of statistically significant intron inclusion revealed by Veridical. In fact, evidence indicates that a significant portion of the genome is transcribed 34 , and it is estimated that 95% of known genes are alternatively spliced 30 . Defective mRNA splicing can lead to multiple alternative transcripts including those with retained introns, cassette exons, alternate promoters/terminators, extended or truncated exons, and reduced exons 35 . In breast cancer, exon skipping and intron retention were observed to be the most common form of alternative splicing in triple negative, non-triple negative, and HER2 positive breast cancer 36 . In normal tissue, intron retention and exon skipping has been predicted to affect 2572 exons in 2127 genes and 50 633 exons in 12 797 genes, respectively 37 . In addition, previous studies suggest that the order of intron removal can influence the final mRNA transcript composition of exons and introns 38 . Intron inclusion observed in normal tissue may result from those introns that are removed from the transcript at the end of mRNA splicing. Given that these splicing events are relatively common in normal tissues, it becomes all the more important to distinguish expression patterns that are clearly due to the effects of splicing mutations -one of the guiding principles of the Veridical method. Veridical is an important analytical resource for unsupervised, thorough validation of splicing mutations through the use of companion RNA-Seq data from the same samples. The approach will be broadly applicable for many types of genetic abnormalities, and should reveal numerous, previously unrecognized, mRNA splicing mutations in exome and complete genome sequences.
Author contributions PKR elaborated the problem and conceived of the analytic solution. CV implemented the algorithm, wrote, and tested the Veridical software and its accompanying Perl and R scripts. CV and SND designed the methods used, generated, and interpreted mutation results. BCS wrote and tested the Shannon pipeline, which provided the predictions evaluated by Veridical. BCS and CV implemented procedures and code for software distribution. CV, SND, and PKR wrote the manuscript, which has been approved by all of the authors.
Competing interests PKR is the inventor of US Patent 5,867,402 and other patents pending, which underlie the prediction and validation of mutations. He is one of the founders of Cytognomix, Inc. which is developing software based on this technology for complete genome or exome splicing mutation analysis. BCS is an employee of Cytognomix, Inc. CV and SND declare that they have no competing interests. Both trial and licensed versions of Veridical are fully functional, however the length of the trial period is limited. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Acknowledgments
We acknowledge the TCGA consortium for providing the data used in this article. The data used in this article is part of TCGA project #988: "Predicting common genetic variants that alter the splicing of human gene transcripts". This work was made possible by the facilities of the Shared Hierarchical Academic Research Computing Network (SHARCNET) and Compute/Calcul Canada.

Supplementary materials Veridical variant input format
This input format most easily accepts formatted output from the Shannon Pipeline. In particular, all variants of interest should be concatenated into a single file. Once a, tab-delimited, concatenated file has been generated, it can easily be formatted correctly by using FilterShannonPipelineResults.pl. One can also manually ensure the following: the header line has no quotation marks or special characters, empty columns have been replaced by a period (.) and each variant line contains only a single gene (comma-delimited gene lists must be split such that there is only one gene per line). If one wishes Veridical to consider variants pertaining to more than one experimental sample, a comma-delimited list of experimental samples, in the form of BAM file names, must be provided as the key column. The key column must always contain at least one file name that is present as the base name of one of the files listed in the BAM file list that must be passed to Veridical.
Alternatively, one can prepare the input format as follows. The header must contain at least the following, case-insensitive, values to which the file's columns must adhere to: chromosome, splice&coordinate, strand, type, gene, location, location_type, heterozygosity, variant, input, key. The column headers need only contain the given text (i.e. a column labeled gene_name would be sufficient to satisfy the above requirement for a "gene" column). Column headers with ampersands (&) denote that all words joined by this symbol must be present for that column (i.e. Splice_site_coordinate satisfies the "splice&coordinate" requirement). The order of the columns is immaterial. The input column can contain any identifier for the variant and need not be unique. The location column specifies if the site is natural or cryptic. For Veridical, all that matters is that cryptic variants contain the word "cryptic" as part of their value in this column and that non-cryptic variants do not. The location_type column is only used for cryptic variants and specifies if the variant is intronic or exonic. It is not currently used by the program. This column must be present but can always be set to null (i.e.).
A few rows from a sample variant file is provided below (text wrapped for readability): Veridical exome annotation input format This input format can be generated via ConvertToExomeAnotation.pl. The file must be tab-delimited, excepting its header, which must be comma-delimited. It must have the following, case-insensitive, header columns, to which its data must adhere: transcript, chromosome, exon chr start, exon chr end, exon rank, gene. The column headers need only contain the given text (i.e. a column labeled gene_name would be sufficient to satisfy the above requirement for a "gene" column). The order of the columns is immaterial.

Veridical output
If a variant contains any validating reads, Veridical outputs the variant in question, along with some summary information and a table specifying the numbers of each validating read type detected for both the experimental and control samples. Within the output of Veridical, the phrase: "Validated (x) variant n times" means that the variant was validated mainly for splicing consequence x and has n validating reads. The variant will only appear within the *.filtered output file if the p-value for either junction-spanning or read-abundance-based reads for splicing consequence x was statistically significant (defined, by default, as: p < 0.05). After the variant being validated is provided, along with its primary predicted splicing consequence, the output is divided into two sections with identical contents: one for the experimental sample(s) and another for control samples. The summary enumerates the number of reads of each splicing consequence, partitioned by evidence type (junction-spanning or read-abundance-based), and by sample type (tumour or normal for control samples, and only tumour for experimental samples). A table describing the number of each read type for every file follows this summary. An example of this output, for the variant within RAD54L, as shown by Figure 7 and the last portion of Table 2, is provided. While Veridical outputs this as plain text, with the table in a tab-delimited format, we provide this output as an Excel document with descriptions of the meaning of each table heading, to clarify the presentation of the data. All input and output files for the five variants presented are provided.

VeridicalOutExample.xls
contains the output for the variant within RAD54L, along with descriptions of the terms used and the output format. all.vin contains the input variant file.
allTumoursBAMFileList.txt and allNormalsBAMFileList.txt are the BAM file lists for tumour and normal samples, respectively. all.vout contains the Veridical output. The exome file can be retrieved using ConvertToExomeAnnotation.pl, available with the other programs at: www.veridical.org. The BAM file lists contain the TCGA file UUID, followed by a slash, followed by the file name. The RNA-Seq data itself can be downloaded from TCGA at: https:// tcga-data.nci.nih.gov/tcga/.  18 March 2014 Referee Report: This paper from the Rogan group presents a methodology for validation of DNA sequencing variants that alter mRNA splicing. While variants of the most conserved splice site nucleotides at the intron-exon boundary can be predicted to cause splice defects with high reliability, it remains difficult to predict whether variants deeper in the intron or those that potentially affect exonic splicing enhancers actually cause splice defects. RNA-seq data, when coupled with variant data, potentially provide a means of correlating variation data with observations of (mis-)splicing patterns.
The program fulfils an important need in the community, the results appear promising and will be of special interest to groups performing RNA-seq analysis in medical settings. I have only some minor suggestions that the authors may like to consider.

Suggestions:
The explanation of the methodology is relatively difficult to follow, and I wonder if it might not be better to simplify Figures 1 and 2 for didactic sake. For instance, in Figure 1A, it is unclear where the location of variant is. Does the curved line mean that it could be anywhere in the middle C exon? Also, I assume that exons are being shown in blue and reads shown in gray? Also, the legend text is overly complicated: > swap and . While aficionados of first order D E D E logic will follow without problems, I would suggest that it would be better for didactic purposes to delete this and to implicitly assume that < for the sake of this figure. Figure 1B is confusing at D E this point in the manuscript because the motivation for switching the variable , , , and is not A B D E yet clear. On the other hand, panel C and panel D are trivial and do not add anything. I would suggest using Figure 1 to provide one concrete example one a simple level, and stating in the text that the variables are to be switched if the candidate mutation is located on the other side of the exon.
Also, the explanations of the method that are couched in first order logic-like notation are difficult to follow, because it is not stated whether the variant can precede the start of the read (in which C case would be negative). The subscripts for in turn have the subscript s but the variable in C-S r 1 S the formula does not.
Although in the end, I think I follow the overall method, the reader is forced to make arbitrary assumptions in order to interpret the formulae being used to explain the method. A similar comment pertains to the flow chart in Figure 2.
Therefore, I would suggest the authors take some pains to improve the clarity of the explanation of F1000Research 1.

Approved with reservations: 07 March 2014
07 March 2014 Referee Report: The authors describe a method and the associated software, Veridical, for assessing the effects on pre-mRNA splicing of predicted splicing-affecting mutations. To do so the program compares splicing effects, measured by the supporting read counts, in variant-containing (disease) samples against a distribution derived from very large numbers of 'normals', either normal tissue from the same individual or samples from healthy individuals.
The idea is ingenious and novel as applied to mutations affecting splicing, although not in general (see [ ], which exploits the availability of large numbers of samples to identify likely VAAST Yandell et al. 2011 deleterious variants; it is also the premise for the and projects). The software is HapMap 1000 Genomes fast and practical, being able to test thousands of variants in hundreds of samples within hours. This is the first software of its kind, and if accurate it will be a very valuable resource for clinical genomics.
That being said, while the article provides proof-of-concept and clearly demonstrates the potential of the tool with specific examples, there are several missing pieces that are needed to provide the readers with a view of its overall performance and limitations and to help them use it effectively.

Major comments:
The article shows numerous positive examples, however there is no indication of the tool's performance in general. The authors should include the results from running the tool on a full data set, to give potential users an idea of the expected outcome.
Also, several other tools (e.g., , , ) have been developed for the related MATS Miso SpliceTrap problem of discovering alternative splicing events and comparing them among samples. MATS in particular, allows differential splicing analyses with multiple replicates. Ideally the paper would include a comparison with MATS on the data set analyzed; this comparison is informative even if MATS is used with only a subset of the samples.
The method uses the YJ-transformed distribution of supporting read counts across the 'normals' to determine a p-value for the variant, and thus judge its significance and impact on splicing. This is an interesting concept that assumes that with large numbers of 'normals' sample and batch effects will even out; hence, large numbers of samples are required to ensure accuracy. Since these are absolute (non-normalized) counts, however, the method may not work if the variant sample is obtained with a different method, e.g. by rRNA depletion of total RNA whereas most normal samples would come from polyA+ libraries. The authors should clearly discuss this and other possible limitations of their approach.
Related to the above, the authors mention on several occasions the difficulty in identifying intron inclusion (II) events, in particular the large number of false positives. Indeed, IIs are generally difficult to predict due to the presence of intronic reads ('noise') from unspliced RNA. The levels can vary from sample to sample and across the genome, depending on the sample preparation, gene expression level, splicing efficiency, etc. By comparing read counts exclusively among samples and without taking into account the gene-or genome-level background, Veridical is likely to produce many false positives.
In particular, the 14 supporting reads in the left intron on Figure 4 seem hardly sufficient to indicate an II event, all the more as there is a larger number of reads in the neighboring intron (not predicted trans-acting splicing factors could account in part or totally for the splicing changes across samples. The statistical tests properly conducted in Veridical are designed to minimize such possibility but do not rule it out. In addition, the inherent noisy nature of RNAseq datasets also prompts for caution in the conclusions. To me, the direct proof that a DNA mutation changes splicing of its pre-mRNA can only be provided using minigenes and cell transfection (or in vitro splicing), in which the substrate sequences and cellular context are under almost absolute control. Indeed, the Veridical method is reminiscent of GWAS (Genome-Wide Association Studies), in which the genotype in the DNA, wild-type or mutant, is associated to its phenotype, such as normal versus disease (or other traits) in GWAS, or normal versus aberrant splicing in this study. Thus, for me Veridical provides strong associations -but not validations -between DNA mutations and their effects on splicing.
As mentioned briefly at the beginning of Discussion, Veridical has built-in prediction tools to prioritize the mutations that are more likely to affect splicing, such as those mapping to splice sites. Even if other sources and tools are cited, a more extensive explanation of these components of Veridical would help the reader/user.

I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
No competing interests were disclosed. Competing Interests:

Stefania Bortoluzzi
Department of Biology, University of Padova, Padova, Italy

Approved with reservations: 27 January 2014
27 January 2014 Referee Report: The paper " " by Validation of predicted mRNA splicing mutations using high-throughput transcriptome data Viner . presents Veridical, a new software for the interpretation and validation of genetic variants et al identified by DNA sequencing that alter mRNA splicing, leveraging RNA-seq data. The method is based on statistical comparisons of the normalized read counts of abnormally-spliced RNA species in mutant versus non-mutant tissues.
Actually, the interpretation of genetic variants is a difficult and key issue in current research. The integration of genomic and transcriptomic data, namely the use of RNA-seq-based transcriptome characterization as a "molecular phenotype" of cells is useful and meaningful.
The software is standalone (not a web-tool) and it is completed by perl scripts, facilitating data management. The manuscript declare that "Veridical and its associated software programs are available at: www.veridical.org".
Actually, Veridical is commercially available to the scientific community. A trial version lasting 30 days can be downloaded by the website, but in order to obtain binaries, the website requests a registration with an institutional email address -they reserve the right to deny access to users who register with third-party mail servers (Gmail, Yahoo, Hotmail, etc.).
No pricing information is included in the manuscript and, more importantly, in the webpages accessible to download the software, either before or after registration.

Methods
2 Par, Line 5 (minor): "Maximising" is used, but probably the meaning is "increasing" (the number nd of). In general, in many cases in the manuscript, the correspondence between legend and figure can be improved, by indicating more clearly the points specific sentences in the legend refer to. Regarding this issue, for instance in Figure 4 I can see indicated neither the "exon 12" nor the "14 reads" mentioned in the legend. Please indicate (using colors, boxes, arrows or overlapping text) key elements in the figure, and revise all figures using the same criterion.
Page 5 (minor): Consider revising the sentence "Furthermore, this transformation allowed us to avoid the use of non-parametric testing, which has its own pitfalls regarding assumptions of the underlying data distribution", since normally it is assumed that parametric tests ground on assumptions on data distributions, but non-parametric tests by definition can be used without information about data distribution.
End of the next paragraph (major): "It is important to realize, therefore, that the p-values given by Veridical are much more robust when the program is provided with a large number of samples." This is a pretty clear concept. Please, indicate a general rule to the user/reader: How many samples are required? Setting a reasonable minimum can be more useful for experimental design than saying the larger the sample size the most robust the result.

I have two important criticisms about the Results section:
The section is not organized in paragraphs, and mixes performance info (run time using different