Splicing Express: a software suite for alternative splicing analysis using next-generation sequencing data

Sequences

The human genome reference sequence (GRCh37/hg19) was downloaded from UCSC Genome Bioinformatics portal (http://genome.ucsc.edu), as well as RefSeq sequences and their respective annotations (files refFlat.txt.gz and refSeqAli.txt.gz). RNA-Seq data from the Illumina Human Body Map Project Version 2 were downloaded from the EMBL-EBI portal (accession E-MTAB-513). All available tissues were used: adrenal, adipose, brain, breast, colon, heart, kidney, liver, lung, lymph, ovary, prostate, skeletal muscle, testes, thyroid and white blood cells. Analyses were restricted only to the 2 × 50 bp paired-end sequences, which represent approximately 80 million read pairs per tissue.

Transcriptome assembly

Reads from all 16 tissues from the Illumina Human Body Map Project were mapped to the genome using Tophat version 2 (Trapnell, Pachter & Salzberg, 2009). A conservative parameter was maintained (–splice-mismatches 0), and a reference GTF file, containing RefSeq aligments from the UCSC portal, was set for helping the mapping process. The resulting BAM files created by Tophat were, afterwards, used as input for Cufflinks version 2 (Roberts et al., 2011). Use of that software was carried out using default parameters and without any reference file as input. Additional processing, such as sequence clustering, was performed by Splicing Express.

Sequencing clustering

The expressed sequences were clustered by Splicing Express, using a strategy that has been used previously by us (Galante et al., 2004; Galante et al., 2007; Kroll et al., 2012). Basically, the genome coordinates of all expressed sequences were compared to the genome coordinates of Reference Sequences (RefSeqs) entries. Sequences presenting multiple exons and sharing at least one exon-intron boundary with a RefSeq sequence were merged. In the specific case of Splicing Express, all remaining sequences were excluded from the alternative splicing analysis, but were taken into account for measuring the total gene expression. The intronless sequences were merged by checking their overlap to RefSeqs, which should be greater than or equal to 30 bp.

Alternative splicing analysis

The GTF files resulting from the transcriptome assembly of the Illumina Human Body Map Project were processed using Splicing Express. The expressed sequences were clustered to a set of RefSeq sequences, which comprehended 26.318 genes and a total of 50.129 entries. The processing steps are better described in the “Implementation.” In-house Perl scripts were developed for all large-scale analyses, and all events described in the “Analysis” were visually identified by using Splicing Express. A basic analysis using Splicing Express was also performed for the Rat Body Map (Yu et al., 2014), which is available at the Splicing Express website.

Statistics

Splicing Express performs statistical analysis aiming to give significance values to the fold differences between the splicing variants. That analysis is only performed when the GTF files have data regarding the FPKM confidence intervals, which are standardly present in any GTF file created by Cufflinks (Roberts et al., 2011). Those values are gathered to calculate the FPKM Standard Error (SE), given by the equation: ${\displaystyle \mathrm{SE}=\left(C_{\mathrm{hi}}-C_{\mathrm{lo}}\right)/\left(2\ast 1.96\right),}$ where C_hi and C_lo are the FPKM values at the lower and higher limit of the 95% FPKM confidence interval, respectively. Moreover, it is known that more than one expressed sequence can share the same splicing event in a sample. Hence, the FPKM of those expressed sequences have to be summed. Its respective SE is then calculated by the equation: ${\displaystyle {\mathrm{SE}}_x=\sqrt{{\mathrm{SE}}_1^2+\cdots +{\mathrm{SE}}_n^2},}$ where x refers to a splicing cluster and n refers to the expressed sequences from the respective cluster x. Finally, a Welch’s t-test is performed between the two clusters of splicing variants, defined here as a and b. The t-test equation is given by: ${\displaystyle t=\left(X_a-X_b\right)/\sqrt{{\mathrm{SE}}_a^2+{\mathrm{SE}}_b^2},}$ where X is the FPKM of the sample. Probabilities are then calculated from t-values through the Perl Statistics::Distributions library using a degree of freedom equal one.

Availability

Splicing Express source code is available through a repository manager platform at http://www.bioinformatics-brazil.org/splicingexpress, where tutorials (written and video), sample files, and a bug track system are provided.

Implementation

Splicing Express is a suite composed of Perl scripts that executes several processing steps, as schematized in Fig. 1B. Prior to running Splicing Express, parameters about the reference genome, the output directory and the input GTF files must be firstly set. After that, GTF files are checked for consistency and converted to a workable file format. GTF files created by Cufflinks are recommended to be used as input due to some annotation standards referring to sequence names and FPKM expression. Similarly, chromosome names and genome positions from the GTF files must be compatible with the respective genome data provided by the UCSC genome portal.

Expressed sequences are then clustered through a known strategy (see “Material & Methods”). Clustering is needed for identifying the gene of origin of each expressed sequence contained in the GTF files, and is also very important to standardize the annotation process using gene symbols according to the HUGO (Gene Nomenclature Committee). Splicing Express, therefore, allows the use of GTF files showing any kind of gene annotation since those information will be posteriorly rewritten.

The ASEs identification is done through an algorithm based on pairwise comparison, a very common approach used for many other similar tools (Florea, Song & Salzberg, 2013; Kroll et al., 2012; Sammeth, Foissac & Guigó, 2008). Splicing Express, moreover, uses an strategy implemented previously by us in Splooce (Kroll et al., 2012), which is based on regular expressions that are capable of identifying all well-known simple and combined ASEs. In that strategy, the expressed sequences are represented as binary sequences, where exons and introns are represented by 1 and 0, respectively. Those sequences are then pairwisely compared in an exhaustive manner, creating numerical patterns which represent their splicing differences. Currently, Splicing Express is able to identify patterns of simple alternative splicing events, such as exon skipping, intron retention and alternative 5′ and 3′ splicing borders.

Afterwards, a graphic representation for the expression level (FPKM) is created for each gene and for each identified ASE. Graphics are created using an external SVG library, which is available at CPAN (Comprehensive Perl Archive Network). They can be opened in all recent releases of internet browsers, zoomed in or out without loss of quality. Finally, those graphics are gathered together with additional information in HTML files. An example of some graphical displays are shown in Figs. 1C and 1D. A main file (index.html) is created at the root of the results directory tree, providing advanced tools for exploring and filtering events (Fig. 1A). A basic statistical summary and a help section are also provided.

Performance

Splicing Express analysis does not aim speed or efficiency since the major computational processing is done previously by tools like Tophat and Cufflinks. Overall, Splicing Express processing time is fast, and it is able to process 15 samples in about half-hour using an average desktop personal computer. Here, we performed a simple performance analysis. No accuracy and sensitivity parameters were measured because the ASE identification algorithm here proposed is based on an exhaustive approach. As a testing platform, an i7 Core processor with 16 GB RAM machine was chosen, although Splicing Express uses only a single core and needs no more than 1 Gb of RAM memory.

Sixteen GTF files representing each tissue from the Illumina Human Body Map Project were analyzed. The performance was measured for each processing step of Splicing Express as a function of the number of analyzed samples (Fig. 2). The processing was done in triplicate and the results shown represent their average. Total processing time increased almost linearly as more samples were added to the analysis. The linear trend is particularly shaped by the processing time of “Data preparation” and “Annotation” steps. The “ASE identification” step, on the other hand, is based on a pairwise strategy, which may consume exponentialy more processing time as more samples are analyzed. The “ASE identification” processing, however, was very fast and the exponential trend could be barely observed. In turn, the “HTML & Graphics” processing time theoretically decreases as the identification of new alternative splicing events becomes saturated (less graphics are created), meaning that it tends to stabilize as more samples are analyzed. However, it could not be observed for the current data. “HTML & Graphics” step also depends on a moderate use of disk I/O due to a large amount of files being written. Hence, slow downs may be observed when outputing Splicing Express results directly to pen drives or other kind of slow devices.

Analysis

Data from the Illumina Human Body Map Project was used for validating the use of Splicing Express in finding specific events among different samples. In total, 16 distinct tissues, showing high sequence coverage, were processed. As a result, 43,307 events were identified, from which 36% comprehended well-known events (Table 1). Meanwhile, partially known and unknown events accounted for 64% of the data, as expected (Florea, Song & Salzberg, 2013). Known events were defined as any event showing both variants supported by a RefSeq. Partially known events, in turn, were defined as events showing only one of both variants supported by a RefSeq. In order of frequency, taking into account only known events, the most frequent event type was exon skipping, followed by alternative splicing borders and intron retention, as observed before (Kroll et al., 2012). A more detailed analysis, considering tissue-specific alternative splicing events (tsASEs), showed that brain, adrenal gland, lymphnode and testis, among all analyzed tissues, had the highest amount of specific variants (t(df = 14) = 9.644117, p < 0.00001) (Table 2). TsASEs were defined as any splicing variant showing a relative expression greater than or equal to 75% of the total expression in all analyzed tissues. Only TsASEs showing a FPKM higher or equal than 1 were selected.

To further explore the use of Splicing Express, we focused on the analysis of brain, since it is known to show one of the highest levels of skipped exons, alternative 5′ and 3′ splice sites (Yeo et al., 2004). Moreover, a high frequency of tsASEs was found for brain (Table 2). A gene ontology (GO) enrichment analysis of biological processes using DAVID (Huang, Sherman & Lempicki, 2009) showed that genes showing tsASEs for brain are more enriched with terms related to brain functions, such as “axonogenesis,” “cell morphogenesis involved in neuron differentiation” and “neuron projection morphogenesis” (p-value <6 × 10⁻³) (Fig. 3). MACF1, for example, is a microtubule actin crosslinking factor known to determine neuronal positioning by regulating microtubule dynamics and mediating GSK-3 signaling during brain development (Ka et al., 2014), and it was observed to be more expressed in 6 other tissues than brain. However, three ASEs were specific to brain (two exon skipping and one 3′ alternative splicing site). In genomic scale analysis, examples like that could only be related to a specific tissue through meticulous alternative splicing analysis, as that provided by Splicing Express.

Among other events, an interesting one was observed for NRXN3, a neurexin related to presynaptic cell adhesion. In brain, the skipping/inclusion of exon SS#4 is known to control the postsynaptic AMPA receptor trafficking, being therefore an important factor for synaptic plasticity (Aoto et al., 2013). The specificity of NRXN3 and of its event to brain could be easily observed (Figs. 1C and 1D). On the other hand, many genes do not show tsASEs, such as the CELF family of proteins, which are involved in cell-specific and developmentally regulated alternative splicing in brain (Ladd, Charlet-B & Cooper, 2001). Splicing Express, however, was able to confirm the specificity of CELF3 and CELF5 to brain through their respective gene expression. For more informations about those and other genes, a full analysis of the Human Body Map and Rat Body map (Yu et al., 2014) using Splicing Express were made available at http://www.bioinformatics-brazil.org/ASE_HBM.