python包-snpflip

snpflip 也是一个已经停止更新维护的包，用于比对bim文件和参考基因组，查找翻转的位点。

目的

snpflip 用 参考基因组 FASTA 把每个 SNP 位点的 REF 碱基查出来，然后判断：

这个 SNP 的 A1/A2 是在 forward 链上、reverse 链上、还是根本分不清。

输出给你两张关键名单：

• 哪些 SNP 要 –flip
• 哪些 SNP 是“ambiguous（分不清链）”应该考虑剔除

其实就是进行碱基的比对而已，自己也能写，逻辑就是用参考基因组REF 与 A1/A2 的关系，推断链状态。

原理

就是用参考基因组REF 与 A1/A2 的关系，推断链状态。

记：

• REF = FASTA 查出来的参考碱基（A/T/C/G）
• REF_rev = 互补(complement(REF))

snpflip 本质上在做 集合/匹配：

情况	结论（你看到的输出）	为什么
A1==REF 或 A2==REF	forward	至少一个 allele 匹配 forward-strand REF ⇒ SNP 已在参考链
A1==REF_rev 或 A2==REF_rev	reverse	只能在 reverse 链上才能“对齐上” ⇒ 需要 flip
两者都不成立	（通常会报 mismatch）	最可能：坐标/版本不一致，或 REF=N，或数据不是同一组装
A/T 或 C/G 且你只靠“等位基因相等”无法区分链	ambiguous	因为 complement(A)==T，正向反向“看起来都对”

更直白地说：它不是一个复杂模型，就是一个“REF 查表 + 等位基因匹配/互补匹配 + 二义性判断”的过程。

什么时候会出现 mismatch？

下面看一下 mismatch 这种异常可能有哪些情况出现

首先是当参考基因组 REF=N 这种未知碱基或不确定碱基时算 mismatch（如果参考基因组上没有这个物理位置，提取不到这个位置的 REF 碱基应该也算）。

当 REF 时正常的 ACGT 碱基时，什么情况会出现 mismatch 呢？

先复习一下基础知识。

碱基：腺嘌呤（A）、鸟嘌呤（G）、胞嘧啶（C）、胸腺嘧啶（T）（DNA）；尿嘧啶（U）（RNA）。

配对：A=T（两条链间两个氢键），G≡C（三个氢键）。在RNA中，A=U。

转换：嘌呤变为嘌呤（A/G）或嘧啶变为嘧啶（C/T）。即，环结构保持不变。因为化学结构接近，出现比例更高。

颠换：嘌呤变为嘧啶，或反之（A/C，A/T，G/C，G/T）。即，环结构发生改变。结构差异大，相对少见。

但是对于这里比对而言，应该分为下面2种情况：

1. bim文件中SNP变异 非 A/T 或 C/G ，例如 A/G ，此时无论参考基因组 REF 是哪一个碱基，必然与 A1 或 A2 中某一个碱基相同或互补。
2. bim文件中SNP变异 是 A/T 或 C/G ，这种情况两个碱基互补，举例 A/T ，此时正常来说，基因组只可能是 A 或 T，而且无法区分链，即无法判断是否翻转；如果参考基因组是 G 或 C ，就不符合常理了，说明基因型和参考基因组二者不匹配，属于异常情况（报 mismatch）。

输出文件

默认输出目录类似 snpflip_output/，里面常见：

• .reverse：写在 reverse strand 的 SNP 列表
  → 你接着用
  

  plink --bfile mydata --flip snpflip_output.reverse --make-bed --out flipped

• .ambiguous：链无法确定的 SNP 列表
  → 通常做法是检查后决定保留/剔除（很多管线会选择剔除 A/T、C/G 的 ambiguous SNP）
• .annotated_bim：把 forward/reverse/ambiguous 标回去的表，方便人工检查

问题-不匹配W/Z染色体

我使用 0.06 版本，其说明如下，这里已经说了，bim 文件和 fasta 中的染色体只能是数字和 X, Y, M （fasta 中的染色体可以有 chr 前缀，分析的时候会自动剔除 chr 前缀）。

$ snpflip --help
snpflip

Report reverse and ambiguous strand SNPs.
(Visit github.com/endrebak/snpflip for examples and help.)

Usage:
    snpflip --fasta-genome=FA --bim-file=BIM [--output-prefix=PRE]
    snpflip --help
    snpflip --version

Arguments:
    -f FA --fasta-genome=FA     fasta genome
    -b BIM --bim-file=BIM       plink bed file (extended variant information)

Options:
    -h --help                   show this message
    -v --version                show the version number
    -o PRE --output-prefix=PRE  the prefix of the output-files
                                (./snpflip_output/<bim_basename> by default)

