CLOVE: classification of genomic fusions into structural variation events

A structural variant (SV) is a rearrangement of the genome caused by at least two double strand DNA breaks followed by DNA repair. Typically, the term SV is used for events that are greater than 1 kb in size [1]. SVs include large insertions, inversions, balanced or unbalanced translocations, and amplifications and large deletions, collectively referred to as copy number variations (CNVs). A precise understanding of SVs is important in the study of population diversity, cancer [2–4] and other diseases (e.g. Charcot-Marie Tooth [5] and autism [6]).

The increasing usage of high throughput sequencing technologies has led to advances in the discovery and genotyping of structural variants in germline and somatic cells [7–9]. Consequently, a variety of methods have been developed to detect SVs from DNA sequencing data. Different approaches can be classified into four distinct categories: read depth (RD), discordant read pair (DR), split reads (SR), and de novo assembly (DN). RD methods involve counting reads in windows and segmenting the counts [10]. They identify only one class of structural variation (CNVs) and provide neither direct evidence for breakpoints, nor information about genomic organization. Their resolution and accuracy is dependent on sequencing coverage and window size, but is typically of the order of kilobases. Examples of RD methods include readDepth [10] and CNVnator [11]. DR methods use pair end-sequenced DNA fragments that span a breakpoint (typically in the un-sequenced region between the reads). These reads map anomalously or discordantly to the reference genome—further apart or closer together than expected based on the selected fragment size, to different chromosomes or with inverted orientation [12]. The signal of a rearrangement is a cluster of anomalous alignments. The resolution of DR methods is related to fragment size distribution and genome coverage. However, in general, single-nucleotide resolution is not possible with DR methods. BreakDancer [13] is an example of a DR method. SR methods rely on individual reads, which span the breakpoint and are capable of single-nucleotide resolution, although micro-homologies at breakpoint sites may introduce uncertainty. Examples of SR methods include Splitread [14] and Socrates [15]. DN methods use some form of assembly following other evidence to locate the locus of a rearrangement. DN methods typically provide single nucleotide resolution, but can be slow. Examples of DN methods are Cortex [16] and SOAPdenovo [17, 18].

Many tools utilize a hybrid approach that combines multiple lines of evidence to predict SVs. For example, Delly [19] and PRISM [20] use DR evidence and incorporate SR evidence through a targeted Smith-Waterman alignment. CNVer [12] uses DR and RD signals to identify potential copy number changes. CREST [21] uses SR then DN to directly map SVs at single nucleotide resolution, while SMUFIN [22] uses DN first and SR subsequently.

Another strategy is the “consensus” caller approach. For example, MetaSV [23] integrates a set of multiple tools into a pipeline. This approach aims to leverage the specific strengths of tools regarding different aspects of SV calling as well as confidence through agreement of multiple tools, and is therefore considered a meta-caller.

Nevertheless, the identification of SVs remains challenging. Existing methods suffer from a variety of issues relating to sensitivity and precision, positional accuracy and error profiles, and classification into one of the various types of SV. The majority of existing methods only identify breakpoints (a pair of connected breakends), also called genomic fusions (henceforth referred to as fusions), but do not classify the rearrangements further. Several SV callers are capable of limited classification of genomic events, such as insertions and deletions, but fail to classify more complex rearrangements. SVs may be simple (involving only a single fusion, such as an deletion) or more complex events (involving two or more fusions, such as a balanced translocation or inversion). Yang et al. [24] noted the lack of more complex events in the output of SV algorithms and introduced complex deletions, as well as inference of underlying DNA repair mechanisms. The study by Sudmant et al. [25] analysed a large cohort of human genomes for SVs including mobile element insertions. Another notable exception is the work by Escaramis et al. [26], which detects more complex events from aligned read data. However, the method is neither widely used nor cited, and failed to run on our data for purposes of comparison.

