Comprehensive structural variation genome map of individuals carrying complex chromosomal rearrangements

Complex chromosomal rearrangements (CCRs) are rearrangements involving more than two chromosomes or more than two breakpoints. Whole genome sequencing (WGS) allows for outstanding high resolution characterization on the nucleotide level in unique sequences of such rearrangements, but problems remain for mapping breakpoints in repetitive regions of the genome, which are known to be prone to rearrangements. Hence, multiple complementary WGS experiments are sometimes needed to solve the structures of CCRs. We have studied three individuals with CCRs: Case 1 and Case 2 presented with de novo karyotypically balanced, complex interchromosomal rearrangements (46,XX,t(2;8;15)(q35;q24.1;q22) and 46,XY,t(1;10;5)(q32;p12;q31)), and Case 3 presented with a de novo, extremely complex intrachromosomal rearrangement on chromosome 1. Molecular cytogenetic investigation revealed cryptic deletions in the breakpoints of chromosome 2 and 8 in Case 1, and on chromosome 10 in Case 2, explaining their clinical symptoms. In Case 3, 26 breakpoints were identified using WGS, disrupting five known disease genes. All rearrangements were subsequently analyzed using optical maps, linked-read WGS, and short-read WGS. In conclusion, we present a case series of three unique de novo CCRs where we by combining the results from the different technologies fully solved the structure of each rearrangement. The power in combining short-read WGS with long-molecule sequencing or optical mapping in these unique de novo CCRs in a clinical setting is demonstrated.


Introduction
Complex chromosomal rearrangements (CCRs) are rearrangements involving more than two chromosomes or more than two breakpoints [1,2].Balanced, de novo, CCRs are extremely rare and more than half are associated with an affected phenotype [3].Traditionally in a diagnostic laboratory, karyotypically detected complex rearrangements are further analyzed using methods such as array comparative genomic hybridization (aCGH) and fluorescence in situ hybridization (FISH).Which method to choose and the number of experiments needed to resolve its genomic structure will be dependent on the characteristics of the rearrangements under scrutiny.These characteristics include the number of potential breakpoints, as well as if the rearrangement includes unbalanced events such as cryptic deletions, which are commonly the cause of the affected phenotype [4,5].Additionally, factors such as cost, patient phenotype, time, availability of tests and interests of the lab (diagnostic versus research) will all be considered in deciding which methodology to choose.
More recently, whole genome sequencing (WGS) has been used to study complex rearrangements.Compared to the standard cytogenetic techniques, WGS offers a wide range of advantages, including high resolution and the ability to detect rearrangements of any size and type in a single experiment, providing important clues for the mechanism of formation [6,7].Furthermore, it may be critical for diagnosis of genetic diseases as well as for genotype-phenotype studies.However, solving the structure of a complex rearrangement is costly and time consuming, mainly due to the sheer amount of breakpoints involved, the large variety in sizes and types of structural variants involved, challenges in interpretation, as well as specific limitation detection of the WGS method of choice.
Today, Illumina sequencing technology is the most commonly used second-generation WGS technology [8].Illumina sequencing is cost efficient and provides high quality data with a wide range of library preparation methods available.Each of these library preparation methods produces data at different cost and with different properties.In particular, the long insertsize mate-pair (MP) sequencing protocols have been useful for analyzing structural variation, mainly due to its high span coverage which potentially could bridge regions that are hard to map, such as repetitive regions [7,9].On the other hand, the higher coverage PCR-free pairedend (PE) libraries enables the use of split reads for exact breakpoint analysis, and highly precise read depth CNV detection [10,11].These different protocols provide a flexible framework for analyzing different types of variations and answering different questions.However, it is not obvious which method to choose when faced with the problem of solving a complex rearrangement.Even though WGS comes with a range of advantages it has important limitations, including low sensitivity in repetitive regions.Furthermore, haplotyping of long genomic segments is still a challenge for WGS even using single-molecule sequencing technologies.Therefore, phasing multiple rearrangement breakpoints that may have been created concomitantly in complex rearrangements including CCRs and chromothripsis events, can be a daunting task.Moreover, the breakpoints found using WGS commonly needs to be validated using orthogonal techniques, most commonly breakpoint PCR and Sanger sequencing.Hence, even when using WGS, the analysis of complex rearrangements is difficult, and a large number of experiments may be needed.
Bionano optical mapping is a technology that enables detection of large structural variants across the entire genome, potentially helpful to resolve long haplotypes including detection of genomic variants in cis.In contrast to the short-read WGS methods commonly used today, Bionano optical mapping utilize long DNA molecules (>100 kb).The usage of long input DNA molecules enables the optical maps to span repetitive and poorly mapped regions of the genome.Further, each optical map may span multiple adjacent breakpoints, providing additional information on how the breakpoints relate to each other.Bionano optical maps show great promise for SV detection and phasing [12], although optical maps are currently limited by lower resolution compared to sequencing technologies and the availability of the technology itself.Linked-read sequencing is a method provided by 10X Genomics, where linked reads are used to detect structural variation and allows for detection of SVs located within repetitive regions using barcoding of long DNA molecules [13,14] but it is also limited by availability of a comprehensive software and high false-positive rate.
In this study, we have characterized three unique CCRs combining multiple technologies including standard cytogenetic techniques (aCGH and FISH), followed by short-and longread approaches in addition to Bionano optical mapping.For short-reads, we used two different sequencing protocols, Illumina 30X PCR-free PE WGS and Illumina Nextera MP WGS, whereas linked-read sequencing using 10X Genomics Chromium was used to obtain synthetic long-reads.By combining next generation sequencing and cytogenetic techniques, we were able to fully characterize the molecular structures of each unique CCR.Lastly, we performed a comprehensive comparison of each methodology applied here concerning the resolution and sensitivity to detect SVs of different sizes throughout the genome, which we discuss in detail.

