perl:从列表中查找宽度最大的区域
我有一张具有以下结构的桌子
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。所以在新的表格中，我只想知道这个转录本的长度。
我以前使用过散列，所以我想把“gene”和“transcript”作为两个不同的键存储：
#! /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;..。
回答 2
Stack Overflow用户
发布于 2015-04-25 19:08:41
考虑到你想做的事，我会这样想：
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 65Stack Overflow用户
发布于 2015-04-25 19:31:38
使用来自nmax_by的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;
}https://stackoverflow.com/questions/29868889
复制相似问题
腾讯云开发者
