hhconsensus seems to be the recommended first step in converting MSA in fasta to a3m that allows to skip reformat.pl and thus would be expected to generate identical output to a workflow that has reformat.pl as a first step.
Current Behavior
hhconsensus -M X generates a3m's that are different from those generated for alignments first converted to a3m with reformat.pl -M X
Steps to Reproduce (for bugs)
Not a minimal example, just an ad hoc alignment of 383 sequences: seq.txt derived from this dataset.
reformat.pl fas a3m seq.txt /dev/stdout -M 50 -v 0 | hhconsensus -maxres 65535 -i stdin -o reformat_hhcons.a3m
# - XX:XX:XX.XXX INFO: Reading HMM / multiple alignment from standard input ...
hhconsensus -maxres 65535 -i seq.txt -o hhcons.a3m -M 50
# - XX:XX:XX.XXX INFO: seq.txt is in A2M, A3M or FASTA format
The a3m file generated by hhconsensus differs from the one that was created by the reformat.pl | hhconsensus pipeline in the distribution of the gaps and consequently in the sequence of the consensus. Most noticeably, hhcons.a3m has a single sequence in the whole alignment that doesn't start with a gap despite the -M 50 argument:
# skip header and consensus sequence, discard deflines and count sequences not starting with a gap
tail -n+4 reformat_hhcons.a3m | grep -v '^>' | grep -vc ^-`
# 305
tail -n+4 hhcons.a3m | grep -v '^>' | grep -vc ^-
# 1
HH-suite Output (for bugs)
No errors, only the above messages
Context
This seems to be a bug, but I might be missing something. This is certainly not an intuitive behavior.
Expected Behavior
hhconsensusseems to be the recommended first step in converting MSA infastatoa3mthat allows to skipreformat.pland thus would be expected to generate identical output to a workflow that hasreformat.plas a first step.Current Behavior
hhconsensus -M Xgeneratesa3m's that are different from those generated for alignments first converted toa3mwithreformat.pl -M XSteps to Reproduce (for bugs)
Not a minimal example, just an ad hoc alignment of 383 sequences: seq.txt derived from this dataset.
The
a3mfile generated byhhconsensusdiffers from the one that was created by thereformat.pl | hhconsensuspipeline in the distribution of the gaps and consequently in the sequence of the consensus. Most noticeably,hhcons.a3mhas a single sequence in the whole alignment that doesn't start with a gap despite the-M 50argument:HH-suite Output (for bugs)
No errors, only the above messages
Context
This seems to be a bug, but I might be missing something. This is certainly not an intuitive behavior.
Your Environment