Results
Results for the cytogenetic studies are presented in S1 Appendix.BAM files containing all supporting reads for the three rearrangements are available at European Nucleotide Archive (ENA), project number PRJEB26322, and Bionano xmap files are available in the S1 Dataset.

Case 1
The t(2;8;15) complex rearrangement was first indicated by regular karyotyping, however, this analysis only showed an abnormal derivative chromosome 15 and the karyotype was reported out as 46,XX,del(15)(q?22)(S1 Fig) .Further delineating of the chromosomes using FISH and aCGH identified the complex rearrangement involving three chromosomes.Deletions in the breakpoints of chromosomes 2 and 8 were detected using FISH and genome-wide aCGH.Of note, although the 14.5 Mb deletion at 8q23 is within the resolution of classical chromosome analysis, it was not detected, likely masked by the presence of multiple segments involved in the complex rearrangement (S1 Fig) .WGS confirmed the aforementioned chromosomal rearrangements and determined the chromosome 2 deletion sizes to be 2.1 Mb and 2.3 Mb, respectively.Additionally, a genomic segment of approximately 970 kb, originally located between the two deleted parts from chromosome 2, was inserted onto chromosome 15 (fragment D, Fig 1A) together with a small inverted 12 kb fragment also originating from chromosome 2 (fragment C, Fig 1A).The first deletion of 2.1 Mb removed eight protein-coding genes and the second deletion of 2.3 Mb deletion removed three protein-coding genes (Table 1).The deletion on chromosome 8 was determined to be 14.5 Mb and removed 49 protein-coding genes (Table 1).The translocation breakpoint on chromosome 15  A total of eight breakpoints and five breakpoint junctions were identified in Case 1.Three out of five of the breakpoint junctions were detected using all four WGS technologies.The 12 kb inversion on chromosome 2 (fragment C, Fig 1A) was detected using both short-read sequencing (PCR-free PE and MP) technologies and the linked-reads technology with the Supernova pipeline, but not by optical mapping likely due to its small size.One of the translocation junctions between chromosome 2 and chromosome 8 was not detected using the linked-reads technology (Supernova pipeline).None of the breakpoint junctions in Case 1 were found with the linked-reads technology using the Long Ranger pipeline.
Finally, eight breakpoints were located within repeats elements.Sequence homology was found in two out of these junctions, however only short stretches of microhomology indicate that no fusion repeat elements were formed (

Case 2
The t(1;5;10) complex rearrangement was first identified using regular karyotyping and FISH.A large deletion of 14.4 Mb located at chromosome 10p12 was not detected by the initial chromosome analysis, possibly due to the fact that multiple chromosomal segments were involved and the deletion was located in the translocation breakpoint region.Instead, the deletion was identified during characterization of the rearrangement using FISH and later confirmed by aCGH.The molecular cytogenetics data was reported previously [15].WGS confirmed the cytogenetic findings (Fig 1B ), including the deletion (segment F, Fig 1B) that involved 75 protein-coding genes (Table 1).A total of five breakpoints and four breakpoint junctions were identified.Breakpoint junctions revealed 0 to 3 nt of microhomology, and no SNVs or insertions neither in the junctions or flanking the junctions (Table 2, S2 Fig) .Three out of five breakpoints of this complex translocation were located within repeat elements, but no evidence for fusion repeats elements were observed in the breakpoint junctions, meaning that the rearrangement produced a truncated repeat (S2 Fig) .The breakpoints of the large deletion on chromosome 10 were both located far from any repeat element (Table 2).Using the WGS data, the molecular karyotype could be determined to be t(1;10;5)(q31.3;p12.31;q23.2) seq[GRCh37] g.[chr1:pter_cen_196997343::chr5:124956736_qter] g.[chr5:pter_cen_124956731::chr10: 20816168_pter] g.[chr10:qter_cen_20816166::4689760_19120882del::chr1:196997343_qter].No other rare structural variants were identified in the microarray or WGS data.All breakpoint junctions indicated by the cytogenetic analysis were found by all four WGS technologies applied, although only PCR-free PE WGS and linked-read sequencing resolved the junctions on the nucleotide level.

Case 3
First, karyotyping identified an unusual banding pattern on the short arm of chromosome 1.Follow-up analysis using FISH and BAC array identified 14 breakpoint junctions including a 0.87 Mb deletion located at 1p36.2, which was previously reported in Lindstrand et al. (2008) [16].Combined analysis of the four different WGS technologies confirmed those 14 breakpoint junctions and unveiled 12 additional ones (Table 3) originating from chromosomal segments translocated, inverted, and deleted involving both the short and the long arm of chromosome 1 (Fig 2).Hence, a total of 33 breakpoints and 26 breakpoint junctions were   identified after WGS analysis.The WGS analysis identified seven deletions at the rearrangement breakpoints (Table 1).In the breakpoint junctions, microhomology was found in 10 junctions (2-6 nt), and one junction contains a 9 bases long non-templated insertion (Table 3, S2 Fig) .Every breakpoint junction contained at least one repeat region, mostly Alu elements (Table 3).The junctions that do not involve Alu elements are located within a wide range of distinct types of repeat regions.The second major group of repeats were found to be LINE elements (11 breakpoints), remaining breakpoints seem to be scattered randomly amongst various repeats including simple repeats, LTR elements, and retroviral elements (Table 3).In total, 26 breakpoint junctions were detected in Case 3 using combined analysis.Two junctions, junction 1 and junction 2, were not detected using any short-read (PCR-free PE and MP) method, likely because they map within low-copy repeats (LCRs) in a 125 kb large intron of TTC34.Breakpoints located within such regions may require longer reads to be correctly mapped, and indeed only optical mapping and linked-read sequencing were able to detect those junctions (Fig 3 , Table 4).Of note, all 26 breakpoint junctions were detected using the linked-read WGS technology, although we were required to use both the Long Ranger and Supernova pipelines as they complement each other (Table 4).PCR-free PE WGS detected 24   4).

