Bioinformatic pipelines for whole transcriptome sequencing data exploitation in leukemia patients with complex structural variants

Data generation and preprocessing

Our pipelines were developed and tested on total RNA-Seq data generated from ten CLL cases (T1-T10) with cSVs that were previously detected using genomic CytoScan™ HD Arrays and analyzed with Chromosome Analysis Suite (Thermo Fisher Scientific). These cases were identified during the long-term clinical research at the University Hospital Brno and were classified as chromothripsis as they showed clustered copy number alterations on a limited number of chromosomes. Genomic breakpoint localization was extracted for all detected copy number variants and losses of heterozygosity. Total RNA-Seq libraries were prepared with TruSeq^® Stranded Total RNA kit (Illumina) with Ribo-Zero ribosomal RNA depletion and sequenced using an Illumina HiSeq 2500 machine producing 125bp long pair-end reads. Sequencing read quality was evaluated in FastQC software. Adapter sequences were trimmed from raw reads using Trimmomatic software (Bolger, Lohse & Usadel, 2014) (version 0.32) according to sequencing facility standards. In parallel, we performed GeneChip^® Human Transcriptome Arrays 2.0 (Thermo Fisher Scientific) to complement expression data from the RNA-Seq experiment. These data were then analyzed by the Transcriptome Analysis Console (Thermo Fisher Scientific) by comparing each consecutive sample to all other samples using a one-way ANOVA statistical test to identify a unique expression pattern in each sample. Identified fusion gene junctions and their respective sequencing reads were visualized in the Integrative Genomic Viewer (Thorvaldsdottir, Robinson & Mesirov, 2013).

The study was approved by the Ethical Committee of the University Hospital Brno under the ref. no. 15-31834A. All patients involved in the study provided their written informed consent to the research use of their samples.

RNA-Seq data processing for differential genes expression analysis

Processed reads were mapped to the hg38 human genome reference using STAR, a splice-aware aligner (Dobin et al., 2013). We chose the genome reference over the transcriptome as recommended by best practices for RNA-Seq data when non-canonical junctions and fusion transcripts are of interest (Conesa et al., 2016). The parameters of the STAR aligner were set to default settings according to best practices.

Gene expression analysis from RNA-Seq data is based on counting reads covering gene regions (defined by reference relative GTF file) and statistical analysis of these read counts. In the first step, we used an htseq-count script (Anders, Pyl & Huber, 2015), which efficiently counts reads that align to or overlap with more than one gene. For this purpose, the software parameter union was applied. The reads were then assigned to gene regions irrespective of DNA strand orientation. A read count table was used as an input for consequent statistical analysis of differential gene expression; low-expressed genes were filtered out from the dataset to decrease dataset complexity and avoid FP. Of the ten CLL cases in our cohort, we expected at least six samples to be covered by at least one read. However, these criteria can be modified in the pipeline according to the actual dataset analyzed. Pre-filtered read count table was then normalized using the rpkm() function from edgeR Bioconductor package (Robinson, McCarthy & Smyth, 2010). The function applies the RPKM (reads per kilobase per million mapped reads) method, which performs data normalization based on a read length and a total number of sequencing reads with respect to gene length. Obtained log₂ scale normalized expressions were subjected to further computational steps, for which we developed a novel statistical pairwise comparison (PComp) approach that allows obtaining a gene expression profile of individual samples.

Design and evaluation of novel statistical approach PComp

For the read count (expression) table resulting from RNA-Seq data processing steps, linear regression model with confidence bands containing expression values for all possible sample combinations (i.e., T1–T2, …, T1–T10) was applied to identify gene outliers in particular sample. Only genes constantly lying outside these bands in all possible sample combination were stored. In the next step, one sample t-test was applied to assess the p.value to every single gene outlier in a given sample according to the gene expression in the rest of the samples. All gene candidates with significant p.value that can be specified manually (default 0.05) were considered as differentially expressed. Graphical representation of PComp algorithm is depicted in Fig. 2 with artificial expression values. Obtained results were compared to the transcriptomic array results representing a gold standard for gene expression analysis.