Note:
    To enable snpflip to match the chromosomes in the `.bim` and `.fa` files, the
    chromosomes in the `.bim` file must be called `1, 2, ... X, Y, M`
    while the chromosomes in the fasta file must be called the same, or
    `chr1, chr2, ... chrX, chrY, chrM`

    If there is whitespace after the chromosome name in the fasta file,
    additional text is allowed, e.g. the below fasta id is accepted:
    >1 dna:chromosome chromosome:GRCh37:1:1:249250621:1

实际测试发现，鸡的W/Z染色体会导致它直接无法分析。

查看代码是 reference_genome.py 其中这两个函数的问题，源代码如下。

• convert_fa_chromosome_names(fasta_dict) 是这三个函数的主函数，调用了另外2个函数，其作用是修改 fasta 染色体名称，并提取需要的染色体。
  • 输入参数 fasta_dict 是一个字典，键为 fasta 文件中的染色体，值为序列。
    • 输出是提取修改后的字典，按照原本代码的意思是只保留提取数字和 X, Y, M 的染色体（如果有 前缀chr 则剔除chr ，有后缀 T 也剔除）
• _is_valid_chromosome(name) : 剔除 fasta 中的非法染色体（存在 _ 或 .的染色体）
• _convert_fa_name(name) : 对常规染色体进行转换，进行标准化。这里使用了一个正则表达式，允许的染色体格式是 chr 开头，数字或 X, Y, M 为中间字符，以 T 结尾的格式，然后提取其中数字或 X, Y, M 作为标准化后的染色体。这里问题就是，由于这里没有对染色体进行进一步的校验，如果染色体既不是上面的非法染色体（存在 _ 或 .的染色体），也不是符合正则表达式的格式，那么这里会直接报错退出。

def convert_fa_chromosome_names(fasta_dict):

    return {_convert_fa_name(k): fasta_dict[k] for k in fasta_dict.keys()
            if _is_valid_chromosome(k)}


def _is_valid_chromosome(name):

    for c in "_.":
        if c in name.split(None, 1)[0]:
            return False

    return True


def _convert_fa_name(name):

    bim_name = re.compile(r"(chr)?(?P<chr_symbol>\d+|X|Y|M)(?:T)?")
    result = bim_name.match(name).group("chr_symbol")

    return result

修改后的代码如下，这里修改了第1个和第3个代码，修改点如下

1. _convert_fa_name() 中的正则表达式加入 W 和 Z 和 MT，加入鸡的染色体。
2. 更稳健，如果出现不符合正则的染色体，程序也不会莫名报错退出，而只是剔除这个染色体继续分析。

def convert_fa_chromosome_names(fasta_dict):

    result = {}
    for k in fasta_dict.keys():
        if not _is_valid_chromosome(k):
            continue
        std_name = _convert_fa_name(k)
        if std_name is None:
            continue
        result[std_name] = fasta_dict[k]
    return result


def _is_valid_chromosome(name):

    for c in "_.":
        if c in name.split(None, 1)[0]:
            return False

    return True


def _convert_fa_name(name):

    # ✅ 关键修改：加入 W 和 Z 和 MT
    # ✅ 去掉历史遗留的 (?:T)? 后缀
    bim_name = re.compile(r"(chr)?(?P<chr_symbol>\d+|X|Y|M|MT|W|Z)")
    m = bim_name.match(name)
    if m is None:
        return None
    return m.group("chr_symbol")

想法：两个参考基因组的碱基互补是否表示翻转？

答案是并不一定，下面是 AI 的回答，意思就是有可能旧版这个位置测错了，又或者这个位置由于存在多态性，两次选择作为 REF 的碱基不一样。

在从旧版参考基因组（如 GRCh37）升级到新版（如 GRCh38）的过程中，科学家不仅会调整染色体结构，还会修复测序错误、补充新的测序数据。

• 情况 A（真翻转）：某个区域的 DNA 在两个版本间发生了倒位（Inversion）。此时，新版的序列是旧版的反向互补。
• 情况 B（假翻转/仅碱基更新）：某个位点没有发生倒位，物理方向没变。但是科学家发现旧版参考基因组在这个位置测错了（比如旧版是 A，新版修正为 T），或者这个位置存在多态性，新版选择了一个不同的碱基作为参考。

那应该怎么看？如果真的要这么看的话，就要看前后这一段序列是否都发生了翻转，才能确定是不是真的翻转。

我们不能只看“单点”的碱基是不是互补，而必须提取这个位点前后几十到上百个碱基（侧翼序列）进行整体比对。只有当一大段序列都变成了反向互补时，我们才能拍板确认：“这里确实发生了物理方向的翻转！”