2012-07-26 23 views
2

我正在解析來自psiblast的輸出報告。我使用COG比對和搜索基因數據庫中的匹配(同系物)。我想要做的一件事是找出哪些基因與多個COG相匹配。我的部分腳本如下。錯誤將值推入陣列的散列

我在創建一個數組時保留了所有分配給多個COG的基因的所有COG的問題。

我收到以下錯誤:「不能使用字符串(」COG0003「)作爲ARRAY ref,而」parse_POG_reports.pl 26行67「使用」strict refs「。

我看過其他張貼有關將元素推入數組散列的內容。但我認爲當一個基因與同一個COG有2個匹配時,可能會發生錯誤,並且它會嘗試將相同的COG推入數組(即樣本輸入的最後2行)。這有意義嗎?如果是這樣,我該如何避免這個問題?

use strict; 
use warnings; 

my %maxBits;my %COGhit_count; 
my $Hohits={};my %COGhits; 

my $COG_psi_report=$ARGV[0]; 
open (IN, $COG_psi_report) or die "cannot open $COG_psi_report\n"; 
while (my $line=<IN>){ 
    next if ($line =~/^#/); 
    chomp $line; 
    my @columns = split(/\t/,$line); 
    my $bits=$columns[11]; 
    my $COG=$columns[0]; 
    my $hit=$columns[1]; 
    my $Eval=$columns[10]; 
    next if ($Eval > 0.00001); # threshold for significant hits set by DK 
    $COGhit_count{$hit}++; # count how many COGs each gene is homologous to 
    $COGhits{$hit}=$COG; 
    if ($COGhit_count{$hit}>1) { 
      push @{$COGhits{$hit}}, $COG; # 
    } 
    ## for those that there are multiple hits we need to select top hit ## 
    if (!exists $maxBits{$hit}){ 
      $maxBits{$hit}=$bits; 
    } 
    elsif (exists $maxBits{$hit} && $bits > $maxBits{$hit}){ 
      $maxBits{$hit}=$bits; 
    } 
    $Hohits->{$hit}->{$bits}=$COG; 
} 
close (IN); 

例如輸入:

POG0002 764184357-stool1_revised_scaffold22981_1_gene47608  23.90 159  112  3  1  156  1  153  2e-06 54.2 
POG0002 764062976-stool2_revised_C999233_1_gene54902 23.63 182  121  5  3  169  2  180  2e-06 53.9 
POG0002 763901136-stool1_revised_scaffold39447_1_gene145241  26.45 155  89  3  3  137  5  154  3e-06 53.9 
POG0002 765701615-stool1_revised_C1349270_1_gene168522 23.53 187  115  5  3  169  2  180  5e-06 53.1 
POG0002 158802708-stool2_revised_C1077267_1_gene26470 22.69 216  158  5  3  213  5  216  5e-06 52.7 
POG0003 160502038-stool1_revised_scaffold47906_2_gene161164  33.00 297  154  6  169  424  334  626  6e-40 157 
POG0003 160502038-stool1_revised_scaffold47906_2_gene161164  16.28 172  128  4  23  192  46  203  1e-06 56.6 
POG0003 158337416-stool1_revised_C1254444_1_gene13533 30.06 346  184  7  133  424  57  398  6e-40 155 
POG0003 158337416-stool1_revised_scaffold29713_1_gene153054  28.61 332  194  8  132  424  272  599  2e-38 152 
POG0003 158337416-stool1_revised_scaffold29713_1_gene153054  24.00 200  131  5  1  193  5  190  9e-11 69.3 
+0

謝謝 - 我真的不知道這個事情接受直到最近,會做 – user1249760 2012-07-26 14:53:09

回答

1

你需要擺脫線24(倒數):

$COGhits{$hit}=$COG; 

在這裏面,你設置$COGhits{$hit}一個標值(值$COG)。後來,在第26行中,您正嘗試將$COGhits{$hit}解除引用作爲數組。這不起作用,因爲那裏有一個標量。

只要刪除if並將這些行更改爲此。現在應該這樣做,因爲所有$hit都存儲在數組引用中。的$COGhits

$COGhit_count{$hit}++; # count how many COGs each gene is homologous to 
push @{$COGhits{$hit}}, $COG; 

輸出:

$VAR4 = { 
     '158802708-stool2_revised_C1077267_1_gene26470' => [ 
                  'POG0002' 
                 ], 
     '764062976-stool2_revised_C999233_1_gene54902' => [ 
                  'POG0002' 
                 ], 
     '764184357-stool1_revised_scaffold22981_1_gene47608' => [ 
                   'POG0002' 
                   ], 
     '765701615-stool1_revised_C1349270_1_gene168522' => [ 
                  'POG0002' 
                  ], 
     '763901136-stool1_revised_scaffold39447_1_gene145241' => [ 
                   'POG0002' 
                   ], 
     '160502038-stool1_revised_scaffold47906_2_gene161164' => [ 
                   'POG0003', 
                   'POG0003' 
                   ] 
    }; 

如果你想但同時標量和數組引用,試試這個代碼。 雖然我不推薦這個

$COGhit_count{$hit}++; # count how many COGs each gene is homologous to 
if ($COGhit_count{$hit} == 1) { 
    $COGhits{$hit}=$COG;    # Save as scalar 
} 
elsif ($COGhit_count{$hit} == 2) { # If we've just found the second hit, 
    my $temp = $COGhits{$hit};  # save the first and convert $COGhits{$hit} 
    $COGhits{$hit} = [];    # to an array ref, then push both the old and 
    push @{$COGhits{$hit}}, $temp, $COG; # the new value in it. 
} elsif ($COGhit_count{$hit} > 2) { 
    push @{$COGhits{$hit}}, $COG; # Just push the new value in 
} 

思想:你可能有$COGhits{$hit}=$COG第一,但隨後注意到,有時可以有多個值,讓你添加的push行,但你沒有意識到,你實際上不得不更換老線。

0

它告訴你究竟是你做錯了什麼。

$COGhits{$hit}=$COG; # <--- scalar 
if ($COGhit_count{$hit}>1) { 
     push @{$COGhits{$hit}}, $COG; # <--- array 
} 

您不能分配的值作爲非引用類型,然後嘗試autovivify它作爲參考類型。 Perl會執行後者,但如果您已經在該位置存儲了衝突的數據類型,則不會。另外,如果這是由一些奇蹟第一次工作(它不會),並且你不止一次運行這個操作,那麼任何你可能通過推送自動生成的數組都會被標量非參考分配。

我不知道你在做什麼,但第一行應該是可能被刪除。


而不是構建的,你要決定是否將永遠是$COG不止一個規範的$hit值。如果可能的話,簡單地用push替換那4條線是要走的路。

我之前做過多用途結構插槽,他們在很大程度上是一個難以維護的問題。但是,如果你想要做這樣的事情,你可以這樣做:

my $ref = \$hashref->{ $key }; # autovivifies slot as simple scalar. 
           # it starts out as undefined. 
if (ref $$ref) {    # ref $ref will always be true 
    push @$$ref, $value; 
} 
else { 
    $$ref = defined($$ref) ? [ $$ref, $value ] : $value; 
} 

但是你必須要訪問一些不同的方式混合樹每次寫分支邏輯。用標量得到的性能節省,有點被測試和分支所吞噬。

所以我不再做太多了。我事先決定關係是1-1還是1-n。像下面這樣的例程可以在一定程度上使它更直接地處理這些類型的表。

sub get_list_from_hash { 
    my ($hash, $key) = @_; 
    my $ref = \$hash->{ $key }; 
    return unless defined($$ref); 
    return ref($$ref) ? @$$ref : $$ref; 
} 

sub store_in_hash { 
    $_[0] = {} unless ref $_[0]; 
    my ($hash, $key, @values) = @_; 
    my @defined = grep {; defined } @values; 
    unless (@defined) { 
     delete $hash->{ $key }; 
     return; 
    } 

    my $ref = \$hash->{ $key }; 
    if (ref $$ref) { 
     push @$$ref, @defined; 
    } 
    elsif (defined $$ref) { 
     $$ref = [ $$ref, @defined ]; 
    } 
    elsif (@values > 1) { 
     @$$ref = @defined; 
    } 
    else { 
     ($$ref) = @defined; 
    } 
} 
+0

我是這麼認爲的,但我以前試過$ COGhits {$命中} [0] = $ COG;首先,然後按照simbabque的建議刪除第24行(這會導致空數組,例如ARARAY(0xaf9b58)) – user1249760 2012-07-26 14:54:39

+1

除非命中映射到COG爲1-1,否則刪除第24行是沒有意義的。刪除賦值是有意義的,特別是如果你期望'$ COGhit_count {$ hit}'的值永遠大於1.但分配給數組的第一個元素沒有任何意義。 – Axeman 2012-07-26 15:08:08

+0

當然你對這個映射是正確的。沒有想到這一點。 – simbabque 2012-07-26 15:14:52