Perl - извлекать повторяющиеся последовательности из файла мульти-fast

Question

У меня есть большой файл fasta input.fasta , который состоит из множества повторяющихся последовательностей. Я хочу ввести имя заголовка и извлечь все последовательности с соответствующим заголовком. Я знаю, что это легко сделать с помощью awk / sed / grep, но мне нужен Perl-код.

input.fasta

>OGH38127_some_organism
PAAALGFSHLARQEDSALTPKHYTWTAPGEGDVRAPCPVLNTLANHEFLPHNGKNITVDK
AITALGDAMNISPALATTFFTGGLKTNPTPNATWFDLDMLHKHNVLEHDGSLSRRDMHFD
TSNKFDAATFANFLSYFDANATVLGVNETADARARHAYDMSKMNPEFTITSSMLPIMVGE
SVMMMLVWGSVEEPGAQRDYFEYFFRNERLPVELGWTPGETEIGVPVVTAMITAMVAASP
TDVP
>ABC14110_some_different_org_name
WWVAPGPGDSRGPCPGLNTLANHGYLPHDGKGITLSILADAMLDGFNIARSDALLLFTQ
AIRTSPQYPATNSFNLHDLGRDQLNRHNVLEHDASLSRADDFFGSNHIFNETVFDESRAY
AMLANSKIARQINSKAFNPQYKFTSKTEQFSLGEIAAPIIAFGNSTSGEVNRTLVEYFFM
NERLPIELGWKKSEDGIALDDILRVTQMISKAASLITPSALSWTAETLTP
>OGH38127_some_organism
LPWSRPGPGAVRAPCPMLNTLANHGFLPHDGKNISEARTVQALGRALNIEKELSQFLFEK
ALTTNPHTNATTFSLNDLSRHNLLEHDASLSRQDAYFGDNHDFNQTIFDETRSYWPHPVI
DIQAAALSRQARVNTSIAKNPTYNMSELGLDFSYGETAAYILILGDKDFGKVNRSWVEYL
FENERLPVELGWTRHNETITSDDLNTMLEKVVN
.
.
.

Я пробовал использовать следующий скрипт, но он не дает никакого вывода.

script.pl

#!/perl/bin/perl -w
use strict;
use warnings;

print "Enter a fasta header to search for:\n";
my $head = <>;

my $file = "input.fasta";
open (READ, "$file") || die "Cannot open $file: $!.\n";
my %seqs;
my $header;

while (my $line = <READ>){
    chomp $line;
    $line =~ s/^>(.*)\n//;
    if ($line =~ m/$head/){
        $header = $1;
    }
}
close (READ);

open( my $out , ">", "out.fasta" ) or die $!;

my @count_seq = keys %seqs;
foreach (@count_seq){
    print $out $header, "\n";
    print $out $seqs{$header}, "\n";
}

exit;

Пожалуйста, помогите мне исправить этот сценарий. Спасибо!

Answer 1

Если вы используете модуль Bioperl Bio :: SeqIO для обработки файлов fasta, это становится действительно простым:

#!/usr/bin/perl
use warnings;
use strict;
use Bio::SeqIO;

my ($file, $name) = @ARGV;
my $in = Bio::SeqIO->new(-file => $file, -format => "fasta");
my $out = Bio::SeqIO->new(-fh => \*STDOUT, -format => "fasta");

while (my $s = $in->next_seq) {
  $out->write_seq($s) if $s->display_id eq $name;
}

запустить с perl grep_fasta.pl input.fasta OGH38127_some_organism

Answer 2

Нет необходимости хранить последовательности в памяти, вы можете распечатать их непосредственно при чтении файла. Используйте переменную-флаг ($inside в примере), которая сообщает вам, читаете ли вы нужную последовательность или нет.

#! /usr/bin/perl
use warnings;
use strict;

my ($file, $header) = @ARGV;

my $inside;
open my $in, '<', $file or die $!;
while (<$in>) {
    $inside = $1 eq $header if /^>(.*)/;
    print if $inside;
}

Запуск от имени

perl script.pl file.fasta OGH38127_some_organism > output.fasta