SVInterpreter is a Python-CGI developed Web application, freely available on https://dgrctools-insa.min-saude.pt/cgi-bin/SVInterpreter.py. The code is accessible at https://github.com/DGRC-PT/SVInterpreter, and can be run locally with an Apache configuration.

TAD data from 10 tissue or cell types, available at YUE Lab website ¹ , were accessed for SVInterpreter. The regions bordering TADs–TAD boundaries—known to potentially restrict interactions of regulatory elements, were predicted using the Dixon pipeline (Dixon et al., 2012), whereas loops were established by Peakachu (Salameh et al., 2020).

For the chromosome Y, the TAD average size was calculated for each tissue or cell line, varying from 815 kb for lymphoblastoid cell line GM12878 to 1.8 Mb for bladder tissue (human genome assembly GRCh38/Hg38), and used as reference (Supplementary Table S1).

Full description of data sources used by SVInterpreter is available in Supplementary Table S2.

SVInterpreter analyzes any type of balanced and unbalanced SVs larger than 1 kb (translocations, inversions, insertions, deletions, and duplications) and retrieves a table of compiled information to assist their interpretation. Complex SVs must be subdivided in distinct SVs and analyzed separately (Supplementary Figure S1). Optionally, the user can apply SVInterpreter to any genomic region, without specifying the SV type.

SVs can be mapped within cell- or tissue-specific TADs, using the breakpoints as signpost. In this case, by default, TADs affected by breakpoints (brTADs) are retrieved, with the possibility of including up to five additional breakpoint flanking TADs (TAD−5 to TAD+5). Alternatively, instead of a TAD based analysis, SVs can be analyzed within a genomic region defined by its genomic position (Supplementary Figure S2).

To run SVInterpreter, a series of general parameters, such as genome version, tissue, or cell line to be used as reference for TAD and loop definition, and SVs or genomic region-specific parameters (Figure 1 and Supplementary Figure S2), are required.

From the selected specific genomic regions, SVInterpreter data were downloaded from public databases (last updated by March 31, 2021) or are automatically retrieved through an Application Programming Interface (API). From the breakpoints, all functional and non-functional genomic elements are retrieved, whereas from the remaining region, only protein-coding genes, lincRNAs, lncRNAs, functional, and non-functional genomic elements with a GTEx expression pattern are selected (Ardlie et al., 2015). Then, associated data are collected, including human disorders, cancer-specific rearrangements, phenotypes reported in animal models, genome-wide association studies (GWAS) data, and bibliography (Supplementary Table S2). The data are organized into a table, with indication of the breakpoint positions following the International System for Human Cytogenomic Nomenclature 2020 (McGowan-Jordan et al., 2020). In addition, to help visualization and interpretation of the SVs within the analyzed genomic regions, links to UCSC genome browser are made available on the output table. In this UCSC genome browser session, the selected genomic region is depicted, highlighting the SV (breakpoint or deleted/duplicated region). Native UCSC genome browser tracks compatible with the output table are shown, together with custom tracks, including the cell line/tissue-specific TADs and chromatin loops. Further description is available in Supplementary Methods.

For CNVs, SVInterpreter offers an option of performing an overlap search between the query CNV and those curated in several public CNV databases and published datasets (Supplementary Figure S2). The overlap specifications are similar to our previously published CNV-Content Tool (David et al., 2020), which retrieves the best hit by database, with the respective overlap percentage and variant frequency. In addition to the SVInterpreter standard output table, a detailed overlap table is available for download on the output web page.

Furthermore, to facilitate the evaluation of CNVs according to the ACMG guidelines, together with the standard output, the scoring parameters, as presented on the CNV pathogenicity calculator ² , are retrieved. SVInterpreter outputs the scores for the parameters that are possible to be established automatically and then performs an automatic calculation of the final score and their respective class (Riggs et al., 2020).

The output table(s) are written in XLSX format and made available for download. Further description of the output, a step-by-step tutorial, and an application example is available in Supplementary Methods and Supplementary Table S3.

Optionally, the proband’s HPO-based phenotypic features can be inputted for phenotype comparison (Köhler et al., 2019).

