【问题】对区间匹配计算

huang6894

各位大神，想继续请教一个问题：

对于下面两个文本，文本1每段是以两个换行符区分的，每段的第一行是该段的总描述，第2行至最后一行是该段所有区间的信息，其中，第一行第1列和第9列是整个文本唯一的，这里作为key值。第2行至最后一行中第1列是区间名称，第2列是区间起始位置，第3列是终止位置。
如果引入文本2：如果文本2第2、5、6列相同时，判断第3列数字所在文本1的区间（第2列等于文本1每段第一行的第2列且该数字位于后面区间信息的某行的起始终止位置间的），如果同一个区间内有66%以上的饱和度（比如长度为100的区间，有66个以上的点在这个区间）即输出文本二所在该区间内的最小的点至最大的点以及对应文本1的key值和区间名称。如果多个区域的坐标数字是连续的，也需要用:合并。比如同样key下同时出现了两个连续区间比如CDS6和intron6需要表示为CDS6：intron6....


文本1：

1. 
   
2. NM_001145717    chr6    +       36238237        36276372        36238237        36275493        8       PNPLA1  cmpl    cmpl
   
3. CDS1    36238437        36238441        0       1       1
   
4. intron1 36238442        36259446        .       205     206
   
5. CDS2    36259447        36259329        1       2       206
   
6. intron2 36259330        36260837        .       438     439
   
7. CDS3    36260838        36260903        0       3       439
   
8. intron3 36260904        36261966        .       504     505
   
9. CDS4    36261967        36262176        0       4       505
   
10. intron4 36262177        36263140        .       714     715
    
11. CDS5    36263141        36263201        0       5       715
    
12. intron5 36263202        36269637        .       775     776
    
13. CDS6    36269638        36270246        1       6       776
    
14. intron6 36270247        36274068        .       1384    1385
    
15. CDS7    36274069        36274153        1       7       1385
    
16. intron7 36274154        36275363        .       1469    1470
    
17. CDS8    36275364        36275493        2       8       1470
    
18. 3-UTR   36275494        36276372        .       8       *1      *879
    
19. 
    
20. NM_002364       chrX    +       2336324        2336356        2336328        2336357        2       MAGEB2  cmpl    cmpl
    
21. 5-UTR   2336324        2336326        .       1       -103    -6
    
22. intron1 2336327       2336329       .       -6      -5
    
23. 5-UTR   2336330        2336337        .       2       -5      -1
    
24. CDS1    2336338        2336347        0       2       1
    
25. 3-UTR   2336348        2336356        .       2       *1      *549
    
26. 
    
27. NM_019056       chrX    -       47001615        47004609        47001716        47004078        3       NDUFB11 cmpl    cmpl
    
28. 5-UTR   47004079        47004609        .       1       -1      -531
    
29. CDS1    47003872        47004078        0       1       1
    
30. intron1 47002144        47003871        .       208     207
    
31. CDS2    47001983        47002143        0       2       208
    
32. intron2 47001840        47001982        .       369     368
    
33. CDS3    47001716        47001839        2       3       369
    
34. 3-UTR   47001615        47001715        .       3       *101    *1

复制代码

文本2：

1. 
   
2. 16642     chr6    36238435 18      HOM     del
   
3. 16642     chr6    36238437 17      HOM     del
   
4. 16642     chr6    36238438 18      HOM     del
   
5. 16642     chr6    36238439 17      HOM     del
   
6. 16642     chr6    36238440 17      HOM     del
   
7. 16642     chr6    36238441 24      HOM     del
   
8. 16642     chr6    36238442 24      HOM     del
   
9. 16642     chr6    36238443 24      HOM     del
   
10. 16642     chr6    36238444 24      HOM     del
    
11. 16642     chr6    36238445 24      HOM     del
    
12. 16642     chr6    36238446 24      HET     del
    
13. 16642     chr6    3623844724      HET     del
    
14. 16642     chrX    2336324 24      HET     del
    
15. 16642     chrX    2336325 24      HOM     del
    
16. 16642     chrX    2336326 24      HOM     del
    
17. 16642     chrX    2336327 24      HOM     del
    

复制代码

结果：

1. 
   
2. 16642     chr6    36238436-36238445     HOM     del  NM_001145717   PNPLA1 CDS1:intron1
   
3. 16642     chrX    2336325-2336326        HOM     del  NM_002364   MAGEB2  5-UTR

复制代码

我试图用perl把文本1存到一个二维hash里面并把每个区间的长度存到另一个hash里面，然后读入文本2，每一行去获取区间，遇到一个+1，最后和区间的长度比较，如果满足条件（66%）的话，输出，可是最后什么结果都没有，而且吃了很多资源和时间（文本1：58M；文本2：12M），想向大神们求个解决方案，感谢感谢了。。。

yestreenstars

Debug by yourself.

1. #!/usr/bin/perl
   
2. use warnings;
   
3. use strict;
   
4. 
   
5. my %hash1;
   
6. open(FILE1, "b") or die;
   
7. push(@{$hash1{(split)[1]}{(split)[4]}{(split)[5]}}, [split]) while <FILE1>;
   