Here, we present a new method, CLOVE, which integrates calls from one or more breakpoint (or SV) detection methods and (re-)classifies the SV. Our method creates a graph data structure from the provided breakpoint information and then looks for patterns that are characteristic of more complex rearrangement types (e.g. balanced translocations). CLOVE is not another SV caller, but integrates multiple independent breakpoint predictions from other tools into a single, more accurate and potentially more complex event. This makes the output of these other tools more interpretable and increases the precision. A better-categorized output allows for better filtering or prioritization of specific events that are most relevant to the biological interpretation or experimental validation. CLOVE is the first meta-SV-caller that (i) can use any set of input (from current or future SV algorithms) to (ii) re-classify the data into more complex SVs than the original call sets.

Our algorithm is capable of handling the breakpoint or SV calls produced by a variety of tools – in fact, it can handle multiple sets of fusion at the same time, benefiting from the increase in sensitivity in raw calls. CLOVE can be thought of as augmenting existing breakpoint callers. It allows for stratification of SVs into (i) sets of simple events that collectively present the signature of a complex event, and (ii) remaining simple events that pass or don’t pass a read depth check. CLOVE presents the stratified SVs in VCF format with additional SV types, statistics about the read depth of events, and the levels of support for events (i.e. how many SV callers support the SV). In the following, we distinguish between the notion of basic and complex SV types. Simple or basic SV types are those that are fully represented by a single fusion. Complex SVs contain more than one fusion.

Basic SV types

We introduce the following terminologies and conventions in order to discuss the different types of SVs that can be identified in sequencing data:

1. (1)

   A fusion is a pair of loci that are adjacent in the donor genome but on different chromosomes or separated on the same chromosomes in the reference. The two separate locations of a fusion are referred to as breakends.

    
2. (2)

   Each breakend has an orientation. We define the orientation to be “+” if the breakend occurs on the 3′ end of the fused region (on its right) in the reference, and “−” if it occurs at the 5’ end (on its left). If a fusion has DR support, the discordant reads flanking the breakend will be aligned “pointing” towards the break and the read strand will coincide with breakend orientation. Similarly, for a fusion with SR support, if the aligner maps the paired reads to the same side (left) of the breakend with the 3′ end of the fragment clipped, then the breakend has a “+” orientation. Figure 1 illustrates an example of a breakend in an intra-chromosomal event and DR support.

    
3. (3)

   Based on the definitions above, we define a fusion as a pair of loci and directions: chromosome₁ (chr₁), position₁, orientation₁, chromosome₂ (chr₂), position₂, orientation₂.

    
4. (4)

   If a fusion is an intra-chromosomal event, i.e. the event occurs in the same chromosome (chr₁==chr₂), by convention, we assume position₁ ≤ position₂.

    

Based on our definitions, a single fusion can be classified into eight categories (Table 1). Furthermore, we have assigned an event type to each category. These definitions are in line with the literature [19]. The labels are motivated by the fact that some of the more basic structural rearrangements in DNA bear the corresponding signatures. Note that the sub-classifications of inter-chromosomal translocation types require a chromosome ordering and assume chr₁ < chr₂. This definition is arbitrary but necessary to distinguish the two non-inverted events. Most SV callers that do provide event classification often follow the naming scheme presented in Table 1, or something similar. However, while the breakend signatures for these types are indeed consistent with the events that give them their names, these are not the only events that can cause such a signature. Overlooking this fact is often a source of classification error and confusion when interpreting fusions.

Chromosomes|	Orientations
	+		−
chr1 = chr2	+	Inversion type 1 (INV1)	Deletion type (DEL)
	−	Tandem duplication (TAN)	Inversion type 2 (INV2)
chr1 ≠ chr2	+	Inter-chromosomal inversion type 1 (INVTX1)	Inter-chromosomal translocation type 1 (ITX1)
	−	Inter-chromosomal translocation type 2 (ITX2)	Inter-chromosomal translocation type 2 (INVTX2)

Rows refer to the orientation of the first breakend and columns to the orientation of the second breakend. Simple events may be combined into complex event types; for example, an inversion (an inverted segment of DNA) is comprised of the two simple events INV1 and INV2

Complex SV types