General method comparison for solving complex chromosomal rearrangements
Comparing the results of the different technologies, it was found that all technologies tested here could detect the majority of the junctions (Tables 4 and 5).PCR-free PE WGS and the linked-read technology present the highest detection rate: PE detected all but two junctions whereas linked-reads detected all but one junction.Moreover, both PCR-free PE WGS and linked-read WGS share the highest resolution (1 bp).MP WGS present a resolution of ~400bp with a false-negative rate of 9% (fails to report three junctions out of the total 35 junctions) (Table 5).Lastly, optical mapping has the lowest detection rate since it fails to report seven junctions of the total 35 junctions (20%) and the lowest resolution of ~7 kb (Table 5).

Delineating and identifying derivative chromosomes from WGS data
The three CCRs presented here were previously identified using chromosome analysis, therefore the WGS analysis was focused on delineating the structure of the derivative chromosomes.However, to simulate a "WGS first" scenario and evaluate the utility of each of the techniques applied here for SV detection without a prior hypothesis, we ran FindSV developed for PCRfree PE and MP WGS data.The FindSV analysis pipeline is described in the WGS methods section.Generated calls were thereafter ranked based on i) frequency in the SweGen SV database [17] ii) amount of discordant read pairs and iii) size in base pairs, or chromosomal position.Based on these filtering criteria, calls pinpointing the breakpoints of the complex rearrangements were ranked high in the generated lists of rare SVs in all cases (S1 Table ).The FindSV pipeline generated 162 calls from the PCR-free PE data in Case 3 from which 34 wellsupported interchromosomal calls were ranked above the intrachromosomal rearrangement on chromosome 1 (S1 Table ).These calls are most likely due to mobile elements [18], therefore we concluded that filtering of PCR-free PE data needs to be optimized to minimize false-  positive interchromosomal calls due to repetitive elements and favor ranking of potentially false-negative intrachromosomal rearrangements.Lastly, not all breakpoint junctions will be detected in the FindSV output data.FindSV may fail to detect breakpoint junctions that are located in repetitive regions as well as regions that are poorly covered due to various technical artifacts or low DNA quality.All highly ranked calls are manually inspected in IGV, which allows for detection of additional breakpoints and small aberrations in the junctions.Analysis of the linked-read and optical map pipelines in a simulated "WGS-first" scenario was more challenging because both technologies produce a large number of calls (>1000).In addition, there is still a lack of frequency databases for these technologies, which makes filtering and ranking based on public database frequencies currently not possible.We concluded that using optical maps or linked-read sequencing technologies as the initial screening techniques for CCRs might not be feasible with current pipelines until more frequency data is available.Furthermore, optical mapping has a lower resolution limitation which hamper its use to detect smaller and more complex junctions, including the small 12 kb inversion of fragment C in Case 1 (Fig 3A -3C).Importantly, however, we found that optical mapping and linked-read WGS performs better than short-read WGS in highly repetitive regions or paralogous regions (such as Junctions 1 and 2, Case 3, located within LCRs).Example of the junctions that could not be detected using any of the short-read WGS protocols is shown in Fig 3D and Fig 3E.

