免费注册 查看新帖 |

Chinaunix

  平台 论坛 博客 文库
12下一页
最近访问板块 发新帖
查看: 4340 | 回复: 11

求助一个perl数据区间的统计! [复制链接]

论坛徽章:
0
发表于 2017-06-16 05:31 |显示全部楼层
本帖最后由 xiayuu119 于 2017-06-16 05:42 编辑

一个区间问题累加。麻烦大神了!

输入数据:
A.txt
---------------------------------------------------------
#ID        site        a        b        c        d        e
chr1        1        0/1        1/2        4/4        -        0/6
chr1        3        0/1        1/2        4/4        4/5        0/6
chr1        7        0/1        1/1        4/4        -        0/6
chr1        10        0/1        1/1        4/4        4/5        0/6
chr1        22        0/2        1/1        4/4        -        0/6
chr1        23        0/2        1/1        4/4        4/5        0/6
chr1        25        0/2        1/1        4/4        4/5        0/6
chr1        26        0/2        1/1        4/4        4/5        0/6
chr1        30        1/1        1/1        4/4        4/5        0/6
chr1        35        1/1        1/2        4/4        4/5        0/6
chr1        41        1/1        1/2        4/4        4/5        0/6
chr1        47        1/1        1/2        4/4        4/5        0/7
chr1        51        0/2        1/2        4/4        4/5        0/7
chr1        57        0/2        1/2        4/4        4/5        0/7
chr1        58        0/2        1/3        4/4        4/5        0/7
chr1        61        0/2        1/3        4/4        4/5        0/7
chr1        63        1/2        1/3        4/4        4/5        0/7
chr1        67        1/2        1/3        4/4        4/5        0/7
chr2        7        1/2        1/3        4/4        4/5        0/7
chr2        11        1/2        1/2        4/4        4/5        0/6
chr2        15        1/2        1/2        4/4        4/5        0/6
chr2        17        1/2        1/2        4/4        4/5        0/6
chr2        18        1/2        1/2        4/4        4/5        0/6
chr2        25        0/2        1/2        4/4        4/5        0/6
chr2        26        0/2        1/2        4/4        4/5        0/6

B.txt
---------------------------------------------------------
ID        ID2        start        end
chr1        g1        1        20
chr1        g2        22        29
chr1        g3        25        50
chr1        g4        45        68
chr2        g5        2        20
chr2        g6        18        31


结果:
在文件B.txt中start-end区间,A.txt文件site的位点进行计算,“/”左右的数据分别累加。(缺失数据为“-”,不计算)

C.txt
---------------------------------------------------------
ID2        ID        a        b        c        d        e
g1        chr1:1-20        0/4        4/6        16/16        8/10        0/24
g2        chr1:22-29        0/8        4/4        16/16        12/15        0/24
g3        chr1:25-50        4/8        6/9        24/24        24/30        0/37
g4        chr1:45-68        3/13        7/18        28/28        28/35        0/49
g5        chr2:2-20        5/10        5/11        20/20        20/25        0/31
g6        chr2:18-31        1/6        3/6        12/12        12/15        0/18        

论坛徽章:
145
技术图书徽章
日期:2013-10-01 15:32:13戌狗
日期:2013-10-25 13:31:35金牛座
日期:2013-11-04 16:22:07子鼠
日期:2013-11-18 18:48:57白羊座
日期:2013-11-29 10:09:11狮子座
日期:2013-12-12 09:57:42白羊座
日期:2013-12-24 16:24:46辰龙
日期:2014-01-08 15:26:12技术图书徽章
日期:2014-01-17 13:24:40巳蛇
日期:2014-02-18 14:32:59未羊
日期:2014-02-20 14:12:13白羊座
日期:2014-02-26 12:06:59
发表于 2017-06-16 19:52 |显示全部楼层
回复 1# xiayuu119

$ awk -f get_chrSE A.txt B.txt
ID2     ID             a       b       c       d       e      
g1      chr1:1-20      0/4     4/6     16/16   8/10    0/24   
g2      chr1:22-29     0/8     4/4     16/16   12/15   0/24   
g3      chr1:25-50     4/8     6/9     24/24   24/30   0/37   
g4      chr1:45-68     3/13    7/18    28/28   28/35   0/49   
g5      chr2:2-20      5/10    5/11    20/20   20/25   0/31   
g6      chr2:18-31     1/6     3/6     12/12   12/15   0/18


