如何用Python计算疟原虫线粒体VCF文件的等位基因频率？

一、INFO列中AF的计算逻辑 #

根据VCF文件头部定义，INFO字段的AF严格遵循公式 AF = AC / AN：

• AC：所有样本中该ALT等位基因的总计数（每个ALT对应一个AC值，与ALT列表顺序一致）
• AN：所有样本中已明确分型（called genotypes）的等位基因总数

对于疟原虫线粒体这类多拷贝基因组，这里的AC通常是所有样本中ALT等位基因的拷贝数/测序深度总和，AN则是所有样本中REF+ALT等位基因的拷贝数/测序深度总和，AF本质是ALT等位基因在群体中的占比。

二、无现成AC/AN时的计算方案 #

由于线粒体是多拷贝基因组，常规二倍体GT（0/0、0/1、1/1）的统计方式不适用，应基于FORMAT字段中的AD（等位基因深度）计算，这是最准确的方式：

• 总ALT计数：所有样本的AD字段中ALT深度的总和
• 总等位基因数：所有样本的AD字段中REF+ALT深度的总和
• 重新计算的AF = 总ALT深度 / 总REF+ALT深度

三、Python实现方案 #

方案1：基于已转换的Pandas DataFrame #

假设你的VCF已转为名为vcf_df的DataFrame，代码如下：

import pandas as pd

# 获取所有样本列（排除前9个固定列：CHROM到FORMAT）
sample_cols = vcf_df.columns[9:]

# 从样本字符串中提取REF和ALT深度
def parse_ad(sample_val):
    fields = sample_val.split(':')
    ref_d, alt_d = map(int, fields[1].split(','))
    return ref_d, alt_d

# 按变异位点计算重新计算的AF
def compute_site_af(row):
    total_ref = 0
    total_alt = 0
    for col in sample_cols:
        ref_d, alt_d = parse_ad(row[col])
        total_ref += ref_d
        total_alt += alt_d
    return total_alt / (total_ref + total_alt) if (total_ref + total_alt) > 0 else 0.0

# 添加新列存储重新计算的AF
vcf_df['recalculated_AF'] = vcf_df.apply(compute_site_af, axis=1)

方案2：基于原始VCF文件（使用pysam） #

直接读取VCF文件处理，无需提前转DataFrame：

import pysam

with pysam.VariantFile("sample.vcf") as vcf:
    for record in vcf:
        total_ref = 0
        total_alt = 0
        # 遍历所有样本的基因型数据
        for sample_data in record.samples.values():
            ad = sample_data.get('AD')
            if ad:
                total_ref += ad[0]
                total_alt += ad[1]
        # 计算AF
        af = total_alt / (total_ref + total_alt) if (total_ref + total_alt) > 0 else 0.0
        # 输出或处理结果
        print(f"位点 {record.chrom}:{record.pos} 重新计算的AF: {af:.6f}")

关键说明 #

疟原虫线粒体为多拷贝单亲遗传，GT字段的二倍体表示（如0/0）仅为caller默认格式，无法准确反映线粒体拷贝的等位基因组成，而AD字段直接提供了测序层面的等位基因深度，是计算群体AF的可靠依据。

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