GASOLINE: detecting germline and somatic structural variants from long-reads data

Abstract

Long-read sequencing allows analyses of single nucleic-acid molecules and produces sequences in the order of tens to hundreds kilobases. Its application to whole-genome analyses allows identification of complex genomic structural-variants (SVs) with unprecedented resolution. SV identification, however, requires complex computational methods, based on either read-depth or intra- and inter-alignment signatures approaches, which are limited by size or type of SVs. Moreover, most currently available tools only detect germline variants, thus requiring separate computation of sample pairs for comparative analyses. To overcome these limits, we developed a novel tool (Germline And SOmatic structuraL varIants detectioN and gEnotyping; GASOLINE) that groups SV signatures using a sophisticated clustering procedure based on a modified reciprocal overlap criterion, and is designed to identify germline SVs, from single samples, and somatic SVs from paired test and control samples. GASOLINE is a collection of Perl, R and Fortran codes, it analyzes aligned data in BAM format and produces VCF files with statistically significant somatic SVs. Germline or somatic analysis of 30[Formula: see text] sequencing coverage experiments requires 4-5 h with 20 threads. GASOLINE outperformed currently available methods in the detection of both germline and somatic SVs in synthetic and real long-reads datasets. Notably, when applied on a pair of metastatic melanoma and matched-normal sample, GASOLINE identified five genuine somatic SVs that were missed using five different sequencing technologies and state-of-the art SV calling approaches. Thus, GASOLINE identifies germline and somatic SVs with unprecedented accuracy and resolution, outperforming currently available state-of-the-art WGS long-reads computational methods.

Figure 1

GASOLINE workflow. Panel (a) shows the steps to calculate NRO for a pair of SV signatures coordinates. Once the gap- and split-alignments coordinates ($S_i$ and $E_i$ for $i,j=[1,2]$) have been extracted from each read (a1), these are used to calculate the size of non-overlapping ($N{O}_{12}^1$ and $N{O}_{12}^2$, blue lines) and overlapping ($O_{12}$) segments for each pair of signatures (a2). $N{O}_{12}^1$, $N{O}_{12}^2$, $O_{12}$ and the total size of the two intervals ($L^1$ and $L^2$) are then used to calculate the $NR{O}_{12}$ coefficient (a3). In panel (b) is reported the workflow followed by GASOLINE for the detection of germline SVs in a sample. After signatures extraction, the tool calculates the $NR{O}_{\mathit{ij}}$ between all the signature pairs and generates an NRO matrix (b1) that is used as adjacency matrix to create an undirected graph by filtering out $NR{O}_{\mathit{ij}}$ values smaller than a predefined threshold (continuous edges represent $NR{O}_{\mathit{ij}}>NR{O}_{\mathit{thr}}$, while dotted edges $NR{O}_{\mathit{ij}}<NR{O}_{\mathit{thr}}$ (b2). The undirected graph is then analyzed with the Eppstein–Löffler–Strash algorithm to extract maximal cliques that represent clusters of SV signatures that can be assumed to be generated from the same SV event (b3). Next, all the SV signatures of a cluster are used to estimate the genomic coordinate (orange segment) of each SV event (b4). Finally, the number of SV signatures of a cluster and the total number of reads aligned in the breakpoints are used for genotyping with a maximum-likelihood Bayesian classification algorithm (b5). In panel (c) are reported the steps that GASOLINE follows for detecting somatic SVs. Somatic SVs are identified by comparing the SV signatures of a test (cancer) sample with a control (normal) sample. The SVs detected in the test sample (c1) are compared with the SVs signatures extracted from the control sample (c2) by calculating the NRO: SV signatures with a $NR{O}_{\mathit{SV}}$ larger than a predefined threshold are considered to be generated from the SV event of the test sample (c3). Statistical significance of each somatic SV is calculated by applying the Fisher’s exact test on the contingency table of (c4): $N_{\mathit{RRT}}$ (number of reads without SV signatures in test sample), $N_{\mathit{VRT}}$ (number of reads with the SV signatures in test sample), $N_{\mathit{RRC}}$ (number of reads without SV signatures in control sample), $N_{\mathit{VRC}}$ (number of reads with the SV signatures in control sample) . SVs with a p-value smaller than a predefined significance threshold are considered somatic.

Figure 2

Global performance of GASOLINE and the other three tools in the detection of synthetic and real germline SVs. Panels a-f report the F1 score obtained by the four tools in the analysis of simulated inversions (a,d), duplications (b,e) and translocations (c,f). Results are reported for ONT (a–c) and PacBio (d,e) synthetic reads aligned with minimap2. Panels (g–j) report precision and recall obtained by the four tools in the analysis of the NA24385 datasets for the detection of small SVs (g,h) and large SVs (i,j) with ONT (g,i) and PacBio (h,j) data. The curves in panels (g–j) were obtained by ordering all the SVs as a function of number of supporting reads and calculating precision and recall including SVs with decreasing number of reads. Panels (k–n) show F1 score obtained by the four tools for the NA24385 datasets in the detection of small (k–l) and large SVs (m–n) with ONT (k,m) and PacBio (l,n) datasets at different sequencing coverages.

Figure 3

Performance of GASOLINE on the detection of somatic SVs. Panels (a–d) show F1 score obtained by the four tools in the detection of simulated small (a,b) and large somatic SVs (c,d) with ONT (a,c) and PacBio (b,d) datasets at different sequencing coverage. Panels (e–o) report the precision-recall obtained by GASOLINE and the other three tools in the detection of deletions (e,j), insertions (f,l), duplications (g,m), inversions (h,n) and translocations (i,o) of the Valle-Inclan et al. true-set for the COLO829 cell lines sequenced with ONT (e–i) and PacBio (j–o) technologies. The results for GASOLINE were reported for different somatic p-value thresholds ($5\times {10}^{-1}$, $1\times {10}^{-1}$, $5\times {10}^{-2}$, $1\times {10}^{-2}$, $5\times {10}^{-3}$, $1\times {10}^{-3}$, $5\times {10}^{-4}$, $1\times {10}^{-4}$). All the results reported in the panels are based on mimimap2 alignment data.