☰

计算每 2 行的总和，如果总和小于特定值，则将其替换为另一个值

tag-icon • tag-icon python perl bioinformatics
计算每 2 行的总和，如果总和小于特定值，则将其替换为另一个值

我有一个基因型矩阵（带有表格空间），有 200 万行和 12 列。列是个体，行是 SNP。我每个人的每个 SNP 有 2 行，一行是参考等位基因的数量，另一行是替代等位基因的数量（每 2 行对应一个 SNP，这意味着第 1 行和第 2 行对应于 SNP 1，第 3 行和第 4 行对应于 SNP 1） SNP2，第 5 行和第 6 行对应于 SNP 3)。

这是一个示例（2 个 SNP 和 8 个个体）：

head genotype
2   3   1   0   0   3   5   3       
18  15  19  18  16  15  13  17      
2   1   0   0   0   1   1   1           
18  19  18  16  20  17  17  23

对于每个 SNP，如果参考和替代等位基因的总和小于 20，我想用 0 替换这两个等位基因，如果它们等于或大于 20，我想保留它们。这是我想要的输出

head (desired_output)
    2   0   1   0   0   0   0   3       
    18  0   19  0   0   0   0   17      
    2   1   0   0   0   0   0   1           
    18  19  0   0   20  0   0   23

知道如何准确地做到这一点吗？

答案1

这个想法是将连续的行保存在两个数组中，然后通过相应的索引比较数组元素。

将其保存到文件中，例如“twenty.awk”

#/usr/bin/env awk

# ref https://www.gnu.org/software/gawk/manual/html_node/Join-Function.html
function join(array, start, end, sep,    result, i)
{
    if (sep == "")
        sep = " "
    else if (sep == SUBSEP) # magic value
        sep = ""
    result = array[start]
    for (i = start + 1; i <= end; i++)
        result = result sep array[i]
    return result
}

{
    split($0, a)
    getline

    for (i=1; i<=NF; i++)
        if (a[i] + $i < 20)
            a[i] = $i = 0

    print join(a, 1, NF)
    print
}

然后运行

awk -f twenty.awk data.file | column -t > data.file.twenty

相关内容