VCF统计文件与MAF匹配的REF/ALT字段适配代码求助

处理逻辑 #

• 首先逐字符比对原REF（第4列）和ALT（第5列）序列，找到最长公共前缀的结束位置
• 插入变异（ALT长度 > REF长度）：截掉ALT序列中与REF公共的前缀部分，REF2保留原REF值
• 缺失变异（REF长度 > ALT长度）：截掉REF序列中与ALT公共的前缀部分，ALT2保留原ALT值
• 单碱基变异等长场景直接保留原REF/ALT值即可

代码实现 #

版本1：纯Python逐行处理（适合大体积VCF文件） #

def get_common_prefix_len(s1, s2):
    min_len = min(len(s1), len(s2))
    for i in range(min_len):
        if s1[i] != s2[i]:
            return i
    return min_len

input_file = "你的输入文件路径.txt"
output_file = "你的输出文件路径.txt"

with open(input_file, 'r') as f_in, open(output_file, 'w') as f_out:
    # 处理表头
    header = f_in.readline().strip().split('\t')
    header.extend(['REF2', 'ALT2'])
    f_out.write('\t'.join(header) + '\n')
    # 逐行处理数据
    for line in f_in:
        line = line.strip()
        if not line:
            continue
        parts = line.split('\t')
        ref = parts[3]
        alt = parts[4]
        prefix_len = get_common_prefix_len(ref, alt)
        if len(alt) > len(ref):
            # 插入变异场景
            ref2 = ref
            alt2 = alt[prefix_len:]
        elif len(ref) > len(alt):
            # 缺失变异场景
            ref2 = ref[prefix_len:]
            alt2 = alt
        else:
            # 等长变异（SNP、MNP）直接保留原值
            ref2 = ref
            alt2 = alt
        parts.extend([ref2, alt2])
        f_out.write('\t'.join(parts) + '\n')

版本2：Pandas处理（适合小文件，方便后续数据分析） #

import pandas as pd

def calc_ref2_alt2(row):
    ref = row['REF']
    alt = row['ALT']
    min_len = min(len(ref), len(alt))
    prefix_len = 0
    for i in range(min_len):
        if ref[i] != alt[i]:
            break
        prefix_len += 1
    if len(alt) > len(ref):
        return pd.Series([ref, alt[prefix_len:]])
    elif len(ref) > len(alt):
        return pd.Series([ref[prefix_len:], alt])
    else:
        return pd.Series([ref, alt])

# 读取文件，若为空格分隔则将sep参数改为'\s+'
df = pd.read_csv("你的输入文件路径.txt", sep='\t')
df[['REF2', 'ALT2']] = df.apply(calc_ref2_alt2, axis=1)
# 输出结果
df.to_csv("你的输出文件路径.txt", sep='\t', index=False)

验证结果 #

输入测试数据运行代码后，输出完全符合预期：
CHR POS ID REF ALT REF2 ALT2
chr11 71579744 rs71049992 A ACAGCAGCTGGACTGGGAGCAGCAGGACCTG A CAGCAGCTGGACTGGGAGCAGCAGGACCTG
chr11 124880551 rs71859853 CCGGAGT C CGGAGT C

注意事项 #

• 若输入文件为空格分隔，修改对应代码的分隔符参数即可
• 代码仅新增REF2、ALT2两列，原有杂合/纯合标记等其他字段全部保留不修改

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