25. 序列比对原理与算法¶

句子A：我今天去超市买了苹果和香蕉
句子B：我昨天去商店买了苹果和橘子

序列1: A G C T
序列2: A C T

第一步：画一个 (m+1) × (n+1) 的表格（m和n是两条序列的长度）

        -    A    G    C    T
   -  [ 0] [-1] [-2] [-3] [-4]    ← 第一行/列用gap罚分填满
   A  [-1] [ 1] [ 0] [-1] [-2]
   C  [-2] [ 0] [ 0] [ 1] [ 0]
   T  [-3] [-1] [-1] [ 0] [ 2]  ← 右下角就是最终得分

第二步：每个格子从三个方向选最大值：
   ↖ 对角线（匹配/错配）：上一格分数 + 匹配分or错配罚分
   ← 左边（序列2插入gap）：左格分数 + gap罚分
   ↑ 上边（序列1插入gap）：上格分数 + gap罚分

第三步：从右下角回溯到左上角，走过的路径就是最优比对结果

匹配（match）：+1 或 +5
错配（mismatch）：-1 或 -4

罚分 = d × L
d = 每个gap的罚分（如 -2）
L = gap长度

例：gap长度5 → 罚分 = -2 × 5 = -10

罚分 = 开gap罚分(Gap Open) + 延伸罚分(Gap Extension) × (L-1)

例：Gap Open = -10, Gap Extension = -1, gap长度5
罚分 = -10 + (-1) × 4 = -14

而线性模型：-2 × 5 = -10（惩罚不够）

步骤：
1. 切词（Word/Seed）：把查询序列切成短片段
   - 蛋白质：默认3个氨基酸一个词（word size=3）
   - DNA：默认11个碱基一个词（word size=11）

2. 查表：把这些短词和数据库中所有序列快速比对
   - 只保留得分超过阈值 T 的匹配（"命中"/hit）

3. 延伸（Extension）：从每个命中位置向两边延伸
   - 延伸时计分，分数下降超过阈值就停止
   - 得到 HSP（High-scoring Segment Pair，高分片段对）

4. 评估：用统计学方法评估 HSP 的显著性（E-value）

原理演示（以字符串 "banana$" 为例）：

第一步：生成所有旋转（cyclic rotations）
banana$    ← 原始
anana$b    ← 向左旋转1位
nana$ba    ← 向左旋转2位
ana$ban    ← 向左旋转3位
na$bana    ← 向左旋转4位
a$banan    ← 向左旋转5位
$banana    ← 向左旋转6位

第二步：按字典序排序
$banana
a$banan
ana$ban
anana$b
banana$
na$bana
nana$ba

第三步：取最后一列 → BWT结果
annb$aa    ← 这就是BWT变换后的字符串

特性：相同字符聚集在一起，极易压缩！

搜索模式串的核心操作（Backward Search）：
- 从模式串的最后一个字符开始
- 每次用 LF-mapping（Last-to-First映射）缩小候选范围
- 最终得到模式串在原文中所有出现位置

时间复杂度：O(m)，m是模式串长度
与参考基因组大小无关！

概念：
- k-mer：长度为k的连续子串（如k=5：ATCGA）
- w-minimizer：在连续w个k-mer中，hash值最小的那个k-mer

示例（k=3, w=3）：
序列：A T C G A C T
k-mers: ATC, TCG, CGA, GAC, ACT
窗口1 (ATC, TCG, CGA) → 选hash最小的
窗口2 (TCG, CGA, GAC) → 选hash最小的
窗口3 (CGA, GAC, ACT) → 选hash最小的

效果：只需存储和比较 ~2/(w+1) 的k-mer
      大大减少了需要处理的种子数量

1. 收集种子（Seeding）
   - 对参考和query都提取minimizer
   - 通过hash表找到共同的minimizer（锚点/anchor）

2. 链式连接（Chaining）
   - 把距离合理的锚点串成链（chain）
   - 用动态规划找最优链（考虑锚点间距和方向一致性）
   - 时间复杂度：O(n log n)，n是锚点数

3. 局部比对（Alignment）
   - 只对链覆盖的区域做精确比对（带gap的SW）
   - 大大减少了需要精确比对的区域

假设有4条序列A,B,C,D
1. 两两比对得到距离矩阵
2. 根据距离建引导树：((A,B),(C,D))
3. 先比对A和B → 得到 AB
4. 再比对C和D → 得到 CD  
5. 最后比对 AB 和 CD → 得到 ABCD
问题：步骤3如果错了，后面全跟着错（error propagation）

# ============ Python 序列比对实操 ============
# 环境：conda activate bioinfo
# 安装：pip install biopython

# --- 示例1：手动实现简化版 Needleman-Wunsch ---
import numpy as np  # 数值计算库

