There's some preprocessing steps that I'm not quite clear on. To clarify: you are doing the following: - 2x101bp illumina paired-end sequencing - adapter trimming which can reduce your reads to as short as 25bp due to highly fragmented source material - merge overlapping read pairs into a longer synthetic read - what level of contamination do you have, and how are you dealing with it?

Are you error correcting the overlapping bases?

You are correct in asserting that there isn't much of signal in your data for SV tools to work with. Since your genome is so small, de novo assembly is not computationally prohibitive so that's probably worth trying but you're going to lose most of reads since they're so short. You'll also assembly a whole lot of junk if you have a high level of contamination. If you're interested in ancient novel sequence that has since been lost, this is probably your only option.

If you use a traditional SV tool, then you're just going to have to deal with the high signal to noise ratio due to your data. Pindel works pretty well with short reads since it was developed in the era of 2x36bp sequencing but can only detect certain events and has a high false discovery rate. Another option is to use GRIDSS (disclaimer: my software) and see how much of a signal it can get out of the noise. GRIDSS augments the traditional read pair and split read signals with a breakend assembly step. Essentially, it attempts an assembly at every genomic locus that there might be an SV event, then does split read analysis on the breakend assembly contigs to identify the other loci involved. This approach allows some of the reads (such as short soft clips) that are ignored in traditional split read analysis to contribute to the variant call.

With such short fragments, there's not much of a signal to start with but it's still better than manually inspecting 3Mb of sequence in IGV.

I have to say most of the approaches on SV-detection depend on paired-end data. No insert-size so already there is a data loss.

Also, since 30 bases (starts at 24-26bp) would be the least you need to get correct mapping information you would lose some data. Even if you need split-reads, if your read is 50bp and 30+20 is the split, the ambiguity of mapping for 20bp is high or almost unmappable. You do get information from longer reads, one developer suggested that you need atleast 70bp to actually trust break-points after various levels of testing (GRIDSS).

I think you might be using the leehom post-processing, usually 10bp for merging overlapping PE-data is quite minimal. FLASH uses 25 and they have an error rate close to 0.25%.

So whatever you might want to do with the existing SV-detection pipelines might be limiting. Also, you have an error-rate which might prove difficult at such small read-lengths. As your colleague suggested, the best way is to go for atleast local-assemblies where you expect break-points or a complete denovo. It might not be that subjective, might be time-consuming - but we are not that close to push-button bioinformatics anyway :)

We've had pretty decent luck with structural variation in bacteria with split read alignments using bwa mem, then generating BigWig files from alignments filtered to only contain split read alignments as a coverage metric. As Rohit mentions you would probably want to use only reads >= 70nt long.

In our case these were also paired-end, but frankly this was only used to corroborate the possible sites; normally read depth with bacterial sequences is pretty high and the splits piled up to a high end degree we just looked for spikes of reads using the coverage information in IGV. From this we were able to identify the cause of a pleiotropic mutation (transposon insertion in a gene).