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

huang6894

回复 9# expert1


    感谢！

cao627

@jason680

1. 
   
2. $ cat file1
   
3. a 2 7
   
4. a 4 9
   
5. A 2 7
   
6. A 4 6
   
7. 
   
8. $ awk '{k=$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;else if($1)print $1,$2,$3;for(m=0;m++<3;)$m=a[N,m]}}' file1
   
9. A 2 6  #A 2 7
   
10. a 2 9
    

复制代码

seesea2517

居然还真有一群夜猫。
本猫可要睡了，终于可以睡了。

jason680

回复 12# cao627

Here you are

$ awk '{k=$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
A 2 7
a 2 9

   

expert1

回复 12# cao627


    你这种情况和

http://bbs.chinaunix.net/thread-3620300-2-1.html

很相似啊，但这里是升序排列的。

expert1

回复 8# jason680


    求解释一下。看的费力啊。

jason680

回复 16# expert1

# awk '{
  k=sprintf("%s_%020d",$5,$2);  # for key
  for(n=0;n++<6;)a[k,n]=$n;     
  c[k]                          
  }
  END{
    $5="";
        t=asorti(c,s);                                # sort all data
        for(n=0;n++<t+1;){
          N=s[n];
          if($5==a[N,5]&&$3+1>=a[N,2])     # check overlay area
            a[N,2]=$2;                     # keep the minimum value
          else if($5)
            print ">"$5"_E"++C,$6,$1,$4,$2,$3;
          if($5!=a[N,5]) C=0;              # reset count
          for(m=0;m++<6;)                  # save last data in $1, $2, $3, ...
            $m=a[N,m]
        }
}' FILE

for example
sorted data
refseq1 860921  861380  +       NM_152486       SAMD11  <== last , $1,$2,...
refseq1 861102  861593  +       NM_152486       SAMD11  <== Now  , a[N,1], a[N,2],...

if($5==a[N,5]&&$3+1>=a[N,2])
  $5==a[N,5] is true   , the $5(NM_152486) are the same
  $3+1>=a[N,2] is true , 861380 > 861102
   the end of last one is bigger than start of next one (Now)
   to keep start value of the "last"
   by a[N,2]=$2;
Note: there hav a bug that didn't check $3 for maximum value


refseq1 879088  880161  +       NM_152486       SAMD11  <== last , $1,$2,...
refseq1 935046  935752  -       NM_021170       HES4    <== Now  , a[N,1], a[N,2],...

if($5==a[N,5]&&$3+1>=a[N,2])
$5==a[N,5] is false
and ouput the "last" data with count
else if($5)
    print ">"$5"_E"++C,$6,$1,$4,$2,$3;

  the $5 are different and reset count
  if($5!=a[N,5]) C=0;              # reset count   
   

expert1

回复 17# jason680


    print ">"$5"_E"++

他这个规则咋算的？ 第5列相等一次，++一次吗，？第一次是1，然后第二次相等就2，以此类推？

还有貌似他的数据没得什么规律，只是起始数据$2<$3，貌似你这个有bug,没有判断正确，

$ awk '{k=sprintf("%s_%020d",$5,$2);for(n=0;n++<6a[k,n]=$n;c[k]}END{$5="";t=asorti(c,s);for(n=0;n++<t+1{N=s[n];if($5==a[N,5]&&$3+1>=a[N,2])a[N,2]=$2;else if($5)print ">"$5"_E"++C,$6,$1,$4,$2,$3;if($5!=a[N,5])C=0;for(m=0;m++<6$m=a[N,m]}}' <<<"refseq2 87813  88508  -       NM_001002919    FAM150B
> refseq2 87383  88092  -       NM_001002919    FAM150B
> refseq2 86090  86543  -       NM_001002919    FAM150B
> refseq2 85923  86403  -       NM_001002919    FAM150B
> refseq2 82911  89375  -       NM_001002919    FAM150B
> refseq2  63  280352  -        NM_001002919    FAM150B"
>NM_001002919_E1 FAM150B refseq2 - 63 86543
>NM_001002919_E2 FAM150B refseq2 - 87383 88508

按理应该是63---280352才对。

jason680

回复 18# expert1

1. ++ is for C variable
   print ">"$5"_E"++C   

2. thank you for notice again with bug
  I said it had the bug in 16L(2013-09-17 09:51:26)
  ca0627 said it had some bug in 11L and I solved issue in 13L
  

expert1

回复 19# jason680


    你这个针对cao627的fix. 好像考虑的还是不全，这种问题是非常麻烦的，想了半天都没什么好的思路。

awk '{k=$1""$2;for(n=0;n++<3a[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]}}' <<<"
> A 10 30
> A 20 40
> A 50 70
> A1 38 168"
A 10 30
A1 38 168
A 20 40
A 50 70


我这里是gawk.

按理10-30,20-40应该合并为10-40才对。

这个可以简化为，三个字段，在第一个字段相等的情况下，找出最小/最大坐标（$2<$3，这个应该可以肯定），但数据都没什么顺序和规律（$2,$3不是有规律的排序过的）