BACKGROUND
genomic structural variants  are defined as alterations in the dna that affect > <dig> bp that may delete, insert, duplicate, invert, or move genomic sequence  <cit> .
svs are shown to be common in human genomes  <cit> , which caused increased interest in the characterization of both normal , and disease-causing large variants  <cit> .
furthermore, svs are known to be one of the driving forces of creation of new haplotypes  <cit>  and evolution  <cit> .
a subset of svs, namely copy number variations , were initially identified using bacterial artificial chromosome  and oligo array comparative genomic hybridization   <cit> , and snp genotyping arrays  <cit> .
a more detailed map of svs was made possible using fosmid end sequencing  <cit> ; however this method was too expensive and time-consuming since it involved creating and plating of fosmid libraries followed by sanger sequencing.
introduction of high throughput sequencing  finally made it possible to screen the genomes of many  to thousands  <cit>  of individuals.
although there are now many algorithms to discover and genotype svs using hts data  <cit> , they mainly focus on cnvs, which change the amount of dna, such as deletions, duplications, insertions, and retrotranspositions.
other types of svs, namely balanced rearrangements such as inversions and translocations are harder to detect due to the fact that their breakpoints usually lie within complex repeats that reduce mappability.
balanced rearrangements also do not alter the read depth, which makes the use of read depth signature  <cit>  irrelevant for their detection.
therefore, very few attempts have been taken to characterize inversions which are reliable only for small inversions  , and exhibit high false discovery rates in translocation call sets  <cit> .
another algorithm, gasvpro  <cit>  is also able to detect inversions with a size limit up to  <dig> kbp; however its sensitivity and specificity for large inversions are yet untested.
characterization of larger genomic inversions using hts remains an open problem.
motivation and approach
most known examples of large inversions have been identified in studies on human disease where inversions have no detectable effect in parents, but increase the risk of a disease-associated rearrangement in the offspring.
in the williams-beuren syndrome, for example, 25–30% of transmitting parents carry a  <dig>  mbp inversion encompassing a commonly deleted region, whereas the same inversion is present in only 6% of the general population  <cit> .
similarly, a polymorphic inversion has been reported at 15q11-q <dig> that gives rise to a de novo deletion resulting with the angelman syndrome  <cit> .
two more striking examples are found in the sotos syndrome  <cit>  and the 17q <dig>  microdeletion syndrome .
in each of these disorders, where a de novo microdeletion arises, every parent studied to date carries an inversion at the same region.
all these inversions are enriched in segmental duplications at their breakpoints, leading to an increased susceptibility to non-allelic homologous recombination , which in turn elevates risk for disease-causing rearrangements in the offspring.
the typical presence of duplicated sequences at the inversion boundaries is also the major challenge for inversion detection.
creation of a map of inversion polymorphisms will provide valuable information regarding their distribution and frequency in the human genome.
such a map will be important for future studies aimed to unravel how inversions and the segmental duplications architecture associated with inverted haplotypes contribute to genomic susceptibility to disease rearrangements.
to fill this gap, the invfest  <cit>  database aims to collect a comprehensive set of inversions reported in the literature.
currently invfest hosts  <dig> validated inversions, of which  <dig> are larger than  <dig> mbp.
the common method to discover inversions is to analyze the read pair signature  <cit> , where the mapping strand of the read pairs spanning the inversion breakpoints will be different from what is expected .
for example, the illumina platform generates read pairs from opposing strands but if the dna frag ment spans an inversion breakpoint, the read pairs will have a discordant size and they will map to the same strand.
when the inversion is large, the real mapping distance between pairs increases, therefore increasing the chance of incorrect mapping due to the common repeats that usually map at the inversion breakpoints and on other chromosomes.
another complication in accurate detection of large inversions arises as other types of svs might occur within inversions or around inversion breakpoints further confusing the sequence signatures.fig.
<dig> sequence signatures used by the valor algorithm.
in the presence of an inverted haplotype in the sequenced genome, we look for both read pair and split clone signatures.
paired-end reads that span the inversion breakpoints will be mapped to the same strand with a large distance between them, instead of the concordant read pairs that map to opposing strands  <cit> .
large insert clones will show mapping properties similar to the split read sequence signature  <cit> , but since we do not have the full clone sequence, or sufficient coverage to assemble clones, we interrogate lengths of contiguous read mapping
large molecule sequencing for long range contiguity
the hts platforms generate data at very high rates with minimal cost.
however, since both the hts reads , and the dna fragments are very short , the mappability of the hts data is dramatically reduced in repeat-rich regions that harbor most of the inversion breakpoints.
recently several methods have been developed to address these issues, and effectively to obtain information over a longer range.
there are substantial differences between different approaches, yet, all have the same basic principles.
first, the dna is broken into very large molecules , then diluted and separated into a number of pools.
each pool is later barcoded and sequenced using the illumina platform.
the first method that followed this procedure, called pooled clone sequencing , used fosmid and bacs to clone the input dna  <cit>  to generate 40-to-  <dig> kb molecules.
the alternative approaches avoid the cloning step and generate different average molecule sizes.
these include the truseq synthetic long reads , 10x genomics linked- reads  <cit>  , dovetail genomics , and cpt-seq  <cit> .
in the remainder of this paper, we collectively refer to these technologies as long range sequencing for simplicity.
our approach to discover large genomic inversions up to tens of thousands base pairs using long range sequencing follows from the observation that, dna molecules  that span the inversion breakpoint will be split into two sections when mapped to the reference genome, also separated by a distance approximately the size of the inversion.
we call this sequence signature as split clones , which is similar to the split read sequence signature used by several sv discovery tools such as delly  <cit>  and pindel  <cit>  but has the advantage that can span over repetitive regions.
in contrast to split reads, incorrect mappings in large inversion regions using split clone signatures are less pronounced.
this is because split clones are identified from many reads where each read pair is mapped concordantly, rather than using shorter alignments of single reads and then encompass larger regions that are longer than the repeats.
combining split-clone signatures with paired-end read signatures, we can distinguish the true paired-end read signatures and here even one pair of paired-end reads would suffice to detect the presence of the inversion and thus this approach allows us to detect very large inversions with low false discovery rate.
based on these observations, we developed a novel combinatorial algorithm and statistical heuristics called valor  .
briefly, valor searches for both read pair and split clone sequence signatures  using the mapping locations of long range sequencing reads, and requires split clones from different pools to cluster at the same putative inversion breakpoints .
ambiguity due to multiple possible pairings of split clones are resolved using an approximation algorithm for the maximal quasi clique problem  <cit> .
other tools such as variationhunter  <cit>  and clever  <cit>  use the set-cover or the equivalent maximum clique formulation  <cit>  to cluster the variants.
this approximation fails in clustering large inversions because it aims to detect complete cliques with low cardinality, which results in identifying a single breakpoint within repetitive regions of the genome .fig.
<dig> the valor algorithm.
we start with mapping paired-end reads for each pool.
we then separate read pairs that map to the same strand, and generate two files for the two strands.
we use the concordantly mapping reads to reconstruct clone locations, and calculate depth and breadth of coverage for the clones.
next, we remove clones with low coverage values, which results in the inferred clone locations.
we identify the split clones from this list, and then find psc, then we remove those with insufficient read pair support.
remaining pscs are used to construct a graph, which we then use to find maximal quasi cliques that signal possible inversion locations.
the clone and read pair support are then recalculated for the merged pscs, those with low support and those that intersect with reference assembly gaps, or intersect with segmental duplications in both breakpoints are discarded
fig.
<dig> clustering split clones to detect inversions.
a we first identify clone locations that are shorter than the expected clone size, but when paired with another short clone found in the same pool, the total length sums up to a full clone length.
we refer to such clones as “split clones”.
b we then cluster pairs of split clones that are mapped to approximately the same breakpoints.
note that due to read mapping errors and our clone reconstruction heuristics, a split clone may be identified as spanning a breakpoint.
c finally we cluster multiples of split clones from different pools if they agree on breakpoint location and the size of the inversion.
gap: size of the region between the start and end locations of split clones from different pools.
overlap: size of the overlapping region of split clones from different pools
valor proves its potential when tested on simulated data, and it is able to discover previously characterized large inversions in the genome of a human individual  using pooled bac sequence data.
additionally we tested valor using 10xg data generated from the same genome  <cit>  and obtained similar results.
kitzman et al.
<cit>  note that large clone  size is required to span segmental duplication blocks, and smaller clones such as fosmids may not be sufficient to detect inversions around segmental duplications.
in contrast we found on simulated data that fosmids perform as well as bac clones, if not better, given that fosmids are still statistically larger than most segmental duplications.
in conclusion, the theoretical minimum inversion size detectable by valor is limited by clone length, i.e.
<dig> kbp when bacs are used.
RESULTS
we designed various simulations to benchmark valor in terms of precision and specificity under ideal conditions, and to compare performance of valor with some of the popular sv discovery tools developed for whole genome shotgun  sequence data, such as lumpy  <cit> , delly  <cit> , and gasvpro  <cit> .
additionally, we tested valor using both pcs data  <cit>  and 10x genomics data generated from the genome of a human individual  .
simulations
using varsim  <cit> , we implanted  <dig> simulated inversions in different sizes  to the human reference genome .
varsim generated a diploid simulated genome, where 252/ <dig> inversions intersected with another simulated inversion.
varsim uses databases of previously validated variation in the human genome and it allows for inserting few novel inversions .
the simulated inversions are representative of realistic inversions.
we then simulated a wgs library at 50x coverage using art  <cit>  to benchmark the wgs-based tools.
the read length was  <dig> bp and the average fragment size was set to  <dig> bp with a standard deviation of  <dig> bp.
in order to test valor, we randomly generated a set consisting of  <dig> pools of simulated fosmid  clones from the same simulated chromosomes at 5x physical coverage.
we then simulated paired-end reads at 10x sequence coverage for each pool using art with a read length of  <dig> bp, average fragment size of  <dig> bp and standard deviation of  <dig> bp.
however, by random chance, no fosmid clones spanned breakpoints of  <dig> of  <dig> inversions larger than  <dig> kbp, and  <dig> of  <dig> inversions greater than  <dig> kbp.
we mapped both libraries to the reference genome using bwa-mem  <cit> , sorted and removed duplicates using samtools  <cit>  and picard  <cit> , and realigned around indels using gatk indelrealigner  <cit>  with default parameters.
we used valor on the simulated clone data set, and three popular sv discovery tools  on the wgs simulation.
in our tests, we required at least 50% reciprocal overlap between inversion intervals using the bedtools suite  <cit> .
here, due to the presence of heterozygous inversions, a predicted inversion may intersect with two simulated inversions.
we classify both such inversions as correctly identified in such cases.
we note that valor can theoretically discover inversions at least as large as the clone size; however, using fosmid clones valor was able to identify smaller inversions , although its performance and accuracy increase for larger variants .
of the  <dig> inversions , valor was able to accurately detect  <dig>  with only  <dig> false positives out of  <dig> calls .
the performance of valor increases slightly for larger inversions : valor discovered  <dig> out of  <dig> inversions , with  <dig> false positive calls .
among the wgs-based algorithms, delly performed the best in terms of recall  but it suffered from a high false discovery rate .
on the contrary, lumpy had the best precision  with low recall rates .
gasvpro performance was the lowest ; however, this is mainly due to the fact that gasvpro was designed for svs smaller than  <dig> kbp  <cit> .
overall, valor showed the best balance of precision and recall rates.table  <dig> simulation results on grch37
valor
valor
we implanted  <dig> inversions to the reference genome  using varsim and simulated two libraries, one pooled fosmid clone sequencing library for valor, and one wgs data set.
<dig> inversions had size > <dig> kbp, and  <dig> were > <dig> kbp.
a <dig> inversions  and  <dig> inversions  had no clone coverage.
bwhen inversions that had no clone coverage at breakpoints are removed.
tp true positive, fp false positive, fn false negative.
found: number of simulated inversions that intersect  with calls.
precision: positive predictive value, calculated as tp/.
recall: sensitivity, calculated as tp/.
note that due to diploid simulated inversions, one call may intersect with multiple implanted inversions
we also performed three more simulation experiments to comprehensively test the valor performance in the presence of other svs and segmental duplications, and its ability to differentiate between inversions and inverted duplications .
characterizing false discovery using a haploid genome
all tools developed for sv discovery suffer from high false discovery rates.
we used a haploid genome  to characterize the false positive calls.
to achieve this, we mapped illumina wgs reads at 100x coverage to the chm <dig>  assembly  <cit> .
following the fact that the chm <dig> cell line is haploid, and both the assembly and the wgs reads are derived from the same dna resource, we do not expect to find any real svs, but assembly errors may present themselves as variation.
we applied the same analysis steps described above.
overall, delly predicted  <dig> , gasvpro found  <dig> , and lumpy called  <dig> inversions, which shows the high false discovery rate of these tools.
similarly, we simulated clones  from the chm <dig> assembly at 5x physical coverage.
we then generated illumina reads using art at 10x sequence coverage per simulated clone.
in this experiment, valor returned no inversions as expected from the data.
pcs data
we tested valor using a real pcs data set generated from the genome of an individual of european descent .
we used genomic dna from na <dig> to construct the library.
high molecular weight dna was isolated, partially ecori digested, and subcloned into pcc1bac vector  to create a ~ <dig> kbp insert library using previously described protocols  <cit> .
we then split a portion of this library to  <dig> sets of  <dig> pools each, with  <dig> clones per pool for set  <dig>   <dig> clones per pool for set  <dig> and  <dig> clones per pool for set  <dig>  each pool was expanded by direct liquid outgrowth after infection.
we next constructed  <dig> barcoded sequencing libraries per each set, for a total of  <dig> sequencing libraries  <cit> .
libraries from each set were indexed with barcodes, combined and sequenced using the illumina hiseq platform .
upon sequencing a total of  <dig>  clones  we obtained  <dig> x expected physical depth of coverage.
after read mapping and clone reconstruction,  <dig> % of the genome was covered by one or more clones.
we mapped the paired-end reads from a total of  <dig> pools to the reference genome using bwa-mem.
average fragment length of the paired-end reads was ~ <dig> bp, with a standard deviation of ~ <dig> bp.
using valor, we reconstructed the clone locations, which showed an average clone length of ~ <dig> kbp and a standard deviation of ~ <dig> kbp.
the mapping quality and coverage of data was very low and almost all the pools in set  <dig> were contaminated and did not map to any chromosome, leaving us with  <dig> sets  and even lower coverage.
using valor, we identified a total of  <dig> inversions larger than  <dig> kbp .
we accurately detected all previously validated large inversions  .
additionally, valor was able to accurately predict a new inversion of  <dig> mb in size at the 16p <dig>  locus that we validated using fluorescent in situ hybridization .table  <dig> summary of validation of inversions predicted in the genome of na <dig> using valor
valor returns four coordinates for each inversion for two breakpoints.
the coordinates above are the inner breakpoint predictions, and are from the grch <dig> reference genome.
the invfest database reports inversions in ncbi reference build  <dig> coordinates, we thus converted the coordinates using the liftover tool
10x genomics linked reads
the linked-read sequencing  technology, developed recently by the 10x genomics incorporation sequencing technology, shows similarities to pcs, the dna is sheared into large molecules that are pooled, barcoded, and sequenced using the illumina platform.
reasoning from the similarities, we tested the valor using 10xg data generated from the same  genome  <cit> , and two other individuals of the same trio .
this data set was provided as bam files where approximately  <dig>  pools were tagged with barcodes per individual and each pool included ~ <dig> mbp of sequence.
valor is not yet trio-aware, thus we analyzed each sample separately.
as in the pcs test, valor could detect all previously known large inversions  and the same novel inversion that we validated using fish .table  <dig>
valor predictions in the ceph trio genomes using 10x genomics data
inversions detected using 10xg.
check marks denote that the inversion is found in the genome of the corresponding individual.
those inversions marked with  are listed as predicted and those marked with  are listed as validated in invfest
whole genome sequencing of na <dig> and chm1
the simulation results above showed that the valor’s performance was comparable to the performance of lumpy and delly, but benchmark on real data is warranted.
we tested lumpy, delly, and gasvpro using a pcr-free wgs data set at 50x coverage.
we downloaded the bam file that corresponds to na <dig> from the european nucleotide archive, sequenced as part of the “illumina platinum genomes” project , and applied the aforementioned sv discovery tools.
delly was able to find all of the known large inversions; however, it returned a total of  <dig>  inversion calls.
we note that the known inversions that delly correctly identified had low quality score, thus simple score-based filtration from the large call set would also remove the true positives.
lumpy returned  <dig> inversion calls, but it failed to discover any of the previously known large inversions.
gasvpro found  <dig> inversions of size > <dig> kbp , and it was able to find one of the large inversions.
such results are expected since these tools are designed for smaller inversions.
in addition we used the chm <dig> resource to test the sensitivity of wgs-based tools and compared the results from the wgs-based tools against a well-curated and validated database of svs generated from the same library  <cit> .
we mapped the chm <dig> illumina wgs reads to the human reference genome , and used the same tools to discover svs.
briefly, out of  <dig> validated inversions reported in  <cit> ), delly identified  <dig> and lumpy found  <dig> inversions.
gasvpro failed to characterized any known and validated inversions.
since there are no real pcs data sets generated from chm <dig>  and the average validated inversion length is ~ <dig> kb, we could not test valor with the same approach.
experimental validation
after the initial split clone clustering and maximal quasi clique approximation , we filtered those inversion clusters without sufficient read pair signature support .
additionally, when pcs data is used, we also fi﻿ltered ﻿﻿those that were within segmental duplications or intersected with gaps.
we found that valor could correctly identify two inversions that were previously validated within the same genome .
these inversions include a  <dig> mbp inversion in 8p <dig>   <cit>  and a  <dig> mbp inversion in 15q <dig>   <cit> .
we then selected  <dig> additional inversions for experimental validation that were not experimentally tested before , but  <dig> of them could not be tested by fish due to the amount of segmental duplications at the breakpoints.
our validation efforts resulted in the discovery of a novel, previously undocumented  <dig> mb inversion at the 16p <dig>  locus .
moreover we compared our calls with previously validated inversions reported in invfest and found a match for one more inversion at the 17q <dig> locus.fig.
<dig> new inversion validated at chr16: <dig> ,093– <dig> , <dig>  in the na <dig> genome.
a read-pair and split-clone signature used by valor to detect the inversion.
a-b split clones span the proximal breakpoint and c-d split clones span the distal.
b experimental validation of the novel inversion discovered using interphase fish
discussion
in this paper, we presented a novel algorithm, valor, to characterize large genomic inversions using two of the new sequencing methods that were initially developed to improve haplotype phasing.
we showed that valor is compatible with both pcs and linked-read sequencing technologies.
although it suffers from high fdr using real data , valor was able to identify all previously validated inversion events, and also discover a novel variant.
furthermore, valor performed better with simulated data suggesting that the relatively lower performance with the na <dig> pcs data set may be improved with higher clone  coverage with more pools with higher sparsity .
we note that valor still performed better than wgs-based tools in terms of large inversion detection sensitivity.
there are multiple directions that we can take to further improve valor.
first, to reduce the fdr, we can incorporate split read sequence signature  <cit> , and we can perform local de novo assembly around the predicted breakpoint intervals with an approach similar to tigra  <cit> .
however, since both of these methods need high and relatively uniform sequence coverage, they might not be suitable to directly apply to the data sets we used.
instead, it will be better to simultaneously use wgs data generated from the genome of the same individual.
since the pcs and 10xg methods also require additional wgs data for haplotype phasing, it can be expected to generate matching pcs/10xg-wgs data sets from the same genomes.
another future research on valor will be testing and improving its abilities to discover smaller, yet still large inversions .
in this paper, we focused on inversions larger than  <dig> kbp, because the upper size limit for gasvpro  <cit>  algorithm is  <dig> kbp, and only such large inversions can be reliably tested using fish.
note that validating smaller inversions is a more difficult task using fiber fish or pcr, if the breakpoints do not lie within unique regions.
in addition, the clone size distribution should be tighter to ensure clone reconstruction method does not artificially “merge” split clones into a single interval.
we still would like to investigate valor’s performance using real fosmid data, but this may require additional algorithmic enhancements especially in the presence of nearby segmental duplications  <cit> .
in this paper we present fosmid simulation experiments, and there is currently only one pooled fosmid sequencing dataset  <cit>  generated from the genome of a gujarati indian individual .
however, this data set has even lower coverage and data quality, since this is the first data ever generated with pcs during its development phase.
we would like to apply valor to a newer fosmid-based dataset and evaluate its per formance with experimental validation.
the chicago method from dovetail genomics, the truseq synthetic long reads  <cit> , and the cpt-seq technology  <cit>  are other candidates for further valor development.
valor can also be extended to characterize other forms of large sv, including deletions, insertions, direct and inverted duplications.
each of these types of svs present themselves with different split clone signatures .
we also note that, determining the location of a segmental duplication event is yet a largely unsolved problem, even when long reads are used  <cit> .
it may also be possible to discover translocations using split clones; however, chance of finding incorrect split clones will also increase, causing a reduction in the performance of maximal quasi clique approximation.
CONCLUSIONS
valor is the first algorithm that can discover large genomic inversions using hts technologies.
our understanding of the phenotypic effects of inversions is still limited, and one of the reasons of this is the lack of reliable and cost effective methods to characterize such events.
this is also true for other complex rearrangements such as duplications and translocations.
improvements in characterization of large complex rearrangements will help us better understand the biological mechanisms that lead to phenotypic difference, disease, and evolution.

 genomic structural variants are defined as alterations in the dna that affect > bp that may delete, insert, duplicate, invert, or move genomic sequence. svs are shown to be common in human genomes, which caused increased interest in the characterization of both normal, and disease-causing large variants. furthermore, svs are known to be one of the driving forces of creation of new haplotypes and evolution. a subset of svs, namely copy number variations, were initially identified using bacterial artificial chromosome and oligo array comparative genomic hybridization, and snp genotyping arrays. a more detailed map of svs was made possible using fosmid end sequencing; however this method was too expensive and time-consuming since it involved creating and plating of fosmid libraries followed by sanger sequencing. introduction of high throughput sequencing finally made it possible to screen the genomes of many to thousands of individuals. although there are now many algorithms to discover and genotype svs using hts data, they mainly focus on cnvs, which change the amount of dna, such as deletions, duplications, insertions, and retrotranspositions. other types of svs, namely balanced rearrangements such as inversions and translocations are harder to detect due to the fact that their breakpoints usually lie within complex repeats that reduce mappability. balanced rearrangements also do not alter the read depth, which makes the use of read depth signature irrelevant for their detection. therefore, very few attempts have been taken to characterize inversions which are reliable only for small inversions, and exhibit high false discovery rates in translocation call sets. another algorithm, gasvpro is also able to detect inversions with a size limit up to kbp; however its sensitivity and specificity for large inversions are yet untested. characterization of larger genomic inversions using hts remains an open problem. motivation and approach most known examples of large inversions have been identified in studies on human disease where inversions have no detectable effect in parents, but increase the risk of a disease-associated rearrangement in the offspring. in the williams-beuren syndrome, for example, 25–30% of transmitting parents carry a mbp inversion encompassing a commonly deleted region, whereas the same inversion is present in only 6% of the general population. similarly, a polymorphic inversion has been reported at 15q11-q that gives rise to a de novo deletion resulting with the angelman syndrome. two more striking examples are found in the sotos syndrome and the 17q microdeletion syndrome. in each of these disorders, where a de novo microdeletion arises, every parent studied to date carries an inversion at the same region. all these inversions are enriched in segmental duplications at their breakpoints, leading to an increased susceptibility to non-allelic homologous recombination, which in turn elevates risk for disease-causing rearrangements in the offspring. the typical presence of duplicated sequences at the inversion boundaries is also the major challenge for inversion detection. creation of a map of inversion polymorphisms will provide valuable information regarding their distribution and frequency in the human genome. such a map will be important for future studies aimed to unravel how inversions and the segmental duplications architecture associated with inverted haplotypes contribute to genomic susceptibility to disease rearrangements. to fill this gap, the invfest database aims to collect a comprehensive set of inversions reported in the literature. currently invfest hosts validated inversions, of which are larger than mbp. the common method to discover inversions is to analyze the read pair signature, where the mapping strand of the read pairs spanning the inversion breakpoints will be different from what is expected. for example, the illumina platform generates read pairs from opposing strands but if the dna frag ment spans an inversion breakpoint, the read pairs will have a discordant size and they will map to the same strand. when the inversion is large, the real mapping distance between pairs increases, therefore increasing the chance of incorrect mapping due to the common repeats that usually map at the inversion breakpoints and on other chromosomes. another complication in accurate detection of large inversions arises as other types of svs might occur within inversions or around inversion breakpoints further confusing the sequence signatures.fig. sequence signatures used by the valor algorithm. in the presence of an inverted haplotype in the sequenced genome, we look for both read pair and split clone signatures. paired-end reads that span the inversion breakpoints will be mapped to the same strand with a large distance between them, instead of the concordant read pairs that map to opposing strands. large insert clones will show mapping properties similar to the split read sequence signature, but since we do not have the full clone sequence, or sufficient coverage to assemble clones, we interrogate lengths of contiguous read mapping large molecule sequencing for long range contiguity the hts platforms generate data at very high rates with minimal cost. however, since both the hts reads, and the dna fragments are very short, the mappability of the hts data is dramatically reduced in repeat-rich regions that harbor most of the inversion breakpoints. recently several methods have been developed to address these issues, and effectively to obtain information over a longer range. there are substantial differences between different approaches, yet, all have the same basic principles. first, the dna is broken into very large molecules, then diluted and separated into a number of pools. each pool is later barcoded and sequenced using the illumina platform. the first method that followed this procedure, called pooled clone sequencing, used fosmid and bacs to clone the input dna to generate 40-to- kb molecules. the alternative approaches avoid the cloning step and generate different average molecule sizes. these include the truseq synthetic long reads, 10x genomics linked- reads, dovetail genomics, and cpt-seq. in the remainder of this paper, we collectively refer to these technologies as long range sequencing for simplicity. our approach to discover large genomic inversions up to tens of thousands base pairs using long range sequencing follows from the observation that, dna molecules that span the inversion breakpoint will be split into two sections when mapped to the reference genome, also separated by a distance approximately the size of the inversion. we call this sequence signature as split clones, which is similar to the split read sequence signature used by several sv discovery tools such as delly and pindel but has the advantage that can span over repetitive regions. in contrast to split reads, incorrect mappings in large inversion regions using split clone signatures are less pronounced. this is because split clones are identified from many reads where each read pair is mapped concordantly, rather than using shorter alignments of single reads and then encompass larger regions that are longer than the repeats. combining split-clone signatures with paired-end read signatures, we can distinguish the true paired-end read signatures and here even one pair of paired-end reads would suffice to detect the presence of the inversion and thus this approach allows us to detect very large inversions with low false discovery rate. based on these observations, we developed a novel combinatorial algorithm and statistical heuristics called valor. briefly, valor searches for both read pair and split clone sequence signatures using the mapping locations of long range sequencing reads, and requires split clones from different pools to cluster at the same putative inversion breakpoints. ambiguity due to multiple possible pairings of split clones are resolved using an approximation algorithm for the maximal quasi clique problem. other tools such as variationhunter and clever use the set-cover or the equivalent maximum clique formulation to cluster the variants. this approximation fails in clustering large inversions because it aims to detect complete cliques with low cardinality, which results in identifying a single breakpoint within repetitive regions of the genome.fig. the valor algorithm. we start with mapping paired-end reads for each pool. we then separate read pairs that map to the same strand, and generate two files for the two strands. we use the concordantly mapping reads to reconstruct clone locations, and calculate depth and breadth of coverage for the clones. next, we remove clones with low coverage values, which results in the inferred clone locations. we identify the split clones from this list, and then find psc, then we remove those with insufficient read pair support. remaining pscs are used to construct a graph, which we then use to find maximal quasi cliques that signal possible inversion locations. the clone and read pair support are then recalculated for the merged pscs, those with low support and those that intersect with reference assembly gaps, or intersect with segmental duplications in both breakpoints are discarded fig. clustering split clones to detect inversions. a we first identify clone locations that are shorter than the expected clone size, but when paired with another short clone found in the same pool, the total length sums up to a full clone length. we refer to such clones as “split clones”. b we then cluster pairs of split clones that are mapped to approximately the same breakpoints. note that due to read mapping errors and our clone reconstruction heuristics, a split clone may be identified as spanning a breakpoint. c finally we cluster multiples of split clones from different pools if they agree on breakpoint location and the size of the inversion. gap: size of the region between the start and end locations of split clones from different pools. overlap: size of the overlapping region of split clones from different pools valor proves its potential when tested on simulated data, and it is able to discover previously characterized large inversions in the genome of a human individual using pooled bac sequence data. additionally we tested valor using 10xg data generated from the same genome and obtained similar results. kitzman et al. note that large clone size is required to span segmental duplication blocks, and smaller clones such as fosmids may not be sufficient to detect inversions around segmental duplications. in contrast we found on simulated data that fosmids perform as well as bac clones, if not better, given that fosmids are still statistically larger than most segmental duplications. in conclusion, the theoretical minimum inversion size detectable by valor is limited by clone length, i.e. kbp when bacs are used.  we designed various simulations to benchmark valor in terms of precision and specificity under ideal conditions, and to compare performance of valor with some of the popular sv discovery tools developed for whole genome shotgun sequence data, such as lumpy, delly, and gasvpro. additionally, we tested valor using both pcs data and 10x genomics data generated from the genome of a human individual. simulations using varsim, we implanted simulated inversions in different sizes to the human reference genome. varsim generated a diploid simulated genome, where 252/ inversions intersected with another simulated inversion. varsim uses databases of previously validated variation in the human genome and it allows for inserting few novel inversions. the simulated inversions are representative of realistic inversions. we then simulated a wgs library at 50x coverage using art to benchmark the wgs-based tools. the read length was bp and the average fragment size was set to bp with a standard deviation of bp. in order to test valor, we randomly generated a set consisting of pools of simulated fosmid clones from the same simulated chromosomes at 5x physical coverage. we then simulated paired-end reads at 10x sequence coverage for each pool using art with a read length of bp, average fragment size of bp and standard deviation of bp. however, by random chance, no fosmid clones spanned breakpoints of of inversions larger than kbp, and of inversions greater than kbp. we mapped both libraries to the reference genome using bwa-mem, sorted and removed duplicates using samtools and picard, and realigned around indels using gatk indelrealigner with default parameters. we used valor on the simulated clone data set, and three popular sv discovery tools on the wgs simulation. in our tests, we required at least 50% reciprocal overlap between inversion intervals using the bedtools suite. here, due to the presence of heterozygous inversions, a predicted inversion may intersect with two simulated inversions. we classify both such inversions as correctly identified in such cases. we note that valor can theoretically discover inversions at least as large as the clone size; however, using fosmid clones valor was able to identify smaller inversions, although its performance and accuracy increase for larger variants. of the inversions, valor was able to accurately detect with only false positives out of calls. the performance of valor increases slightly for larger inversions: valor discovered out of inversions, with false positive calls. among the wgs-based algorithms, delly performed the best in terms of recall but it suffered from a high false discovery rate. on the contrary, lumpy had the best precision with low recall rates. gasvpro performance was the lowest; however, this is mainly due to the fact that gasvpro was designed for svs smaller than kbp. overall, valor showed the best balance of precision and recall rates.table simulation results on grch37  valor  valor we implanted inversions to the reference genome using varsim and simulated two libraries, one pooled fosmid clone sequencing library for valor, and one wgs data set. inversions had size > kbp, and were > kbp. a inversions and inversions had no clone coverage. bwhen inversions that had no clone coverage at breakpoints are removed. tp true positive, fp false positive, fn false negative. found: number of simulated inversions that intersect with calls. precision: positive predictive value, calculated as tp/. recall: sensitivity, calculated as tp/. note that due to diploid simulated inversions, one call may intersect with multiple implanted inversions we also performed three more simulation experiments to comprehensively test the valor performance in the presence of other svs and segmental duplications, and its ability to differentiate between inversions and inverted duplications. characterizing false discovery using a haploid genome all tools developed for sv discovery suffer from high false discovery rates. we used a haploid genome to characterize the false positive calls. to achieve this, we mapped illumina wgs reads at 100x coverage to the chm assembly. following the fact that the chm cell line is haploid, and both the assembly and the wgs reads are derived from the same dna resource, we do not expect to find any real svs, but assembly errors may present themselves as variation. we applied the same analysis steps described above. overall, delly predicted, gasvpro found, and lumpy called inversions, which shows the high false discovery rate of these tools. similarly, we simulated clones from the chm assembly at 5x physical coverage. we then generated illumina reads using art at 10x sequence coverage per simulated clone. in this experiment, valor returned no inversions as expected from the data. pcs data we tested valor using a real pcs data set generated from the genome of an individual of european descent. we used genomic dna from na to construct the library. high molecular weight dna was isolated, partially ecori digested, and subcloned into pcc1bac vector to create a ~ kbp insert library using previously described protocols. we then split a portion of this library to sets of pools each, with clones per pool for set clones per pool for set and clones per pool for set each pool was expanded by direct liquid outgrowth after infection. we next constructed barcoded sequencing libraries per each set, for a total of sequencing libraries. libraries from each set were indexed with barcodes, combined and sequenced using the illumina hiseq platform. upon sequencing a total of clones we obtained x expected physical depth of coverage. after read mapping and clone reconstruction, % of the genome was covered by one or more clones. we mapped the paired-end reads from a total of pools to the reference genome using bwa-mem. average fragment length of the paired-end reads was ~ bp, with a standard deviation of ~ bp. using valor, we reconstructed the clone locations, which showed an average clone length of ~ kbp and a standard deviation of ~ kbp. the mapping quality and coverage of data was very low and almost all the pools in set were contaminated and did not map to any chromosome, leaving us with sets and even lower coverage. using valor, we identified a total of inversions larger than kbp. we accurately detected all previously validated large inversions. additionally, valor was able to accurately predict a new inversion of mb in size at the 16p locus that we validated using fluorescent in situ hybridization.table summary of validation of inversions predicted in the genome of na using valor valor returns four coordinates for each inversion for two breakpoints. the coordinates above are the inner breakpoint predictions, and are from the grch reference genome. the invfest database reports inversions in ncbi reference build coordinates, we thus converted the coordinates using the liftover tool 10x genomics linked reads the linked-read sequencing technology, developed recently by the 10x genomics incorporation sequencing technology, shows similarities to pcs, the dna is sheared into large molecules that are pooled, barcoded, and sequenced using the illumina platform. reasoning from the similarities, we tested the valor using 10xg data generated from the same genome, and two other individuals of the same trio. this data set was provided as bam files where approximately pools were tagged with barcodes per individual and each pool included ~ mbp of sequence. valor is not yet trio-aware, thus we analyzed each sample separately. as in the pcs test, valor could detect all previously known large inversions and the same novel inversion that we validated using fish.table valor predictions in the ceph trio genomes using 10x genomics data inversions detected using 10xg. check marks denote that the inversion is found in the genome of the corresponding individual. those inversions marked with are listed as predicted and those marked with are listed as validated in invfest whole genome sequencing of na and chm1 the simulation results above showed that the valor’s performance was comparable to the performance of lumpy and delly, but benchmark on real data is warranted. we tested lumpy, delly, and gasvpro using a pcr-free wgs data set at 50x coverage. we downloaded the bam file that corresponds to na from the european nucleotide archive, sequenced as part of the “illumina platinum genomes” project, and applied the aforementioned sv discovery tools. delly was able to find all of the known large inversions; however, it returned a total of inversion calls. we note that the known inversions that delly correctly identified had low quality score, thus simple score-based filtration from the large call set would also remove the true positives. lumpy returned inversion calls, but it failed to discover any of the previously known large inversions. gasvpro found inversions of size > kbp, and it was able to find one of the large inversions. such results are expected since these tools are designed for smaller inversions. in addition we used the chm resource to test the sensitivity of wgs-based tools and compared the results from the wgs-based tools against a well-curated and validated database of svs generated from the same library. we mapped the chm illumina wgs reads to the human reference genome, and used the same tools to discover svs. briefly, out of validated inversions reported in ), delly identified and lumpy found inversions. gasvpro failed to characterized any known and validated inversions. since there are no real pcs data sets generated from chm and the average validated inversion length is ~ kb, we could not test valor with the same approach. experimental validation after the initial split clone clustering and maximal quasi clique approximation, we filtered those inversion clusters without sufficient read pair signature support. additionally, when pcs data is used, we also fi﻿ltered ﻿﻿those that were within segmental duplications or intersected with gaps. we found that valor could correctly identify two inversions that were previously validated within the same genome. these inversions include a mbp inversion in 8p and a mbp inversion in 15q. we then selected additional inversions for experimental validation that were not experimentally tested before, but of them could not be tested by fish due to the amount of segmental duplications at the breakpoints. our validation efforts resulted in the discovery of a novel, previously undocumented mb inversion at the 16p locus. moreover we compared our calls with previously validated inversions reported in invfest and found a match for one more inversion at the 17q locus.fig. new inversion validated at chr16:,093–, in the na genome. a read-pair and split-clone signature used by valor to detect the inversion. a-b split clones span the proximal breakpoint and c-d split clones span the distal. b experimental validation of the novel inversion discovered using interphase fish discussion in this paper, we presented a novel algorithm, valor, to characterize large genomic inversions using two of the new sequencing methods that were initially developed to improve haplotype phasing. we showed that valor is compatible with both pcs and linked-read sequencing technologies. although it suffers from high fdr using real data, valor was able to identify all previously validated inversion events, and also discover a novel variant. furthermore, valor performed better with simulated data suggesting that the relatively lower performance with the na pcs data set may be improved with higher clone coverage with more pools with higher sparsity. we note that valor still performed better than wgs-based tools in terms of large inversion detection sensitivity. there are multiple directions that we can take to further improve valor. first, to reduce the fdr, we can incorporate split read sequence signature, and we can perform local de novo assembly around the predicted breakpoint intervals with an approach similar to tigra. however, since both of these methods need high and relatively uniform sequence coverage, they might not be suitable to directly apply to the data sets we used. instead, it will be better to simultaneously use wgs data generated from the genome of the same individual. since the pcs and 10xg methods also require additional wgs data for haplotype phasing, it can be expected to generate matching pcs/10xg-wgs data sets from the same genomes. another future research on valor will be testing and improving its abilities to discover smaller, yet still large inversions. in this paper, we focused on inversions larger than kbp, because the upper size limit for gasvpro algorithm is kbp, and only such large inversions can be reliably tested using fish. note that validating smaller inversions is a more difficult task using fiber fish or pcr, if the breakpoints do not lie within unique regions. in addition, the clone size distribution should be tighter to ensure clone reconstruction method does not artificially “merge” split clones into a single interval. we still would like to investigate valor’s performance using real fosmid data, but this may require additional algorithmic enhancements especially in the presence of nearby segmental duplications. in this paper we present fosmid simulation experiments, and there is currently only one pooled fosmid sequencing dataset generated from the genome of a gujarati indian individual. however, this data set has even lower coverage and data quality, since this is the first data ever generated with pcs during its development phase. we would like to apply valor to a newer fosmid-based dataset and evaluate its per formance with experimental validation. the chicago method from dovetail genomics, the truseq synthetic long reads, and the cpt-seq technology are other candidates for further valor development. valor can also be extended to characterize other forms of large sv, including deletions, insertions, direct and inverted duplications. each of these types of svs present themselves with different split clone signatures. we also note that, determining the location of a segmental duplication event is yet a largely unsolved problem, even when long reads are used. it may also be possible to discover translocations using split clones; however, chance of finding incorrect split clones will also increase, causing a reduction in the performance of maximal quasi clique approximation.  valor is the first algorithm that can discover large genomic inversions using hts technologies. our understanding of the phenotypic effects of inversions is still limited, and one of the reasons of this is the lack of reliable and cost effective methods to characterize such events. this is also true for other complex rearrangements such as duplications and translocations. improvements in characterization of large complex rearrangements will help us better understand the biological mechanisms that lead to phenotypic difference, disease, and evolution. 