8. close(FILE1);
   
9. 
   
10. my($chr, %hash2, @result);
    
11. open(FILE2, "a") or die;
    
12. while(<FILE2>){
    
13.         my @field = split;
    
14.         $chr = $field[1] and $hash2{$chr} = join("\t", (split)[0, 8]) and next if $#field == 10;
    
15.         next if /^\s*$/;
    
16.         for my $i (keys %{$hash1{$chr}}){
    
17.                 for my $j (keys %{$hash1{$chr}{$i}}){
    
18.                         my @num;
    
19.                         while(@{$hash1{$chr}{$i}{$j}}){
    
20.                                 my $temp = $hash1{$chr}{$i}{$j}[0][2];
    
21.                                 shift @{$hash1{$chr}{$i}{$j}} and next if $temp < $field[1];
    
22.                                 last if $temp > $field[2];
    
23.                                 push(@num, $temp);
    
24.                                 shift @{$hash1{$chr}{$i}{$j}};
    
25.                         }
    
26.                         push(@result, [16642, $chr, @num[0, -1], $i, $j, $hash2{$chr}, $field[0]]) if @num && ($num[-1] - $num[0] + 1) / ($field[2] - $field[1] + 1) > 0.66;
    
27.                 }
    
28.         }
    
29. }
    
30. close(FILE2);
    
31. 
    
32. while(@result){
    
33.         if($result[1]){
    
34.                 if(join("\t", @{$result[0]}[1, 4, 5]) eq join("\t", @{$result[1]}[1, 4, 5])){
    
35.                         if($result[0][3] + 1 == $result[1][2]){
    
36.                                 $result[1][2] = $result[0][2];
    
37.                                 $result[1][-1] = "$result[0][-1]:$result[1][-1]";
    
38.                         }else{
    
39.                                 print join("\t", join("\t", @{$result[0]}[0, 1]), "$result[0][2]-$result[0][3]", join("\t", @{$result[0]}[4..7])), $/;
    
40.                         }
    
41.                         shift @result;
    
42.                         next;
    
43.                 }
    
44.         }
    
45.         print join("\t", join("\t", @{$result[0]}[0, 1]), "$result[0][2]-$result[0][3]", join("\t", @{$result[0]}[4..7])), $/;
    
46.         shift @result;
    
47. }

复制代码

1. [root@localhost ~]# cat a
   
2. NM_001145717    chr6    +       36238237        36276372        36238237        36275493        8       PNPLA1  cmpl    cmpl
   
3. CDS1    36238437        36238441        0       1       1
   
4. intron1 36238442        36238446        .       205     206
   
5. CDS2    36259447        36259329        1       2       206
   
6. intron2 36259330        36260837        .       438     439
   
7. CDS3    36260838        36260903        0       3       439
   
8. intron3 36260904        36261966        .       504     505
   
9. CDS4    36261967        36262176        0       4       505
   
10. intron4 36262177        36263140        .       714     715
    
11. CDS5    36263141        36263201        0       5       715
    
12. intron5 36263202        36269637        .       775     776
    
13. CDS6    36269638        36270246        1       6       776
    
14. intron6 36270247        36274068        .       1384    1385
    
15. CDS7    36274069        36274153        1       7       1385
    
16. intron7 36274154        36275363        .       1469    1470
    
17. CDS8    36275364        36275493        2       8       1470
    
18. 3-UTR   36275494        36276372        .       8       *1      *879
    
19. 
    
20. NM_002364       chrX    +       2336324        2336356        2336328        2336357        2       MAGEB2  cmpl    cmpl
    
21. 5-UTR   2336324        2336326        .       1       -103    -6
    
22. intron1 2336327       2336329       .       -6      -5
    
23. 5-UTR   2336330        2336337        .       2       -5      -1
    
24. CDS1    2336338        2336347        0       2       1
    
25. 3-UTR   2336348        2336356        .       2       *1      *549
    
26. 
    
27. NM_019056       chrX    -       47001615        47004609        47001716        47004078        3       NDUFB11 cmpl    cmpl
    
28. 5-UTR   47004079        47004609        .       1       -1      -531
    
29. CDS1    47003872        47004078        0       1       1
    
30. intron1 47002144        47003871        .       208     207
    
31. CDS2    47001983        47002143        0       2       208
    
32. intron2 47001840        47001982        .       369     368
    
33. CDS3    47001716        47001839        2       3       369
    
34. 3-UTR   47001615        47001715        .       3       *101    *1
    
35. [root@localhost ~]# cat b
    
36. 16642     chr6    36238435 18      HOM     del
    
37. 16642     chr6    36238437 17      HOM     del
    
38. 16642     chr6    36238438 18      HOM     del
    
39. 16642     chr6    36238439 17      HOM     del
    
40. 16642     chr6    36238440 17      HOM     del
    
41. 16642     chr6    36238441 24      HOM     del
    
42. 16642     chr6    36238442 24      HOM     del
    
43. 16642     chr6    36238443 24      HOM     del
    
44. 16642     chr6    36238444 24      HOM     del
    
45. 16642     chr6    36238445 24      HOM     del
    