Sensitivity of WGS SV calling compared to CNVs detected by aCGH
To compare the utility of all WGS technologies used here, we also performed a comparison of all polymorphic CNVs present in the three probands of the present study detected with aCGH (S2 Table ).In brief, CNV calls obtained by WGS were compared to those from two different aCGH platforms: a custom designed 400K genome-wide array with ~2000 targeted high-resolution genes and a commercial 1M medical exon array from OGT with ~5000 targeted highresolution genes.The array calls that were also found with at least one of the WGS methods were considered confirmed and hence more likely to be true variants.The sensitivity for detecting such confirmed CNVs in each case and method is given in S2 Table and the summarized sensitivity per method in Table 6.The data suggest that PCR-free PE WGS has a high detection rate on both small (<0.1 kb) and large (>10 kb) SVs and offers the highest overall detection rate.Linked-read WGS technology has the second highest detection rate, mainly because of the Supernova pipeline and the Long Ranger small deletions algorithm.Lastly, MP WGS and optical mapping perform similarly: both technologies perform well on the larger CNVs, but fail to detect a large number of smaller variants (S2 Table ).Overall, it is clear that the high resolution of linked-reads WGS and PCR-free PE WGS allows for the detection of small SVs that MP WGS and optical mapping does not have the required resolution to detect.

