0

我遇到了一个相当独特的问题。我有 2 个正在阅读的文件。这 2 个文件的小版本如下所示：

文件 1

chr1    9873    12227   11873   2354    +   NR_046018   DDX11L1
chr1    760970  763155  762970  2185    +   NR_047520   LOC643837

文件2

chr1    9871    0   
chr1    9872    1
chr1    9873    1
chr1    9874    2
chr1    9875    1
chr1    9876    3
chr1    9877    3
chr1    760970  1
chr1    760971  1
chr1    760972  1
chr1    760973  2
chr1    760974  3
chr1    760975  3
chr1    760976  4
chr1    760977  5
chr1    760978  6
chr1    760979  7
chr1    760980  6
chr1    760981  7
chr1    760982  8
chr1    760983  9
chr1    760984  10
chr1    760985  11
chr1    760986  12
chr1    760987  10
chr1    760988  9
chr1    760989  6

问题

  1. 从第一个文件中，我必须从每一行中提取第二个元素并将其作为$start. 结束位置由 确定$end = $start + 10。

  2. 基于$start，我现在必须获取第二个文件，并查看每行的第二个元素。找到$start后，我需要将 5 个一组的第三个元素的下 5 个对应值相加，最多为$end.

因此，我以 5 个为一组进行求和，将获得 2 个求和值 $end。$start + 10

如果某些值 upto$end不存在于第二个文件的第二个元素中，代码不应停止，它应继续执行求和并将总和显示为 0（如果不存在连续的 5 个元素组）。

以此处的文件为例，来自File1，第二个元素 = 9873，它被分配给$start. 因此$end将是$start+109883。

从File2中，一旦$start在该行的第 2 个元素中找到，接下来 5 行的第 3 个元素必须作为 1 组求和，接下来的 5 个值作为第 2 组求和直到$end。

笔记

这里可以在File2中看到，$end即 9883 不存在。因此，从 9879 到 9883 的值的总和必须是zero。它不能将 760970 的值相加......

期望的输出

chr1    9873    12227   11873   2354    +   NR_046018   DDX11L1      10   0
chr1    760970  763155  762970  2185    +   NR_047520   LOC643837    8   25

注意事项

  1. 在处理实际文件时， $end = $start+10,000(而不是 $end = $start+10 )
  2. 此外，在同一个注释中，将对 25 个值的组求和（而不是 5 个），在处理实际文件时总共获得 400 个值。
  3. 如果 $file2 的第二个元素中不存在一系列值，则求和应正常进行，如果不存在连续的 25 个值对，0则应打印输出。
  4. 每个文件包含 > 100 万行。

代码

到目前为止，我编写的代码能够做到以下几点：

  1. 从文件中读取。
  2. 分配$start和$end来自file1
  3. 从file2，将所有第二个元素推入 array @c_posn；所有第三个元素到数组@peak中。
  4. 检查是否$start存在于@c_posn

我无法弄清楚如何进行求和部分。我曾想过创建一个哈希，其中 2nd 文件的所有 2nd 元素进入keys和 3rd 元素进入values。但是哈希是无序的。@c_posn所以我为第二个元素@peaks和第三个元素创建了 2 个数组。但现在我不知道如何同时比较 2 个数组（以确保 760970 的值不会相加）

use 5.012;
use warnings;
use List::Util qw/first/;

my $file1 = 'chr1trialS.out';
my $file2 = 'b1.wig';

open my $fh1,'<',$file1 or die qw /Can't_open_file_$file1/;
open my $fh2,'<',$file2 or die qw /Can't_open_file_$file2/;

my($start, $end);
while(<$fh1>){
    my @val1 = split;
    $start = $val1[1]; #Assign start value
    $end = $start + 10; #Assign end value
    say $start,"->",$end; #Can be commented out
}

my @c_posn;
my @peak;

while(<$fh2>){
    my @val2 = split;   
    push @c_posn,$val2[1]; #Push all 2nd elements 
    push @peak, $val2[2];  #Push all 3rd elements        
}           

if (first { $_ eq $start} @c_posn) { say "I found it! " } #To check if $start is present in @c_posn

say "@c_posn"; #just to check all 2nd elements are obtained
say "@peak"; #just to check all 3rd elements are obtained

感谢您花时间解决我的问题。如果需要任何澄清，请务必问我。我将不胜感激任何评论/回答。

4

3 回答 3

2

你对哈希有正确的想法。是否已订购并不特别相关，因为如果我理解正确，您正在寻找 11 个特定值（9873、9874、9875... 9883），而不是文件中的起始值和下一个 10（9873 ,... 9877, 760970,... 760975)。

根据您的描述，我将采取以下措施：

#!/usr/bin/env perl

use strict;
use warnings;

my $sum_interval = 5;   # number of lines to group into each sum
my $sum_count = 2;      # number of sums to generate
my @sums;               # final results of the operation

my %lookup;
open my $fh2, '<', 'file2.txt' or die "Can't open file 2: $!";
while (<$fh2>) { 
  my @data = split;
  $lookup{$data[1]} = $data[2];
}
close $fh2;

open my $fh1, '<', 'file1.txt' or die "Can't open file 1: $!";
while (my $line = <$fh1>) { 
  my @line_sums;
  my $start = (split /\s+/, $line)[1];
  for my $interval_num (0 .. $sum_count - 1) {
    my $cur_sum = 0;
    my $interval_start = $start + ($sum_interval * $interval_num);
    for (0 .. $sum_interval - 1) {
      # use || instead of // for Perl older than 5.10
      $cur_sum += $lookup{$interval_start + $_} // 0;
    }
    push @line_sums, $cur_sum;
  }
  push @sums, \@line_sums;
}
use Data::Dumper; print Dumper(\@sums);

