基于后验平均基因型(PMG)计算等位基因频率的实现方法问询

核心代码 #

{
    # 初始化单行统计计数器
    num_denom = 0  # 对应公式分子项累计值
    valid_indiv = 0 # 有效参与统计的个体总数（排除NA）
    # 从第4列开始遍历所有个体的PMG值，跳过前3列非PMG信息
    for (i=4; i<=NF; i++) {
        pmg = $i
        # 跳过当前NA字段，继续处理后续列，不跳转行
        if (pmg == "NA") continue
        pmg_val = pmg + 0 # 转为数值避免字符串比较偏差
        # 严格按给定规则做PMG基因型分箱
        if (pmg_val >=0 && pmg_val <0.5) {
            # 纯合参考型
            num_denom += 1
            valid_indiv +=1
        } else if (pmg_val >=0.5 && pmg_val <1.5) {
            # 杂合型
            num_denom +=1
            valid_indiv +=1
        } else if (pmg_val >=1.5 && pmg_val <=2) {
            # 纯合变异型
            num_denom +=2
            valid_indiv +=1
        }
        # 超出PMG取值范围的异常值自动跳过，不纳入统计
    }
    total_allele = 2 * valid_indiv # 对应公式分母的总有效等位基因数
    # 处理无有效样本的除零异常
    if (total_allele == 0) {
        printf "\"%s\", \"NA\", NA\n", $1
        next
    }
    # 按要求格式输出结果：第一列值、有效等位数/总有效等位数、频率值
    printf "\"%s\", \"%d/%d\", %.6f\n", $1, num_denom, total_allele, num_denom/total_allele
}

功能适配说明 #

• 自动适配可变列数：循环遍历到行尾总字段数位置，后续增减个体列无需修改代码
• 自动识别并跳过所有值为NA的字段，不受NA数量影响
• 严格执行PMG区间分箱规则，替代之前直接使用原始PMG值统计的偏差逻辑
• 无任何中间结果输出，仅输出要求格式的最终统计值
• 内置除零异常处理：当某行无有效统计个体时自动输出NA，避免程序报错

使用方法 #

将上述代码保存为pmg_freq.awk，执行以下命令即可得到结果：

awk -f pmg_freq.awk input_file.txt > result.txt

注意：代码默认使用awk原生空白分隔符适配大部分文本格式，如果输入文件为逗号、分号等固定分隔格式，只需在代码最开头添加BEGIN {FS="对应分隔符"}即可快速适配。

内容的提问来源于stack exchange，提问作者Bot75