Finding incorrectly formatted needle in apparently mostly formatted properly haystack

[Nicholas_Bokulich]

@devonorourke

I’m assuming “old” was Mangan_bplot_vsearch.qzv, while “new” is " Mangan_bplot_vsearch_p90m50k.qzv".

Yes, I was using those same definitions.

@devonorourke

…odd. and yeah, I’m moving on from this for now. I’ll let you know if it jumps up again.

Thanks! Do let me know.

@devonorourke

What’s been frustrating me is that I’ve found vsearch OTUs being assigned as an arthropod, but blast pointing to a mammal. Not just any mammal, but a bat. And not just any bat, but the exact bat species I know should be in our sample. The blast search is absolutely correct in this one OTU instance.

Is this NCBI blast or classify-consensus-blast or standalone blastn?

@devonorourke

Okay, with that context, hopefully this explains my stupid parameters. I thought what I was doing with my vsearch perc-identity parameter was making it more stringent for a “hit” in the global alignment search to be retained. I think that’s right.

Yes, that will definitely make it more stringent (though percent-id is better for toggling stringency)... but since taxonomy is based on consensus here that is also causing problems for you! Here's what is happening:

250k unique seqs (let's say 1% == 2.5k of these are bat, rest are arthropod)
You grab 50k of these
2.5k will be bat (presumably should be the top hits!) and the remainder will be bug (provided they are above the similarity threshold).
your query will be classified as some type of arthropod because 47.5k/50k = 0.95 >> 0.51 (default consensus threshold)

So this is really a quirk of your database and the characteristics of the sequences and their similarities... but a really cool "edge case" for testing the limits of this classifier! Also makes a great test for optimization since you have very distinct species that are evidently somewhat similar at this locus...

@devonorourke

You could also use --maxaccepts 0 --maxrejects 0 , which will cause vsearch to consider all database sequences. It will take much longer as all the heuristics are bypassed

We use maxrejects=0, so the whole database is searched for matches before ranking the top hits, so using maxaccepts=1 will find you the true top hit (unless if I am missing something!).

But that is a great question — and actually blast maxaccepts works in that non-intuitive way (takes first N hits that exceed percent id, then skip the rest, instead of finding the true top hits) so great sleuthing!!!

@devonorourke

Error: Invalid value for "--p-maxaccepts": 0 is smaller than the minimum valid value 1.

Oops, that's just a very minor bug in setting the parameter limits, and I have raised an issue to fix it here. You solved it with the correct workaround — setting a very high number for maxaccepts.

@devonorourke

What I want is the most stringent terms possible, at least for one pass of the data, to determine which OTUs are absolutely bat-associated sequences.

Presumably bat sequences will be very similar to the ref.... so you can take a first pass of the classifier with a low maxaccepts and high percent-id to match bat seqs. Better yet, since you probably just want to remove the bat seqs since these are host DNA, just use qiime quality-control exclude-seqs to remove host DNA before classification!

Then do the real round of classification. You should probably use parameters similar to the defaults... use a smaller maxaccepts and use percent-id to limit matches to relevant matches

@devonorourke

I think the thing I should mention in all of this is that I just want to get the OTUs as correct as possible

percent-id will be the parameter to focus on for correctness. If you think that the samples contain many species that are not in the ref (or notice that empirically via high numbers of unassigned), then you need to dial back the percent-id and increase maxaccepts to get a broader consensus.

@devonorourke

Regarding benchmarking, I do conveniently have some mock communities that have been spiked into some of my datasets. These consist of about 20 known arthropod sequences and I could likely use them to run the benchmarking tools your team has so wonderfully created.

Awesome! That would be a great dataset to add to our mock community repository if you would like to share (probably after publishing).

@devonorourke

Feel free to follow up with me if there are more things to know than what has been provided in the link you’ve sent.

This thread discusses "getting started" — tax-credit was a sprawling project and I've done my best to make it easy to navigate but there is still a lot there. Follow the instructions in the READMEs instead of trying to navigate the notebooks on your own!

Feel free to get in touch with me over the forum or directly via email if you have any questions or issues, I would be very glad to help you get this running through tax-credit, and am personally interested in the optimal settings for COI.

Good luck!