第2章 等位基因频率的基本计算公式

好，咱们进入正题。等位基因频率——说白了就是某个等位基因占该基因座上所有等位基因总数的比例。你手头有一堆样本的基因型数据，怎么算？公式本身不复杂，但细节容易翻车。

我刚开始做群体遗传分析时，遇到过一个坑。当时用PLINK算出的频率和我手动核对的数值差了零点几，排查了半天才发现是缺失数据计数方式搞混了。后来我养成一个习惯：凡是涉及频率计算，先拿一个小数据集手算一遍，再交给软件跑。

2.1 基础公式：从基因型计数到等位基因频率

假设一个二倍体群体，某基因座上有两个等位基因：A 和 a。你的数据是每个样本的基因型：AA、Aa、aa。记下它们的个体数：

• N_AA：纯合野生型个体数
• N_Aa：杂合个体数
• N_aa：纯合突变型个体数
• N_total = N_AA + N_Aa + N_aa：总样本数

那么等位基因 A 的频率 p 和等位基因 a 的频率 q 就是：

p = (2 × N_AA + 1 × N_Aa) / (2 × N_total)
q = (2 × N_aa + 1 × N_Aa) / (2 × N_total)
p + q = 1

为什么分母乘以2？因为每个二倍体个体贡献两个等位基因。这儿容易犯的一个错误是直接用个体数当分母。我见过有人写成 p = N_AA / N_total，那算出来的是基因型频率，不是等位基因频率。两码事。

2.2 实战：手工算一遍

咱们拿一个简单例子练手。假设你测了100个人，基因型分布如下：

套公式：

p = (2×36 + 48) / (2×100) = (72 + 48) / 200 = 120 / 200 = 0.6
q = (2×16 + 48) / (2×100) = (32 + 48) / 200 = 80 / 200 = 0.4

验证一下：0.6 + 0.4 = 1。没问题。

2.3 多等位基因情况

现实中很多基因座不止两个等位基因，比如微卫星位点或者HLA系统。公式需要扩展一下。

设某基因座有 k 个等位基因：A1, A2, ..., Ak。对于等位基因 Ai，它的频率 pi 等于：

pi = (纯合 AiAi 个体数×2 + 所有包含 Ai 的杂合个体数之和) / (总个体数 × 2)

写成计数形式：

令 n_i = 样本中出现的等位基因 Ai 的总个数
pi = n_i / (2 × N_total)

所有等位基因频率之和必须为1：∑ pi = 1。

我曾经处理过一个STR位点有12个等位基因的数据集。手动算每个等位基因的计数时，用Excel的COUNTIF函数加乘2，再除以总条数——但一定要注意杂合子里的两条等位基因要分别计入对应等位基因的计数里。这里我栽过一次：批量算的时候忘了给杂合子计数翻倍，结果频率和等于0.5而不是1，排查了好久才发现。

2.4 缺失数据怎么处理？

现实数据永远不会干干净净。一定会遇到部分样本某个基因型缺失（比如测序深度不够、分型失败）。这时候分母要不要变？

我个人推荐两种常见策略：

1. 完全剔除（Complete Deletion）：对于该基因座，只保留所有样本都有有效基因型的数据。简单粗暴，但如果缺失率很高会损失样本量。
2. 部分保留（Available Deletion）：分子只统计有效等位基因数，分母也只算有效样本贡献的等位基因总数。即去掉缺失样本再算频率。

假设你有100个样本，其中2个样本的基因型缺失（记为NA）。剩下98个样本的计数如下：

按部分保留策略：

有效总样本数 = 50 + 30 + 18 = 98
有效等位基因总数 = 98 × 2 = 196

p = (2×50 + 30) / 196 = (100 + 30) / 196 = 130 / 196 ≈ 0.6633
q = (2×18 + 30) / 196 = (36 + 30) / 196 = 66 / 196 ≈ 0.3367

2.5 用代码验证你的手工计算

我建议你在正式跑全基因组数据之前，先写个小脚本验证自己的逻辑。下面是一个简单的Python片段，可以自己改数据测试：

def allele_freq(genotype_counts, alleles=('A','a')):
    """
    genotype_counts: dict, 格式如 {'AA':36, 'Aa':48, 'aa':16}
    """
    total_individuals = sum(genotype_counts.values())
    total_alleles = total_individuals * 2

    freq = {}
    for allele in alleles:
        count = 0
        for gt, n in genotype_counts.items():
            # 计算该基因型中当前等位基因出现的次数
            allele_occurrences = gt.count(allele)
            count += allele_occurrences * n
        freq[allele] = count / total_alleles
    return freq

# 测试
data = {'AA':36, 'Aa':48, 'aa':16}
print(allele_freq(data))

输出应该得到 {'A': 0.6, 'a': 0.4}。如果你拿真实数据跑，发现结果奇怪，建议先检查缺失标记，再检查计数逻辑。

2.6 章节知识图谱：频率计算的核心流程

下面这张SVG图总结了整个计算逻辑，从原始基因型到频率输出，中间经过计数、比对、校验三个步骤。