Aligning to a genome with an "assembly gap" (NNNNN) calls a deletion

[dannagifford]

Hello,

Something I noticed recently working with some new de novo assemblies I made. For context, we made the assembly using shovill, and as (it turns out) our new strain was very close to an existing reference genome, we used RagTag to help scaffold the assemblies. (And then prokka for annotation.)

RagTag introduces N's at assembly gaps.

We then used breseq to map the original Illumina reads back to the ragtag assembly. Where there are N's/assembly gaps, breseq calls these as a deletion, but:

• the MC plots show reads spanning this region---they just don't match the N's in the reference genome
• the JC evidence suggests there's a new junction at either side of the N's

I'm not sure if this logic makes sense---I think the fact that there are reads that span the assembly gap means that, possibly, those reads could be used to correct the gap? I don't think it's a "deletion"?

I've attached the html files for the JC/MC evidence. I could possibly share more if needed, though I wonder if you could replicate this just by replacing arbitrary genome with N's...

Ns-called-as-deletion.zip

Shovill and RagTag for info:
https://github.com/tseemann/shovill
https://github.com/malonge/RagTag
https://github.com/tseemann/prokka

[jeffreybarrick]

At least for the one you are showing, breseq is telling you to cleanly delete all 100 N's here.*

When breseq shows read bases in lowercase, it indicates that part of the read does not match the reference sequence shown above. This occurs when it splits the reads that are supporting a JC. If you look at ....TCAATGAGAA on the left side in the reference, that exactly matches the lowercase part on the right side of every read.

Probably, you expected there to be a big gap here and the reads to align to both sides (that would make things more obvious!), but what breseq does internally for JC is that it makes two copies of the read, each with one "half" aligned. Normally these are distant from each other in the chromosome, so keeping it as one and putting gaps into the CIGAR string wouldn't make sense. Here though, it would be nicer if it did that. If you click on the other side of the JC, you'll get uppercase and lowercase reversed for each read.

If there had been different bases where the N's are in the reference sequence, the junction would say something like +ATGGATTAGGC (whatever the novel sequence was).

Does this make sense?

(*We can't entirely rule out something more complicated going on b/c you have a few reads with the consistent gap and T substitution on the right side, which could indicate something more complicated involving near-repeats that are close in size to your read length. Probably, this disagreement is why RagTag refused to join these two sides and instead put in the scaffolding.)