【问题标题】:need to optimize perl script需要优化perl脚本
【发布时间】:2015-10-12 23:51:08
【问题描述】:

我有两个文件,每个大约 700,000 行,看起来像:

文件1

1 rs58108140 0 10583 A G
1 rs189107123 0 10611 G C
1 rs180734498 0 13302 T C
1 rs144762171 0 13327 C G
22 rs1001008 0 44928369 G A

文件2

hg19chrc snpid a1 a2 bp info or se p ngt
1 rs4951859 C G 729679 0.631 0.97853 0.0173 0.2083 0
1 rs142557973 T C 731718 0.665 1.01949 0.0198 0.3298 0
22 rs1001008 A G 44928369 0.969 0.98649 0.0107 0.2023 0

我有一个脚本,我在其中查找文件 1 中的字段 1 和 4 以及文件 2 中的字段 1 和 5 匹配的行。然后,如果文件 1 中的字段 5 和 6 与文件 2 中的字段 3 和 4 匹配,我只是打印出文件 2 中的行。但是,如果字段 3 和 4(文件 2)相对于字段 5 和 6 ( File1),我从 File2 中取字段 7 的倒数,并打印带有调整后字段 7 的行:

#! perl -w                                                                                                           
use strict;
use warnings;

my @kgloci;
open( my $loci_in, "<", "File1" ) or die $!;
while (<$loci_in>) {
    my ($chr, $snpID, $dist, $bp, $A1, $A2) = split;
    next if m/Chromosome/;
    push @kgloci, [$chr, $snpID, $dist, $bp, $A1, $A2];
}
close $loci_in;

my $filename = shift @ARGV;
open( my $input, "<", "File2" ) or die $!;
while (<$input>) {
    next if m/hg19chrc/;
    my ($chr, $snpID, $A1, $A2, $bp, $info, $or, $se, $p, $ngt) = split;
    foreach my $kglocus (@kgloci) {
        if (    $chr == $kglocus->[0]
                and $bp == $kglocus->[3]
                and $A1 eq $kglocus->[4] ){
        print "$chr $snpID $A1 $A2 $bp $info $or $se $p $ngt\n";
            next;
        }
        elsif ( $chr == $kglocus->[0]
            and $bp == $kglocus->[3]
            and $A1 eq $kglocus->[5]){
        my $new_or = 1/$or;
        print "$chr $snpID $A1 $A2 $bp $info $new_or $se $p $ngt\n";  
        next;
        }
    }
}
close($input);

照原样,脚本将运行数天。有人能指出提高效率的方法吗?

【问题讨论】:

  • 将文件导入数据库?

标签: perl optimization


【解决方案1】:

使用哈希而不是数组。

以下脚本的输出应该与您的相同(除非 File1 中的 A1 和 A2 相同):

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

my %lookup;
open my $LOCI, '<', 'File1' or die $!;
while (<$LOCI>) {
    next if /Chromosome/;
    my ($chr, $snpID, $dist, $bp, $A1, $A2) = split;
    $lookup{"$chr:$bp:$A1"} = 1;
    $lookup{"$chr:$bp:$A2"} = 2;
}

open my $IN, '<', 'File2' or die $!;
while (<$IN>) {
    next if m/hg19chrc/;
    my ($chr, $snpID, $A1, $A2, $bp, $info, $or, $se, $p, $ngt) = split;
    my $new_or = ( sub {},
                   sub { shift },
                   sub { 1 / shift },
                 )[ $lookup{"$chr:$bp:$A1"} || 0 ]->($or);
    print "$chr $snpID $A1 $A2 $bp $info $new_or $se $p $ngt\n"
        if defined $new_or;
}

这就是我创建测试数据的方式:

perl -MList::Util=shuffle -wE '
    $i = 1;
    $pos = 1;
    for (1 .. 10000) {
        say join " ", $i, "rs$pos", int rand 5, 1000 + int rand 10000,
            (shuffle(qw(A C T G)))[0,1];
        ++$i if rand 1 > .99;
        $pos += int rand 20;
    }' > File1

perl -wE '
    $i = 1;
    $pos = 1;
    for (1 .. 10000) {
        say join " ", $i, "rs$pos", (map qw(A C T G)[rand 4],1,2),
            1000 + int rand 10000, rand 1, rand 5, rand 1, rand 1,
            int rand 3;
        ++$i if rand 1 > .99;
        $pos += int rand 20;
    }' > File2

结果:旧脚本耗时 16 秒,新脚本不到 0.1 秒。对于较大的文件(700_000 行),新脚本需要 4 秒。

【讨论】:

  • 这个脚本很快,谢谢。但是我确实注意到,在某些情况下,您的脚本为字段 7 输出了不正确的值。例如,其中 File1 为:22 rs1296744 0 17968155 A G,File2 为 22 rs1296744 A G 17968155 0.956 0.9977 0.0108 0.831 0,输出为 22 rs1296744 A G 17968155 0.956 1.00230530219505 0.0108 0.831 0,其中字段 7 实际上应保持为 0.9977。奇怪的是,我的脚本为这一行输出了两行,一行带有正确的字段 7,另一行带有调整后的字段 7
  • @theo4786:可能还有另一行影响输出。只有这两行,我从两个脚本中得到了相同的结果。
  • 是的,你是对的,对于这些实例,有两个输入行包含 File1 中的重复字段 1 和 4,但其他字段中的值不同。
  • @theo4786:你能显示这样的行并指出正确的输出吗?
  • 例如,File1 中的这些行:22 rs1296744 0 17968155 A G22 rs201055523 0 17968155 AGG A。正确的输出是22 rs1296744 A G 17968155 0.956 0.9977 0.0108 0.831 0。我相信我可以通过$lookup{"$chr:$snpID:$bp:$A1"} 实现这一点,因为我想要 $snpID 也匹配的输出。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 2012-01-06
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2015-07-17
  • 1970-01-01
相关资源
最近更新 更多