# Produce Multiple Sequence Alignments using Profile HMMs

I have a few multiple sequence alignments and a few groups (let's say N) of unaligned sequences. I would like to learn the emission and transition probabilities from the multiple sequence alignments, to build with these distributions a profile HMM and then to produce with the Viterbi algorithm N multiple sequence alignments from the groups of unaligned sequences.

I have a known profile HMM structure, I have a Viterbi implementation but I am having trouble figuring out how to estimate emission distributions and transition distributions from the available multiple sequence alignments. Any ideas?

At the risk of plugging my own package you can derive a profile HMM using the R package aphid. It isn't on CRAN yet but you can download it from github using the devtools package:

devtools::install_github("shaunpwilkinson/aphid")
library(aphid)

I recommend importing your sequences as a DNAbin object first using the ape package if possible, since aphid will recognize the residues as nucleotides and treat the ambiguities appropriately. Once you have your DNAbin object (say x) which can either be a matrix of aligned sequences or a list of non-aligned sequences, you can derive a profile HMM by running

phmm <- derivePHMM(x)

From there, your transition and emission probabilities can be returned by calling phmm$A and phmm$E, respectively.

Otherwise if you don't want to use R you could try either the HMMER (www.hmmer.org/) or SAM (https://compbio.soe.ucsc.edu/sam.html) package.

• Maybe it was not clear from the question but I'm rather interested in a theoretical answer than a concrete implementation. As soon as I know what exactly I am searching I can decide upon using tools or implementing myself. Thanks, though. :) – Eszter May 24 '17 at 15:55
• You can't go past this book: Durbin R, Eddy SR, Krogh A, Mitchison G (1998) Biological sequence analysis: probabilistic models of proteins and nucleic acids. Cambridge University Press, Cambridge, United Kingdom. Chapter 5.3 is all about deriving profile HMMs from multiple sequence alignments. – Shaun Wilkinson May 25 '17 at 3:32