$ cat get_chrSE.awk
BEGIN{
  S=SUBSEP;
  split("8,15,8,8,8,8,8",f,",")
}
function P(n,s){
  s="";
  for(n=1;n<=7;++n){
    s=s sprintf("%-"f[n]"s",d[n]);
    d[n]="";
  }
  print s
}
function A(s1,s2,nt,dt){
  split(s1,b,"/");
  nt=b[1];
  dt=b[2];
  split(s2,b,"/");
  nt+=b[1];
  dt+=b[2];
  return(nt"/"dt);
}

function X(g,c,s,e,n,m,t){
  d[1] = g;
  if(c=="ID"){
    d[2] = c;
    for(n=3;n<=7;++n)
      d[n]=a["#"c,"site",n];
    P();
    next
  }
  d[2] = c":"s"-"e
  for(n=s;n<=e;++n){
    if(!a[c,n])continue;
    for(m=3;m<=7;++m)
      d[m] = A(d[m],a[c,n,m]);
  }
  P();
}
FNR==NR{
  a[$1,$2]=1;
  for(n=3;n<=NF;++n)
  a[$1,$2,n]=$n;
  next
}
{
  X($2,$1,$3,$4)
}

论坛徽章:
307
程序设计版块每周发帖之星
日期:2016-04-08 00:41:33操作系统版块每日发帖之星
日期:2015-09-02 06:20:00每日论坛发贴之星
日期:2015-09-02 06:20:00程序设计版块每日发帖之星
日期:2015-09-04 06:20:00每日论坛发贴之星
日期:2015-09-04 06:20:00每周论坛发贴之星
日期:2015-09-06 22:22:00程序设计版块每日发帖之星
日期:2015-09-09 06:20:00程序设计版块每日发帖之星
日期:2015-09-19 06:20:00程序设计版块每日发帖之星
日期:2015-09-20 06:20:00每日论坛发贴之星
日期:2015-09-20 06:20:00程序设计版块每日发帖之星
日期:2015-09-22 06:20:00程序设计版块每日发帖之星
日期:2015-09-24 06:20:00
发表于 2017-06-17 01:18 |显示全部楼层
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;

  4. sub Load{
  5.         local @ARGV = @_;
  6.         my %hData = ();
  7.         my %hToPos = ();
  8.         my $cnt = 0;
  9.         while(<>){
  10.                 next if(/\A[A-Z]/);
  11.                 my ($id, $id2, $start, $end) = split;
  12.                 push(@{$hData{$id}}, [$start, $id2, $end]);
  13.                 $hToPos{$id2} = $cnt++;
  14.         }
  15.         \%hData, \%hToPos;
  16. }

  17. my ($rhData, $rhToPos) = Load(pop);
  18. my @aOut = ();
  19. while(<>){
  20.         next if(/\A#/);
  21.         my ($id, $site, @aT) = split;
  22.         @aT = map {/-/ ? [0, 0] : [split('/', $_)]} @aT;
  23.         foreach my $it (grep{$site >= $_->[0] and $site < $_->[-1]} @{$rhData->{$id}}){
  24.                 my $pos = $rhToPos->{$it->[1]};
  25.                 if(!defined($aOut[$pos])){
  26.                         $aOut[$pos] = [$it->[1], "$id:$it->[0]-$it->[-1]", @aT];
  27.                         next;
  28.                 }
  29.                 foreach(2 .. $#{$aOut[$pos]}){
  30.                         $aOut[$pos][$_] = [$aOut[$pos][$_][0] + $aT[$_ - 2][0], $aOut[$pos][$_][-1] + $aT[$_ - 2][-1]];
  31.                 }
  32.         }
  33. }

  34. foreach my $it (@aOut){
  35.         print(join("\t", @{$it}[0, 1], map {join('/', @{$it->[$_]})} 2 .. $#{$it}), "\n");
  36. }
复制代码

perl abc.pl a b
----------------------------------------------------------------------------
g1      chr1:1-20       0/4     4/6     16/16   8/10    0/24
g2      chr1:22-29      0/8     4/4     16/16   12/15   0/24
g3      chr1:25-50      4/8     6/9     24/24   24/30   0/37
g4      chr1:45-68      3/13    7/18    28/28   28/35   0/49
g5      chr2:2-20       5/10    5/11    20/20   20/25   0/31
g6      chr2:18-31      1/6     3/6     12/12   12/15   0/18

论坛徽章:
0
发表于 2017-06-17 01:45 |显示全部楼层
sunzhiguolu 发表于 2017-06-17 01:18
perl abc.pl a b
----------------------------------------------------------------------------
g1  ...

太感谢了!! 献上膝盖!简洁实用!

论坛徽章:
0
发表于 2017-06-20 02:35 |显示全部楼层
回复 3# sunzhiguolu

您好! 感谢您的帮助! 少量数据可以完成,但是我的数据非常巨大,我无法使用这个成功跑通(似乎都存入内存了)。有没有更好的解决方法呢?

论坛徽章:
307
程序设计版块每周发帖之星
日期:2016-04-08 00:41:33操作系统版块每日发帖之星
日期:2015-09-02 06:20:00每日论坛发贴之星
日期:2015-09-02 06:20:00程序设计版块每日发帖之星
日期:2015-09-04 06:20:00每日论坛发贴之星
日期:2015-09-04 06:20:00每周论坛发贴之星
日期:2015-09-06 22:22:00程序设计版块每日发帖之星
日期:2015-09-09 06:20:00程序设计版块每日发帖之星
日期:2015-09-19 06:20:00程序设计版块每日发帖之星
日期:2015-09-20 06:20:00每日论坛发贴之星
日期:2015-09-20 06:20:00程序设计版块每日发帖之星
日期:2015-09-22 06:20:00程序设计版块每日发帖之星
日期:2015-09-24 06:20:00
发表于 2017-06-20 10:13 |显示全部楼层
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;

  4. sub Load{
  5.         local @ARGV = @_;
  6.         my %hData = ();
  7.         while(<>){
  8.                 next if(/\A[A-Z]/);
  9.                 my ($id, $id2, $start, $end) = split;
  10.                 push(@{$hData{$id}}, [$start, $id2, $end]);
  11.         }
  12.         \%hData;
  13. }

  14. my $rhData = Load(pop);
  15. my $iFile = 'urPath'; #Your intermediate file
  16. open(my $FHw, '>', $iFile); #Create Intermediate File, DataFormat:{$id2, $id, "$start-$end", $tail}
  17. while(<>){
  18.     next if(/\A#/);
  19.     my ($id, $site, $tail) = split(' ', $_, 3);
  20.     foreach my $it (grep {$site >= $_->[0] and $site < $_->[-1]} @{$rhData->{$id}}){
  21.         print $FHw join("\t", $it->[1], $id, "$it->[0]-$it->[-1]", $tail);
  22.     }
  23. }
  24. close($FHw);

  25. $rhData = {};
  26. open(my $FHr, '<', $iFile);
  27. while(<$FHr>){
  28.     my ($id2, $id, $section, @aT) = split;
  29.     @aT = map {/-/ ? [0, 0] : [split('/', $_)]} @aT;
  30.     my $this = "$id:$section";
  31.     if(!exists($rhData->{$this})){
  32.         $rhData->{$this} = [$id2, @aT];
  33.         next;
  34.     }
  35.     my $who = $rhData->{$this};
  36.     foreach my $pos (0 .. $#aT){
  37.         $who->[$pos + 1][0] += $aT[$pos][0];
  38.         $who->[$pos + 1][-1] += $aT[$pos][-1];
  39.     }
  40. }
  41. close($FHr);

  42. foreach(keys %{$rhData}){
  43.     my ($id2, @aT) = @{$rhData->{$_}};
  44.     print(join("\t", $id2, $_, map {join('/', @{$_})} @aT), "\n");
  45. }

  46. unlink($iFile); #Del intermediate file
复制代码

perl abc.pl a b
----------------------------------------------------------------
g4      chr1:45-68      3/13    7/18    28/28   28/35   0/49
g1      chr1:1-20       0/4     4/6     16/16   8/10    0/24
g2      chr1:22-29      0/8     4/4     16/16   12/15   0/24
g3      chr1:25-50      4/8     6/9     24/24   24/30   0/37
g6      chr2:18-31      1/6     3/6     12/12   12/15   0/18
g5      chr2:2-20       5/10    5/11    20/20   20/25   0/31

论坛徽章:
307
程序设计版块每周发帖之星
日期:2016-04-08 00:41:33操作系统版块每日发帖之星
日期:2015-09-02 06:20:00每日论坛发贴之星
日期:2015-09-02 06:20:00程序设计版块每日发帖之星
日期:2015-09-04 06:20:00每日论坛发贴之星
日期:2015-09-04 06:20:00每周论坛发贴之星
日期:2015-09-06 22:22:00程序设计版块每日发帖之星
日期:2015-09-09 06:20:00程序设计版块每日发帖之星
日期:2015-09-19 06:20:00程序设计版块每日发帖之星
日期:2015-09-20 06:20:00每日论坛发贴之星
日期:2015-09-20 06:20:00程序设计版块每日发帖之星
日期:2015-09-22 06:20:00程序设计版块每日发帖之星
日期:2015-09-24 06:20:00
发表于 2017-06-21 00:51 |显示全部楼层
Jason680 大神,要不您老人家来一个 让我等小白学习下处理超大文件的技巧???!!!

论坛徽章:
307
程序设计版块每周发帖之星
日期:2016-04-08 00:41:33操作系统版块每日发帖之星
日期:2015-09-02 06:20:00每日论坛发贴之星
日期:2015-09-02 06:20:00程序设计版块每日发帖之星
日期:2015-09-04 06:20:00每日论坛发贴之星
日期:2015-09-04 06:20:00每周论坛发贴之星
日期:2015-09-06 22:22:00程序设计版块每日发帖之星
日期:2015-09-09 06:20:00程序设计版块每日发帖之星
日期:2015-09-19 06:20:00程序设计版块每日发帖之星
日期:2015-09-20 06:20:00每日论坛发贴之星
日期:2015-09-20 06:20:00程序设计版块每日发帖之星
日期:2015-09-22 06:20:00程序设计版块每日发帖之星
日期:2015-09-24 06:20:00
发表于 2017-06-21 00:53 |显示全部楼层
rubyish 大神,来难题了。要不您老人家出手,让我等小白学习下处理超大文件的技巧???!!!

论坛徽章:
0
发表于 2017-06-21 02:10 |显示全部楼层
回复 8# sunzhiguolu sunzhiguolu Jason680  rubyish 先叩谢各位大神!!~!让我小小白能够有机会学习学习!!!!


论坛徽章:
307
程序设计版块每周发帖之星
日期:2016-04-08 00:41:33操作系统版块每日发帖之星
日期:2015-09-02 06:20:00每日论坛发贴之星
日期:2015-09-02 06:20:00程序设计版块每日发帖之星
日期:2015-09-04 06:20:00每日论坛发贴之星
日期:2015-09-04 06:20:00每周论坛发贴之星
日期:2015-09-06 22:22:00程序设计版块每日发帖之星
日期:2015-09-09 06:20:00程序设计版块每日发帖之星
日期:2015-09-19 06:20:00程序设计版块每日发帖之星
日期:2015-09-20 06:20:00每日论坛发贴之星
日期:2015-09-20 06:20:00程序设计版块每日发帖之星
日期:2015-09-22 06:20:00程序设计版块每日发帖之星
日期:2015-09-24 06:20:00
发表于 2017-06-22 00:28 |显示全部楼层
回复 9# xiayuu119
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;

  4. sub Load{
  5.     local @ARGV = pop;
  6.     my @aData = ();
  7.     while(<>){
  8.         push(@aData, [split]) if(!/\A[A-Z]/);
  9.     }
  10.     \@aData;
  11. }

  12. my $raData = Load(pop);
  13. my $bFile = pop;
  14. while(my $it = shift(@{$raData})){
  15.     my @aOut = ([0, 0], [0, 0], [0, 0], [0, 0], [0, 0]);
  16.     open(my $FHr, '<', $bFile);
  17.     while(<$FHr>){
  18.         next if(/\A#/);
  19.         my ($id, $site, @aT) = split;
  20.         @aT = map {/-/ ? [0, 0] : [split('/', $_)]} @aT;
  21.         if($it->[0] eq $id and $site >= $it->[-2] and $site < $it->[-1]){
  22.             foreach my $pos (0 .. $#aT){
  23.                 $aOut[$pos][0] += $aT[$pos][0];
  24.                 $aOut[$pos][-1] += $aT[$pos][-1];
  25.             }
  26.         }
  27.     }
  28.     close($FHr);
  29.     print(join("\t", $it->[1], "$it->[0]:$it->[-2]-$it->[-1]", map {join('/', @$_)} @aOut), "\n");
  30. }
复制代码
perl abc.pl a b
------------------------------------------------------------
g1      chr1:1-20       0/4     4/6     16/16   8/10    0/24
g2      chr1:22-29      0/8     4/4     16/16   12/15   0/24
g3      chr1:25-50      4/8     6/9     24/24   24/30   0/37
g4      chr1:45-68      3/13    7/18    28/28   28/35   0/49
g5      chr2:2-20       5/10    5/11    20/20   20/25   0/31
g6      chr2:18-31      1/6     3/6     12/12   12/15   0/18


您需要登录后才可以回帖 登录 | 注册

本版积分规则 发表回复

  

北京盛拓优讯信息技术有限公司. 版权所有 京ICP备16024965号-6 北京市公安局海淀分局网监中心备案编号:11010802020122 niuxiaotong@pcpop.com 17352615567
未成年举报专区
中国互联网协会会员  联系我们:huangweiwei@itpub.net
感谢所有关心和支持过ChinaUnix的朋友们 转载本站内容请注明原作者名及出处

清除 Cookies - ChinaUnix - Archiver - WAP - TOP