按一列分组，根据该列查找最小值和最大值并计算滑动窗口 awk 的平均值

Question

我正在对一些基因组数据进行滑动窗口分析，并使用 awk 来完成。

• 第 1 列：表示染色体名称。
• 第 2 列和第 3 列：表示位置（基于 0 的符号）。
• 第 4 列和第 5 列：感兴趣的值。

我需要做什么来获取 X 行（窗口大小为 X）并且：

1. 获取第 1 列的所有不同值（即此窗口内的所有不同染色体）
2. 为第 1 列中的每个不同值获取第 2 列的最小值（即每个染色体的最小值）
3. 为第 1 列中的每个不同值获取第 3 列的最大值（即，对于每个染色体，最大值）
4. 第 4 列的平均值。
5. 第 5 列的平均值。

然后开始相同的过程（取 X 行），但在下面开始 Y 行（即 window step=Y）。

例如，假设我使用 10 (X=10) 的窗口大小和 5 (Y=5) 的窗口步长。我的数据示例可能是：

A   85  86  .043    .021                
A   86  87  .031    .014                
A   87  88  .035    .016                
A   88  89  .033    .015                
A   89  90  .031    .014                
A   90  91  .031    .014                
A   91  92  .032    .015                
B   1   2   .030    .013                
B   2   3   .038    .018                
B   3   4   .032    .014                
B   4   5   .030    .013                
B   5   6   .034    .016                
B   6   7   .032    .015                
B   10  11  .033    .015                
B   11  12  .045    .022                
B   12  13  .055    .029                
B   13  14  .059    .032                
B   20  21  .058    .031                
B   22  23  .059    .031                
B   24  25  .064    .035                
B   26  27  .063    .034                
B   28  29  .058    .031                
B   30  31  .063    .034                
B   108 109 .063    .034                
B   109 110 .067    .037                
B   110 111 .066    .037                
B   111 112 .061    .033                
B   112 113 .061    .033                
B   113 114 .056    .029                
B   114 115 .058    .031

我想要的输出是：

10  A_85_92;B_1_4   0.0336  0.0154
15  A_90_91;B_1_12  0.0337  0.0155
20  B_4_25  0.0469  0.0239
25  B_12_110    0.0609  0.0328
30  B_26_115    0.0616  0.0333 

这个输出有：

• 第一列：窗口标识符。可以是窗口的第一个或最后一个 NR。在这种情况下是最后一个 NR，即 10、15 等。-
• 该窗口内的不同染色体，以及每个染色体内的开始（最小）和结束（最大）位置。在此示例中，对于第一个窗口（从第 1 行到第 10 行），我有 A 和 B 染色体，A 从位置 85 到位置 92，B 从位置 1 到 4。使用“_”和染色体分隔位置使用“;”，但这不是强制性的。
• 在 $4 和 $5 列上计算的平均值。在我们的示例中，分别为 0.0336 和 0.0154。

下一个窗口从第 5 行转到第 15 行，然后从 10 转到第 20 行，然后从 15 转到第 25 行，以此类推。

到现在为止：

    awk -v OFS="\t" 'BEGIN{window=10;step=5}  
{
mod=NR%window; if(NR<=window){count++}
else
{ N[$1]++;{min=$2}{if ($2 < min) min = $2};{max=$3}{if ($3 > max) max = $3}; sum1-=array1[mod]; sum2-=array2[mod]}
sum1+=$4;
sum2+=$5;
array1[mod]=$4;
array2[mod]=$5;
} 
(NR%slide)==0{for (p in N) print NR,p, max, min, sum1/count, sum2/count}'  toy

结果：

15  B   11  12  0.0337  0.0155
20  B   24  25  0.0469  0.0239
25  B   109 110 0.0609  0.0328
30  B   114 115 0.0616  0.0333

因此，我无法正确获取窗口内的所有染色体以及每个染色体的最小值和最大值。另外我正在使用应该是 10 的第一个窗口，但我不知道为什么。

任何输入？提前致谢

Answer 1

$ cat tst.awk
BEGIN {
    winSize = 10
    winStep = 5
    OFS = "\t"
}
{ buf[NR % winSize] = $0 }
(NR >= winSize) && ((NR % winStep) == 0) { prt() }

function prt(   sum,f,i,idx,beg,end,prev,ranges) {
    for (i=1; i<=winSize; i++) {
        idx = (NR+i) % winSize
        split(buf[idx],f)

        if ( f[1] != prev ) {
            ranges = (i > 1 ? ranges end ";" : "") f[1] "_" f[2] "_"
            prev = f[1]
        }
        end = f[3]

        sum[4] += f[4]
        sum[5] += f[5]
    }

    print NR, ranges end, sum[4] / winSize, sum[5] / winSize
}

.

$ awk -f tst.awk file
10  A_85_92;B_1_4   0.0336  0.0154
15  A_90_92;B_1_12  0.0337  0.0155
20  B_4_25  0.0469  0.0239
25  B_12_110    0.0609  0.0328
30  B_26_115    0.0616  0.0333

Answer 2

尽管 Ed Morton 已经给出了一个公认的解决方案，但我还是想分享我的解决方案。

不同之处在于，最小/最大/平均值是针对第一列中的值计算的，而不是针对最后 x 行的组。

输出：

5 A_85_90   0.0346 0.016
10 A_90_92   0.0126 0.0058
10 A_85_92   0.00063 0.00109
10 B_1_4   0.02 0.009
15 B_4_12   0.0348 0.0162
20 B_12_25   0.059 0.0316
20 B_4_25   0.00295 0.00239
25 B_26_110   0.0628 0.034
30 B_110_115   0.0604 0.0326
30 B_26_115   0.00302 0.00333

脚本：

function init(x) {
        pmin[x]=min[x];
        pmax[x]=max[x];
        ps4[x]=s4[s];
        ps5[x]=s5[x];
        min[x]=maxvalue;
        max[x]=-maxvalue;
        s4[x]=0;
        s5[x]=0
}
function calcWindow(x) {
        if ((pmin[x]!=maxvalue) && pmin[x]!="") {
                cmin=pmin[x]<min[x] ? pmin[x] : min[x];
                cmax=pmax[x]>max[x] ? pmax[x] : max[x];
                cs4=(ps4[x]+s4[x])/window;
                cs5=(ps5[x]+s5[x])/window;
                print NR, x "_" cmin "_" cmax, " ", cs4/window, cs5/window ;
        }
}
BEGIN {
        maxvalue=999999;
        window=10;
        windowstep=5;
}
{
        if (!($1 in min)) { init($1) }
        if ($2<min[$1]) { min[$1]=$2 }
        if ($3>max[$1]) { max[$1]=$3 }
        s4[$1]+=$4;
        s5[$1]+=$5;
}
NR%windowstep==0{
        for (i in min) {
                if (min[i]!=maxvalue) {
                        print NR, i "_" min[i] "_" max[i], " ", s4[i]/windowstep, s5[i]/windowstep ;
                        if (NR%window==0) calcWindow(i);
                        init(i);
                }
        }
}