时间戳

首页 > 解决方案 > 用awk计算滑动窗口的中位数

bash - 用awk计算滑动窗口的中位数

问题描述

我需要生成一个包含数百万行的滑动窗口并计算第 3 列的中值。我的数据看起来像这样，第 1 列始终相同，第 2 列等于行号，第 3 列是我需要中值的信息为了：

HiC_scaffold_1  1   34
HiC_scaffold_1  2   34
HiC_scaffold_1  3   36
HiC_scaffold_1  4   37
HiC_scaffold_1  5   38
HiC_scaffold_1  6   39
HiC_scaffold_1  7   40
HiC_scaffold_1  8   40
HiC_scaffold_1  9   40
HiC_scaffold_1  10  41
HiC_scaffold_1  11  41
HiC_scaffold_1  12  41
HiC_scaffold_1  13  44
HiC_scaffold_1  14  44
HiC_scaffold_1  15  55

我需要这样的结果，假设滑动窗口为 4 并四舍五入到最接近的整数。在真实数据集中，我可能会使用 1000 的滑动窗口：

HiC_scaffold_1  4   35
HiC_scaffold_1  5   37
HiC_scaffold_1  6   38
HiC_scaffold_1  7   39
HiC_scaffold_1  8   40
HiC_scaffold_1  9   40
HiC_scaffold_1  10  40
HiC_scaffold_1  11  41
HiC_scaffold_1  12  41
HiC_scaffold_1  13  41
HiC_scaffold_1  14  43
HiC_scaffold_1  15  44

我在这里找到了下面的脚本来做我想做的事，但是是平均的，而不是中位数：

awk -v OFS="\t" 'BEGIN {
        window = 4
        slide = 1
}

{
        mod = NR % window
        if (NR <= window) {
                count++
        } else {
                sum -= array[mod]
        }
        sum += $3
        array[mod] = $3
}

(NR % slide) == 0 {
        print $1, NR, sum / count
}
' file.txt

和这个用 awk 从这里计算中位数的脚本：

sort -n -k3 file.txt |
awk '{
        arr[NR] = $3
}

END {
        if (NR % 2 == 1) {
                print arr[(NR + 1) / 2]
        } else {
                print $1 "\t" $2 "\t" (arr[NR / 2] + arr[NR / 2 + 1]) / 2
        }
}
'

但我不能让他们一起工作。另一个问题是中位数计算需要排序输入。我还找到了这个datamash解决方案，但我不知道如何使用滑动窗口有效地工作。

标签： bashawkmediansliding-window

解决方案

下面假设功能的可用性，asort由 GNU awk (gawk) 提供。程序由 wsize 参数化，窗口大小——这里是 4：

gawk -v wsize=4 '
   BEGIN { 
    if (wsize % 2 == 0) { m1=wsize/2; m2=m1+1; } else { m1 = m2 = (wsize+1)/2; } 
   }
   function roundedmedian() {
     asort(window, a);
     return (m1==m2) ? a[m1] : int(0.5 + ((a[m1] + a[m2]) / 2));
   }
   function push(value) {
     window[NR % wsize] = value;
   }
   NR < wsize { window[NR]=$3; next; }
   { push($3);
     $3 = roundedmedian();
     print $0;
   }'

推荐阅读