For this, the HPO ontology provided by the HPO.db package ³ and the links between genes, diseases, and terms provided by R data file (RDA) are a prerequisite. Since these were deprecated, we developed in-house scripts and used the June 2021 HPO release data ⁴ (Köhler et al., 2019) to create state-of-the-art HPO.db and RDA files. The scripts and guidelines are available at https://github.com/DGRC-PT/HPOSim_Helper_Scripts.

The phenotype similarity is evaluated based on phenotype similarity score (PhenSSc), maximum similarity score (MaxSSc), and p-value (p), which are calculated for each combination of inputted phenotype and Online Mendelian Inheritance in Man (OMIM) ⁵ phenotype associated with functional genomic elements within the analyzed region. This is performed by HPOSim—getTermListSim function that calculates pairwise similarities between HPO terms, using the information content (IC) of the most informative ancestor shared by both terms (Deng et al., 2015). The IC is a numeric value associated to each term, which inversely reflects the number of diseases annotated by the term, or any of its descendent terms. That is, terms with higher ICs annotate fewer diseases, being more specific, whereas lower ICs are associated to most common terms. When comparing groups of HPO terms, the getTermListSim result is the mean of the ICs of the pairwise comparisons, and reflects the similarity between the said groups, where higher scores represent higher similarities.

For PhenSSc, the inputted clinical features and the ones associated to a disorder are compared. This score reflects the similarity between the inputted phenotypic traits and the ones used to describe the disorder.

For MaxSSc, the inputted clinical features are compared with themselves, which means that MaxSSc consists of the mean of the ICs of the inputted terms. This metric was developed by us to reflect the maximum similarity score that can be obtained from the inputted terms, and to be used in comparison with PhenSSc.

The p-value, which reflects the probability of obtaining the PhenSSc by random chance, was adapted from Redin et al. (2017). In sum, for each disorder that PhenSSc and MaxSSc was previously calculated, a random set of HPO terms is selected. Most importantly, this set must have the same number of terms as the input, to limit the bias. The similarity score is then calculated between this set of terms and a disorder-associated phenotype (simulated score). Then, this is repeated 100 (n) times, where each time a different set of HPO terms is selected, and a new simulated score is obtained. Finally, the disorder specific p-value is calculated as: P = ∑ i = 1 n [ s i m u l a t e d s c o r e i ≥ P h e n S S c ] n (1)

Phenotypes mainly composed of terms common in a wide range of disorders, as global developmental delay, or intellectual disability (with 1,386 and 1,619 associated OMIM disorders⁴, respectively), can present high PhenSSc, close to MaxSSc, and a high p-value as well. In these cases, the high p-value reflects the high probability of the phenotype to overlap by chance, warning for the limited significance of PhenSSc. Hence, ideally, the PhenSSc should be close to MaxSSc and present a p-value as close to zero as possible.

For retrospective analysis, 97 and 125 previously published and unpublished balanced (translocations, inversions) and unbalanced (insertions, deletions and duplications) SVs were selected, respectively (Table 1; Supplementary Table S5) (David et al., 2003, 2009, 2015, 2018, 2020; Redin et al., 2017; Riggs et al., 2020). Of note, about half of those published by Redin et al. (2017) were previously analyzed by Zepeda-Mendoza et al. (2017), with the position_effect ⁶ tool, for identification of additional candidate genes.

For effectiveness evaluation in clinical setting, nine prenatal cases (three without associated ultrasound abnormalities, four with isolated increased nuchal translucency, one with limb abnormalities, and another with multisystemic traits) and 11 postnatal cases (with isolated organ-specific or complex multisystem disorders) were used (Table 1; Supplementary Table S5). These were randomly selected among those referral for clinical diagnosis, of which genomic variants were identified by CytoScan 750K (nine cases), CytoScan HD microarrays (six cases), and long-insert genomic sequencing (liGS) (five cases). Microarray and liGS analysis was carried out as previously described (David et al., 2018, 2020).

The microarray data were processed using Chromosome Analysis Suite 4.2.0.80 with NetAffx 20200828 (GRCh37/Hg19) and with the detection criteria of, at least, 15 probes within 35 kb for gains and losses. Selected SVs were manually interpreted based on the following criteria: absence/presence of OMIM genes, their association with autosomal dominant (AD) or recessive (AR) disorders, disruption of genes by the breakpoints, haploinsufficiency/triplosensitivity, and genotype/phenotype correlation (Silva et al., 2019). For this, data available at UCSC genome browser ⁷ , Decipher ⁸ , ClinGen ⁹ , ClinVar ¹⁰ , OMIM, DGV ¹¹ , Unique ¹² , and Orphanet ¹³ databases were used.

