时间戳

首页 > 解决方案 > perl- 从多 fasta 文件中提取重复序列

perl - perl- 从多 fasta 文件中提取重复序列

问题描述

我有一个很大的 fasta 文件input.fasta，其中包含许多重复的序列。我想输入一个标题名称并提取出所有具有匹配标题的序列。我知道这可以用 awk/sed/grep 轻松完成，但我需要一个 Perl 代码。

输入法

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

我已经尝试使用以下脚本，但它没有提供任何输出。

脚本.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;

请帮我纠正这个脚本。谢谢！

标签： perl

解决方案

如果你使用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

推荐阅读