从GVF文件提取引物区域SNP数据异常问题求助

问题背景 #

我需要从1000GENOMES-phase_3.gvf.gz文件中提取引物对应区域的SNP数据（MAF、rs编号、参考/替代等位基因），输入的引物信息为BLAST格式：

Primer_0|CYP2C19    NC_000010.11    100.000 23  0   0   1   23  94781749    94781771    2.65e-04    46.1    23 CCAGAGCTTGGCATATTATCT Within specified genomic region

使用以下Python脚本提取时，29条引物仅能找到20条对应的SNP数据，但通过命令行查询无结果的引物却能定位到数据：

import gzip

def fetch_snps_for_primer(chrom, start, end, sequence, gvf_file):
    """Fetch SNPs from the compressed GVF file by processing it line by line in Python."""
    results = []
    
    # Open the compressed GVF file using gzip
    with gzip.open(gvf_file, 'rt') as gvf:  # 'rt' means read text mode
        for line in gvf:
            if line.startswith("#"):  # Skip header lines
                continue
            
            parts = line.strip().split("\t")
            if len(parts) < 9:
                continue  # Skip malformed lines
            
            chrom_gvf = parts[0]
            pos_start = int(parts[3])  # SNP start position
            info_field = parts[8]

            # Check if this SNP is in the current primer region
            if chrom_gvf == chrom and start <= pos_start <= end:
                snp_id = None
                ref = None
                alts = []

                # Extract SNP ID, reference, and variant sequences
                for field in info_field.split(";"):
                    if field.startswith("Dbxref=dbSNP_"):
                        snp_id = field.split(":")[1]  # Extract rsID from Dbxref
                    elif field.startswith("Reference_seq="):
                        ref = field.split("=")[1]
                    elif field.startswith("Variant_seq="):
                        alts = field.split("=")[1].split(",")

                if snp_id:
                    results.append({
                        'Chromosome': chrom_gvf,
                        'Start': pos_start,
                        'End': pos_start,  # SNPs typically are a single base
                        'Primer Sequence': sequence,
                        'SNP ID': snp_id,
                        'Reference Allele': ref,
                        'Alternate Alleles': ','.join(alts)
                    })

    return results

排查方向及修复方案 #

1. 染色体命名不匹配 #

GVF文件中的染色体名称格式可能与引物数据不一致（比如引物用NC_000010.11，GVF用chr10或10），脚本中严格的字符串匹配会导致遗漏。

修复：添加染色体名称归一化函数，统一两者的命名格式：

def normalize_chromosome(chrom):
    # 转换NC_*格式为chrXX格式
    if chrom.startswith("NC_"):
        chrom_num = chrom.split("_")[2].split(".")[0]
        return f"chr{chrom_num}"
    # 处理其他可能的格式（如纯数字转chr前缀）
    elif chrom.isdigit():
        return f"chr{chrom}"
    return chrom

然后修改脚本中的匹配逻辑：

chrom_gvf_normalized = normalize_chromosome(chrom_gvf)
chrom_normalized = normalize_chromosome(chrom)
# 同时处理引物位置可能的反向情况（负链BLAST结果）
primer_start = min(start, end)
primer_end = max(start, end)

if chrom_gvf_normalized == chrom_normalized and primer_start <= pos_start <= primer_end:

2. SNP ID提取逻辑不全 #

脚本仅提取Dbxref=dbSNP_开头的rs编号，但GVF中SNP ID可能存在其他格式（比如Dbxref=rsXXXXXX、ID=rsXXXXXX），导致无匹配的SNP被过滤。

修复：扩展SNP ID提取逻辑，覆盖更多格式：

snp_id = None
ref = None
alts = []
maf = None

# 先从第2列（ID字段）尝试提取
if parts[1].startswith("rs"):
    snp_id = parts[1]