For liGS, SVs larger than 1 kb, and CNVs identified by discordant pair clustering and coverage analysis, were selected. Then, among these, novel variants and SVs overlapping a reported variant (Collins et al., 2017; Chaisson et al., 2019, Gnomad ¹⁴ ) with a database frequency <1%, and affecting loss-of-function (LoF) sensitive genes [with an expected vs. observed ratio (oe) of LoF variants of <0.35] and/or associated to AD disorders, were indicated for clinical evaluation. On average, per individual, 11 SVs were selected for analysis.

Three evaluators classified the Riggs dataset (Riggs et al., 2020) of unbalanced SVs. Therefore, based on the following criteria: (1) a classification equal by at least two of the evaluators, or (2) a median classification that reflects dissimilar evaluations, we merged them into a consensus classification (Supplementary Table S5).

To allow the comparison between published predicted outcomes of SVs and SVInterpreter-based prediction, criteria for translocations, inversions, and insertions were adapted from the previously described ones (Table 2) (Redin et al., 2017; David et al., 2020). For CNVs, ACMG guidelines were applied (Riggs et al., 2020). In addition, the same genome version and reference cell line as in the original publications was used. If available, the proband’s phenotype was inputted. For variants without pre-set of reference cell line, the human embryonic stem cell was used. By default, for all types of variants, the brTAD was used as reference, with rare exceptions. For CNVs, the overlap search against all available databases with a minimum mutual overlap of 70% was applied. The full set of variants and parameters used is available at Supplementary Table S5.

For retrospective analysis, 97 balanced and 125 unbalanced previously published SVs were reevaluated (Table 1; Supplementary Table S5). With the exception of chromosome 21, SVs are distributed regularly along the genome, with an average of 12 rearrangements per chromosome. Nevertheless, the larger number of translocations (n = 15), inversions (n = 5), and insertions (n = 3) involved chromosome 1, chromosomes 2 and X, and chromosome 3, respectively (Supplementary Table S4).

These variants were reevaluated by SVInterpreter, and based on its retrieved data, their clinical outcome was predicted according to the established parameters (Figure 2; Supplementary Table S5).

The first level of analysis involves functional and non-functional genomic elements localized within the brTADs and their annotation data, which is usually sufficient to evaluate a SV. For clinical outcome prediction of a gene disruption, SVInterpreter retrieves gene-specific annotation data such as the LoF sensitivity, Genomics England PanelApp ¹⁵ data, its association with disorders and respective phenotypic overlap, animal model data, gene expression patters, and GWAS data.

Concomitantly, the disruption of major genes by de novo BCA breakpoints leading to major AD developmental disorders, as retrieved by SVInterpreter, indicated the pathogenicity of ANKRD11 (OMIM *611192; proband DGRC0016) and WDR26 (OMIM *617424; proband DGRC0025) (David et al., 2020). In the abovementioned cases, the calculated similarity between the inputted phenotypes and of gene-associated disorders localized within the analyzed regions played a major part on the interpretation, where ANKRD11 PhenSSc was 2.64 (p = 0.02; MaxSSc = 4.01) and WDR26 PhenSSc was 2.31 (p = 0.02; MaxSSc = 2.91).

If the full extent of the clinical features cannot be explained by disruption or misregulation of a candidate gene, or the breakpoint is within an intergenic region, in search for potential position effect, annotation data of all genomic elements within a brTAD must be evaluated. Several data retrieved by SVInterpreter can suggest position effect events. In addition to the phenotypic overlap and expression pattern, disruption of chromatin loops and GeneHancer clusters of interactions are important signs for possible position effect (Gloss and Dinger, 2018).

DGAP131 t(1;5)(p31;q33)dn, was originally classified by Redin et al. (2017) as variant of uncertain significance (VUS). SVInterpreter showed MEF2C’s (OMIM *600662) GeneHancer cluster of interactions and 10 of its 14 chromatin loops disrupted by the chromosome 5 breakpoint. The PhenSSc of 1.82 (p = 0.02; MaxSSc = 3.1) corroborated the proposed position effect. Likewise, MEF2C was indicated as a potential candidate gene by Zepeda-Mendoza et al. (2017).

