Glimmer Parse Error

[nhartwic]

Currently operating on funannotate version 1.8.2 and using docker/singularity to isolate my install. Link below...

docker://nolanhartwick/salk-containers:fun

...I'm running into an intermittent bug with glimmer. Some assemblies hit this error, some don't with no apparent rhyme or reason. Glimmer version is 3.0.4. Stack Trace below...

GlimmerHMM parsing error: 1     GlimmerHMM      CDS     643381  643647  .       +       0       ID=1.cds15##sequence-region 2 1 25532557
Traceback (most recent call last):
  File "/opt/conda/envs/myenv/bin/funannotate", line 706, in <module>
    main()
  File "/opt/conda/envs/myenv/bin/funannotate", line 696, in main
    mod.main(arguments)
  File "/opt/conda/envs/myenv/lib/python3.7/site-packages/funannotate/predict.py", line 1553, in main
    args.out, 'predict_misc'), GlimmerPredictions)
  File "/opt/conda/envs/myenv/lib/python3.7/site-packages/funannotate/library.py", line 6002, in runGlimmerHMM
    glimmer2gff3(glimmerRaw, output)
  File "/opt/conda/envs/myenv/lib/python3.7/site-packages/funannotate/library.py", line 5931, in glimmer2gff3
    '\t')
ValueError: not enough values to unpack (expected 9, got 4)

I suppose I'll work on getting the glimmer logs and raw output. If this is a known error of some kind, let me know. Thanks in advance for any assistance

[nhartwic]

Update. Glimmer logs and output files for a failed run: glimmer_results.tar.gz

Seems like issue is on line 151753 of "glimmerHMM.output.raw".

 10      GlimmerHMM      CDS     6900000 6900084 .       +       1       ID=10.cds1103.2;Parent=10.path1.gene1103;Name=10.path1.gene1103;Note=final-exon
 1
 ##sequence-region 11 1 14171033
 11      GlimmerHMM      mRNA    2309    2458    .       +       .       ID=11.path1.gene1;Name=11.path1.gene1
 11      GlimmerHMM      CDS     2309    2458    .       +       0       ID=11.cds1.1;Parent=11.path1.gene1;Name=11.path1.gene1;Note=final-exon

Can someone confirm that glimmhmm.pl/glimmerhmm is responsible for producing this file and not funannotate?

Assuming the above is true, there must be an error in the perl script or in the glimmerhmm executable. I suppose I'll try to recompile glimmerhmm and see if that fixes anything.

[nextgenusfs]

Hi @nhartwic, so is that 1 on its own line? It might just be the ##sequence-region line. Let me look at the code quick.

[nhartwic]

Yes, its on its own line, You can verify by looking at "glimmerHMM.output.raw" after you unpack the tar archive. It may be the sequence region line, though I don't know why that line would get mangled.

[nextgenusfs]

Here is the function generating that output, its indeed the raw output from glimmer.

funannotate/funannotate/library.py

Lines 5994 to 5999 in bf4a348

	# now run GlimmerHMM prediction
	# glimmhmm.pl <glimmerhmm_program> <fasta_file> <train_dir> <options>
	glimmerRaw = os.path.abspath(os.path.join(dir, 'glimmerHMM.output.raw'))
	cmd = ['perl', which_path('glimmhmm.pl'), which_path(
	'glimmerhmm'), os.path.abspath(fasta), os.path.abspath(tmpdir), '-g']
	runSubprocess2(cmd, dir, log, glimmerRaw)

[nextgenusfs]

I can have the parser ignore lines that don't have the proper length/num of tabs. Hopefully that won't cause some other unintended issues? I have not seen this error before, possible its a glimmer compilation issue I guess.

[nextgenusfs]

The parser is here:

funannotate/funannotate/library.py

