的Perl：從列表

Question

尋找具有最大寬度的區域我有具有以下結構

gene transcript exon length 
A  NM_1   1  10 
A  NM_1   2  5 
A  NM_1   3  20 
A  NM_2   1  10 
A  NM_2   2  5 
A  NM_2   3  50 
B  NM_5   1  10 
...  ...   ...  ... 

所以基本上，該表包括與所有人類基因列的表。第二列包含抄本名稱。同一個基因可以有多個轉錄本。第三列包含一個外顯子編號。每個基因由多個外顯子組成。第四列包含每個外顯子的長度。

現在我想創建一個新的表看起來像這樣：

gene transcript length 
A  NM_2   65 
B  NM_5   10 
... ...   ... 

所以我基本上想要做的就是找到每個基因的轉錄時間最長。 這意味着當每個基因（列基因）有多個轉錄本（列轉錄本）時，我需要對該基因轉錄本的所有外顯子的長度列中的值進行求和。

所以在這個例子中有兩個基因A的轉錄本：NM_1和NM_2。每個都有三個外顯子。 NM_1 = 10 + 5 + 20 = 35這三個值的總和，NM_2是10 + 5 + 50 = 65。因此對於基因A，NM_2是最長的轉錄本，所以我想把它放在新表格中。對於基因B，只有一個轉錄本，一個長度爲10的外顯子。因此，在新表格中，我只想報告這個轉錄本的長度。

我和哈希工作過，所以我想存儲「基因」和「成績單」是兩個不同的密鑰：

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

open(my $test,'<',"test.txt") || die ("Could not open file $!"); 
open(my $output, '+>', "output.txt") || die ("Can't write new file: $!"); 

# skip the header of $test # I know how to do this 

my %hash =(); 
while(<$test>){ 
    chomp; 
    my @cols = split(/\t/); 
    my $keyfield = $cols[0]; #gene name 
    my $keyfield2 = $cols[1]; # transcript name 
    push @{ $hash{$keyfield} }, $keyfield2; 

...

Answer 1

考慮你想要什麼這樣做，我就會想是這樣的：

use strict; 
use warnings; 

my %genes; 

my $header_line = <DATA>; 

#read the data 
while (<DATA>) { 
    my ($gene, $transcript, $exon, $length) = split; 
    $genes{$gene}{$transcript} += $length; 
} 

print join("\t", "gene", "transcript", "length_sum"), "\n"; 

foreach my $gene (keys %genes) { 
    #sort by length_sum, and 'pop' the top of the list. 
    my ($longest_transcript) = 
     (sort { $genes{$gene}{$b} <=> $genes{$gene}{$a} or $a cmp $b } 
      keys %{ $genes{$gene} }); 
    print join("\t", 
     $gene, $longest_transcript, $genes{$gene}{$longest_transcript}), 
     "\n"; 
} 


__DATA__ 
gene transcript exon length 
A  NM_1   1  10 
A  NM_1   2  5 
A  NM_1   3  20 
A  NM_2   1  10 
A  NM_2   2  5 
A  NM_2   3  50 
B  NM_5   1  10 

輸出

gene transcript length_sum 
B NM_5 10 
A NM_2 65 

Answer 2

從List::UtilsBy使用nmax_by（數字最大值）功能可以大大減少不必要的操作。該程序將總長度累加到散列中，然後使用nmax_by爲每個基因挑選最長的轉錄本。

我推測你可以在$fh上打開輸入文件，而不是使用DATA句柄？或者你可以通過命令行將路徑傳遞給輸入文件，只需使用<>而不是<$fh>而不顯式打開任何內容。

use strict; 
use warnings; 

use List::UtilsBy qw/ nmax_by /; 

my $fh = \*DATA; 

<$fh>; # Drop header line 

my %genes; 

while (<$fh>) { 
    my ($gene, $trans, $exon, $len) = split; 
    $genes{$gene}{$trans} += $len; 
} 

my $fmt = "%-7s%-14s%-s\n"; 
printf $fmt, qw/ gene transcript length /; 
for my $gene (sort keys %genes) { 
    my $trans = nmax_by { $genes{$gene}{$_} } keys %{ $genes{$gene} }; 
    printf ' '.$fmt, $gene, $trans, $genes{$gene}{$trans}; 
} 


__DATA__ 
gene transcript exon length 
A  NM_1   1  10 
A  NM_1   2  5 
A  NM_1   3  20 
A  NM_2   1  10 
A  NM_2   2  5 
A  NM_2   3  50 
B  NM_5   1  10 

輸出

gene transcript length 
A  NM_2   65 
B  NM_5   10 

---

更新

這裏是nmax_by大大縮短的版本，將工作爲您測試。您可以在程序的頂部添加此，或者如果你寧願把它放在結尾，那麼你需要與sub nmax_by(&@);在頂部，因爲它有一個原型

sub nmax_by(&@) { 
    my $code = shift; 
    my ($max, $maxval); 
    for (@_) { 
    my $val = $code->($_); 
    ($max, $maxval) = ($_, $val) unless defined $maxval and $maxval >= $val; 
    } 
    $max; 
}