Resources chevron_right Protocols chevron_right HMMER
Protocol // Annotation

HMMER

Protein homology search using profile hidden Markov models to detect remote sequence relationships.

Align Your FASTA Sequences

mafft --maxiterate 1000 --localpair --thread 10 AA_seq.fasta > AA_seq.mafft

The output of the alignment is also in FASTA format, but to create an HMMER profile, the program prefers input in Stockholm format. To convert, use the following script:

#!/usr/bin/env perl -w

BEGIN {
    use FindBin;
    use lib $FindBin::Bin;
    push @INC, $FindBin::Bin;
}

my $usage = "Usage: $0 <gapped FASTA alignment file(s)>\n";

my @argv;
while (@ARGV) {
    my $arg = shift;
    if ($arg =~ /^-/) {
        if ($arg eq "-h") { print $usage; exit }
        else { die $usage }
    } else {
        push @argv, $arg;
    }
}
push @argv, "-" unless @argv;

# loop through FASTA files
foreach my $fasta (@argv) {
# read FASTA file
    my %seq;
    my @name;
    my $name;
    open FASTA, "<$fasta" or die "Couldn't open '$fasta': $!";
    while (<FASTA>) {
        if (/^\s*>\s*(\S+)/) {
            $name = $1;
            die "Duplicate name: $name" if defined $seq{$name};
            push @name, $name;
        } else {
            if (/\S/ && !defined $name) {
                warn "Ignoring: $_";
           } else {
                s/\s//g;
                $seq{$name} .= $_;
            }
        }
    }
    close FASTA;

# check all seqs are same length
    my $length;
    my $lname;
    foreach my $name (@name) {
        my $l = length $seq{$name};
        if (defined $length) {
            die "Sequences not all same length ($lname is $length, $name is $l)" unless $length == $l;
        } else {
            $length = length $seq{$name};
            $lname = $name;
        }
    }

# print Stockholm output
    print "# STOCKHOLM 1.0\n";
    foreach my $name (@name) {
        print $name, " ", $seq{$name}, "\n";
    }
    print "//\n";
}

To run this script, create a text file with it:

nano fasta2stockholm.pl

Turn the script into an executable file:

chmod +x fasta2stockholm.pl

Now you are ready to run it:

perl ./fasta2stockholm.pl jukA.mafft > jukA.stock

Creating a HMMER Profile

hmmbuild AA_seq.hmm AA_seq.stock

Run Your HMMER Search

hmmsearch AA_seq.hmm Sequences2search.fasta > Sequences2search.out