Lines 5907 to 5968 in bf4a348

	def glimmer2gff3(input, output):
	'''
	scaffold_39 GlimmerHMM mRNA 23692 25015 . + . ID=scaffold_39.path1.gene12;Name=scaffold_39.path1.gene12
	scaffold_39 GlimmerHMM CDS 23692 23886 . + 0 ID=scaffold_39.cds12.1;Parent=scaffold_39.path1.gene12;Name=scaffold_39.path1.gene12;Note=initial-exon
	scaffold_39 GlimmerHMM CDS 24282 24624 . + 0 ID=scaffold_39.cds12.2;Parent=scaffold_39.path1.gene12;Name=scaffold_39.path1.gene12;Note=internal-exon
	scaffold_39 GlimmerHMM CDS 24711 25015 . + 2 ID=scaffold_39.cds12.3;Parent=scaffold_39.path1.gene12;Name=scaffold_39.path1.gene12;Note=final-exon
	scaffold_39 GlimmerHMM mRNA 25874 27899 . - . ID=scaffold_39.path1.gene13;Name=scaffold_39.path1.gene13
	scaffold_39 GlimmerHMM CDS 25874 26973 . - 2 ID=scaffold_39.cds13.1;Parent=scaffold_39.path1.gene13;Name=scaffold_39.path1.gene13;Note=final-exon
	scaffold_39 GlimmerHMM CDS 27257 27899 . - 0 ID=scaffold_39.cds13.2;Parent=scaffold_39.path1.gene13;Name=scaffold_39.path1.gene13;Note=initial-exon
	'''
	with open(output, 'w') as outfile:
	outfile.write(("##gff-version 3\n"))
	exonCounts = {}
	GeneCount = 1
	skipList = []
	idsSeen = {}
	with open(input, 'r') as infile:
	for line in infile:
	if line.startswith('##sequence-region'):
	idsSeen = {}
	if line.startswith('#') or line.startswith('\n'):
	continue
	line = line.strip()
	contig, source, feature, start, end, score, strand, phase, attributes = line.split(
	'\t')
	ID, Parent, Name = (None,)*3
	info = attributes.split(';')
	for x in info:
	if x.startswith('ID='):
	ID = x.replace('ID=', '')
	elif x.startswith('Parent='):
	Parent = x.replace('Parent=', '')
	if Parent and Parent in skipList:
	continue
	if feature == 'mRNA':
	genelen = int(end) - int(start)
	if genelen < 150:
	if not ID in skipList:
	skipList.append(ID)
	continue
	geneID = 'glimmerG_'+str(GeneCount)
	transID = 'glimmerT_'+str(GeneCount)+'-T1'
	idsSeen[ID] = (geneID, transID)
	outfile.write('{:}\t{:}\t{:}\t{:}\t{:}\t{:}\t{:}\t{:}\tID={:};Alias={:};\n'.format(
	contig, source, 'gene', start, end, score, strand, phase, geneID, ID))
	outfile.write('{:}\t{:}\t{:}\t{:}\t{:}\t{:}\t{:}\t{:}\tID={:};Parent={:};Alias={:};\n'.format(
	contig, source, 'mRNA', start, end, '.', strand, '.', transID, geneID, ID))
	GeneCount += 1
	elif feature == 'CDS':
	if Parent in idsSeen:
	geneID, transID = idsSeen.get(Parent)
	if not transID in exonCounts:
	exonCounts[transID] = 1
	else:
	exonCounts[transID] += 1
	num = exonCounts.get(transID)
	outfile.write('{:}\t{:}\t{:}\t{:}\t{:}\t{:}\t{:}\t{:}\tID={:}.exon{:};Parent={:};\n'.format(
	contig, source, 'exon', start, end, '.', strand, '.', transID, num, transID))
	outfile.write('{:}\t{:}\t{:}\t{:}\t{:}\t{:}\t{:}\t{:}\tID={:}.cds;Parent={:};\n'.format(
	contig, source, feature, start, end, score, strand, phase, transID, transID))
	else:
	print(('GlimmerHMM parsing error: {:}'.format(line)))

So would be something like adding if line.count('\t') < 8: continue

[nhartwic]

I would rather figure out why glimmer is outputting junk than try to make funannotate handle junk. This may not be possible though. I'm working on building a new container with a manual glimmer install and will try to test again.

I also am certain that these glimmer errors occur when writing out the models, not when writing the sequence region. There are a handful of errors that look like this when I dug in deeper...

12	GlimmerHMM	mRNA	8008355	8009100	.	+	.	ID=12.path1.gene1019;Name=12.path1.gene1019
12	GlimmerHMM	CDS	8008355	800840
##sequence-region 13 1 13541320
13	GlimmerHMM	mRNA	4280	4816	.	+	.	ID=13.path1.gene1;Name=13.path1.gene1
13	GlimmerHMM	CDS	4280	4816	.	+	0	ID=13.cds1.1;Parent=13.path1.gene1;Name=13.path1.gene1;Note=final-exon

[nextgenusfs]

Yes, good point. I just glanced at the file and noticed what seemed like an extra config name before every ##sequence-region, but from what you have here its more of a problem. Seems like output lines are being truncated somewhat randomly?

[nextgenusfs]

I'll push this code, so it will complete but output the offending lines to stdout:

>>> import funannotate.library as lib
>>> lib.glimmer2gff3('predict_misc/glimmerHMM.output.raw', 'test.glimmer.out')
ERROR parsing GlimmerHMM Raw output in line 151753:
   1
ERROR parsing GlimmerHMM Raw output in line 168669:
   12	GlimmerHMM	CDS	8008355	800840
ERROR parsing GlimmerHMM Raw output in line 273366:
   31	GlimmerHMM	CDS	1182602	1182691	.	+	0	ID=31.cds1
ERROR parsing GlimmerHMM Raw output in line 281298:
   35	GlimmerHMM	CDS	653391	653512	.	-	0	ID=
ERROR parsing GlimmerHMM Raw output in line 286182:
   4
ERROR parsing GlimmerHMM Raw output in line 288006:
   42	GlimmerHMM	CDS	823638	8237

[nhartwic]

Well, good news is your commits seem to be handling Glimmers errors. Bad news is recompiling glimmer didn't accomplish really accomplish anything. I'm getting different errors now, but ultimately they look like the same kind of error, entries are cut off, always immediately prior to a "sequence region" tag.

I also have been looking more closely at the regions that seem to be causing the errors. Perhaps unsurprisingly, these regions have large N gaps. I think these are what is causing the issues with glimmer, though I don't know exactly why. My next update will probably come after I have shrunk these gaps to a standard 500 bp or something.