46. 16642     chr6    36238446 24      HET     del
    
47. 16642     chr6    36238447 24      HET     del
    
48. 16642     chrX    2336324 24      HET     del
    
49. 16642     chrX    2336325 24      HOM     del
    
50. 16642     chrX    2336326 24      HOM     del
    
51. 16642     chrX    2336327 24      HOM     del
    
52. [root@localhost ~]# ./p | column -t
    
53. 16642  chr6  36238437-36238445  HOM  del  NM_001145717  PNPLA1  CDS1:intron1
    
54. 16642  chrX  2336325-2336326    HOM  del  NM_002364     MAGEB2  5-UTR
    
55. [root@localhost ~]#

复制代码

yestreenstars

解释一下第一条结果的每一列是怎么来的吧

pitonas

现在问题来了喔。     
问题已经难很了, 然而，意想不到的事发生了. .

理解问题比回答这个问题更困难 ...

1. 问题本困难，你又不举例，
   
2. 描述又乱写，让人费疑猜？

复制代码

huang6894

回复 2# yestreenstars


    嗯嗯~结果的1、2、4、5就是文件2的1、2、4、5、6，第三列是文件2第三列数字所在该区域的最小值-最大值，第7列是该区域的key值也就是文本1上该区域的第一行第一列，第8列是第二个key值也就是文本1上该区域的第一行第9列，第十列是落在的区域，如果多个区域，用：合并~

比如：6642     chr6    36238436-36238445     HOM     del  NM_001145717   PNPLA1 CDS1:intron1

因为文本2上：
16642     chr6    36238435 18      HOM     del    第二列是chr6，排除掉文本1段落第一行第二列不是chr6的，然后36238435未落在任何区域，此时36238435上$h{chr6}{HON}{del}的覆盖度为1
16642     chr6    36238437 17      HOM     del    第二列是chr6，排除掉文本1段落第一行第二列不是chr6的，然后36238437落在CDS1，此时CDS1上$h{chr6}{HON}{del}覆盖度为1
16642     chr6    36238438 18      HOM     del    第二列是chr6，排除掉文本1段落第一行第二列不是chr6的，然后36238438落在CDS1，此时CDS1上$h{chr6}{HON}{del}覆盖度为2
16642     chr6    36238439 17      HOM     del    第二列是chr6，排除掉文本1段落第一行第二列不是chr6的，然后36238439落在CDS1，此时CDS1上$h{chr6}{HON}{del}覆盖度为3
16642     chr6    36238440 17      HOM     del    第二列是chr6，排除掉文本1段落第一行第二列不是chr6的，然后36238440落在CDS1，此时CDS1上$h{chr6}{HON}{del}覆盖度为4
16642     chr6    36238441 24      HOM     del    第二列是chr6，排除掉文本1段落第一行第二列不是chr6的，然后36238441落在CDS1，此时CDS1上$h{chr6}{HON}{del}覆盖度为5
16642     chr6    36238442 24      HOM     del     第二列是chr6，排除掉文本1段落第一行第二列不是chr6的，然后36238442落在intron1，此时CDS1上$h{chr6}{HON}{del}覆盖度为1
16642     chr6    36238443 24      HOM     del     第二列是chr6，排除掉文本1段落第一行第二列不是chr6的，然后36238442落在intron1，此时CDS1上$h{chr6}{HON}{del}覆盖度为2
16642     chr6    36238444 24      HOM     del     第二列是chr6，排除掉文本1段落第一行第二列不是chr6的，然后36238442落在intron1，此时CDS1上$h{chr6}{HON}{del}覆盖度为3
16642     chr6    36238445 24      HOM     del    第二列是chr6，排除掉文本1段落第一行第二列不是chr6的，然后36238445落在intron1，此时CDS1上$h{chr6}{HON}{del}覆盖度为4
16642     chr6    36238446 24      HET     del    第二列是chr6，排除掉文本1段落第一行第二列不是chr6的，然后36238446落在intron1，此时CDS1上$h{chr6}{HON}{del}覆盖度为1
16642     chr6    36238447 24      HET     del    第二列是chr6，排除掉文本1段落第一行第二列不是chr6的，然后36238447落在CDS2，此时CDS2上$h{chr6}{HON}{del}覆盖度为1
--------------------------------------------------------------------------------
CDS1的长度为36238437-36238441 ：5，$h{chr6}{HON}{del}覆盖度为5；符合条件
intron1的长度36238442-36259446:5,$h{chr6}{HON}{del}覆盖度为4；符号条件
考虑到CDS1和intron1是相邻的，合并之
所以输出6642     chr6    36238436-36238445     HOM     del  NM_001145717   PNPLA1 CDS1:intron1

huang6894

回复 3# pitonas


    对不起。。。迷迷糊糊的没说清楚

yestreenstars

擦，文本2被改了，早上不是这样的！

huang6894

回复 6# yestreenstars


    哈哈。。一贯作风，大神原谅我

yestreenstars

请问36238436是怎么来的？

huang6894

回复 8# yestreenstars


    36238436===>36238437...失误失误。嘿嘿