awk 查找重叠

Question

我有一个包含如下所示列的文件。

Group   Start        End
chr1    117132092    118875009
chr1    117027758    119458215
chr1    103756473    104864582
chr1    105093795    106219211
chr1    103354114    104747251
chr1    102741437    105235140
chr1    100090254    101094139
chr1    100426977    101614730
chr2    86644663     87767193
chr2    82473711     83636545
chr2    83896702     85079032
chr2    83876122     85091910
chr2    82943211     84350917
chr3    89410051     90485635
chr3    89405753     90485635
chr3    86491492     87593215
chr3    82507157     83738004
chr3    85059618     86362254

我想找到每组中这些坐标之间的重叠（按 chr1、chr2、chr3.. 分组）。

如果与同一组中的其他人至少有 50% 的重叠，则必须检查开始和结束坐标。如果至少有 50% 的重叠，则必须在第 3 列和第 4 列（这是重叠区域的范围）中报告新的开始和结束坐标。如果它们不重叠，则必须在第 3 列和第 4 列中报告原始开始和结束。

为了更清楚，让我们取前两行

                 117132092..........118875009
         117027758...........................119458215

由于它们彼此重叠至少 50%，因此重叠的范围在输出中报告为新开始和新结束。并且第 3 行和第 4 行不与其他行重叠，因此原始坐标在第 3 列和第 4 列中报告为新起点和新终点。同样，由于第 5 行和第 6 行彼此重叠 50%，因此它们的范围报告为新坐标第 3 列和第 4 列中的开始和新结束。 这是预期的输出：

Group   Start     End         NewStart   NewEnd   
chr1 117132092 118875009  117027758   119458215
chr1 117027758 119458215  117027758   119458215
chr1 103756473 104864582  103354114   104864582
chr1 105093795 106219211  105093795   106219211
chr1 103354114 104747251  102741437   105235140
chr1 102741437 105235140  102741437   105235140
chr1 100090254 101094139  100090254   101614730
chr1 100426977 101614730  100090254   101614730
chr2 86644663 87767193    86644663    87767193
chr2 82473711 83636545    82473711    83636545 
chr2 83896702 85079032    83876122    85091910
chr2 83876122 85091910    83876122    85091910
chr2 82943211 84350917    82943211    84350917
chr3 89410051 90485635    89405753    90485635
chr3 89405753 90485635    89405753    90485635
chr3 86491492 87593215    86491492    87593215
chr3 82507157 83738004    82507157    83738004
chr3 85059618 86362254    85059618    86362254

我已经用 R 编程语言实现了这一点，但是原始文件太大并且需要很长时间才能运行。有人可以帮助在 awk 中做到这一点。

Answer 1

使用 Gnu Awk 版本 4，您可以尝试：

gawk -f a.awk file file

a.awk 在哪里：

NR==FNR {
    if (FNR>1) {
        a[$1][++i]=$2
        b[$1][i]=$3
    }
    next
}
FNR==1 {
    fmt="%-7s%-10s%-10s%-10s%-10s\n"
    printf fmt,"Group","Start","End","NewStart","NewEnd" 
}
FNR>1{
    $4=$2; $5=$3
    n=checkInside($1,$2,$3)
    if (n>0) {
        ff=0; x=$2; y=$3
        for (i=1; i<=n; i++) {
            ar=a[$1][R[i]]; br=b[$1][R[i]];
            getIntersect($2,$3,ar,br)
            getLargest($2,$3,ar,br)
            ovl=((i2-i1)/($3-$2))*100;
            ovr=((i2-i1)/(br-ar))*100;
            if (ovl>50 && ovr>50) {
                if (r1<x) x=r1
                if (r2>y) y=r2
                ff=1
            }
        }
        if (ff) {
            $4=x; $5=y
        }
    }
    printf fmt,$1,$2,$3,$4,$5
}

function getLargest(x1,y1,x2,y2) {
    r1=(x1<=x2)?x1:x2
    r2=(y1>=y2)?y1:y2
}

function getIntersect(x1,y1,x2,y2) {
    if (x1>=x2 && x1<=y2) {
        i1=x1;
    } else {
        i1=x2;
    }
    i2=(y1<=y2)?y1:y2
}

function checkInside(g,x,y,i,j,x1,y1) {
    R["x"]=0
    for (i in a[g]) {
        x1=a[g][i]; y1=b[g][i];
        if ((x>=x1 && x<=y1) || (y>=x1 && y<=y1)) {
            if (!(x==x1 && y==y1))
                R[++j]=i
        }
    }
    return j
}

输出：

Group  Start     End       NewStart  NewEnd    
chr1   117132092 118875009 117027758 119458215 
chr1   117027758 119458215 117027758 119458215 
chr1   103756473 104864582 103354114 104864582 
chr1   105093795 106219211 105093795 106219211 
chr1   103354114 104747251 102741437 105235140 
chr1   102741437 105235140 102741437 105235140 
chr1   100090254 101094139 100090254 101614730 
chr1   100426977 101614730 100090254 101614730 
chr2   86644663  87767193  86644663  87767193  
chr2   82473711  83636545  82473711  83636545  
chr2   83896702  85079032  83876122  85091910  
chr2   83876122  85091910  83876122  85091910  
chr2   82943211  84350917  82943211  84350917  
chr3   89410051  90485635  89405753  90485635  
chr3   89405753  90485635  89405753  90485635  
chr3   86491492  87593215  86491492  87593215  
chr3   82507157  83738004  82507157  83738004  
chr3   85059618  86362254  85059618  86362254