高度嵌套生物信息学处理

Question

我有一个解析BAM的脚本

E00489:44:HNNVYCCXX:3:1101:24890:5616   99      22      16052014        150M
E00489:44:HNNVYCCXX:1:1110:21704:27345  99      22      16052044        150M
E00489:44:HNNVYCCXX:1:1217:2372:69519   163     22      16052044        150M
E00489:44:HNNVYCCXX:3:2123:8948:16779   99      22      16052044        150M
E00489:44:HNNVYCCXX:2:2213:2920:25534   147     22      16052054        146M4S
E00489:44:HNNVYCCXX:2:2206:5020:71717   83      22      16052055        145M5S
E00489:44:HNNVYCCXX:4:2206:12642:40829  99      22      16052056        144M6S

(在输入此脚本之前，BAM文件实际上是通过cut -f1-4,6运行的，因此这只是字段的子集。)

第一列是一种唯一的ID，第二列是位标志。第三和第四部分描述了人类基因组中的染色体和位置。第五个是雪茄串，它主要是脚本所解析的。

我很少使用Perl，但它似乎效率很低。它确实按预期工作，但速度很慢。我想加快它的速度，如果可能的话，也会使它更易读和更容易理解。

#!/bin/perl

#initialize hashes
my %softhash;
my %IDhash;
my $file1Name = $ARGV[0] . 'all_sc_positions.txt';
my $file2Name = $ARGV[0] . 'edge_sc_positions.txt';
open(my $fh_all, '>', $file1Name);
open(my $fh_edge, '>', $file2Name);

#for each line
while (my $line = ) {

    #read in line and parse
    chomp($line);
    my @a = split("\t", $line);
    my $start = $a[3];
    my $cigar = $a[4];
    my @b = split(/[A-Z]/, $cigar);# keeps track of numbers
    my @c = split(/[0-9]*/, $cigar); #keeps track of letters
    my $loc = $start; #distance from start of read
    my $var_start;
    my $var_end;

    #loop over type of alignment in cigar string, build hashes of candidate locations
    for (my $i = 0; $i <= $#c; $i++) {

        #if there is softclipping or an indel, find the location and store in hash
        if ($c[$i] eq "S" || $c[$i] eq "I" || $c[$i] eq "D") {

            #find softclipping location
            if ($i == 1) {

                $var_start = $start - $b[0];
                $var_end = $start;

                if ($c[$i] eq "S") {

                    for (my $pos = $var_start; $pos < $var_end; $pos++) {
                        $softhash{$a[2]}{$pos}++;
                    }

                    $var_start = $var_end - 1;
                    $var_end = $var_start; 
                }

            } else {

                for (my $j = $i-2; $j >= 0; $j--) {# subtract 2 from i, because the first value from c will always be empty
                    $loc = $loc + $b[$j]; 
                }

                $var_start = $loc;
                $var_end = $loc + $b[($i-1)];

                if ($c[$i] eq "S") { 

                      for (my $pos=$var_start; $pos<$var_end; $pos++) {
                           $softhash{$a[2]}{$pos}++;
                      }

                      $var_end = $var_start;
                }
            }

            $IDhash{$a[2]}{$var_start}{$var_end}{$c[$i]}++;         
        }       
    }
}

#write out edge features
foreach my $key_chr (sort(keys %IDhash)) {

        foreach my $key_start (sort { $a <=> $b } (keys %{  $IDhash{$key_chr}  })) {

                foreach my $key_end (sort { $a <=> $b } (keys %{  $IDhash{$key_chr}{$key_start}  })){

                    print $fh_edge "$key_chr\t$key_start\t$key_end\t";
                    my $sum = $IDhash{$key_chr}{$key_start}{$key_end}{I} + $IDhash{$key_chr}{$key_start}{$key_end}{D} + $IDhash{$key_chr}{$key_start}{$key_end}{S};
                    print $fh_edge "$sum,";

                    for my $l ('S','I','D') {

                    if (defined($IDhash{$key_chr}{$key_start}{$key_end}{$l})) {
                         print $fh_edge  "$IDhash{$key_chr}{$key_start}{$key_end}{$l},";
                    } else {
                        print $fh_edge "0,";
                    }
            }   

            print $fh_edge "\n";        
        }
    }
}

#write out "all" features
foreach my $key_chr (sort(keys %softhash)) {

        foreach my $key_pos (sort { $a <=> $b } (keys %{  $softhash{$key_chr}  })) {

             print $fh_all "$key_chr\t$key_pos\t";
             print $fh_all "$softhash{$key_chr}{$key_pos}\n";
        }
}

