Fix "Invalid alphabet type" error in nhmmer - iss#271

[traviswheeler]

This addresses https://github.com/EddyRivasLab/hmmer/issues/271.

In iss#271, running the command

% nhmmer --cpu 0 --dna query.fasta target.fasta

produces the following error:

Parse failed (sequence file assembly.fasta):
Line 2: unexpected char A; expected FASTA to start with >

The error arises when two things are true of the first sequence in the target file: (1) It does not contain all four nucleotide characters (e.g. it contains only A,C, and Ts, and has no Gs) (2) It is longer than the buffer length used when reading a sequence in nhmmer (4000)

What's happening is: (1) nhmmer calls esl_sqfile_GuessAlphabet() for the target sequence, even if the --dna (or --rna) flag is set (around line 704 of nhmmer.c).

• This is bad. It should respect the --dna/--rna flag.
• That's fixed in this commit.
• But supposing --dna/--rna weren't used, we would like this to be done (and done correctly), so let's dig deeper into what's happening: (2) That function calls esl_sq_GuessAlphabet(), which returns the status eslENOALPHABET, causing esl_sqfile_GuessAlphabet() to return, and skipping the code that usually resets the buffer pointer. That failure to reset the buffer pointer shouldn't be a problem (the function is returning an error, after all!), but it explains the error message we'll see later. (3) esl_sqfile_GuessAlphabet() immediately returns to nhmmer, returning the eslENOALPHABET status. (4) nhmmer fails to check the returned status, and continues processing. (i) It should respond to the status by alerting the user with an appropriate error message. It doesn't do that, which is bad. (ii) Because the buffer pointer wasn't reset, and processing continues, the next ReadWindow call in nhmmer starts reading from the buffer at the 4001st character of the target sequence. When it tries to read the header starting at the buffer position, it sees some character other than '>', and throws an error, with a confusing message. This is solved by capturing the esl_sqfile_GuessAlphabet() status and throwing an error if appropriate (duh). That now happens in this commit.

So overall, the commit (i) obeys the --dna/--rna flag, so that it only bothers to call esl_sqfile_GuessAlphabet() on the target if it hasn't been told. (ii) checks the status returned by esl_sqfile_GuessAlphabet(), and reacts appropriately.

[ptrebert]

@traviswheeler can you estimate when this will be merged? Sorry for being annoying about this... :-)

[traviswheeler]

@ptrebert: it's out of my hands now, and will depend on @npcarter or @cryptogenomicon (or someone else in the group?) finding time to validate the PR. In the meantime, you can clone the nhmmer-invalid-alph branch from my repo, and use it until the PR has been merged.

[ptrebert]

ok, thanks for the info