Searching FASTA file for motif and returning title line for each sequence containing the motif
Below is the code I have for searching a FASTA file entered at the command line for a user-provided motif. When I run it and enter a motif that I know is in the file it returns 'Motif not found'. I'm only a beginner in Perl, and I can't fugure out how to get it to print motif found, let alone return the title line. I would appreciate any help in resolving this.
Thanks.
use warnings;
use strict;
my $motif;
my $filename;
my @seq;
#my $m开发者_运维百科otif_found;
my $scalar;
$filename = $ARGV[0];
open (DNAFILE,$filename) || die "Cannot open file\n";
@seq = split(/[>]/, $filename);
print "Enter a motif to search for; ";
$motif = <STDIN>;
chomp $motif;
foreach $scalar(@seq) {
if ($scalar =~ m/$motif/ig) {
print "Motif found in following sequences\n";
print $scalar;
} else {
print "Motif was not found\n";
}
}
close DNAFILE;
There is no point in "rolling your own" Fasta parser. BioPerl has spent years developing one, and it would be silly not to use it.
use strict;
use Bio::SeqIO;
my $usage = "perl dnamotif.pl <fasta file> <motif>";
my $fasta_filename = shift(@ARGV) or die("Usage: $usage $!");
my $motif = shift(@ARGV) or die("Usage: $usage $!");
my $fasta_parser = Bio::SeqIO->new(-file => $fasta_filename, -format => 'Fasta');
while(my $seq_obj = $fasta_parser->next_seq())
{
printf("Searching sequence '%s'...", $seq_obj->id);
if((my $pos = index($seq_obj->seq(), $motif)) != -1)
{
printf("motif found at position %d!\n", $pos + 1);
}
else
{
printf("motif not found.\n");
}
}
This program only finds the (1-based) position of the first motif match in each sequence. It can easily be edited to find the position of each match. It also may not print things exactly in the format you want/need. I'll leave these issues as "an exercise for the reader." :)
If you need to download BioPerl, try this link. Let me know if you have any issues.
For bioinformatics questions like this, I've found the BioStar forum very helpful.
You're trying to read from the filename, not the file handle.
Replace
@seq = split(/[>]/, $filename);
by
@seq = <DNAFILE>
(or split it if you need to - I don't know what your split /[>]/ is supposed to be doing: there's no point in putting a single character in []).
精彩评论