变量名称可能可以改进，但您可以将$sum_intervaland更改$sum_count为 25 和 400，它应该在您的实际应用程序中同样工作。

如果您提供的样本数据被放入file1.txtand file2.txt，这将产生输出：

$VAR1 = [
          [
            10,
            0
          ],
          [
            8,
            25
          ]
        ];

如果我手动进行总和，此输出与我得出的结果相匹配。

请注意，我与您的规范略有不同，因为它总和从$startto$start + 9而不是$start + 10因为您说它应该为两组五个和$startto求和$start + 1011 个项目。

编辑：将初始伪代码修改为完整的可运行程序。

于 2013-01-02T11:04:16.957 回答
2

这很简单，如果b1.wig足够小可以读入内存中的哈希，从第 2 列获取键，从第 3 列获取值。然后必须做的就是访问每个序列中的每个键，如果 a相应的哈希元素不存在（因此访问它返回undef）。

你还没有说你想如何将新总数与现有数据分开，chr1trialS.out所以我使用了空格。当然，如果需要，这很容易改变。

use strict;
use warnings;

use constant SAMPLE_SIZE => 10;
use constant CHUNK_SIZE => 5;

my $file1 = 'chr1trialS.out';
my $file2 = 'b1.wig';

my %data2;
{
  open my $fh, '<', $file2 or die $!;

  while (<$fh>) {
    my ($key, $val) = (split)[1,2];
    $data2{$key} = $val;
  }
}

open my $fh, '<', $file1 or die $!;

while (<$fh>) {
  chomp;
  my $key = (split)[1];
  my @totals;
  my $n = 0;
  while ($n < SAMPLE_SIZE) {
    push @totals, 0 if $n++ % CHUNK_SIZE == 0;
    $totals[-1] += $data2{$key++} // 0;
  }
  print "$_ @totals\n";
}

输出

chr1    9873    12227   11873   2354    +   NR_046018   DDX11L1 10 0
chr1    760970  763155  762970  2185    +   NR_047520   LOC643837 8 25
于 2013-01-02T18:52:00.200 回答
1

这是我目前的解决方案：

#!/usr/bin/perl

use 5.012; use warnings;

my $file1 = Reader->open("<", "filename1");
my $file2 = Reader->open("<", "filename2");

my $groupsize = 5;
my $step = 10;
my $sum_number = int($step / $groupsize) + ($step % $groupsize ? 1 : 0); # ceil($step/$groupsize)

use constant DEBUG_FLAG => 0;
sub DEBUG (@)   { say STDERR "DEBUG: ", @_ if DEBUG_FLAG }

LINE1:
while (my $line1 = $file1->readline) {
    my (undef, $start) = split ' ', $line1, 3;
    my $end = $start + $step;
    my @sums = (0) x $sum_number; # initialize all fields to zero
    my $i = 0;
    my $last;
    LINE2:
    while (my $line2 = $file2->readline) {
        my (undef, $key, $val) = split ' ', $line2, 4;
        if ($start > $key) { # throw away all keys that are too small
            DEBUG "key $key too small for start $start";
        } elsif ($key >= $end) { # termination condition
            DEBUG "key $key too large for end $end";
            $file2->pushback($line2);
            last LINE2;
        } else {
            $last = $key unless defined $last;
            $i += $key - $last; # get interval. This may be set to "1" as an optimization
            DEBUG "counting ($i): $sums[$i/$groupsize] + $val at $key";
            $sums[$i/$groupsize] += $val;
            $last = $key;
        }
    }
    DEBUG "inner loop broken";
    say join "\t", $line1, @sums; # assuming tab-seperated output
}

{
    package Reader;
    # There is probably a CPAN module for this ... :/
    use Carp;
    use constant DEBUG_FLAG => 0;
    sub open :method {
        my ($class, $mode, $filename) = @_;
        open my $fh, $mode, $filename or die qq(Can't open "$filename": $!);
        bless [$fh, []] => $class;
    }
    sub readline :method {
        my $self = shift;
        return shift @{ $self->[1] } if @{ $self->[1] };
        my $line = scalar readline $self->[0];
        chomp $line if defined $line;
        carp "readline: " . ($line // "undef") if DEBUG_FLAG;
        return $line;
    }
    sub pushback {
        my ($self, $line) = @_;
        carp "pushback: " . ($line // "undef") if DEBUG_FLAG;
        unshift @{ $self->[1] }, $line;
        return $self;
    }
    sub eof :method {
        my $self = shift;
        eof $self->[0];
    }
}

输出：

chr1    9873    12227   11873   2354    +   NR_046018   DDX11L1         10      0
chr1    760970  763155  762970  2185    +   NR_047520   LOC643837       8       25

此解决方案假定两个输入文件都按第二个字段按升序排序，并且不会请求重叠序列。如果可以满足这些条件，它将在恒定内存和线性时间中执行。如果不是，它会产生垃圾，您可能有兴趣使用其他答案（线性内存、线性时间、无限制）。事实上，Dave Sherohman 的答案总体上不那么脆弱，并且可能在大多数输入上执行得更快。

根据您的系统，如果您放弃所有面向对象，并内联缓冲行（或者更确切地说，一行）的代码，则可能会提高速度。

关于$i = $key - $last：如果跳过键，代码将继续工作，并且仍将数字添加到正确的存储桶中。如果您可以断言不会跳过任何键，或者正确的总和无关紧要（ID 小于 的前五行$end，而不是应该添加接下来的五个 ID），那么删除$last变量并简单地递增$i1 即可。

于 2013-01-02T12:44:21.593 回答