编辑：吉斯特！

Answer 1

非代码/非常高级别的注意事项.

为了提高性能，你可以做的第一件事就是得到一台更好的电脑。我的计算机已经有几年的历史了，但是它在100000行的示例中运行未经修改的代码不到1秒，而您说它只需一分钟。(当然，仍然值得记住的是，更好的算法放大了更好的硬件的性能优势，所以我也将关注硬件)。

我在生物信息学方面的经验非常有限，而且我以前从未使用过这种特殊的格式。我只看了一下医生。但是，在链接到雪茄烟字符串的描述之后，以及从该页面到https://samtools.github.io/hts-specs/SAMv1.pdf的另一个链接之后，我观察到该语句。

S可能只会在他们和雪茄串的末端之间做H手术。

然而，100000行示例中的1368行违反了该限制.我建议您检查代码，以检查它不依赖于它，因为您是领域专家，并且比我们大多数人都更了解代码的含义。

我认为值得指出的是你的观察

我很少使用Perl，但它似乎效率很低。

这里的密码不多。也许你应该考虑把它移植到一种你知道得更好、效率更高的语言上？(我猜想使用PyPy或Cython运行Python会更快，而且Python知识在您的领域中似乎足够普遍，不会为以后从您那里继承它的人造成问题)。

代码注意事项

可能是虫子？

    my $loc = $start; #distance from start of read
    ...
    for (my $i = 0; $i <= $#c; $i++) {
        if ($c[$i] eq "S" || $c[$i] eq "I" || $c[$i] eq "D") {
            ...
            #find softclipping location
            if ($i == 1) {
                ...
            } else {
                for (my $j = $i-2; $j >= 0; $j--) {# subtract 2 from i, because the first value from c will always be empty
                    $loc = $loc + $b[$j]; 
                }
                ...
            }
            ...
        }       
    }

如果我正确地解释了事情，$loc = $start不应该就在$j循环之前吗？否则，带有几个I/D/S的雪茄串将重复计算一些$b[$j]。

此外，我真的不明白I和D是如何被同等对待的。难道他们中的一个不应该让$var_start倒下吗？

注意：由于这个例子太快，我无法获得很多有用的分析信息，这是基于常识的。速度优化常常违背常识。我建议一个一个地测试这些建议，看看哪些在大型数据集中起作用，哪些不起作用。

最明显的优化涉及到softhash，特别是循环。

for (my $pos = $var\_start; $pos < $var\_end; $pos++) { $softhash{$a[2]}{$pos}++; }

如果范围趋向于适度宽，那么这就做了大量的工作，优化的时机已经成熟。具体来说，您可以用

                    $softhash{$a[2]}{$var_start}++;
                    $softhash{$a[2]}{$var_end}--;

，然后将末尾的输出循环修改为

foreach my $key_chr (sort(keys %softhash)) {
    my $accum = 0;
    my $prev_pos = -1;
    foreach my $key_pos (sort { $a <=> $b } (keys %{  $softhash{$key_chr}  })) {
         if ($accum > 0) {
             for (my $i = $prev_pos; $i < $key_pos; $i++) {
                 print $fh_all "$key_chr\t$i\t$accum\n";
             }
         }
         $accum += $softhash{$key_chr}{$key_pos};
         $prev_pos = $key_pos;
    }

主处理应该比以前快得多，输出循环应该稍微快一些，因为(keys %{ $softhash{$key_chr} })需要排序的项较少。

如果我以上提到的$loc是正确的，那么可能值得消除$j上的循环，代之以无条件更新$var_start和$var_end。这可能无助于性能，但它可能会使代码更具可读性。

为了便于阅读，我要做的其他事情主要是与名字有关。我认为提取$chromosome = $a[2]是有帮助的；我想知道$var_start和$var_end是否比$range_start和$range_end更好；而且我很确定key_前缀并不能传递多少有用的信息。例如，重用$chromosome而不是$key_chr。

另一个可读性问题:我对循环for (my $i = 0; $i <= $#c; $i++)感到困惑，因为它显然运行了太多次。要么在1开始迭代，要么找到某种方法使索引从0开始。如果您无条件地实现更新范围的建议，我认为您不再需要引用$b的旧值，可能可以用while ($cigar =~ /([0-9]+)([A-Z])/g)或类似的东西替换循环。(未经测试)。