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

Posted 2023-02-14

【中文标题】用awk计算滑动窗口的中位数

【英文标题】：Calculate median of a sliding window with awk

【发布时间】：2020-07-05 00:31:06

【问题描述】：

我需要生成一个包含数百万行的滑动窗口并计算第 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

我找到了以下脚本 here 来做我想做的事，但只是为了平均，而不是中位数：

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

这个脚本用于计算来自here的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 解决方案，但我不知道如何使用滑动窗口有效地工作。

【问题讨论】：

为什么第二列在输出中移动了？第一个中位数和最后一个中位数是如何计算的？哦，前 3 行被忽略了？

datamash 很棒，但它没有窗口；全有或全无。

输出从第 4 行开始，因为第 4 行是可以生成滑动窗口 4 的第一行。最后的中位数是最后一行和前三行的中位数。

【参考方案1】：

将 GNU awk 用于asort()：

$ cat tst.awk
BEGIN 
    OFS = "\t"
    window = 4
    befMid = int(window / 2)
    aftMid = befMid + (window % 2 ? 0 : 1)

 array[NR % window] = $3 
NR >= window 
    asort(array,vals)
    print $1, $2, int( (vals[befMid] + vals[aftMid]) / 2 + 0.5 )

.

$ awk -f tst.awk file
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

【参考方案2】：

以下假设函数 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;
   ' 

【参考方案3】：

以下带有 GNU awk 的脚本似乎生成了您提供的输出：

awk -v OFS='\t' -v window=4 '

    # I store the numbers in an array `nums` indexed with `1 ... window`
    mod = NR % window + 1;
    nums[mod] = $3;


# If the count of numbers is greater or equal the window,
# we can start calculating the median.
NR >= window 

    # Copy the array nums, cause we need to sort it.
    for (i = 1; i <= window; ++i) 
        copy[i] = nums[i];
    

    # Sort the copy.
    # asort is a GNU extension if I remember.
    # For non-gnu, write a sorting function yourself.
    asort(copy);

    # Calculate the median.
    # I hope that is ok.
    half = int( (window + 1) / 2 );
    if (window % 2 == 0) 
        # You seem to want to round 0.5 up.
        # Just add 1 and round down.
        median = int( (copy[half] + copy[half + 1] + 1) / 2 );
     else 
        median = copy[half];
    

    # Output
    print $1, $2, median 
'