Comparison of genome wide SV calls from different WGS technologies
Next, we used the Case 2 WGS data, which were of the best quality among the three cases, to compare the sensitivity, number of calls, reported variant sizes, and reported variant types of the WGS technologies to one-another.Using SVDB (https://github.com/J35P312/SVDB)[10], we found that the technologies differ regarding the overlap of the variant calls (Table 7).In particular, optical mapping and the Illumina based technologies detect only a few hundred (6.6%) overlapping calls (Table 7).In contrast, MP and PCR-free PE WGS produce the largest fraction of overlapping calls with nearly 71% of the MP calls also detected by PCR-free PE WGS.The overlap between PCR-free PE WGS and linked-reads WGS was 42% (Table 7).
There are many reasons why the number of overlapping calls may differ between technologies, including artifacts and features specific to each technology and calling pipeline that can influence the output (S2 Table ).Bias towards variant types is also observed, for instance optical maps report no duplications and 6733 insertions while PE WGS reports 454 duplications but no insertions (Table 8).Moreover, due to distinct resolutions, each methodology also present a Next, we compared the overlap between the variant calls and three public datasets: the Database of Genomic Variants (DGV) [19], the HG002 integrated call-set compiled by the Genome In A Bottle consortium (GIAB) (https://github.com/genome-in-a-bottle)[20], as well as the CNV list published in a previous paper by Conrad et al. (2010) [21] (Table 9).Notably, the performance between the technologies differs depending on which database they are compared to.Comparing the methodologies used here to the CNVs reported in Conrad et al. [21], we found that Bionano and PCR-free PE WGS achieve similar numbers.MP WGS reports fewer variants matching the CNVs listed in Conrad et al, and linked-read WGS detect only a few of those CNVs (Table 9).
In contrast, linked-read WGS performs better on the DGV database, and detects the largest number of deletions (Table 9).Despite the high detection rate of these deletions, linked-read WGS performs relatively poor on all other variant types.Similar to the comparison of the Conrad et al. dataset [21], optical maps and PCR-free PE WGS report nearly the same amount of variants.However, optical maps detect a greater number of insertions, while PCR-free PE WGS detects a greater number of deletions.
MP WGS detects the smallest number of variants, however, most of the MP variant calls do match a DGV variant indicating that the precision of MP WGS is high.In contrast, a smaller fraction of optical maps and Long Ranger calls are similar to DGV variants, indicating a lower precision of these technologies.
Lastly, the technologies were compared to the HG002 (GIAB) integrated call-set, which consists of mainly small SVs (<0.1 kb) (S4 Table) detected using both long-and short-read sequencing methods (ftp://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/AshkenazimTrio/analysis/NIST_SVs_Integration_v0.6/README_SV_v0.6.txt).It was found that linked-read WGS and PCR-free PE WGS detect the largest number of those variants, while optical maps and MP WGS produce relatively few hits.These results are not surprising given the resolution and size bias that those methodologies present (Table 5).

Discussion
It has previously been demonstrated that cytogenetically balanced CCRs often contain cryptic deletions in the breakpoints, explaining the clinical phenotype of the carrier patients in many  cases [5].Here we report two individuals with complex interchromosomal translocations (Case 1 and Case 2) and one individual with a very complex intrachromosomal rearrangement (Case 3) for whom we used multiple combined cytogenetic and molecular methodologies to refine the alterations in their genome content which aided clinical assessment.Del: deletion, dup: duplication, inv: inversion, ins: insertion, SV: structural variant, MP: mate-pair, PE: paired-end.
The Total row indicates the total number of variants of each type in the three datasets: DGV [19], the HG002 integrated call-set [22] and Conrad et al. (2010) [21].BssSI and BspQI are the two restriction enzymes used in the Bionano optical mapping experiment, Long Ranger: the 10X Genomics mapping-assembly based pipeline, Supernova: the results of a custom pipeline utilizing the Supernova de novo assembler. https://doi.org/10.1371/journal.pgen.1007858.t009 In Case 1, a 5 Mb deletion on chromosome 2 and a 14.5 Mb deletion on chromosome 8 was identified using FISH mapping.Further analysis using aCGH and WGS, showed that the 5 Mb deletion actually consisted of two deletions, 2.1 Mb and 2.3 Mb, respectively.WGS could also demonstrate that a 970 kb genomic segment, originally located in-between the two deletions, had been translocated onto chromosome 15.The deletion on chromosome 8 in Case 1 covers the TRPSII (Trichorhinophalangeal syndrome type II, TRPSII, Langer Gideon Syndrome, LGS, MIM# 150230) locus and most of the characteristic symptoms of TRPSII were present in Case 1.
In Case 2, the translocated piece of chromosome 10 contained a 14.5 Mb deletion, encompassing the typical region for hypoparathyroidism, sensorineural deafness, and renal disease syndrome (HDRS, Barakat syndrome, MIM# 146255), characterized by the triad hypoparathyroidism, renal dysplasia and hearing loss.The common cause of HDR syndrome is mutations in GATA3 [23].Deletions of 10p are recurrent, and GATA3 has been pinpointed as the causative gene for HDR syndrome seen in 10p deletions.The size and location of 10p deletions vary, as well as the clinical picture [15].
Case 3 presented with mild dysmorphic features, psychomotor delay, ectopic left kidney, minor hearing disability, dysphasia, feeding difficulties, mild short stature, and developmental delay.Thorough analysis of the 26 breakpoints on chromosome 1 revealed six known diseaserelated OMIM disease genes to be disrupted or affected by a deletion (RYR2, MFN2, CAMTA1, SLC9A1, PRDM16 and PLOD1).Two of the OMIM genes (CAMTA1 and PRDM16) were identified using FISH in a previous study of the same case, while remaining genes were novel findings using WGS [16].CAMTA1 is known to cause autosomal dominant cerebellar ataxia (non-progressive) with mild intellectual disability (MIM# 614756), with phenotypes including delayed psychomotor development, cerebellar ataxia, intellectual disability, neonatal hypotonia, and variable dysmorphic features, some of them consistent with the phenotype in Case 3. Remaining genes have not been associated with any phenotypes present in Case 3 and the full phenotype (ectopic kidney, hearing disability) in this individual could not fully be explained by the WGS analysis, but needs further investigation of the genes affected by the rearrangement.
There are a number of events that could lead to complex rearrangements, including replication-based mechanisms (fork-stalling and template-switching (FoSTeS) model) [24] and chromothripsis [25].WGS allowed for detailed analysis of all breakpoint junctions on the nucleotide level in Case 1 and Case 2, and all junctions but one in the two complex translocations were blunt, with maximum of two nucleotides microhomology.In one breakpoint junction, there were two small insertions of three and four nucleotides, respectively.They did not seem to be duplicated or templated from nearby sequences, which would have indicated a replicative error mechanism, instead the mutational signatures in both Case 1 and Case 2 indicate non-homologous end-joining (NHEJ), characteristic of chromothripsis rearrangement junctions.The characteristics of the breakpoint junctions in Case 3 makes it likely that the complex rearrangement of chromosome 1 formed through a single catastrophic event.These characteristics include randomness of the DNA fragment joins, the DNA fragments appear to be randomly joined in inverted/non-inverted orientation, and the ability to walk along the derivative chromosome [26].However, the rearrangement does not include a regularity of oscillating copy-number states [26], but only involves a small number of randomly spread deletions.In addition, the breakpoints of the q and p arms are separated by a 171 Mb DNA fragment, hence the breakpoints are not clustered as would be expected [26].Further, the p and q arms of chromosome 1 seem to have been brought close together in a ring-like formation.Possibly, the p and q arm were brought together before the scattering of the chromosome, otherwise, the fragments of one arm would either be lost, or the fragments would be less prone to translocate between the arms.Hence, the rearrangement of chromosome 1 in Case 3 could have arisen either from the halted formation of a ring chromosome, or even through a chromothripsis event of a ring chromosome.
The four WGS technologies performed were utilized in three different settings: i) for solving the derivative chromosome structure of the three CCRs, ii) for a comparison of detection rate of polymorphic CNVs first detected by aCGH in the three cases and iii) for a general assessment of genome wide SV calls from the three cases as compared to calls present in the public datasets [19,21,22].
First, we found that no technology provides a significantly higher detection rate than the others regarding the ability to detect and solve the structure of the de novo CCRs presented herein.This is partly a result of the relatively small number of CCR junctions in this study (35), but also due to the high detection rate of all four WGS technologies.Moreover, although those CCRs are complex in nature, the majority of the breakpoint junctions could be uniquely mapped to regions that do not present complex repetitive genomic patterns such as LCRs, satellites, centromeric or telomeric repeats.This is exemplified for LCR-containing junctions 1 and 2 of Case 3 that could only be resolved by optical mapping and linked-read WGS.Those technologies require longer DNA molecules and therefore are more appropriated to resolve repeats than short reads.These regions are also more prone to rearrangements [27], and hence the ability to resolve junctions mapping within those structures is of great importance.
Although the detection rates of these methods are similar, the resolutions differ.Both PCRfree PE WGS and linked-read sequencing can resolve most breakpoints to base-pair resolution.In contrast, optical mapping provides the lowest resolution (~7 kb), which likely explains why it failed to detect breakpoints involving the smaller fragments of Case 1 and Case 3.
To assess whether the CCRs would have been detected in a "WGS-first" scenario the MP and PE WGS data was filtered based on allele frequency, but filtered variant lists of the linkedread and optical mapping calls could not be obtained as there are no frequency databases available for those technologies.Given the low similarity between PCR-free PE WGS and these two methods, the databases such as the SweGen cohort [17] or 1000 Genomes [28] would be of very limited use.Hence, in order to make linked-read WGS and optical mapping usable in a clinical setting, large populations would need to be sequenced using these methods, and the data made available through frequency databases.
Finally, comparing the four technologies to polymorphic CNVs, it was found that PCR-free PE WGS provides the highest sensitivity, closely followed by linked-read sequencing.MP WGS and optical mapping performed similarly, with almost half the detection rate of PCRfree PE WGS.Notably, both these two methods failed in detecting smaller CNVs (<10 kb).Furthermore, PCR-free PE WGS did find a significant number of CNVs that linked-read WGS failed to detect: these CNVs did either belong to Case 1, having partially degraded DNA and too short molecules to be sequenced by linked-read WGS, or the CNVs were subsequently detected using CNVnator instead, which is a read-depth caller and able to detect variants as small as 2 kb [29].These variants exemplify that the large amounts of high performing Illumina WGS callers provide an edge over these more recent methods, whose pipelines are less mature, and still undergoing rapid development.
Lastly, the technologies were evaluated through a general assesment of the calls, as well as through comparison to the DGV [19], HG002 integrated call-set [22] and Conrad et al. (2010) datasets [21].This comparison provided two valuable insights not found through the previous analyses.First, the technologies produce significantly different amounts of variants: both linkedread WGS and optical mapping produce more than 10,000 calls on a single individual.In contrast, MP WGS generates less than 1000 calls on a single individual.Given the relatively similar detection rate on the CCRs and the fact that nearly all MP WGS calls are found in DGV, the precision of MP WGS is likely to be high compared to optical mapping and linked-read WGS.Second, it is also clear that each method report variants of different types and sizes and that the reported variant not always reflect the nature of the rearrangement in a given sample.This is particularly observed in the optical mapping results, which do not report any duplications at all, even though optical mapping clearly do detect duplications.Similarly, the Long Ranger SV caller for the linked-read WGS data is the only caller to report variants of "unknown" type.
In aggregate, through this comparison, we found that the four WGS methods produce variants of different sizes and these findings are in accordance to the resolution estimates.Notably, there is only a small overlap between the technologies and we only present two orthogonal methods (Bionano optical mapping and Illumina sequencing).Furthermore, many of the calls presented in DGV [19] or Conrad et al. [21] as structural variation are not validated.The results are nevertheless consistent with the previously shown results: MP WGS provides lower detection rate, but higher precision.PCR-free PE WGS, linked-read WGS, and optical mapping detect similar numbers of variants overall, however there are some biases toward certain variant types (for example, optical mapping reports a large number of insertions).Given the similar numbers of detected variants but different numbers of reported variants, the precision of the technologies are likely to be different: we observed that MP WGS provides the highest precision, PCR-free PE WGS the second most precision, and linked-read WGS or optical mapping the worst precision, depending on how the pipelines are combined.
The present study is limited by the fact that we have not compared any third-generation sequencing technologies, which have shown great promise for detection of structural variation in complex repeat regions [30,31].However, we found that short-read WGS combined with optical mapping is a powerful combination for analyzing CCRs.Combined, these two technologies would enable detection and validation of most breakpoints in two experiments, at maximum resolution.Given the current high cost of single molecule sequencing, a combined approach could be the most cost efficient.In this study we were able to detect all breakpoint junctions except one using linked-read WGS.The only missing breakpoint junction was in Case 1, where the DNA was partly degraded and no replacement DNA was available due to the individual being deceased.Taken together, and looking at the linked-read WGS result for Case 2 and Case 3, we are confident that the missing junction would have been found with better quality input DNA.Hence, in the cases reported here the linked-read sequencing identified all rearrangement breakpoints including those located in repetitive regions and is a valid WGS method of choice to detail complex rearrangements that often have breakpoints in repetitive regions.However, before this can be used in a clinical setting more user-friendly analysis software, as well as better reference data for filtering is desirable.
In conclusion, these findings demonstrate how different high throughput genomic methods can add clinically relevant information to conventional molecular analysis methods and enable characterization of the true nature of de novo CCRs.Finally, Case 3 demonstrates the need for long-molecule sequencing or complementary optical mapping to short-read sequencing to be able to map the structure of a highly complex rearrangement with breakpoints in repetitive regions.

Materials and methods
Materials and methods for the cytogenetic studies are presented in S1 Appendix and quality control (QC) data of the WGS methods are presented in S5 Table.

Cases
Three cases were studied following referral to the Clinical Genetics at Karolinska University Hospital, Stockholm, Sweden, due to clinical symptoms indicating genetic testing.Parental samples from Case 2 and Case 3 were analyzed using karyotyping and FISH, showing that the rearrangement had occurred de novo on the paternal allele [15,16].Parental samples from Case 1 and Case 3 were sequenced using linked-read sequencing, and the rearrangement in Case 1 was also originating from the paternal allele.Clinical parameters and phenotypes of the included cases are presented in Table 10.

Whole genome sequencing
Short-read whole genome sequencing.Genomic DNA from three individuals were sequenced using two separate Illumina WGS protocols, a 30X PCR free protocol, as well as a 3X 2.5 kb insert-size mate-pair protocol, at NGI (National Genomics Infrastructure), Stockholm, Sweden (https://ngisweden.scilifelab.se/).The data was processed using the NGI-piper and structural variants were analyzed using the FindSV pipeline (https://github.com/J35P312/FindSV), a pipeline combining CNVnator V0.3.2 [29] and TIDDIT V2.2.4 [10].Briefly, CNVnator detects CNVs based on read depth signatures, while TIDDIT detects a wide range of SVs by searching for clusters of discordant pairs and split-reads.Hence these callers are complementary, and together they make use of all SV signatures within the WGS data.The outputs of these two callers were combined into one single Variant Calling Format (VCF) file, which was annotated by variant effect predictor (VEP) 89 [32] and filtered based on the VCF file quality flag.The VCF file was subsequently sorted based on a local structural variant frequency database.The PCR-free PE calls were ranked using a database built from the SweGen cohort [17], while the MP WGS dataset was filtered using an in-house database consisting of The variant calls of Case 2 were compared to the DGV, an integrated SV call-set produced by GIAB (ftp://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/AshkenazimTrio/analysis/NIST_SVs_Integration_v0.6/README_SV_v0.6.txt), and the list of CNVs presented in Conrad et al. from 2010 [21].The DGV and the Conrad et al. datasets were converted into VCF files using the DGV2vcf and Conrad2vcf scripts (https://github.com/J35P312/convert2vcf). The resulting VCF files were then split into one VCF per variant type, and compared to the WGS technologies using SVDB merge V1.2.2, which was run using the compare_conrad and compare_dgv scripts (https://github.com/J35P312/convert2vcf). The integrated GIAB call-set was downloaded from the GIAB FTP (ftp://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/AshkenazimTrio/analysis/NIST_SVs_Integration_v0.6/HG002_SVs_Tier1_v0.6.vcf.gz), and filtered for variants detected using PacBio or Complete Genomics data.The filtering was performed by reading the CGcalls and PBcalls entries in the info column of the VCF file, and any variant having non-zero CGcalls or PBcalls value were kept while all other variants were filtered out.
Finally, the four technologies were compared by computing intersect of each pairwise technology-combination, and by calculating the sizes and number of variant calls.These scripts are also made available through (https://github.com/J35P312/convert2vcf). Across all comparisons, two variants were considered similar if their overlap exceed a Jaccard index of 0.4, and if the distance between the breakpoints is less than 100 kb.Any variants not meeting these criteria were considered dis-similar.Throughout these comparisons, only high quality calls were considered: a call was considered to be of high quality if the VCF filter flag was set to PASS.

Breakpoint junction PCR
Primers flanking all junctions except junction 1 and junction 2 in Case 3 were designed approximately 500 base pairs away from the estimated breakpoints.Same primers were subsequently used for sequencing using the Sanger method.Primer sequences are available on request.Breakpoint PCR was performed by standard methods using Phusion High-Fidelity DNA Polymerase (ThermoFisher Scientific, Waltham, MA, USA).Sequences obtained were aligned using BLAT (UCSC Genome Browser) [34] and visualized in CodonCode Aligner (CodonCode Corp., Dedham, MA, USA).

Fig 1 .
Fig 1.Molecular characterization of two complex interchromosomal rearrangements.(A) The complex chromosomal rearrangement identified in Case 1, t(2;8;15) (q34;q23.3;q21.3).A schematic illustration of the derivative chromosomes 2, 8 and 15 is shown to the left.The 11 genomic fragments are labeled from A-K and colored according to the parental chromosomal origin with chromosome 2 in purple, chromosome 8 in green and chromosome 15 in orange.Three fragments from chromosome 2 have been translocated to both derivative 8 and derivative 15 with one fragment inverted (Fragment C, indicated by the orientation of the letter inside the box).To the right, two linear representations of the whole genome sequencing results are shown.On top, the fragments are outlined on the parental chromosomes with links between fragments illustrated as dashed lines.The bottom diagrams show the final derivative chromosomes.For both fragment copy number status is indicated as black (normal) or red (deleted) and inverted orientation is highlighted in blue and marked by an arrow.(B) The complex chromosomal rearrangement identified in Case 2, 46,XY,t (1;10;5)(q31.3;p12.31;q23.2),as in A. The eight genomic fragments are labeled from A-H and chromosome 1 origin is shown in blue, chromosome 5 in green and chromosome 10 in pink.https://doi.org/10.1371/journal.pgen.1007858.g001 given in kb from the p telomere (Hg19).The microhomology (MH) column indicates the presence of microhomology between position A and B, while the insertion column describes any inserted sequence.The repeat column presents any repeat found within 100 bp of the breakpoints that could have mediated the formation of the derivative chromosome.The UCSC repeat masker was used to determine the position of repeats throughout the genome.Chr: chromosome, Pos: position, SV: structural variant, MH: microhomology https://doi.org/10.1371/journal.pgen.1007858.t002

Fig 2 .Fig 3 .
Fig 2. Molecular characterization of an extremely complex intrachromosomal rearrangement of chromosome 1.Whole genome sequencing enabled mapping of the rearranged chromosome 1.Chromosome 1 had been disrupted in 33 positions, and puzzled back together with deletions or inversions of some fragments.Junction numbers and characteristics for each junction are listed in Table 3.On the left side is a schematic of human chromosome 5 with red boxes indicating regions containing breakpoints.The 34 chromosomal fragments labeled from A-i are shown as a vertical linear diagram with fragment size given in Mb.In the middle the whole genome sequencing results are illustrated.Links between fragments are shown as dashed lines and the junction numbers are given.The final rearranged chromosome is shown on the right with deleted fragments in red and inversions indicated by an arrow.https://doi.org/10.1371/journal.pgen.1007858.g002 junctions, and on the other side are several poorly annotated reads with an unexpectedly high coverage, indicative of repetitive regions.In F, the optical maps are split at positions chr1:2581232 and chr1:2684269 (left), and the same maps continue at chr1:14176001 and chr1:27452000, respectively (right).The short-read WGS data in A, B, D and E is visualized in IGV and the optical maps data in C and F is from Bionano access.https://doi.org/10.1371/journal.pgen.1007858.g003

Fig 4 .
Fig 4. Percentages of the amount and type of structural variant (SV) calls generated by each technology and analysis pipeline, based on size.Approximately 90% of calls from Bionano optical mapping are between 10-100 kb in size and few calls are larger than 100 kb or smaller than 1 kb.Calls from paired-end (PE) WGS data is quite evenly spread, with slightly more calls <1 kb.The structural variant calls from the mate-pair (MP) WGS are sized between 10-100 kb.The linked-read WGS data was analyzed using three algorithms with different strengths, the Supernova de novo assembler, the Long Ranger large SV algorithm and the Long ranger deletion (Del) algorithm.Combining all callers, the linked-read WGS data produce a very high number of calls <1 kb.https://doi.org/10.1371/journal.pgen.1007858.g004

Table 3 . Overview of breakpoints and genomic regions involved in Case 3.
Genomic coordinates are given in Hg19.The microhomology (MH) column indicates the presence of any microhomology between the start and end positions, while the insertion column describes any inserted sequence.The repeat column presents any repeat found within 100 bp of the breakpoints.OMIM disease related genes are shown in bold text.Jct: junction, MH: microhomology, N.i.: no information � = estimated breakpoint https://doi.org/10.1371/journal.pgen.1007858.t003

Table 5 . General comparison of the sequencing technologies.
related to the complex chromosomal rearrangement that were detected from the total amount of structural variants for each method.The total number of breakpoint junctions supported by split reads is indicated within parentheses.Resolution is the median distance between the generated calls and the exact breakpoint position.https://doi.org/10.1371/journal.pgen.1007858.t005

Table 6 . Comparison of the four WGS technologies versus array comparative genomic hybridization. Method Case 1 custom 400K array Case 2 custom 400K array Case 3 custom 400K array
� Long Ranger-large SV �� Long Ranger-deletions ��� Supernova.Del: deletion, Dup: duplication, Sens: sensitivity.The aCGH row (confirmed variants) indicates the number of aCGH variants that were confirmed with any of the four WGS technologies.Each sample was analyzed using two arrays: a custom 400K array, as well as a commercial medical exome 1M array.https://doi.org/10.1371/journal.pgen.1007858.t006

Table 7 . Matrix table comparing the amount of genome wide overlapping SV calls detected by four WGS technol- ogies in Case 2.
WGS: whole genome sequencing, MP: mate-pair, PE: paired-end, within parentheses: each pairwise comparison given as percentages of overlapping SV calls.https://doi.org/10.1371/journal.pgen.1007858.t007biastowardsthe SV sizes it detects.For instance, linked-read WGS and PCR-free PE WGS produce a large fraction of small variant calls (<0.1 kb), in contrast, most of the optical mapping and MP calls are sized 10-100 kb (Fig4).