The computation and output of more complex SVs, such as translocations or inversions, is typically not addressed by existing SV calling algorithms (or only parts thereof as discussed below). Figure 2a and b illustrate the rearrangement patterns created by interspersed duplications and translocations on a single chromosome. Both patterns consist of two to three lower order events of the tandem duplication or the deletion type. This observation motivates our approach to further classify SV calls: the SV output referring to deletions and duplications at sites where a more complex event has occurred is often confusing to the user, and not easily identifiable from the list of breakends. Note that duplications and translocations upstream of the insertion site have a slightly different pattern with the deletion event pointing at the insertion site. Figure 2c shows the signature of an inversion event. Some SV algorithms, e.g. CREST, can classify this type of rearrangement. Figure 2d shows the signature of inverted interspersed duplications on the same chromosome. Similar to the classifications in Fig. 2a-d, parts e-g show the signatures for complex inter-chromosomal events. The list of rearrangements presented here is probably not complete as far as all potential complex events go, but we believe covers the majority of relevant classes. The signal of complex events shown above is consistent with existing discussion of structural variants (such as by [26]). To our knowledge CLOVE is the first meta caller that categorises such a comprehensive list of SVs from fusion calls.

The complex SV events described above share a common principle: there is at least one locus in the rearrangement pattern where two or more breakpoints (or simple SV events) have a common breakend coordinate. This observation guides the design of CLOVE to discover rearrangement patterns in SV data.

With one exception, the patterns for each of the rearrangement events are also unique and specific. Due to the symmetric nature of the intra-chromosomal translocation event, it cannot be determined whether the block has moved as shown in Fig. 2b, or indeed the block marked with an “A” has translocated to just after the second deletion event. Notice that both alternatives result in the same string of DNA. We define the convention that the smaller of the two alternative blocks is considered to have translocated.

Breakpoint graph construction and analysis

CLOVE has two major components, i.e. complex rearrangement pattern matching and read depth validation. Before these stages commence, the breakend calls are parsed into the algorithm. At the time of publication, the method supports outputs from Socrates, Delly, CREST, Gustaf [27], and BedPE (a standard format for fusion data). Furthermore, any SV caller can be re-classified, as long as it provides the sextet of information for each fusion introduced in Section 2.1. Any number of call sets can be input to CLOVE at the same time, which increases sensitivity if this adds non-redundant information.

The main internal data structure of CLOVE is a graph where nodes represent genomic coordinates (or coordinate intervals) and edges the fusions predicted by the SV callers. Fusions are added one at the time, and their coordinates are compared to the existing set of nodes. If a coordinate is within a user defined distance of a node, it is added to that node (modifying its coordinate interval if necessary). Otherwise, a new node (or pair of nodes) is created and added to the graph. Furthermore, one edge per fusion is added to the appropriate node(s). The edges (fusions) are labelled with one of the eight variant types introduced in Table 1. This data structure is similar, but not identical, to the breakpoint graph proposed by Bafna and Pevzner [28].

After the graph is constructed, CLOVE refines the existing edges in the graph. CLOVE scans the graph for redundant edges, i.e. edges between the same pair of nodes with identical SV type and merges them.

The refined graph is subsequently analysed to identify complex rearrangement patterns. Afterwards, the classified SV events are validated using read depth information from the original short read dataset. Figure 3 illustrates the workflow of CLOVE to classify provided fusion data into improved SV calls.

The two main steps are now explained in more detail.

1. (1)

   Complex rearrangement pattern matching: The initial classification step analyses the breakpoint graph in order to search for patterns described in Section 2.2. We iterate over the nodes in ascending coordinate order on each chromosome. If a node is adjacent with at least two edges, every edge-pair is investigated upon matching one of the complex SV rearrangement patterns (illustrated in Fig. 2). Intra-chromosomal duplication types are investigated first by matching tandem duplication and deletion types (or inversion1/2) accordingly. A further search for a deletion of the duplicated interval may then change the classification to a translocation. Any identified complex event is added to the graph and the (up to three) contributing events are removed from the graph. This step concludes once every node has been visited.

    
