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