I have been trying to align a set of amino acid sequences using HMM profiles and the HMMER package. In my set, each sequence consists of two domains, for which two distinct HMM profiles are available. If I try aligning the sequences to one of the two profiles the sequence parts that correspond to that profile are aligned well, but the parts corresponding to the second profile are misaligned.
Apart from the obvious option of aligning a number of these sequences manually and creating a new profile containing both domains, is there a way to either perform an alignment using two profiles instead of one, or somehow combine or merge the two profiles into one, so that both domains are aligned correctly?
1 answer
hmmalign performs global alignments to your sequences, not local alignments. If you're trying to get hmmalign to align a complete sequence (containing domains A and B) to a profile of domain A, then the sequence of domain B will either be left as an overhang or misaligned, because there is nothing for it to align to: i.e. the sequence AAAABBBB does not align globally to model AAAA in the first place, so you've set your problem up wrong.
If you want to get an alignment of domain A subsequences, you have to first extract the matching domain A subsequences (using coordinates from hmmsearch, and the esl-sfetch tool, perhaps), then align those subsequences to your domain A profile.
If you want to get an alignment of complete sequences (including both domains A and B), you have to build a profile of complete sequences including both domains A and B.
For a bit more info, here's a blog post on using esl-sfetch to extract subsequences: http://cryptogenomicon.org/?p=394
Log in to answer this question.