2. (2)

   Read depth validation: The second stage of the classification procedure again traverses the graph in order to analyse the breakpoint consistency with the read depth within the according intervals. Deletion events are expected to include a relatively low read depth in-between the two breakends, while tandem duplications are expected to include an increased coverage. The read depth is established for each individual interval by querying the input BAM file. This value is compared to the expected coverage and standard deviation. A single deviation value is used in this step independent of the size of the interval. Although this is not the most rigorous approach, it works well in practice. The expected deviation is a user-supplied parameter – to be chosen reflecting the coverage spread of the analysed read data. It is key for this step to take place after the classification step. This way deletion types and tandem duplication type events have already been merged into higher order events that do not change the read depth (for example, translocations). Events that do not fulfil the expected coverage response are rejected from the classification, and instead output as events of secondary interest due to read depth inconsistency.

    

After the classification stages, CLOVE presents the results in a new output file in VCF format. The read depth is presented along with the levels of support for each event (informing about the number of tools that support an event and the number of events within the call sets).

Clove is implemented in java and makes use of the htsjdk samtools implementation to handle genomic intervals.

To evaluate our method, we investigate results on simulated and real genomic data. The motivation behind the usage of simulated data is to allow full control over the sequence content, making an assessment of the sensitivity and precision of different tools possible. In the ideal case every simulated structural variant should be found and classified in the sequencing data, and no additional ones. Application to real data is important, but our knowledge about the variants present is incomplete.

We demonstrate that CLOVE is able to assign meaningful labels to rearrangement events. Furthermore, it increases the accuracy of the output of existing SV tools by removing a large proportion of false positive events. In the following two subsections, we compare the initial call set of different SV algorithms with the re-classified calls produced by CLOVE. For these comparisons, we investigate the sensitivity, precision, and accuracy statistics, which calls for the classic contingency tables of true positives (TP), false positives (FP), false negatives (FN) (and true negatives). Furthermore, we modify the standard approach as follows: we introduce a fifth contingency called half-true-positives (HTP). HTPs are defined as SV calls at the correct location but with the wrong event label (and potentially with too few fusions). This type allows us to make a fair comparison between the raw output (which often has few actual TPs) and the classified calls, without being overly harsh on existing methods – meaning, in the following HTPs are going to be counted as TPs for the calculations of recall, precision, and accuracy.

Simulated data results

For this experiment, we use chromosomes 21 and 22 of the human reference (hg19) as the underlying genome. The simulation workflow is as follows:

1. (1)

   The genetic material is divided into non-overlapping bins of length 500kbp. The binning strategy prevents events from overlapping and thus generating unresolvable complexities.

    
2. (2)

   A random event type is chosen for each bin. The types we use here are deletion, tandem duplication, inversion, intra- and inter-chromosomal duplications, intra- and inter-chromosomal translocations, intra- and inter-chromosomal inverted duplications, and intra- and inter-chromosomal inverted translocations. For each event a random starting or insertion position is chosen within the bin, as well as a random length or duplication/translocation interval, respectively.

    
3. (3)

   A new reference is generated that contains the above modifications.

    
4. (4)

   Reads are simulated from the alternative reference (using SimSeq [29]) at a sequencing depth of 30× and 100 bp read length.

    
5. (5)

   The reads are aligned to the original reference containing chromosomes 21 and 22 with Bowtie2 [30] in local alignment mode.

    
6. (6)

   SVs are called with a variety of tools from the aligned read data. All tools are run with standard parameters.

    
7. (7)

   CLOVE re-classifies the output of the SV methods generated in the previous step, and the results are compared to the list of variants created in Step 2). Clove is run without any parameters, except the coverage option, which is set to “-c 30 8” for this experiment.

    

