求大神帮写一个去除交集的perl脚本！！

清泉一边

智能匹配！回复 14# xiumu2280


   

expert1

有点月经问题，awk的关联数组就能搞掂，在shell/perl版我看到此类问题不下3次了。

Monox

今天描了一眼这个帖子的回复，满眼“大神”，Perl 版大神真多。
上面是题外话了。楼主这样的需求我写过好多次了，不过每次用完就扔了（其实应该有备份，文件太多懒得翻了，而且可能还需要修改才能完全符合楼主的需求）。
其实这个问题有一个大的方面要考虑，就是生物信息方面的数据一般数据量都很大，用一般的方法处理都会需要很多内存。当然使用集群的话内存问题可以解决，但是你需要大内存的时候任务排不上队，算上排队的时间不一定比用低内存 footprint 的方法有效。而另一方面，生物信息数据处理的输入往往是上流处理的输出，一般都是排好序的（如果不是的话可以调整上流分析流程，当然中间加一步排序操作应该也是可以接受的）。然后处理的时候用几个全局变量保存当前处理到哪个 scafforld 了，当前处理到的坐标范围是哪，就可以两个文件按行读取，且把一段 segment 的范围降到最低了就可以直接输出这个区段，就可以继续处理下个区段了。这样，实际的内存的 footprint 在 250M 以内（大概），时间也还可以接受（相当于把两个文件按行读取的时间消耗）。
不需要使用 hash，甚至数组都不需要，使用几个标量就够了。

pitonas

{:2_172:}
不好意思啊, 师傅可以提供一个 perl 代码? 我来研究研究。
谢谢了。

回复 53# Monox


   

Monox

回复 54# pitonas

好吧，其实我一般很懒得帖现成代码的，下面是一个demo我刚随便写的，里面有很多重复代码，感兴趣的人可以对代码进行重构：

1. 
   
2. 
   
3.     #!/usr/bin/env perl
   
4. use v5.18;
   
5. use autodie;
   
6. 
   
7. my $file1 = shift;
   
8. my $file2 = shift;
   
9. my $out_file = shift;
   
10. 
    
11. open my $fh_in1, "<", $file1;
    
12. open my $fh_in2, "<", $file2;
    
13. open my $fh_out, ">", $out_file;
    
14. 
    
15. my $line = <$fh_in1>;
    
16. my ($current_scaffold1, $current_start1, $current_end1) = split /\s+/, $line;
    
17. 
    
18. $line = <$fh_in2>;
    
19. my ($current_scaffold2, $current_start2, $current_end2) = split /\s+/, $line;
    
20. 
    
21. while (1) {
    
22.     unless ($current_scaffold1 eq $current_scaffold2) {
    
23.         say $fh_out join "\t", $current_scaffold1, $current_start1, $current_end1;
    
24. 
    
25.         last if eof $fh_in1;
    
26. 
    
27.         $line = <$fh_in1>;
    
28.         ($current_scaffold1, $current_start1, $current_end1) = split /\s+/, $line;
    
29.         next;
    
30.     }
    
31. 
    
32.     if ($current_start2 > $current_end1) {
    
33.         say $fh_out join "\t", $current_scaffold1, $current_start1, $current_end1;
    
34. 
    
35.         last if eof $fh_in1;
    
36. 
    
37.         $line = <$fh_in1>;
    
38.         ($current_scaffold1, $current_start1, $current_end1) = split /\s+/, $line;
    
39. 
    
40.         next;
    
41.     }
    
42. 
    
43.     if ($current_start1 > $current_end2) {
    
44.         if (eof $fh_in2) {
    
45.             say $fh_out join "\t", $current_scaffold1, $current_start1, $current_end1;
    
46.             last;
    
47.         }
    
48. 
    
49.         $line = <$fh_in2>;
    
50.         ($current_scaffold2, $current_start2, $current_end2) = split /\s+/, $line;
    
51.         next;
    
52.     }
    
53. 
    
54.     say $fh_out join "\t", $current_scaffold1, $current_start1, $current_start2 - 1 unless $current_start1 == $current_start2;
    
55. 
    
56.     if ($current_end1 == $current_end2) {
    
57.         last if eof $fh_in1;
    
58. 
    
59.         $line = <$fh_in1>;
    
60.         ($current_scaffold1, $current_start1, $current_end1) = split /\s+/, $line;
    
61. 
    
62.         # the following code is unneccesary, but you can still uncomment it if you
    
63.         # so desire
    
64.         # if (eof $fh_in2) {
    
65.         #     say $fh_out join "\t", $current_scaffold1, $current_start1, $current_end1;
    
66.         #     last;
    
67.         # }
    
68.         #
    
69.         # $line = <$fh_in2>;
    
70.         # ($current_scaffold2, $current_start2, $current_end2) = split /\s+/, $line;
    
71.         next;
    
72.     }
    
73. 
    
74.     if ($current_end1 < $current_end2) {
    
75.         last if eof $fh_in1;
    
76.         # If you so desire, you can uncomment the following line.
    
77.         # $current_start2 = $current_end1 + 1;
    
78. 
    
79.         $line = <$fh_in1>;
    
80.         ($current_scaffold1, $current_start1, $current_end1) = split /\s+/, $line;
    
81.         next;
    
82.     }
    
83. 
    
84.     $current_start1 = $current_end2 + 1;
    
85.     # the following code is unneccesary, but you can still uncomment it if you
    
86.     # so desire
    
87.     # if (eof $fh_in2) {
    
88.     #     say $fh_out join "\t", $current_scaffold1, $current_start1, $current_end1;
    
89.     #     last;
    
90.     # }
    
91.     #
    
92.     # $line = <$fh_in2>;
    
93.     # ($current_scaffold2, $current_start2, $current_end2) = split /\s+/, $line;
    
94. }
    
