序列比对算法原理详解¶

# === Needleman-Wunsch 全局比对算法实现 ===
import numpy as np  # 导入numpy用于矩阵操作

def needleman_wunsch(seq1, seq2, match=1, mismatch=-1, gap=-2):
    """
    全局比对算法
    seq1, seq2: 输入序列
    match: 匹配得分（默认+1）
    mismatch: 错配罚分（默认-1）
    gap: 空位罚分（默认-2）
    """
    m, n = len(seq1), len(seq2)  # 获取两条序列长度

    # 初始化得分矩阵，大小为(m+1) x (n+1)
    score = np.zeros((m + 1, n + 1), dtype=int)

    # 初始化第一行和第一列（累积gap罚分）
    for i in range(m + 1):
        score[i][0] = i * gap  # 第一列：seq1前i个字符vs空序列
    for j in range(n + 1):
        score[0][j] = j * gap  # 第一行：空序列vs seq2前j个字符

    # 填充得分矩阵（核心步骤）
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            # 判断当前位置是匹配还是错配
            if seq1[i-1] == seq2[j-1]:
                diag = score[i-1][j-1] + match      # 对角线：匹配
            else:
                diag = score[i-1][j-1] + mismatch    # 对角线：错配

            up = score[i-1][j] + gap    # 上方：seq1插入gap
            left = score[i][j-1] + gap  # 左方：seq2插入gap

            # 取三个方向的最大值
            score[i][j] = max(diag, up, left)

    # 回溯找最优比对路径
    align1, align2 = "", ""  # 存储比对结果
    i, j = m, n  # 从右下角开始回溯

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

        if i > 0 and j > 0 and score[i][j] == score[i-1][j-1] + current_score:
            align1 = seq1[i-1] + align1  # 对角线回溯：两条序列都取字符
            align2 = seq2[j-1] + align2
            i -= 1
            j -= 1
        elif i > 0 and score[i][j] == score[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 align1, align2, score[m][n]  # 返回比对结果和最终得分

# 测试
seq1 = "ACGTACGT"  # 序列1
seq2 = "ACGACG"    # 序列2
a1, a2, s = needleman_wunsch(seq1, seq2)
print(f"序列1: {a1}")   # 打印比对后的序列1
print(f"序列2: {a2}")   # 打印比对后的序列2
print(f"得分:  {s}")     # 打印最终得分

# === Smith-Waterman 局部比对核心差异 ===
def smith_waterman(seq1, seq2, match=2, mismatch=-1, gap=-1):
    """局部比对：只找最相似的片段"""
    m, n = len(seq1), len(seq2)
    score = np.zeros((m + 1, n + 1), dtype=int)  # 初始化全0（关键区别1）

    max_score = 0    # 记录最大得分
    max_pos = (0, 0) # 记录最大得分位置

    for i in range(1, m + 1):
        for j in range(1, n + 1):
            if seq1[i-1] == seq2[j-1]:
                diag = score[i-1][j-1] + match
            else:
                diag = score[i-1][j-1] + mismatch

            up = score[i-1][j] + gap
            left = score[i][j-1] + gap

            # 关键区别2：加入0选项，负分变0
            score[i][j] = max(0, diag, up, left)

            # 记录全局最大值位置（回溯起点）
            if score[i][j] > max_score:
                max_score = score[i][j]
                max_pos = (i, j)

    return max_score, max_pos  # 返回最高分和位置

步骤1: 建立查询序列的k-mer索引（种子表）
        ACGTACGT → ACG, CGT, GTA, TAC, ACG, CGT

步骤2: 扫描数据库，找精确匹配的种子（命中）
        数据库序列中找到 ACG 出现在位置 15, 203, 1567...

步骤3: 种子延伸（Smith-Waterman）
        从种子位置向两端延伸，计算比对得分
        得分低于阈值时停止延伸

步骤4: 统计显著性评估（E值计算）
        E = K * m * n * e^(-λS)
        K, λ: 统计参数  m,n: 序列长度  S: 比对得分

# === BLAST 基本使用 ===

# 1. 建库（把参考序列变成可搜索的数据库）
makeblastdb \
    -in reference.fasta \       # 输入：参考序列文件
    -dbtype nucl \              # 数据库类型：nucl=核酸, prot=蛋白质
    -out my_db \                # 输出：数据库名称前缀
    -parse_seqids               # 解析序列ID，方便后续提取

# 2. 核酸vs核酸比对（blastn）
blastn \
    -query query.fasta \        # 查询序列
    -db my_db \                 # 目标数据库
    -out results.txt \          # 输出文件
    -outfmt 6 \                 # 输出格式6=表格格式（最常用）
    -evalue 1e-5 \              # E值阈值（越小越严格）
    -num_threads 8 \            # 使用8个线程
    -max_target_seqs 5          # 每条query最多报5个hit

# 3. 蛋白vs蛋白比对（blastp）
blastp \
    -query protein.fasta \      # 蛋白查询序列
    -db nr \                    # NR蛋白数据库
    -out blastp_results.txt \   # 输出文件
    -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" \
    -evalue 1e-10               # 更严格的阈值

# 4. 核酸翻译后vs蛋白库（blastx）—— 用于功能注释
blastx \
    -query contigs.fasta \      # 核酸序列（如组装的contigs）
    -db uniprot_sprot \         # UniProt/SwissProt蛋白库
    -out blastx_results.txt \   # 输出
    -outfmt 6 \                 # 表格格式
    -evalue 1e-5                # E值阈值

# 5. 解读outfmt 6的12列含义
# qseqid   查询序列ID
# sseqid   目标序列ID  
# pident   一致性百分比
# length   比对长度
# mismatch 错配数
# gapopen  gap打开数
# qstart   查询起始位置
# qend     查询结束位置
# sstart   目标起始位置
# send     目标结束位置
# evalue   E值
# bitscore Bit得分

# === MAFFT 多序列比对 ===
# MAFFT是目前最推荐的MSA工具，速度快精度高

# 自动选择最优策略（推荐）
mafft --auto input.fasta > aligned.fasta  # --auto让MAFFT根据序列数量自动选策略

# 高精度模式（序列<200条时推荐）
mafft --localpair --maxiterate 1000 input.fasta > aligned.fasta
# --localpair: 使用Smith-Waterman做配对比对
# --maxiterate 1000: 迭代优化1000次

# 快速模式（序列>10000条时）
mafft --retree 1 input.fasta > aligned.fasta  # 只建一次引导树，速度最快

# === MUSCLE5 多序列比对（2022新版） ===
muscle -align input.fasta -output aligned.fasta  # 新版命令更简洁

# === 查看比对质量 ===
# 用TrimAl修剪低质量区域
trimal \
    -in aligned.fasta \         # 输入：比对后的序列
    -out trimmed.fasta \        # 输出：修剪后的序列
    -automated1                 # 自动选择最佳修剪参数

选择决策树：

比对的是什么？
├── 核酸序列
│   └── 使用默认 NUC44 矩阵（+1/-1 或 +2/-3）
│
└── 蛋白质序列
    ├── 序列相似度 > 80% → BLOSUM80
    ├── 序列相似度 60-80% → BLOSUM62（默认，最常用）
    ├── 序列相似度 35-60% → BLOSUM45
    └── 序列相似度 < 35% → PAM250（远缘同源）

白话解释：
- BLOSUM数字越大 → 对高相似序列越敏感（近亲检测）
- BLOSUM数字越小 → 对低相似序列越敏感（远亲检测）
- PAM系列：基于进化模型，PAM1=1%突变率
- BLOSUM系列：基于保守区块，更适合数据库搜索

# BLAST五大程序
blastn   核酸 vs 核酸库     最快，敏感度一般
blastp   蛋白 vs 蛋白库     蛋白同源搜索标准
blastx   核酸→蛋白 vs 蛋白库  功能注释推荐
tblastn  蛋白 vs 核酸→蛋白库  在基因组中找蛋白
tblastx  核酸→蛋白 vs 核酸→蛋白库  最慢，最敏感

# E值参考
E < 1e-50   极高置信度同源
E < 1e-10   高置信度同源
E < 1e-5    中等置信度
E < 0.01    弱同源，需谨慎
E > 1       可能是随机匹配

# 常用参数
-evalue     E值阈值（默认10，建议设1e-5）
-outfmt 6   表格输出（最方便解析）
-max_target_seqs  每条query最多报几个hit
-num_threads  线程数
-word_size   种子长度（blastn默认11，越小越敏感越慢）