记录一下从fasta提取ref碱基和比对提取alt碱基的过程

记录一下这个过程，后续可以复用。

需求

写个脚本，提取2个参考基因组 6a 和 7b 的 ref 和 alt 碱基

输入文件1格式：前两列的SNP名称均为 chrom_position ，分别是 6a 和 7b 的SNP (chrom_pos, 1_index)，隔，空字符串表示缺失。

SNP_6a	SNP_7b	geno_alleles_6a	geno_alleles_7b
10_11127142	10_10700428	T/C	T/C
10_13026587	10_12598321	A/G	A/G
10_12625800	10_12196013		A/G
10_12678847	10_12249356	A/G	A/G

输入文件2：6a参考基因组 fasta ，格式 fa.gz

输入文件3：7b参考基因组 fasta ，格式 fa.gz

输出文件： 输入文件1后面增加4列 ref_6a， alt_6a 和 ref_7b， alt_7b

脚本逻辑，以 6a 为例，7b同理

1. 如果 SNP_6a 为空，则 ref_6a 和 val_6a 均为空，continue

2. 如果 SNP_6a 不为空，提取6a参考基因组相应位点(SNP_6a)的碱基，即 ref_6a，下面开始提取 val_6a

   2.1 优先比较相同参考基因组的2个碱基，如果 geno_alleles_6a 不为空字符串，提取geno_alleles_6a的两个碱基，与 ref_6a 比较。如果 ref_6a 在这两个碱基中，则另一个为 val_6a ；否则 val_6a 为空字符串

   2.2 如果 geno_alleles_6a 为空但 geno_alleles_7b 不为空：

   ​ 2.2.1 提取 geno_alleles_7b 的两个碱基，与 ref_6a 比较。如果 ref_6a 在这两个碱基中，则另一个为 val_6a

   ​ 2.2.2 不然对 geno_alleles_7b 的两个碱基进行翻转，再次比较，如果 ref_6a 在这两个碱基中，则另一个为 val_6a

   ​ 2.2.3 还不行则 val_6a 为空

注意：fasta 文件较大，请使用专业包来从fasta文件中提取ref碱基，提高效率！

AI答复

下面给你一个用 pysam 的 Python 脚本。pysam.FastaFile 底层是 htslib/samtools faidx，支持对大参考基因组做随机访问，比反复读整个 fasta 快很多。

---

一、坐标说明

snp_id = chrom_position，1-based。

pysam 里：

fa.fetch(chrom, pos-1, pos)   # 取 1-based pos 这一位

---

二、脚本：add_ref_alt.py

我略微调整了一点

#!/usr/bin/env python3
import sys
import pysam

def parse_alleles(s):
    """'T/C' -> ['T','C']；空/非法返回 []"""
    if s is None:
        return []
    s = s.strip()
    if s == "":
        return []
    parts = s.split("/")
    if len(parts) != 2:
        return []
    a, b = parts[0].strip(), parts[1].strip()
    if a == "" or b == "":
        return []
    return [a, b]


def pick_alt(ref, alleles):
    """
    给定 ref 和两个等位基因 alleles=[a,b]
    如果 ref 在 alleles 中，返回另一个；否则返回 ''
    """
    if ref in alleles:
        for a in alleles:
            if a != ref:
                return a
    return ""


def revcomp(s):
    tbl = str.maketrans("ACGTacgt", "TGCAtgca")
    return s.translate(tbl)[::-1]


def get_ref(fa, snp_id):
    """
    snp_id = chrom_pos (1-based)
    返回大写 ref 碱基；解析失败返回 ''
    """
    if snp_id is None:
        return ""
    snp_id = snp_id.strip()
    if snp_id == "":
        return ""

    chrom, _, pos = snp_id.partition("_")
    if pos == "":
        return ""
    try:
        pos = int(pos)
    except ValueError:
        return ""

    try:
        base = fa.fetch(chrom, pos - 1, pos) # 0-based 半开区间 [start, end)
    except Exception:
        return ""
    return base.upper()