95. 
    
96. close $fh_in1;
    
97. close $fh_in2;
    
98. close $fh_out;
    

复制代码

Monox

回复 55# Monox


    用前面提供的file1和file2得到的输出：

1. scaffold1        245        290
   
2. scaffold1        676        718
   
3. scaffold1        1027        1112
   
4. scaffold1        1330        1331
   
5. scaffold1        1647        1652
   
6. scaffold1        1911        2033
   
7. scaffold1        2285        2510
   
8. scaffold1        2672        3005
   
9. scaffold1        3127        3844
   
10. scaffold1        3945        3946
    
11. scaffold1        4080        4172
    
12. scaffold1        4353        5019
    
13. scaffold1        5126        5129
    
14. scaffold1        5230        5967
    
15. scaffold1        6133        7598
    
16. scaffold1        7703        8313
    
17. scaffold1        8581        8656
    
18. scaffold1        8814        9155
    
19. scaffold1        9377        9724
    
20. scaffold1        11200        11205
    
21. scaffold1        13849        13880
    
22. scaffold1        13981        13985
    
23. scaffold1        14149        14172
    
24. scaffold1        14476        14758
    
25. scaffold1        15305        15428
    
26. scaffold1        15607        15691
    
27. scaffold1        15798        15897
    
28. scaffold1        16172        16439
    
29. scaffold1        18075        18655
    
30. scaffold1        22506        22525
    
31. scaffold1        53109        55597
    

复制代码

pitonas

{:2_172:}  先收藏，太感谢了。

回复 56# Monox


   

expert1

回复 1# felix0608

1. 
   
2. awk 'NR==FNR{
   
3.         a[$1] = a[$1] > $3 ? a[$1] : $3
   
4.         if(b[$1]) b[$1] = b[$1] < $2 ? b[$1] : $2
   
5.         else b[$1] = $2
   
6.         for(i = $2 ; i <= $3 ; i ++) m[$1" "i] = 1
   
7. 
   
8. }NR>FNR{
   
9.         for(i = $2 ; i <= $3 ; i ++) delete m[$1" "i]
   
10.         }
    
11. END{
    
12.         for(i in a)
    
13.                 {
    
14.                         p = 0
    
15.                         k = 0
    
16. 
    
17.                      for(j = b[i] ; j <= a[i] ; j ++)
    
18. 
    
19.                    {
    
20.                         if(!(i" "j in m)) continue
    
21.                         if(j == p + k)
    
22.                                 k ++
    
23.                         else
    
24.                            {
    
25.                                 if(k)
    
26.                                         print i ,p ,p + k - 1
    
27.                                 p = j
    
28.                                 k = 1
    
29.                             }
    
30.                   }
    
31.                 print i ,p ,p + k - 1
    
32. 
    
33.    }
    
34. }' f1 f2
    
35. 
    

复制代码

expert1

回复 57# pitonas

回复 1# felix0608