def needleman_wunsch(seq1, seq2, match=1, mismatch=-1, gap=-2):
    """
    简化版全局比对（Needleman-Wunsch算法）
    seq1, seq2: 输入的两条序列（字符串）
    match: 匹配得分
    mismatch: 错配罚分
    gap: 空位罚分
    返回：最优得分和比对结果
    """
    m, n = len(seq1), len(seq2)  # 获取两条序列的长度

    # 第一步：初始化得分矩阵（m+1行 × n+1列，多一行一列给gap）
    score_matrix = np.zeros((m + 1, n + 1), dtype=int)

    # 第一行和第一列用gap罚分递增填充
    for i in range(1, m + 1):
        score_matrix[i][0] = score_matrix[i-1][0] + gap  # 第一列
    for j in range(1, n + 1):
        score_matrix[0][j] = score_matrix[0][j-1] + gap  # 第一行

    # 第二步：填表（动态规划核心）
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            # 判断当前位置是匹配还是错配
            if seq1[i-1] == seq2[j-1]:
                diag_score = score_matrix[i-1][j-1] + match    # 对角线：匹配
            else:
                diag_score = score_matrix[i-1][j-1] + mismatch # 对角线：错配

            up_score = score_matrix[i-1][j] + gap     # 从上方来：seq2加gap
            left_score = score_matrix[i][j-1] + gap   # 从左方来：seq1加gap

            # 取三个方向的最大值
            score_matrix[i][j] = max(diag_score, up_score, left_score)

    # 第三步：回溯（从右下角走到左上角）
    align1, align2 = "", ""  # 存储比对结果
    i, j = m, n  # 从右下角开始

    while i > 0 or j > 0:
        current = score_matrix[i][j]

        if i > 0 and j > 0:
            if seq1[i-1] == seq2[j-1]:
                s = match
            else:
                s = mismatch

        if i > 0 and j > 0 and current == score_matrix[i-1][j-1] + s:
            align1 = seq1[i-1] + align1  # 对角线走：两序列都出一个字符
            align2 = seq2[j-1] + align2
            i -= 1
            j -= 1
        elif i > 0 and current == score_matrix[i-1][j] + gap:
            align1 = seq1[i-1] + align1  # 向上走：seq2加gap
            align2 = "-" + align2
            i -= 1
        else:
            align1 = "-" + align1        # 向左走：seq1加gap
            align2 = seq2[j-1] + align2
            j -= 1

    return score_matrix[m][n], align1, align2

# 测试
seq1 = "AGCTGA"
seq2 = "ACGA"
score, a1, a2 = needleman_wunsch(seq1, seq2)
print(f"得分: {score}")    # 输出最优比对得分
print(f"序列1: {a1}")      # 输出比对后的序列1
print(f"序列2: {a2}")      # 输出比对后的序列2
# 预期输出：
# 得分: 0
# 序列1: AGCTGA
# 序列2: A-CG-A


# --- 示例2：用 Biopython 调用 pairwise2 模块做比对 ---
from Bio import pairwise2                        # Biopython的两序列比对模块
from Bio.pairwise2 import format_alignment       # 格式化输出比对结果

seq1 = "AGCTGATCG"
seq2 = "ACGATCG"

# globalxx：全局比对，只计匹配(+1)和错配(0)，不罚gap
alignments = pairwise2.align.globalms(
    seq1, seq2,
    2,    # match得分（匹配+2）
    -1,   # mismatch罚分（错配-1）
    -3,   # gap open罚分（开gap-3）
    -1    # gap extend罚分（延伸gap-1）
)

# 打印最优比对（可能有多个并列最优）
print("=== Biopython 全局比对结果 ===")
for align in alignments[:3]:  # 只看前3个结果
    print(format_alignment(*align))  # 格式化打印


# --- 示例3：用 Biopython 做 BLAST 搜索（在线） ---
from Bio.Blast import NCBIWWW, NCBIXML  # NCBI BLAST接口

# 一条示例16S rRNA序列片段
query_seq = "AGAGTTTGATCCTGGCTCAGGATGAACGCTGGCGGCGTGCCTAATACATGCAAGTCGAGCGG"

# 提交到NCBI BLAST服务器（blastn搜索核酸库）
# 注意：这是在线搜索，需要网络，可能需要几分钟
print("正在提交BLAST搜索...")
result_handle = NCBIWWW.qblast(
    "blastn",       # 程序名：核酸搜核酸
    "nt",           # 数据库：nt（核酸总库）
    query_seq,      # 查询序列
    expect=0.01     # E-value阈值
)

# 解析结果
blast_record = NCBIXML.read(result_handle)  # 解析XML格式结果

print(f"\n找到 {len(blast_record.alignments)} 条比对结果")
# 打印前3条结果
for i, alignment in enumerate(blast_record.alignments[:3]):
    hsp = alignment.hsps[0]  # 取第一个HSP
    print(f"\n--- 结果 {i+1} ---")
    print(f"目标: {alignment.title[:80]}...")  # 目标序列名称
    print(f"得分: {hsp.score}")                # 比对得分
    print(f"E-value: {hsp.expect}")            # E值
    print(f"一致性: {hsp.identities}/{hsp.align_length} "
          f"({hsp.identities/hsp.align_length*100:.1f}%)")  # 一致性百分比

测序数据是什么类型？
├── Illumina 短读（100-300bp）
│   ├── 比对到参考基因组 → BWA-MEM2 / Bowtie2
│   ├── RNA-seq → HISAT2 / STAR（剪接比对）
│   └── 去宿主 → Bowtie2（--un-conc 参数方便）
├── PacBio / Nanopore 长读（>1kb）
│   └── minimap2（-ax map-pb / map-ont）
├── 蛋白质/基因搜索数据库
│   └── BLAST（blastp/blastn/blastx）
└── 多条序列对齐
    ├── <200条 → MUSCLE / MAFFT
    └── >200条 → MAFFT（速度优势）

全局NW头对尾，局部SW找片段
BLAST种子加延伸，E值越小越靠谱
BWT压缩建索引，短读比对秒级处
minimizer挑代表，长读比对用minimap
MSA多条一起排，MAFFT快准第一选
BLOSUM大数看近亲，PAM大数看远亲