Searching FASTA file for motif and returning title line for each sequence containing the motif

Question

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 $motif_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;

Please don't write comments that describe what a single line of code does. They don't add anything. — Fred Foo, Dec 01 '10 at 13:50
Does this handle Fasta files where the DNA sequence is split over multiple lines? I would recommend using BioPerl's Fasta parser. — Daniel Standage, Dec 01 '10 at 15:29

score 4 · Answer 1 · edited Feb 02 '22 at 19:46

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.

score 1 · Answer 2 · answered Dec 01 '10 at 14:02

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 []).

Searching FASTA file for motif and returning title line for each sequence containing the motif

2 Answers2