[root@221 ~]# cat file1
scaffold1       1       16439
scaffold1       18075   18656
scaffold1       22506   22691
scaffold1       53109   55597
scaffold1       64277   76351
scaffold1       84997   85240
scaffold1       91210   98127
scaffold1       98569   99022
scaffold1       101188  102851
scaffold1       105232  105249
scaffold1       113602  116343
scaffold1       124999  125184
scaffold1       125930  126558
scaffold1       134339  134668
scaffold1       141085  141482
scaffold1       141576  142066
scaffold1       142183  142295
scaffold1       142368  142817
scaffold1       143613  145919
scaffold1       146106  146119
scaffold1       146456  146501
scaffold1       147426  147715
scaffold1       148846  149512
scaffold1       149776  149813
scaffold1       150361  150747
scaffold1       158353  160761
scaffold1       165280  165383
scaffold1       165486  165737
scaffold1       166810  173128
scaffold1       174606  175266
scaffold1       176688  179222
scaffold1       183911  184183
scaffold1       186787  188238
scaffold1       190183  191560
scaffold1       192206  195941
scaffold1       202089  202265
scaffold1       202944  203048
scaffold1       208786  208979
scaffold1       211069  211420
scaffold1       212044  212069
scaffold1       217631  217813
scaffold1       218313  219423
scaffold1       222785  226870
a 10 25
a 22 30
a 32 40
c 100 2300
d 1 200
[root@221 ~]#


[root@221 ~]#
[root@221 ~]# cat file2
scaffold1       1       244
scaffold1       291     675
scaffold1       719     1026
scaffold1       1113    1329
scaffold1       1332    1646
scaffold1       1653    1910
scaffold1       2034    2284
scaffold1       2511    2671
scaffold1       3006    3126
scaffold1       3845    3944
scaffold1       3947    4079
scaffold1       4173    4352
scaffold1       5020    5125
scaffold1       5130    5229
scaffold1       5968    6132
scaffold1       7599    7702
scaffold1       8314    8580
scaffold1       8657    8813
scaffold1       9156    9376
scaffold1       9725    11199
scaffold1       11206   13848
scaffold1       13881   13980
scaffold1       13986   14148
scaffold1       14173   14475
scaffold1       14759   15304
scaffold1       15429   15606
scaffold1       15692   15797
scaffold1       15898   16171
scaffold1       16863   16981
scaffold1       18656   18756
scaffold1       18845   18949
scaffold1       19167   19272
scaffold1       20314   20413
scaffold1       20431   20542
scaffold1       20692   20791
scaffold1       20810   20918
scaffold1       21789   21906
scaffold1       21945   22461
scaffold1       22526   22776
scaffold1       22846   23126
scaffold1       23162   23755
scaffold1       23882   24047
scaffold1       24714   24813
a 15 25
a 26 28
c 150 350
c 200 300
c 202 348
d 10 25
d 5 25
a 247 249
a 1024 1078
[root@221 ~]#

根据5,7楼的数据，外加一些测试数据，结果:

a 10 14
a 29 30
a 32 40
c 100 149
c 351 2300
scaffold1 245 290
scaffold1 676 718
scaffold1 1027 1112
scaffold1 1330 1331
scaffold1 1647 1652
scaffold1 1911 2033
scaffold1 2285 2510
scaffold1 2672 3005
scaffold1 3127 3844
scaffold1 3945 3946
scaffold1 4080 4172
scaffold1 4353 5019
scaffold1 5126 5129
scaffold1 5230 5967
scaffold1 6133 7598
scaffold1 7703 8313
scaffold1 8581 8656
scaffold1 8814 9155
scaffold1 9377 9724
scaffold1 11200 11205
scaffold1 13849 13880
scaffold1 13981 13985
scaffold1 14149 14172
scaffold1 14476 14758
scaffold1 15305 15428
scaffold1 15607 15691
scaffold1 15798 15897
scaffold1 16172 16439
scaffold1 18075 18655
scaffold1 22506 22525
scaffold1 53109 55597
scaffold1 64277 76351
scaffold1 84997 85240
scaffold1 91210 98127
scaffold1 98569 99022
scaffold1 101188 102851
scaffold1 105232 105249
scaffold1 113602 116343
scaffold1 124999 125184
scaffold1 125930 126558
scaffold1 134339 134668
scaffold1 141085 141482
scaffold1 141576 142066
scaffold1 142183 142295
scaffold1 142368 142817
scaffold1 143613 145919
scaffold1 146106 146119
scaffold1 146456 146501
scaffold1 147426 147715
scaffold1 148846 149512
scaffold1 149776 149813
scaffold1 150361 150747
scaffold1 158353 160761
scaffold1 165280 165383
scaffold1 165486 165737
scaffold1 166810 173128
scaffold1 174606 175266
scaffold1 176688 179222
scaffold1 183911 184183
scaffold1 186787 188238
scaffold1 190183 191560
scaffold1 192206 195941
scaffold1 202089 202265
scaffold1 202944 203048
scaffold1 208786 208979
scaffold1 211069 211420
scaffold1 212044 212069
scaffold1 217631 217813
scaffold1 218313 219423
scaffold1 222785 226870
d 1 4
d 26 200
[root@221 ~]#

   
   

expert1

再免费送一个区间合并的，懒得贴了，感兴趣可以去看这个，
http://blog.chinaunix.net/uid-20778583-id-3916255.html