RNA-Seq data were also subjected to a differential gene expression analysis using widely used and well-established algorithm limma (Ritchie et al., 2015) to gain an insight on overall performance of our PComp tool. Using limma we compared gene expression in a single sample (group 1) with all other samples (group 2) in all possible combinations and correlated the results to the PComp output.

Identification of gene fusions in cSVs cases from RNA-Seq

Based on a literature search (Liu et al., 2016), four state of the art tools for gene fusion detection—EricScript (Benelli et al., 2012), JAFFA (Davidson, Majewski & Oshlack, 2015), FusionCatcher (Nicorici et al., 2014), and TopHat-fusion (Kim & Salzberg, 2011)—were tested in our approach. TopHat-fusion was excluded due to a high rate of FP results and time-consuming computation in comparison to other methods. Thus, three tools—EricScript, JAFFA, and FusionCatcher—were run in parallel in our bioinformatic workflow for fusion gene identification; default settings were used for all of them. Although read alignment algorithms differ among the selected tools (Kent, 2002; Li & Durbin, 2009; Langmead & Salzberg, 2012; Dobin et al., 2013), the utilized computational algorithms follow similar concept (Kumar et al., 2016). In general, the mapping step consists of two alignments, one to the hg38 reference and the other to putative junction reference for each potential fusion. During these steps, unmapped and discordant pair-end reads that were mapped uniquely to different loci of the genome were identified, and a library of putative fusion junctions was derived. Reads were then realigned to the putative fusion junction sequences and annotated, a split-read signature (reads spanning fusion junctions) was recognized, and potential fusion genes were reported.

For efficient fusion gene detection, we developed an in-house meta-caller, which combines results from the selected tools into a consensus call. Each piece of software applies various metrics enabling classification of fusions to confidence subsets. In the first meta-caller step, detected fusions were filtered in individual callers as follows: for EricScript fusions with an Eric score (ES) of over 0.90, for JAFFA only fusions with the “HighConfidence” tag, and for FusionCatcher only high confidence fusions. In the consequent consensus call only fusions that passed at least two callers were kept which allowed increasing TP fusion rate and removal of FP fusions arising as specific artifacts of an individual caller algorithm (i.e., read alignments errors). The whole fusion gene detection pipeline is schematically visualized in the Fig. 3.

Validation of gene fusions by overlap with genomic arrays

Finally, to estimate TP and FP rates of the individual callers and the meta-caller, coordinates of fusions detected in RNA-Seq data were compared with the coordinates of genomic breakpoints detected by the genomic arrays. The interval of ± 100 kb around the breakpoints was applied in order to adjust to the array resolution and to the fact that the breakpoints can be located in introns, whereas fusions in RNA-Seq data appear in exon-exon boundaries. All potentially TP fusions identified by this approach were inspected visually in the array results.

Implementation of analytical algorithms

