【问题】连续区间合并

huang6894

额，是我上一个帖子的简化版~

现在得到一些这样的文本：

1. 
   
2. 16642     chr1    2336225 18      HOM     del
   
3. 16642     chr1    2336226 17      HOM     del
   
4. 16642     chr1    2336227 17      HOM     del
   
5. 16642     chr1    2336228 18      HOM     del
   
6. 16642     chr1    2336236 17      HOM     del
   
7. 16642     chr1    2336237 17      HOM     del
   
8. 16642     chr1    2336309 24      HOM     del
   
9. 16642     chr1    2336320 24      HOM     del
   
10. 16642     chr1    2336321 24      HOM     del
    
11. 16642     chr1    2336322 24      HET     del
    
12. 16642     chr1    2336323 24      HET     del
    
13. 16642     chr1    2336324 24      HET     del
    
14. 16642     chr1    2336325 24      HOM     del
    

复制代码

我想对上面文本进行处理，如果第2、5、6列相同时，判断第3列是否连续，如果连续的数字达到3行以上，合并之~最后得到：

1. 
   
2. 16642     chr1    2336225-2336228      HOM     del
   
3. 16642     chr1    2336322-2336324      HET     del
   

复制代码

我记得之前有人回答过相关的问题，可是找不到原贴了，只好厚颜请求大神帮帮忙了，谢谢谢谢



----------------------------update--------------------------------------

问题解决，非常感谢两位大师@stanley_tam@yestreenstars
因为只能给一个最佳答案，选了一个时间相对短一点的。。。
stanley_tam大神的答案见：9楼

yestreenstars

1. #!/usr/bin/perl
   
2. use warnings;
   
3. use strict;
   
4. my $first_line = <DATA>;
   
5. my $str = join(" ", (split /\s+/, $first_line)[0, 1, 4, 5]);
   
6. my($min, $max);
   
7. $min = $max = (split /\s+/, $first_line)[2];
   
8. while(<DATA>){
   
9.         my @fields = split;
   
10.         if(join(" ", (split)[0, 1, 4, 5]) eq $str){
    
11.                 if($fields[2] - $max == 1){
    
12.                         $max = $fields[2];
    
13.                         next;
    
14.                 }else{
    
15.                         print join("\t", (split /\s+/, $str)[0, 1], "$min-$max", (split /\s+/, $str)[2, 3]), $/ if $max - $min >= 2;
    
16.                 }
    
17.         }else{
    
18.                 print join("\t", (split /\s+/, $str)[0, 1], "$min-$max", (split /\s+/, $str)[2, 3]), $/ if $max - $min >= 2;
    
19.         }
    
20.         $str = join(" ", (split)[0, 1, 4, 5]);
    
21.         $min = $max = $fields[2];
    
22. }
    
23. print join("\t", (split /\s+/, $str)[0, 1], "$min-$max", (split /\s+/, $str)[2, 3]), $/ if $max - $min >= 2;
    
24. __DATA__
    
25. 16642     chr1    2336225 18      HOM     del
    
26. 16642     chr1    2336226 17      HOM     del
    
27. 16642     chr1    2336227 17      HOM     del
    
28. 16642     chr1    2336228 18      HOM     del
    
29. 16642     chr1    2336236 17      HOM     del
    
30. 16642     chr1    2336237 17      HOM     del
    
31. 16642     chr1    2336309 24      HOM     del
    
32. 16642     chr1    2336320 24      HOM     del
    
33. 16642     chr1    2336321 24      HOM     del
    
34. 16642     chr1    2336322 24      HET     del
    
35. 16642     chr1    2336323 24      HET     del
    
36. 16642     chr1    2336324 24      HET     del
    
37. 16642     chr1    2336325 24      HOM     del

复制代码

q1208c

以前是写过一个. 但解决不了楼主的问题.

个人以为, 楼主自己要的未必只是合并这么简单. 可否先对源文件进行排序, 然后再处理. 这样, 由于只是处理连续读进来的数据, 所以, 不会消耗太多的内存. 也更简单.

huang6894

回复 2# q1208c


    是的，我也是这么做的，可是突然不怎么记得如何合并了，他们确实排序了

q1208c

回复 3# huang6894


排序了还不简单, 把你要的几个变量合并成一个. 存为 $aa

如果 $aa 为 空, 说明是第一行, 读下一行, 然后比较 $3, 如果 $3 跟前一个 $3 是连续的, 那就再读一行. 满足条件就输出. 不然就重置.     

huang6894

回复 4# q1208c


    满足条件就输出. 不然就重置.就是这个不懂 而且实际数据时我是需要满足200行。。。

yinyuemi

回复 1# huang6894


    看你自己发表的帖子就能找到原帖

huang6894

回复 6# yinyuemi


    谢谢SS兄。。。我应该没问过，不然我应该记得~谢谢~

stanley_tam

貌似设置几个变量记录上一行的信息和start number就可以了，这样行不行{:3_193:}

1. #!perl
   
2. use Modern::Perl;
   
3. use FileHandle;
   
4. 
   
5. my $your_data_file = 'a.txt';
   
6. my $fh = FileHandle->new($your_data_file, 'r');
   
7. 
   
8. my $start_num = q{};
   
9. my $last_line = [ q{}, q{}, 0, q{}, q{}, q{} ];
   
10. 
    
11. my $count = 0;
    
12. while (my $line = $fh->getline) {
    
13.     # remove trailing newline
    
14.     chomp $line;
    
15. 
    
16.     # get columns of current line
    
17.     my $current_line = [(split /\s+/, $line)];
    
18. 
    
19.     # compare currentline with last line
    
20.     if ($current_line->[2] == $last_line->[2] + 1) { # number increase by one
    
21.         
    
22.         if ($last_line->[1] eq $current_line->[1] and
    
23.             $last_line->[4] eq $current_line->[4] and
    
24.             $last_line->[5] eq $current_line->[5]) { # column 2, 5, 6 are equal
    
25.             # valid consecutive row, note it
    
26.             ++$count;
    
27.         }
    
28.         else {
    
29.             # column 2, 5, 6 not equal
    
30.             # the first line doesn't count, $count equal 2 suffice
    
31.             if ($count >= 2) {
    
32.                 my @copy = @{ $last_line };
    
33.                 $copy[2] = "$start_num-$copy[2]";
    
34.                 say join "\t", @copy;
    
35.             }
    
36. 
    
37.             # reset
    
38.             $start_num = $current_line->[2];
    
39.             $count     = 0;
    
40.         }
    
41.     }
    
42.     else {
    
43.         # number not consecutive
    
44.         if ($count >= 2) {
    
45.             my @copy = @{ $last_line };
    
46.             $copy[2] = "$start_num-$copy[2]";
    
47.             say join "\t", @copy;
    
48.         }
    
49. 
    
50.         # reset
    
51.         $start_num = $current_line->[2];
    
52.         $count     = 0;
    
53.     }
    
54. 
    
55.     $last_line = $current_line;
    
56. }
    
57. 
    
58. $fh->close;
    
59. 
    
60. __END__

复制代码

chenhao392

貌似可以用BEDtools, mergeBED和multiIntersectBed..自己google之。
遇到生信的问题先想别人已经写过的binary，多上Seqanswers和Biostar查查看。就像写代码的查查CPAN，去stackoverflow一样..
而不是直接写Perl。
BEDtools是C写的，代码效率也不错。

https://www.biostars.org/p/49581/