如何将.vep.vcf文件多等位基因行拆分为双等位基因行

整体处理流程 #

逐行遍历输入VCF，区分头行、常规双等位行、多等位行三类分别处理，所有双等位行直接原样输出，多等位行按ALT数量复制为对应行数后，逐字段完成拆分匹配，最后添加溯源标记写入输出。

分模块实现逻辑 #

• 头行处理
  1. 所有##开头的元信息行直接原样写入输出
  2. 在#CHROM列名行之前，新增溯源字段的元信息声明：##INFO=<ID=OLD_MULTIALLELIC,Number=1,Type=String,Description="Original multiallelic locus before split">
  3. #CHROM开头的列名行直接原样写入
• 常规双等位行处理
  ALT字段无逗号分隔的行不需要修改，直接写入输出
• 多等位行核心拆分逻辑（ALT含逗号分隔的多个等位基因）
  1. 行结构拆分：按制表符切分整行，得到固定8列（CHROM/POS/ID/REF/ALT/QUAL/FILTER/INFO）、FORMAT列、后续所有样本列；将ALT字段按逗号切分为单个ALT的列表，列表长度即为需要拆分出的行数
  2. 溯源标记添加：每个拆分后的行，在INFO字段首位添加OLD_MULTIALLELIC=原CHROM:原POS:原REF:原ALT串标记，方便后续回溯原多等位位点
  3. INFO字段拆分
     • 把原INFO按分号切分为独立的键值对/标志位
     • 无值的标志位（比如INDEL、变异过滤相关标记）所有拆分行直接复用
     • 单值字段（比如AN、DP、MQ这类全局统计值，值本身无逗号分隔，或逗号分隔长度和ALT数量不匹配）所有拆分行直接复用
     • 多值等位基因对应字段（比如AC、AF、VEP输出的CSQ注释字段，逗号分隔长度和ALT数量完全一致），按逗号切分后取当前ALT对应索引位置的值即可，VCF规范明确要求这类字段的顺序和ALT列顺序一一对应，VEP输出完全遵循该规则
  4. FORMAT与样本列拆分
     • 将FORMAT按冒号切分为字段列表，逐字段处理每个样本的对应值
     • GT字段：将原多等位编码转换为双等位编码，原GT中0（REF）保留，对应当前ALT的编号替换为1，对应其他ALT的编号标记为缺失，避免基因型编码错误
     • AD类字段（长度为REF+ALT总数的逗号分隔值，比如REF+2个ALT对应长度为3）：拆分后仅保留REF对应值、当前ALT对应值
     • 单值字段（比如DP、GQ）所有拆分行直接复用
     • PL等基因型似然值字段：常规场景下取对应基因型位置的值即可，高精度分析场景可以后续用标准化工具做二次校正

可直接复用的Python实现代码 #

import sys

def split_multiallelic_vcf(in_path, out_path):
    with open(in_path, 'r') as fin, open(out_path, 'w') as fout:
        for line in fin:
            line = line.rstrip('\n')
            # 处理元信息头行
            if line.startswith('##'):
                fout.write(line + '\n')
                continue
            # 处理列名行，先写入新增的溯源字段元信息
            if line.startswith('#CHROM'):
                fout.write('##INFO=<ID=OLD_MULTIALLELIC,Number=1,Type=String,Description="Original multiallelic locus before split">\n')
                fout.write(line + '\n')
                continue
            # 切分变异行
            cols = line.split('\t')
            chrom, pos, var_id, ref, alt, qual, filt, info, fmt = cols[:9]
            samples = cols[9:]
            alt_list = alt.split(',')
            # 单ALT直接输出
            if len(alt_list) == 1:
                fout.write(line + '\n')
                continue
            # 预处理原INFO字段
            info_items = []
            for item in info.split(';'):
                if '=' in item:
                    k, v = item.split('=', 1)
                    info_items.append((k, v))
                else:
                    info_items.append((item, None))
            fmt_list = fmt.split(':')
            old_tag = f"OLD_MULTIALLELIC={chrom}:{pos}:{ref}:{alt}"
            # 逐ALT生成新行
            for alt_idx, single_alt in enumerate(alt_list):
                # 构建新INFO
                new_info = [old_tag]
                for k, v in info_items:
                    if v is None:
                        new_info.append(k)
                        continue
                    v_split = v.split(',')
                    if len(v_split) != len(alt_list):
                        new_info.append(f"{k}={v}")
                    else:
                        new_info.append(f"{k}={v_split[alt_idx]}")
                new_info_str = ';'.join(new_info)
                # 构建新样本列
                new_samples = []
                for samp_val in samples:
                    samp_split = samp_val.split(':')
                    new_samp = []
                    for f_key, f_val in zip(fmt_list, samp_split):
                        val_parts = f_val.split(',')
                        # 特殊处理GT字段
                        if f_key == 'GT':
                            gt_alleles = f_val.replace('|', '/').split('/')
                            new_gt = []
                            for a in gt_alleles:
                                if a == '.':
                                    new_gt.append('.')
                                elif int(a) == 0:
                                    new_gt.append('0')
                                elif int(a) == alt_idx + 1:
                                    new_gt.append('1')
                                else:
                                    new_gt.append('.')
                            new_samp.append('/'.join(new_gt))
                        # 处理REF+ALT长度的字段（如AD）
                        elif len(val_parts) == len(alt_list) + 1:
                            new_samp.append(f"{val_parts[0]},{val_parts[alt_idx+1]}")
                        # 处理仅ALT长度的字段
                        elif len(val_parts) == len(alt_list):
                            new_samp.append(val_parts[alt_idx])
                        # 单值字段直接复用
                        else:
                            new_samp.append(f_val)
                    new_samples.append(':'.join(new_samp))
                # 拼接新行写入
                new_line = '\t'.join([
                    chrom, pos, var_id, ref, single_alt, qual, filt, new_info_str, fmt
                ] + new_samples)
                fout.write(new_line + '\n')

if __name__ == "__main__":
    split_multiallelic_vcf(sys.argv[1], sys.argv[2])

注意：代码运行方式为python split_vcf.py 输入文件路径 输出文件路径，如果对PL等复杂基因型字段的准确性要求极高，拆分完成后可以用常规VCF标准化工具做一次校验，避免自定义逻辑出现索引错位。

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