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 查表 + 等位基因匹配/互补匹配 + 二义性判断”的过程。

---

输出文件

默认输出目录类似 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")