Steps 1–3 of the simulations is done by using tools introduced in [31]. This simulation workflow is repeated five times, generating new random events and reads each time. The classified SV calls of CLOVE (v0.14) are based on the outputs of Socrates (v1.13.1), Delly (v0.76) and CREST (v0.0.1) for each simulated dataset. Tests are conducted using both the output of each individual tool as input and outputs combined from two or all three tools as input. CLOVE results are then compared to the raw results of the corresponding SV algorithm(s) that have been used as input. Additionally, we include information about runs of MetaSV (using version 0.5.4 and calls provided by Pindel, Lumpy, and BreakDancer) by itself. The performance is measured in terms of sensitivity, precision, and accuracy. Figure 4 summarises the results for the accuracy metric. Note that Fig. 4 is based on calculations treating HTPs as TPs for the purpose of calculating accuracy – specifically, accuracy ≔ (TP + HTP)/(TP + HTP + FP + FN). The detailed individual results and performance per SV type are listed in the supplement (Additional file 1: Table S1). These results include the performance of Lumpy – not shown in Fig. 4. Lumpy performs with an average raw accuracy of 0.71, which improves to 0.82 when applying CLOVE. Interestingly, Lumpy’s performance is better than MetaSV’s, which uses its data as input.

The results show the following trends. First, there is increased accuracy due to the classification of events. Although some tools already call complex variants, there is not a single tool capable of the same level of accuracy as clove when run on multiple inputs. Consolidating the trend with more specific results in the supplementary material, we can conclude that this effect is caused by the increase in precision of re-classified breakpoints. Conversely, the recall is decreased after classification. Obviously, the recall cannot increase, since CLOVE does not create any additional events in the SV input – and a decrease is inevitable if an event does not pass the read depth check, or a complex event is missing an edge leaving the remaining one(s) to be discarded, whereas it is counted as HTP in the raw input. A summary of the sensitivity for each SV tool and CLOVE can be found seen in Additional file 2: Figure S1. It shows that there is no clear strengths and weaknesses for individual types (except for MetaSV struggling with anything but inversions, deletions, and tandem duplications), and that CLOVE’s sensitivity is competitive, but not (greatly) superior to that of individual tools.

Do these results allow us to answer the question if the genotyping is improving a given set of SV calls? This depends on the application. It is obvious that the precision is greatly improved, which is often desirable, sometimes critical to an experiment. This improvement effect is in fact stronger than the loss of sensitivity caused by the genotyping step, leading to overall increased accuracy. To demonstrate this, we can compare accuracy values before and after classification and apply a pairwise t-test. The change in means is statistically significant for each of the analysed methods. Finally, CLOVE allows usage of separate call sets of SVs in a single run of the classification algorithm. The graph cleaning steps discussed above facilitate the removal of redundant edges (events that have been called by two or more of the input methods). Taking the union of all input events, CLOVE has (potentially) more data to classify, which reduces the drop-off in recall after classification. We can indeed observe that the maximum average accuracy is achieved by CLOVE when using input from Crest, Delly, and Socrates at the same time. This is not quite the case for MetaSV, which shows overall low performance values (for both precision and sensitivity) despite having access to the output of multiple tools.

The simulated data and predicted SVs are available from http://bioinf.wehi.edu.au/clove/.

CLOVE has been developed for and experimented with on germline data. In other types of data, such as cancer genomes, types of variants other than those currently identified by CLOVE may be present. For example, there are large scale chromosomal rearrangements, such as chromothripsis and breakage fusion bridge, common among some cancers. Further, compound events, where two or more of the events described in this work co-occur at the same locus, have also been observed in tumour genomes. It would be desirable to add rules and functionality to CLOVE to deal with such events. However, this is met with technical challenges of potentially incomplete fusion signatures (false negative calls), and an explosion of the rule set (for all potential compound events).

We have presented a new method for classifying complex rearrangement from breakpoint calls generated by different algorithms. We have demonstrated CLOVE’s ability to improve the output of standard SV methods by highlighting biologically relevant features, prioritizing, enhancing the precision of the calls, and generally improving accuracy. CLOVE’s independence of the input algorithm makes it a flexible tool to utilize in any SV calling pipeline. Its ability to process joint inputs from multiple methods is an attractive feature, which often leads to even better rearrangement classification, as has been indicated by our results. CLOVE is the first meta-caller that can use the input of any SV algorithm (provided it outputs sufficient information).