this is mainly due to the fact that breakpoints of such events typically lie within segmental duplications or common repeats, which reduces the mappability of short reads. the algorithms developed within the genomes project to identify inversions are limited to relatively short inversions, and there are currently no available algorithms to discover large inversions using high throughput sequencing technologies. although many algorithms are now available that aim to characterize different classes of structural variation, discovery of balanced rearrangements such as inversions remains an open problem. using valor, we also predicted a novel inversion, which we validated using fluorescent in situ hybridization. here we propose a novel algorithm, valor, to discover large inversions using new sequencing methods that provide long range information such as 10x genomics linked-read sequencing, pooled clone sequencing, or other similar technologies that we commonly refer to as long range sequencing. structural variationlong range sequencinglinked-readsinversionread cloudshttp://dx.doi.org/ /501100004963seventh framework programmemcirg 303772alkan can http://dx.doi.org/ /501100003043emboig-2521alkan can http://dx.doi.org/ /100000051national human genome research institutehg004120eichler evan e. firb futuro in ricercarbfr103ce3ventura mario issue-copyright-statement© the author 2017 in this paper, we show that valor is able to accurately discover all previously identified and experimentally validated large inversions in the same genome with a low false discovery rate.
 although many algorithms are now available that aim to characterize different classes of structural variation, discovery of balanced rearrangements such as inversions remains an open problem. this is mainly due to the fact that breakpoints of such events typically lie within segmental duplications or common repeats, which reduces the mappability of short reads. the algorithms developed within the genomes project to identify inversions are limited to relatively short inversions, and there are currently no available algorithms to discover large inversions using high throughput sequencing technologies.  here we propose a novel algorithm, valor, to discover large inversions using new sequencing methods that provide long range information such as 10x genomics linked-read sequencing, pooled clone sequencing, or other similar technologies that we commonly refer to as long range sequencing. we demonstrate the utility of valor using both pooled clone sequencing and 10x genomics linked-read sequencing generated from the genome of an individual from the hapmap project. we also provide a comprehensive comparison of valor against several state-of-the-art structural variation discovery algorithms that use whole genome shotgun sequencing data.  in this paper, we show that valor is able to accurately discover all previously identified and experimentally validated large inversions in the same genome with a low false discovery rate. using valor, we also predicted a novel inversion, which we validated using fluorescent in situ hybridization. valor is available at https://github.com/bilkentcompgen/valor electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords structural variationlong range sequencinglinked-readsinversionread cloudshttp://dx.doi.org/ /501100004963seventh framework programmemcirg 303772alkan can http://dx.doi.org/ /501100003043emboig-2521alkan can http://dx.doi.org/ /100000051national human genome research institutehg004120eichler evan e. firb futuro in ricercarbfr103ce3ventura mario issue-copyright-statement© the author 2017 