# Align sequences to HMM using HMMer
check_call(("hmmalign", "-o", hmmalign,
            "--informat", "fasta",
            hmm, sequences))

# convert alignment from Stockholm to CLUSTAL format
# and clip to just the aligned subsequences
from ClustalTools import MultipleAlignment
hmmaln = ".".join((jobname,"aln"))
alignment = MultipleAlignment.fromStockholm(open(hmmalign))
rf = alignment.colAnnotations["RF"]
start = rf.find("x")
assert(start > -1)
stop = rf.rfind("x")+1
motifaln = alignment[start:stop]
motifaln.writeClustal(open(hmmaln,"w"))

# Generate bootstrapped NJ tree using CLUSTALW
check_call(("clustalw",
            "-infile=%s" % hmmaln,
            "-type=PROTEIN",
            "-bootstrap"))