Hi all
Over the years, I noticed that align.seqs fails to be consistent in treating gaps. In V4, for example, I run align gaps against SILVA (at it some or trimmed to V4) and one can get, on a conserved region, something like
AGCTCC-AATAGCGTATATTAA-AGTT-GT-T-GTGGTTAAAAA-GCTCG
AGCTCC-AATAGCGTATACTAA-TGTT-GT-T-GCAGTTAAAAA-GCTCG
AGCTCC-AATAGCGTATATTAA-AGTT-GT-T-GCGGTTAAAAA-GCTCG
AGCTCC-AATAGCGTATATTAA-AGTT-GT-T-GCGGTTAAAAC-GCTCG
AGCTCC-AATAGCGTATACTAA-TGTT-GT-T-GCAGTTAAAAA-GCTCG
AGCTCC-AATAGCGTATATTAA-GGTT-G–TTGCGGTTAAAAA-GTCCG
AGCTCC-AATAGCGTATATTAA-AGTT-GC-T-GCGGTTAAAAA-GCTCG
AGCTCC-AATAGCGTATACTAA-CGTT-GT-T-GCAGTTAAAAA-GCTCG
AGCTCC-AATAGCGTATATTAA-AGTT-GT-T-GCGGTTAAAAA-GCTCG
AGCTCC-AATAGCGTATACTAA-CGTT-GT-T-GCAGTT-AAAA-GCTCG
AGCTCCAA-TAGCGTATATTAA-AGTT-GT-T-GCGGTTAAAAA-GCTCG
AGCTCC-AATAGCGTATATTAA-AA-TTGT-T-GCGGTTAAAAC-GCTCG
AGCTCC-AATAGCGTATATTAA-AGTT-GT-T-GCGGTT-AAAAAGCTCG
The bold is a place where the aligner decided that the gap moved a bit was better (there are more, but the text editor is having issues accepting more than one bold section… THere are three other events in other gap columns). But looking at all, is clearly a misplacement. Happens both with needleman and gotoh.
Is there any way to mitigate this behavior within mothur? I know this alignment system doesn’t compare between reads, but those mistakes add gaps where they are not, and remove others where there are, and I am worried about the future unique and denoising.
Hi there,
What reference alignment are you using? This is largely dictated by the reference and that is the Achilles heal of this method. You might look at the align.report file for these sequences and see which reference sequence it and another sequence are matching to. If you wanted, you could go in to ARB and fix the alignment and then try again. My sense is that there are enough differences between this and other sequences that the very occasional one nt difference due to alignment doesn’t matter compared to the other differences.
Here’s the sequence with the bolded positions and the sequence from above it…
AGCTCC-AATAGCGTATACTAA-TGTT-GT-T-GCAGTTAAAAA-GCTCG
AGCTCC-AATAGCGTATATTAA-GGTT-G--TTGCGGTTAAAAA-GTCCG
x x x x x xx
I count ~7 differences between the one with the bolded letters an the sequence above it out of ~45 nucleotides. In this case, the 2 differences due to alignment isn’t going to make a difference since 5/45 or 7/45 are both over 10%. Comparing to the sequence below it, there are 5 differences and it isn’t clear that there are alignment errors. Even if we counted 2 alignment errors, the distance would be 3/45 ~ 6.7%
Pat
Hi Pat - on the chunk above you can wee hoe it moved two gaps. Maybe the differences are huge down the sequence, it is true, but I was surprised how often this happens. The ones above, despite looking identical there, then they are going to have difference of 4nt to start (that on a 350 bp region - 18S V4) starts to be a bit worrying.
Am aligning to SILVA, both complete and trimmed to V4.
Am going to give a shot to SINA (taking the alignment degapped and throwing it to SINA) to see if it handles that better. Maybe not, but since it compared to multiple hits might improve?
About fixing it in ARB - this is a massive V4 dataset, with about 16M unique by now, so this is not doable. And this happens relatively frequently (I just do sanity checks getting 100-500 sequences and I see them all the time), even if just 1% would be not doable by hand.
I suspect the quality of the 18S reference alingment is worse than the 16S. Are you using the seed or the nr version of the database? Again, the seed will be better than the nr. As for my comment about fixing it in ARB, that was fixing the reference alingment, not your data.
Pat
I use the nr version, not the seed. I will give it a shot then (nr). SINA will finish soon too, so I will update here with both results
well, SINA was a much better approach and, although a bit more time-consuming, not a big deal. Trying to solve the fact that it changed my T by U 
but I think I will follow up this step. The gaps were all aligned on the 5’ end, but failed more to do so on the 3’ end. Playing a bit with the settings to see if I can solve this! but feeling much better now.
Did you try our seed? FWIW, their reference alignment is somewhat proprietary - they don’t say what’s in it and there are alignments in it that are not generally available.
Pat