对重叠区间交集进行合并处理

jason680

回复 20# expert1

that's your fault not my with learning by analogy ...

  my key is k=sprintf("%s_%020d",$5,$2);

and get result as below:
# cat file1
A 10 30
A 20 40
A 50 70
A1 38 168

# awk '{k=sprintf("%s_%020d",$1,$2);for(n=0;n++<3;)a[k,n]=$n;c[k]}END{$1="";t=asorti(c,s);for(n=0;n++<t+1;){N=s[n];if($1==a[N,1]&&$3>=a[N,2]){a[N,2]=$2;if($3>a[N,3])a[N,3]=$3}else if($1)print $1,$2,$3;for(m=0;m++<3;)$m=a[N,m]}}' file1
A1 38 168
A 10 40
A 50 70

expert1

回复 21# jason680


    学习了，阁下思路奇特，自愧不如啊

blackold

还有一种情况没有考虑——包含，比如[1,40] [10,20]。

expert1

回复 23# blackold


    啥意思黑哥？

blackold

回复 24# expert1


    没仔细看，后面已经修正了。

expert1

回复 25# blackold


    针对19楼我自己的问题，如下，先sort, 由于坐标$2<$3是永远成立的，那么只需要

1. awk 'NR==1{tag=$1;start=$2;end=$3}NR>1{
   
2. 
   
3.         if($1==tag){
   
4.                                 if($2<end && $3>=end)
   
5. 
   
6.                                 end=$3
   
7. 
   
8.                                else if($3<end)next
   
9. 
   
10.                               else{                                                      
    
11.                         
    
12.                                 print tag,start,end;start=$2;end=$3;tag=$1
    
13.                                        
    
14.                                 }
    
15.                         
    
16.                         
    
17.                 }
    
18. 
    
19. 
    
20.         else {                                                   
    
21.                 print tag,start,end;tag=$1;start=$2;end=$3  
    
22.                      
    
23.                         }
    
24.                         
    
25.         }END{print tag,start,end}'

复制代码

expert1

回复 7# huang6894


    首先sort -k5 -k2 file ,使$2（坐标起点）升序排列，由于$2<$3(起点<终点）这个是恒成立的。所以按照我的这个思路第N+1行和第N行比较，然后合并。

1. # awk -v OFS="\t" 'NR==1{tag=$5;start=$2;end=$3;a=$6;b=$4;c=$1}NR>1{
   
2. 
   
3.         if($5==tag){
   
4.                                 if($2<end && $3>end)
   
5. 
   
6.                                       end=$3
   
7. 
   
8.                                else if($3<end) next
   
9. 
   
10.                                else{                                                      
    
11.                         
    
12.                                    print ">"tag"_"++count,a,c,b,start,end;
    
13. 
    
14.                                    start=$2;end=$3;tag=$5;a=$6;b=$4;c=$1
    
15.                                        
    
16.                                 }
    
17.                         
    
18.                         
    
19.                 }
    
20. 
    
21. 
    
22.         else {                                             
    
23.                  print ">"tag"_"++count,a,c,b,start,end;
    
24. 
    
25. tag=$5;start=$2;end=$3;a=$6;b=$4;c=$1
    
26. count=0      
    
27.                      
    
28.                         }
    
29.                         
    
30.         }END{print ">"tag"_"++count,a,c,b,start,end}' <<<"refseq2 279363  280352  -       NM_001002919    FAM150B
    
31. refseq2 282911  283375  -       NM_001002919    FAM150B
    
32. refseq2 285923  286403  -       NM_001002919    FAM150B
    
33. refseq2 286090  286543  -       NM_001002919    FAM150B
    
34. refseq2 287383  288092  -       NM_001002919    FAM150B
    
35. refseq2 287813  288508  -       NM_001002919    FAM150B
    
36. refseq1 934142  935012  -       NM_021170       HES4
    
37. refseq1 934706  935193  -       NM_021170       HES4
    
38. refseq1 934872  935367  -       NM_021170       HES4
    
39. refseq1 935046  935752  -       NM_021170       HES4
    
40. refseq1 860921  861380  +       NM_152486       SAMD11
    
41. refseq1 861102  861593  +       NM_152486       SAMD11
    
42. refseq1 865335  865916  +       NM_152486       SAMD11
    
43. refseq1 866219  866669  +       NM_152486       SAMD11
    
44. refseq1 870952  871476  +       NM_152486       SAMD11
    
45. refseq1 874220  874709  +       NM_152486       SAMD11
    
46. refseq1 874455  875040  +       NM_152486       SAMD11
    
47. refseq1 879088  880161  +       NM_152486       SAMD11"
    

复制代码

结果如下：

