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 имеется только 1 транскрипт с одним экзоном длиной 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

Это делается намного менее неопрятно, используя функцию nmax_by (числовое значение максимального) от List::UtilsBy. Эта программа накапливает общую длину в хеше, а затем выбирает самый длинный транскрипт для каждого гена, используя 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 

---

Update

Вот очень сокращенный вариант 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; 
}