使用id文件的子集fasta文件

Question

我被一些剧本卡住了！

嗯，我在2008年做了这个脚本，现在我正在使用一些修改，我得到了错误！

#!/usr/bin/perl -w

use strict;
use Data::Dumper;

sub getSequences ($) {

    my $file = $_[0];
    open (INFILE, "<$file") || die "Error1 in opening in file: $file. $!\n";

        my @lines = <INFILE>;
        my $header; my %header2seq = ();

        foreach my $line (@lines) {
            chomp $line;
            if ($line =~ /^(>.+)$/o) { 
                $header = $1; 
            }
            else {$header2seq {$header} .= $line; }
        }
#print %header2seq;
    close (INFILE);
    return (\%header2seq);
}

sub MakeSpList ($) {

    my $sp_list = $_[0]; my %sp_names = ();
    open (INFILE2, "<$sp_list") || die "Error2 in opening in file: $sp_list. $!\n";
    my @sps = <INFILE2>;
    foreach my $line (@sps) { chomp $line; $sp_names {$line} = 1; }
    close (INFILE2);
#print Dumper (%sp_names);
    return (\%sp_names);
}



 sub CompareSpList2Sequences ($$$) {
     my $ref_header2seq = $_[0] ; my $ref_sp_names = $_[1]; my $file = $_[2];
     open (OUTFILE, ">$file.subdata") || die ("Error3 in opening out file: $file.subdata. $!\n");
     foreach my $key (keys %$ref_header2seq) {
         $key =~ m/^>([A-Z]+[0-9]+[A-Z+]).+$/o;
        #print "$1\n";
         my $header_sub = $1;
         #print $header_sub, "\n";
         #print $ref_sp_names, "\n";
         if (exists $ref_sp_names -> {$header_sub}) {
             my $seq = $ref_header2seq -> {$key};
             print OUTFILE ">$key\n$seq\n";
         }
     }  
     close (OUTFILE);
     return "42";
}


my $fasta_seqs = $ARGV[0]; my $sp_list = $ARGV[1];


my $ref_header2seq = getSequences ($fasta_seqs);
my $ref_sp_names = MakeSpList ($sp_list);
CompareSpList2Sequences ($ref_header2seq , $ref_sp_names, $fasta_seqs);

exit;

我想做的是：

我有一个包含序列的fasta文件：

：YAL004W YAL004W SGDID:S 000002136，Chr i来自140760-141407，基因组释放64-2-1，可疑的ORF，“可疑的开放阅读框架；根据现有的实验和比较序列数据，不太可能编码一个功能蛋白；完全重叠已证实的基因SSA1 1/yal005c”

S000000004，Chr i，141431-139503，基因组释放64-2-1，反向补体，经验证，“参与蛋白质折叠和NLS-定向核转运的ATPase；HSP70家族成员；与Ydj1p形成伴侣复合体；定位于细胞核、细胞质和细胞壁；98%与模拟Ssa2p一致，但这两种蛋白质之间的细微差异为酵母URE3启动子的繁殖和葡萄糖凝聚酶的空泡降解提供了功能特异性；Hsp104p对启动纤维的一般靶向因子”Ssa2p“

我还有另一份有身份证的文件：

YAL005C YAL012W

当与ID的文件匹配时，我想检索序列和所有标题。

，但我发现了一个错误:不要打印任何东西！

你能帮帮我吗？

提前谢谢。

• 我已经搜索了其他方法(我也无法得到结果)，但我真的想知道这个错误！
• 请勿使用生物操作！

Answer 1

好的，所以第45行是：

     if (exists $ref_sp_names -> {$header_sub}) {

您的错误告诉您$header_sub是未定义的。它的背景是：

my $header_sub = $1;

具体如下：

     $key =~ m/^(>[A-Z])\s.+$/o;

所以-这意味着正则表达式不匹配。我在你的样本数据中没有看到任何>，所以它不能匹配它。当匹配失败时，$1是未定义的，因此您的错误。你从你的print $key语句中得到了什么？

我也要指出- .+$很可能是多余的。同样-- o标志--您可能也不想这样做。http://perldoc.perl.org/perlre.html#Modifiers

Answer 2

你试过使用Bioperl吗？这里有一些代码可以让你开始。

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

my $fasta = shift; #this will just push whatever in cli in.
my $seqio_obj = Bio::SeqIO->(-file => $fasta, -format => 'fasta');

while ( my $seq = $seqio_obj->next_seq){
     print $seq->id . ' = ' . $seq->seq() . "\n";
     #in here you can do your fasta handling with the seq obj
}