>NM_001002919_1 FAM150B refseq2 -       279363  280352
>NM_001002919_2 FAM150B refseq2 -       282911  283375
>NM_001002919_3 FAM150B refseq2 -       285923  286543
>NM_001002919_4 FAM150B refseq2 -       287383  288508
>NM_021170_1    HES4    refseq1 -       934142  935752
>NM_152486_1    SAMD11  refseq1 +       860921  861593
>NM_152486_2    SAMD11  refseq1 +       865335  865916
>NM_152486_3    SAMD11  refseq1 +       866219  866669
>NM_152486_4    SAMD11  refseq1 +       870952  871476
>NM_152486_5    SAMD11  refseq1 +       874220  875040
>NM_152486_6    SAMD11  refseq1 +       879088  880161

2楼也指出来了，为什么那个有重合的没合并呢，我这里是是完全合并了。

expert1

回复 23# blackold


     awk -v OFS="\t" 'NR==1{tag=$5;start=$2;end=$3;a=$6;b=$4;c=$1}NR>1{

        if($5==tag){
                                if($2<end && $3>end)

                                end=$3

                               else if($3<end)next

                               else{                                                      
                        
                                print ">"tag"_"++count,a,c,b,start,end;

start=$2;end=$3;tag=$5;a=$6;b=$4;c=$1
                                       
                                }
                        
                        
                }


        else {                                             
                 print ">"tag"_"++count,a,c,b,start,end;

tag=$5;start=$2;end=$3;a=$6;b=$4;c=$1
count=0      
                     
                        }
                        
        }END{print ">"tag"_"++count,a,c,b,start,end}' <<<"refseq2 279363  280352  -       NM_001002919    FAM150B
refseq2 282911  283375  -       NM_001002919    FAM150B
refseq2 285923  286403  -       NM_001002919    FAM150B
refseq2 286090  286543  -       NM_001002919    FAM150B
refseq2 287383  288092  -       NM_001002919    FAM150B
refseq2 287813  288508  -       NM_001002919    FAM150B
refseq1 934142  935012  -       NM_021170       HES4
refseq1 934706  935193  -       NM_021170       HES4
refseq1 934872  935367  -       NM_021170       HES4
refseq1 935046  935752  -       NM_021170       HES4
refseq1 860921  861380  +       NM_152486       SAMD11
refseq1 861102  861593  +       NM_152486       SAMD11
refseq1 865335  865916  +       NM_152486       SAMD11
refseq1 866219  866669  +       NM_152486       SAMD11
refseq1 870952  871476  +       NM_152486       SAMD11
refseq1 874220  874709  +       NM_152486       SAMD11
refseq1 874455  875040  +       NM_152486       SAMD11
refseq1 879088  880161  +       NM_152486       SAMD11
mnscx 200     240     x       no           abc"c
>NM_001002919_1 FAM150B refseq2 -       279363  280352
>NM_001002919_2 FAM150B refseq2 -       282911  283375
>NM_001002919_3 FAM150B refseq2 -       285923  286543
>NM_001002919_4 FAM150B refseq2 -       287383  288508
>NM_021170_1    HES4    refseq1 -       934142  935752
>NM_152486_1    SAMD11  refseq1 +       860921  861593
>NM_152486_2    SAMD11  refseq1 +       865335  865916
>NM_152486_3    SAMD11  refseq1 +       866219  866669
>NM_152486_4    SAMD11  refseq1 +       870952  871476
>NM_152486_5    SAMD11  refseq1 +       874220  875040
>NM_152486_6    SAMD11  refseq1 +       879088  880161
>equal_1        abc     mnscx   x       1       40
>no_1   abc     mnscx   x       10      40
>no_2   abc     mnscx   x       200     240

测试完全通过

huang6894

回复 27# expert1


    厉害!!!学习了

expert1

哦，拷贝样例的时候忘了这个 refseq1 876324  876886  +       NM_152486       SAMD11了，难怪如此，还有个E。
数据文件加上之后这些时候。结果应该是

>NM_001002919_E1        FAM150B refseq2 -       279363  280352
>NM_001002919_E2        FAM150B refseq2 -       282911  283375
>NM_001002919_E3        FAM150B refseq2 -       285923  286543
>NM_001002919_E4        FAM150B refseq2 -       287383  288508
>NM_021170_E1   HES4    refseq1 -       934142  935752
>NM_152486_E1   SAMD11  refseq1 +       860921  861593
>NM_152486_E2   SAMD11  refseq1 +       865335  865916
>NM_152486_E3   SAMD11  refseq1 +       866219  866669
>NM_152486_E4   SAMD11  refseq1 +       870952  871476
>NM_152486_E5   SAMD11  refseq1 +       874220  875040
>NM_152486_E6   SAMD11  refseq1 +       876324  876886
>NM_152486_E7   SAMD11  refseq1 +       879088  880161