Then, if the protein coding genes or functional genomic elements localized within the brTADs are insufficient to explain the observed phenotype, additional upstream (–1 to –5) and downstream (+1 to +5) flanking TADs are analyzed.

Accordingly, the t(2;11)(q14.2;q14.2) breakpoints reported in proband DGRC0001 were located in intergenic regions, and no gene at the brTADs, capable of explaining the verified phenotype, was found. At TAD+1, SVInterpreter shows that the GeneHancer cluster of interactions of the proposed candidate gene GLI2 (OMIM * 165230) was disrupted by the 2q14.2 breakpoint, confirming the previously proposed position effect (David et al., 2009). Furthermore, the involvement of GLI2 was reinforced by its PhenSSc of 1.33 (p = 0.3; MaxSSc 4.2), with disorders OMIM#615849 and OMIM#610829. Ergo, the translocation was predicted to be likely pathogenic and confirmed the published assertion of the involvement of GLI2 (David et al., 2009).

Furthermore, DGAP107 t(Y;3)(p11.2;p12.3)dn, reported by Redin et al. (2017), presents, among others, neurological defects, urinary tract, and genital abnormalities. They originally classified the SV as potentially pathogenic, due to the disruption of ROBO2 (OMIM *602431). By assessing the associated disorder (OMIM #610878), we realized that ROB O 2 only explained the urinary tract defects (PhenSSc = 1.12; p = 0.08; MaxSSc 1.48). However, SVInterpreter brTAD analysis suggested a position effect on PCDH11Y (OMIM * 400022), which had its GeneHancer cluster of interactions disrupted. The gene function, expression pattern, and animal model data suggest its role in the development of the nervous system, and therefore may explain the neurological defects observed in the proband. Besides, Zepeda-Mendoza et al. (2017) also indicate SRY (OMIM *480000), located at TAD-3, as a candidate gene due to the overlap with the genital abnormalities.

The overlap search of query CNVs in public database data and the automatic ACMG scoring showed to be of utmost utility, since it can clarify immediately the potential significance of deletions and duplications, even in cases where the genomic data are scarce. As such, a 374-kb deletion, arr[GRCh37]10q22.3(81,603,169_81,976,925)x1, in case CK without associated phenotype, was classified by Riggs et al. (2020) as VUS. According to SVInterpreter, the CNV deleted five genes that were not associated to phenotype or reported to be haploinsufficient. The CNV had 100% overlap with a likely benign ClinGen deletion (nsv3896137), and according to its ACMG CNV score of −0.9, the deletion was classified as likely benign.

Overall, more than half (57.2%) of the reevaluated SVs (45 translocations, 8 inversions, 32 deletions, and 42 duplications) were originally classified as VUS, whereas only 10.4% (23) were classified of benign and likely benign (Figures 2A–C). SVInterpreter-based reevaluation of published SVs provided a consistent finding with the original studies on 62.6% of all SVs (39 translocations, 9 inversions, 44 deletions, 45 duplications, and 2 insertions) (Supplementary Table S5). Comparatively with the original classification, the number of variants predicted as VUS decreased by 40% (from 127 to 76) (Figure 2A). For balanced SVs, SVInterpreter-based interpretation led to the reevaluation of 81.1% of the original VUS as pathogenic (9.4%), likely pathogenic (32.1%), and benign (39.6%) (Figures 2B,D,E). In addition, position effect events identified by SVInterpreter sustained the categorization of 30.2% of the potentially pathogenic balanced SVs (Supplementary Table S5). For deletions, the differences between published and tool-based prediction were minor, with similar results obtained by both (Figure 2F). Differently, 19% of the duplications categorized the VUS transited to another category, whereas only three insertions were reclassified from benign to likely pathogenic (Figures 2G,H).

To assess the position effect on distal genes and their contribution on the observed phenotypes, from the 87 balanced SVs published by Redin et al. (2017) and reevaluated by us, Zepeda-Mendoza et al. (2017) also analyzed 44 (Figure 3A). Similar candidate genes were identified in 11 of the SVs (Figure 3B), whereas in 5, neither of them proposed a candidate gene (Figure 3A). SVInterpreter and position_effect identified the same candidate genes for two originally classified VUS and two pathogenic SVs (Figure 3B; Supplementary Table S5). The position_effect tool uniquely identified 24 candidate genes in 19 SVs, where, in 14 of them, the genes were located outside the brTAD (Figure 3C; Supplementary Table S5). Based on expression, phenotypic overlap, and animal model data, SVInterpreter predicted six candidate genes not foreseen by the other two approaches, in five SVs (Figure 3C; Supplementary Table S5).

The effectiveness of this bioinformatic tool in a clinical setting was evaluated by comparative (manual vs. SVInterpreter-based), clinical outcome prediction of SVs identified by Cytoscan 750K microarray in nine prenatal cases, by Cytoscan HD microarray in six postnatal cases, and by liGS in five postnatal cases. Altogether, 145 variants (SVs, CNVs, and cxSVs) were analyzed (Table 1). The average number of SVs per individual, identified by genomic array, was 6, whereas for liGS, it was 289, with 44 balanced and 244 unbalanced variants. From the latter, on average, only 11 SVs (3 balanced and 8 unbalanced) were recognized to be potentially disease causing or pathogenic, and consequently selected for clinical outcome prediction (Supplementary Table S5).

Proband DGRC0004 presented a severe phenotype characterized by global developmental delay, facial dysmorphisms, and heart defects. Among other, the liGS data analysis identified a 67.3-Mb inversion inv(2)(p16.1q14.3) and a 589-kb duplication dup(2)(q21.1). Since, based on the SVInterpreter data, none of the inversion breakpoints disrupted a gene, nor any gene localized within the brTAD supported the verified phenotype, the inversion was classified as benign. Concerning the dup(2)(q21.1), although SVInterpreter identified an identical CNV in a cohort of patients with developmental delay (nsv999864) (Coe et al., 2014), the duplication has the same gene content as reported in benign SVs and did not affect triplosensitive genes, which led to its likely benign classification (ACMG CNV score = −0.9). Furthermore, none of the remaining eight clinically evaluated SVs was predicted to be likely pathogenic or pathogenic; therefore, genomic disorder was excluded in this case. Indeed, exome sequencing identified a pathogenic single-nucleotide variant within KAT6A (OMIM *616268) exon 18, causing AD Arboleda-Tham syndrome (OMIM #616268) (data not shown). The clinical features of this syndrome overlap that of the proband.

As CMA is the technique of choice for identification of CNVs in a clinical setting, the automatic mutual overlap search with CNV public databases and the inclusion of the ACMG scoring system is especially valuable for faster and more informed clinical outcome prediction of these.

A female in her 40s presented a dichorionic diamniotic pregnancy with an elevated risk for aneuploidy following first trimester combined screening test and normal ultrasound examination. Microarray analysis of chorionic villus sample DNA (CS750K07) identified five deletions and two duplications. By manual analysis, due to the absence of genes within the five deleted regions, these were classified as benign, whereas one of the two duplications, encompassing only a non-morbid gene, was classified as likely benign. SVInterpreter confirmed the benign and likely benign classifications, and the absence of overlapping CNVs and triplosensitive genes. In contrast, the remaining 1.1 Mb duplication at 16p13.11, arr[GRCh37]16p13.11(15,416,498_16,527,659)x3, was classified as VUS, since the CNV was overlapped by the 16p13.11 microduplication syndrome, which likely presents an incomplete penetrance and phenotypic variability. SVInterpreter identified four overlapping disorder-associated genes, NDE1 (OMIM *609949), MYH11 (OMIM *160745), ABCC1 (OMIM *158343), and ABCC6 (OMIM *603234). Although these genes are associated to AD or AR disorders, neither of them is triplosensitive or is disrupted by the breakpoints. SVInterpreter identified overlapping duplications that were reportedly classified as pathogenic (nssv15605791), likely pathogenic (nssv15149610), likely benign (nssv15159627), and VUS (nssv15159626). In addition, automatic bibliography search identified publications that described the 16p13.11 microduplication syndrome (PMID: 30287593, PMID: 23637818). Hence, in the absence of prenatal phenotype–genotype correlation, the contradictory classifications of similar duplications, and the overlap with the microduplication syndrome, we maintained the original classification of VUS.

We confirmed the reported manual clinical prediction of SVs identified in 20 individuals analyzed in a clinical setting, with marginal variability between these two approaches (Supplementary Table S5).