Inconsistent reslults with dbCAN

[caonetto]

I am getting some inconsistent results on the annotation of cazymes. Ive manually annotated using the dbCAN server a protein that corresponds to the GH31 group, and the e-value of the hit is 2.2e-122. Ive noticed from the source code of funannotate that the cut-off during hmmscan is 1e-17, however this protein is not getting annotated as part of GH31.
I have manually annotated the protein using the hmmer database that is downloaded for funannotate and the hit is actually there however is not getting passed through to the final output.
This is the score output when running hmmscan manually with the same commands as funannotate,
GH31 - 427 wild87red_001881-T1 - 987 2.3e-151 499.9 3.6 1 2 0.1 67 -2.5 0.0 50 89 237 275 220 287 0.75 -
GH31 - 427 wild87red_001881-T1 - 987 2.3e-151 499.9 3.6 2 2 3.6e-154 2.3e-151 499.9 3.6 3 427 418 859 416 859 0.97 -

The txt file "annotations.dbCAN.txt" located in the annotate_misc folder does not contain this hit.

Any ideas on what could be going on?

Thanks.

[nextgenusfs]

Originally with dbCAN there was post evalue filtering, which was if coverage of the hit was less than 45% than to drop it (Note I wrote this probably 5 years ago, so completely plausible that dbCAN post filtering should be re-evaluated). The code running hmmer-dbCAN is here:

funannotate/funannotate/aux_scripts/hmmer_parallel.py

Lines 104 to 161 in fea594c

	def dbCANsearch(inputList, cpus, evalue, tmpdir, output):
	# run hmmerscan
	dbCAN_out = os.path.join(tmpdir, 'dbCAN.txt')
	dbCAN_filtered = os.path.join(tmpdir, 'dbCAN.filtered.txt')
	lib.runMultiProgress(safe_run2, inputList, cpus, progress=False)
	# now grab results
	resultList = [os.path.join(tmpdir, f) for f in os.listdir(
	tmpdir) if os.path.isfile(os.path.join(tmpdir, f)) and f.endswith('.dbcan.txt')]
	combineHmmerOutputs(resultList, dbCAN_out)
	# now parse results
	Results = {}
	with open(dbCAN_filtered, 'w') as filtered:
	filtered.write(
	"#HMM_family\tHMM_len\tQuery_ID\tQuery_len\tE-value\tHMM_start\tHMM_end\tQuery_start\tQuery_end\tCoverage\n")
	with open(dbCAN_out, 'r') as results:
	for qresult in SearchIO.parse(results, "hmmscan3-domtab"):
	query_length = qresult.seq_len
	hits = qresult.hits
	num_hits = len(hits)
	if num_hits > 0:
	for i in range(0, num_hits):
	hit_evalue = hits[i].evalue
	if hit_evalue > evalue:
	continue
	hit = hits[i].id
	hmmLen = hits[i].seq_len
	hmm_aln = int(hits[i].hsps[0].hit_end) - \
	int(hits[i].hsps[0].hit_start)
	coverage = hmm_aln / float(hmmLen)
	if coverage < 0.45:
	continue
	query = hits[i].query_id
	filtered.write("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%f\n" % (
	hit, hmmLen, query, query_length, hit_evalue,
	hits[i].hsps[0].hit_start,
	hits[i].hsps[0].hit_end,
	hits[i].hsps[0].query_start,
	hits[i].hsps[0].query_end,
	coverage))
	if query not in Results:
	Results[query] = [hit]
	else:
	Results[query].append(hit)
	# run through results and simplify subdomain hits
	with open(output, 'w') as out:
	for k, v in natsorted(Results.items()):
	simplified = []
	for x in v:
	if '_' in x:
	cazy, subdomain = x.rsplit('_', 1)
	if cazy not in simplified:
	simplified.append(cazy)
	else:
	if not x in simplified:
	simplified.append(x)
	for hit in simplified:
	out.write("{}\tnote\tCAZy:{}\n".format(k, hit))

Okay so I had a look here: https://bcb.unl.edu/dbCAN2/help.php. Seems that we should now be running evalue of 1e-15 and filtering at 35% coverage. I will make the changes.

[hyphaltip]

this is fixed with commit 036c9d6