Algorithms and procedures of our novel solution were implemented in the R programming language. The meta-caller and the tool for differential expression are both freely available at GitHub (https://github.com/Hynst/WTS_cSVs_analysis), where the source codes are accessible for download. No installation is needed, as all dependencies (i.e., R packages and libraries) are installed automatically to the R environment.

PComp performance in differential gene expression analysis

We used PComp approach to identify differentially expressed genes on the level of individual samples. We followed a linear regression model and estimated a linear regression line with 99.9% confidence bands (1 − p = 0.999, i.e.,  p = 0.001) for each pair of samples. All genes lying outside the confidence bands were suggested as potentially differentially expressed candidates. We repeated this procedure for each sample pair combination and obtained nine sets of potential candidates per sample. Genes occurring in all nine sets were considered as differentially expressed for a given sample. By applying one sample t-test we assessed the differential expression of a given gene in a given sample according to the rest of the samples. We considered all gene candidates in particular sample with a p < 0.001 by one sample t-test to be differentially expressed. Using this approach, we found 15,200 differentially expressed genes in total (Table 1).

We inspected the results from RNA-Seq obtained using the PComp method and compared them to the results of transcriptomic arrays analyzed by one-way ANOVA. In the ten CLL samples tested, we identified 12,492 deregulated transcripts (including also unannotated ones) using transcriptomics arrays and extracted only annotated gene transcripts which resulted in 6,678 deregulated genes. We created an overlap between PComp and transcriptomic array results and found 1,808 significantly upregulated and 236 significantly downregulated gene candidates (Table 1). We also studied whether the results were concordant between the methods in terms of assigning genes as up-/downregulated and observed good concordance of the methods ranging 97.44–100% for individual samples (Table S1).

Similarly, in RNA-Seq data we identified 14,803 deregulated genes using limma. When we overlapped array data with limma results we found 1168 upregulated and 414 downregulated genes (Table S2). We compared the overlap rate of transcriptomic array with PComp and limma data, respectively, in terms of up-/downregulated genes. Both methods performed differently, we observed 14.7% and 28.9% upregulated genes, and 8.2% and 4.1% downregulated genes using PComp and limma, respectively. Finally, we observed larger overlap between transcriptomic array and PComp approach (Fig. 4), where 2044 deregulated genes were identified in total (compared to 1,585 identified by limma).

Meta-caller approach efficiency of fusion gene identification

In RNA-Seq data, we pre-filtered results of the individual pieces of software in the first meta-caller step to create a fusion gene subset with high confidence. No explicit pre-filtering of lowly expressed fusions was done in order to allow for the capture of lowly expressed fusion genes (i.e., the number of reads spanning fusion junction ≥1). In the next step of the meta-caller, results obtained by ≥2 callers were combined into a consensus call. Following these steps, our pipeline was capable of identifying 40 fusion genes which were potentially TP with high confidence across the ten samples (Table S3). Further, we used breakpoint coordinates generated using genomic arrays and compared them with those identified by RNA-Seq to test TP rates of our pipeline. Only fusions where at least one fusion gene partner was located in a genomic breakpoint were considered as potentially TP. Genomic array data corresponded to RNA-Seq breakpoints in 19 of the 40 cases. These uniquely expressed TP fusions were identified in eight of the ten CLL patients with 1 to 7 fusions per patient.

Considering candidates that were not supported by genomic arrays (21 of 40 features), potential fusions ZMYM5-PSPC1, HACL1-COLQ, MFSD7-ATP5I, and CYTIP-ERMN, occurring recurrently in six, three, three and two instances, respectively, were considered FP because they represented neighboring genes (read-through). Similarly, YWHAZ-ZNF706, LINCPINT-MKLN1, and NFATC3-PLA2G15 were also considered FP due to their adjacent localization in the genome, although not observed recurrently. Altogether, 17 of 21 fusions were considered FP. The remaining four fusions were not confirmed by the independent method; in all of them, at least one fusion partner appeared recurrently in the data, however, with a different partner. All of the four fusions could still represent TP resulting from balanced translocations or inversions that were not detectable by genomic arrays.

To assess the ability of the meta-caller to increase TP rate, we inspected overlap between individual callers and genomic arrays (Table 2). The meta-caller was superior to a single caller approach with the overall 47.5% TP rate, while JAFFA, FusionCatcher, and EricScript showed 9.4%, 8.2%, and 0.5% TP rate, respectively. Still individual callers identified TP fusions (as designated based on the overlap with genomic array results) which were missed by the meta-caller because they were supported only by one caller. Altogether 29 TP fusions were detected by all callers and meta-caller with the highest number of 22 TP fusions by FusionCatcher (Table S4). FusionCatcher and the meta-caller did not show any FP when overlapped with genomic arrays. Detailed results of this comparison are depicted in Table 3. Taking overall TP rates of the callers into account, the meta-caller appears as an efficient tool for fusion gene detection.