# 再从info字段提取
for field in info_field.split(";"):
    if not snp_id:
        if field.startswith("Dbxref="):
            # 处理Dbxref中的多条目（如Dbxref=dbSNP:rs1234,Ensembl:ENST...）
            for entry in field.split("=")[1].split(","):
                if "rs" in entry:
                    snp_id = entry.split(":")[-1]
                    break
        elif field.startswith("ID="):
            if field.split("=")[1].startswith("rs"):
                snp_id = field.split("=")[1]
    if field.startswith("Reference_seq="):
        ref = field.split("=")[1]
    elif field.startswith("Variant_seq="):
        alts = field.split("=")[1].split(",")
    elif field.startswith("MAF="):
        maf = field.split("=")[1]  # 补充提取用户需要的MAF数据

# 只要有rsID就保留结果（即使ref/alts缺失，可后续补全）
if snp_id:
    results.append({
        'Chromosome': chrom_gvf,
        'Start': pos_start,
        'End': int(parts[4]),  # 用GVF中的原始结束位置，兼容多碱基变异
        'Primer Sequence': sequence,
        'SNP ID': snp_id,
        'Reference Allele': ref or "N/A",
        'Alternate Alleles': ','.join(alts) or "N/A",
        'MAF': maf or "N/A"
    })

3. 位置范围判断忽略负链情况 #

BLAST结果中如果引物匹配的是负链，基因组的起始位置会大于结束位置，脚本中start <= pos_start <= end的判断会失效。

修复：使用min()和max()统一引物的位置范围，如上述染色体归一化部分的代码所示。

4. 日志排查异常行 #

添加日志打印，确认是否有合法行被误判为格式错误：

if len(parts) < 9:
    print(f"[WARN] Skipping malformed line: {line.strip()[:100]}...")
    continue

修改后的完整脚本 #

import gzip

def normalize_chromosome(chrom):
    """统一染色体命名格式"""
    if chrom.startswith("NC_"):
        chrom_num = chrom.split("_")[2].split(".")[0]
        return f"chr{chrom_num}"
    elif chrom.isdigit():
        return f"chr{chrom}"
    return chrom

def fetch_snps_for_primer(chrom, start, end, sequence, gvf_file):
    """Fetch SNPs from the compressed GVF file by processing it line by line in Python."""
    results = []
    
    with gzip.open(gvf_file, 'rt') as gvf:
        for line in gvf:
            if line.startswith("#"):
                continue
            
            parts = line.strip().split("\t")
            if len(parts) < 9:
                print(f"[WARN] Skipping malformed line: {line.strip()[:100]}...")
                continue
            
            chrom_gvf = parts[0]
            pos_start = int(parts[3])
            pos_end = int(parts[4])
            info_field = parts[8]

            # 归一化染色体并处理引物位置范围
            chrom_gvf_normalized = normalize_chromosome(chrom_gvf)
            chrom_normalized = normalize_chromosome(chrom)
            primer_start = min(start, end)
            primer_end = max(start, end)

            # 检查SNP是否在引物区域内（覆盖SNP跨位置的情况）
            if chrom_gvf_normalized == chrom_normalized and not (pos_end < primer_start or pos_start > primer_end):
                snp_id = None
                ref = None
                alts = []
                maf = None

                # 优先从ID字段提取rsID
                if parts[1].startswith("rs"):
                    snp_id = parts[1]

                # 从info字段提取详细信息
                for field in info_field.split(";"):
                    if not snp_id:
                        if field.startswith("Dbxref="):
                            for entry in field.split("=")[1].split(","):
                                if "rs" in entry:
                                    snp_id = entry.split(":")[-1]
                                    break
                        elif field.startswith("ID="):
                            entry_val = field.split("=")[1]
                            if entry_val.startswith("rs"):
                                snp_id = entry_val
                    if field.startswith("Reference_seq="):
                        ref = field.split("=")[1]
                    elif field.startswith("Variant_seq="):
                        alts = field.split("=")[1].split(",")
                    elif field.startswith("MAF="):
                        maf = field.split("=")[1]

                if snp_id:
                    results.append({
                        'Chromosome': chrom_gvf,
                        'Start': pos_start,
                        'End': pos_end,
                        'Primer Sequence': sequence,
                        'SNP ID': snp_id,
                        'Reference Allele': ref or "N/A",
                        'Alternate Alleles': ','.join(alts) or "N/A",
                        'MAF': maf or "N/A"
                    })

    return results

内容的提问来源于stack exchange，提问作者Elisa M.M.