def resolve_alt(fa, snp_id, geno_same, geno_other):
    """
    以 6a 为例：
    - ref 来自 6a 参考基因组
    - 优先用 geno_alleles_6a
    - 没有则用 geno_alleles_7b，原方向不行就翻转后再比
    """
    ref = get_ref(fa, snp_id)
    if ref == "":
        return "", ""

    alleles_same = parse_alleles(geno_same)
    alleles_other = parse_alleles(geno_other)

    # 2.1 同参考基因组有 alleles
    if alleles_same:
        alt = pick_alt(ref, alleles_same)
        return ref, alt

    # 2.2 同参考基因组为空，用另一个基因组的 alleles
    if alleles_other:
        alt = pick_alt(ref, alleles_other)
        if alt != "":
            return ref, alt
        # 2.2.2 翻转后再比
        alt = pick_alt(ref, [revcomp(a) for a in alleles_other])
        return ref, alt

    # 都没有
    return ref, ""


def main(in_file, fa6a_path, fa7b_path, out_file):
    fa6a = pysam.FastaFile(fa6a_path)
    fa7b = pysam.FastaFile(fa7b_path)

    with open(in_file) as fin, open(out_file, "w") as fout:
        header = fin.readline().rstrip("\n").split("\t")
        # 原表头 + 4 列
        fout.write("\t".join(header + ["ref_6a", "alt_6a", "ref_7b", "alt_7b"]) + "\n")
        
        j = 1
        for line in fin:
            j += 1
            line = line.rstrip("\n")
            if line == "":
                continue
            cols = line.split("\t")

            # 防止列数不足
            if len(cols) < 4:
                print(f"Error: {j} row in {in_file} has not 4 columns", file=sys.stderr)
                sys.exit(1)

            snp_6a, snp_7b, geno_6a, geno_7b = cols[0], cols[1], cols[2], cols[3]

            # ---- 6a ----
            if snp_6a.strip() == "":
                ref_6a, alt_6a = "", ""
            else:
                ref_6a, alt_6a = resolve_alt(fa6a, snp_6a, geno_6a, geno_7b)

            # ---- 7b ----
            if snp_7b.strip() == "":
                ref_7b, alt_7b = "", ""
            else:
                ref_7b, alt_7b = resolve_alt(fa7b, snp_7b, geno_7b, geno_6a)

            fout.write(
                "\t".join(cols)
                + "\t" + ref_6a
                + "\t" + alt_6a
                + "\t" + ref_7b
                + "\t" + alt_7b
                + "\n"
            )

    fa6a.close()
    fa7b.close()


if __name__ == "__main__":
    if len(sys.argv) != 5:
        print("Usage: python add_ref_alt.py input.tsv ref6a.fa.gz ref7b.fa.gz out.tsv")
        sys.exit(1)

    main(sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4])
    

---

三、运行前准备 - fasta 处理

fasta 如果是压缩格式，必须是 BGZF 压缩 ，不能是普通 gzip压缩格式。

判断是不是 BGZF 。下面代码是 htslib 自带判断，最靠谱。

htsfile 6a.fa.gz

输出示例：

• BGZF：

  6a.fa.gz: BGZF compressed data

• 普通 gzip：

  6a.fa.gz: gzip compressed data

转换代码如下

gunzip ref6a.fa.gz          # 如果是普通 gzip
bgzip ref6a.fa
samtools faidx ref6a.fa.gz

pysam 需要 fasta 有 .fai 索引；第一次运行会自动建，也可以手动建：

samtools faidx ref6a.fa.gz
samtools faidx ref7b.fa.gz

安装一下 pysam 包

mamba install pysam

运行：

python add_ref_alt.py snps.tsv ref6a.fa.gz ref7b.fa.gz snps_with_refalt.tsv

---

四、你给的例子逻辑验证

输入：

10_11127142	10_10700428	T/C	T/C
10_13026587	10_12598321	A/G	A/G
10_12625800	10_12196013		A/G
10_12678847	10_12249356	A/G	A/G

处理规则完全按你写的：

1. SNP_6a 空 → ref_6a/alt_6a 空
2. SNP_6a 非空 → 取参考碱基
3. geno_alleles_6a 非空 → ref 在 {a,b} 中，另一个就是 alt
4. geno_alleles_6a 空但 geno_alleles_7b 非空 → 先用 7b alleles 比；不行就 reverse-complement 再比；还不行 alt 空

---

五、可选增强（如果你后面要扩展）

• 支持多等位 A/C/G
• ref 和 alleles 大小写统一 + N 处理
• 用 pyfaidx 代替 pysam（接口更像 dict）
• 输出 alt 为空时打 reason：no_alleles / ref_not_in_alleles / rc_fail

如果你愿意，我可以再帮你写一个 bedtools getfasta + awk 版本，或者把 alt 为空的原因也写进输出里。