Рассчитать медиану скользящего окна с помощью awk - PullRequest
Новая
       80

Рассчитать медиану скользящего окна с помощью awk

5 голосов
/ 24 марта 2020

Мне нужно создать скользящее окно из миллионов строк и рассчитать медиану столбца 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
        }
}
'

, но я не могу заставить их работать вместе. Еще одна проблема заключается в том, что медианный расчет требует отсортированного ввода. Я также нашел это решение datama sh, но я не знаю, как заставить работать эффективно с помощью скользящего окна.

Ответы [ 3 ]

4 голосов
/ 24 марта 2020

Следующее предполагает наличие функции 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;
   }'
4 голосов
/ 24 марта 2020

Использование 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
3 голосов
/ 24 марта 2020

Следующий скрипт с 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 
}'
...