Поиск файла FASTA для мотива и возврат строки заголовка для каждой последовательности, содержащей мотив

Question

Ниже приведен код для поиска файла FASTA, введенного в командной строке для предоставленного пользователем мотива. Когда я запускаю его и ввожу мотив, который, как мне известно, находится в файле, он возвращает «Мотив не найден». Я только начинающий в Perl, и я не могу понять, как заставить его печатать найденный мотив, не говоря уже о том, чтобы вернуть строку заголовка. Буду признателен за любую помощь в решении этой проблемы.

Спасибо.

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;

Answer 1

Нет смысла «катать свой собственный» парсер Fasta.BioPerl потратил годы на его разработку, и было бы глупо не использовать его.

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");
  }
}

Эта программа находит только (основанную на 1) позицию первого совпадения мотива в каждой последовательности.Это может быть легко отредактировано, чтобы найти положение каждого матча.Он также может не печатать вещи точно в том формате, который вы хотите / нужен.Я оставлю эти вопросы как «упражнение для читателя».:)

Если вам нужно скачать BioPerl, попробуйте по этой ссылке .Дайте мне знать, если у вас есть какие-либо проблемы.

Для подобных вопросов по биоинформатике я нашел форум BioStar очень полезным.

Answer 2

Вы пытаетесь прочитать имя файла, а не дескриптор файла.

Заменить

@seq = split(/[>]/, $filename);

от

@seq = <DNAFILE>

(или разделить его, если вам нужно - я не знаю, что должен делать ваш разделитель / [>] /: нет смысла помещать один символ в []).