生信面试算法与数据结构¶

# ========== 哈希表：生信中最重要的数据结构 ==========

# 场景1：基因ID映射（O(1)查找）
gene_map = {}                                   # 空字典
gene_map["ENSG00000141510"] = "TP53"           # 建立映射
gene_map["ENSG00000012048"] = "BRCA1"
print(gene_map.get("ENSG00000141510", "Unknown"))  # O(1)查找

# 场景2：K-mer计数
from collections import Counter

def count_kmers(seq, k=21):
    """
    统计序列中所有k-mer的频率
    时间复杂度: O(n), n=序列长度
    空间复杂度: O(4^k), k-mer种类数
    """
    kmers = [seq[i:i+k] for i in range(len(seq)-k+1)]  # 生成所有k-mer
    return Counter(kmers)                               # 哈希计数

# 应用：通过k-mer频率判断测序错误
# 正确的k-mer出现频率高，测序错误导致的k-mer出现频率低（1-2次）

# 场景3：两个基因列表求交集
set_a = set(["TP53", "BRCA1", "EGFR", "KRAS"])   # 集合A
set_b = set(["BRCA1", "EGFR", "MYC", "ALK"])     # 集合B
intersection = set_a & set_b                       # 交集：O(min(m,n))
union = set_a | set_b                              # 并集
diff = set_a - set_b                               # 差集（A有B没有）
print(f"交集: {intersection}")                     # {'BRCA1', 'EGFR'}

# 面试要点：set底层也是哈希表，in操作O(1)
# list的in操作是O(n)，查找大量元素时set比list快很多

# ========== 排序：基础但重要 ==========

# Python内置排序：Timsort，O(n log n)
genes = [("TP53", 0.001), ("BRCA1", 0.05), ("EGFR", 0.0001)]

# 按p值排序
sorted_genes = sorted(genes, key=lambda x: x[1])    # 按第2个元素升序
print(sorted_genes)                                   # EGFR, TP53, BRCA1

# 按p值降序
sorted_genes = sorted(genes, key=lambda x: x[1], reverse=True)

# 多级排序：先按染色体，再按位置
variants = [
    ("chr1", 1000, "A", "T"),
    ("chr2", 500, "G", "C"),
    ("chr1", 500, "C", "G"),
]
sorted_variants = sorted(variants, key=lambda x: (x[0], x[1]))  # 元组排序
# 结果：chr1:500, chr1:1000, chr2:500

# 面试题：找出表达量最高的Top 10基因
import heapq
expression = {"TP53": 1500, "BRCA1": 800, "EGFR": 3000, "MYC": 2500}
# 方法1：全排序 O(n log n)
top10 = sorted(expression.items(), key=lambda x: x[1], reverse=True)[:10]
# 方法2：堆 O(n log k)，k=10，大数据更快
top10 = heapq.nlargest(10, expression.items(), key=lambda x: x[1])

# 面试要点：
# 稳定排序(stable sort): 相等元素保持原始顺序 → Python的sorted()是稳定的
# 不稳定排序: 相等元素顺序不确定 → 快速排序

# ========== 二分查找：有序数据的快速查找 ==========
import bisect

def binary_search(arr, target):
    """
    在有序数组中查找target
    时间复杂度: O(log n)
    """
    left, right = 0, len(arr) - 1                # 左右边界
    while left <= right:
        mid = (left + right) // 2                 # 中间位置
        if arr[mid] == target:
            return mid                             # 找到了
        elif arr[mid] < target:
            left = mid + 1                         # 目标在右半边
        else:
            right = mid - 1                        # 目标在左半边
    return -1                                      # 没找到

# 生信应用：BED文件区间查询
# 问题：给定一个基因组位置，快速找到它在哪个基因/外显子区间
def find_region(positions, query):
    """
    在排序的位置列表中找到query应该插入的位置
    应用：判断变异是否在目标区间内
    """
    idx = bisect.bisect_left(positions, query)    # 二分查找插入位置
    return idx

# 示例：exon起始位置已排序
exon_starts = [100, 500, 1000, 2000, 5000, 8000]
query_pos = 750                                    # 变异位置
idx = bisect.bisect_right(exon_starts, query_pos)
print(f"位置{query_pos}在第{idx}个区间之后")       # 在第2个之后（500-1000之间）

# 面试要点：
# BAM文件的索引(.bai)使用类似二分查找来快速定位基因组区域
# tabix索引也是基于排序+二分查找

# ========== 动态规划：序列比对的数学基础 ==========

# 全局比对（Needleman-Wunsch）
def needleman_wunsch(seq1, seq2, match=1, mismatch=-1, gap=-2):
    """
    全局序列比对（Needleman-Wunsch算法）
    时间复杂度: O(m*n)
    空间复杂度: O(m*n)
    """
    m, n = len(seq1), len(seq2)

    # 初始化DP矩阵
    dp = [[0]*(n+1) for _ in range(m+1)]         # (m+1) x (n+1) 矩阵
    for i in range(m+1):
        dp[i][0] = i * gap                        # 第一列：gap累加
    for j in range(n+1):
        dp[0][j] = j * gap                        # 第一行：gap累加

    # 填充DP矩阵
    for i in range(1, m+1):
        for j in range(1, n+1):
            if seq1[i-1] == seq2[j-1]:            # 碱基匹配
                diag = dp[i-1][j-1] + match
            else:                                  # 碱基不匹配
                diag = dp[i-1][j-1] + mismatch
            up = dp[i-1][j] + gap                  # 从上方（seq1加gap）
            left = dp[i][j-1] + gap                # 从左方（seq2加gap）
            dp[i][j] = max(diag, up, left)         # 取最大值

    # 回溯
    align1, align2 = [], []
    i, j = m, n
    while i > 0 or j > 0:
        if i > 0 and j > 0:
            score = match if seq1[i-1] == seq2[j-1] else mismatch
            if dp[i][j] == dp[i-1][j-1] + score:  # 对角线
                align1.append(seq1[i-1])
                align2.append(seq2[j-1])
                i -= 1
                j -= 1
                continue
        if i > 0 and dp[i][j] == dp[i-1][j] + gap:  # 上方
            align1.append(seq1[i-1])
            align2.append('-')
            i -= 1
        else:                                          # 左方
            align1.append('-')
            align2.append(seq2[j-1])
            j -= 1

    return ''.join(reversed(align1)), ''.join(reversed(align2)), dp[m][n]

# 测试
a1, a2, score = needleman_wunsch("ATCGTAC", "ATGTAC")
print(f"Seq1: {a1}")                              # ATC-GTAC 或 ATCGTAC
print(f"Seq2: {a2}")                              # AT-GTAC  或 AT-GTAC
print(f"Score: {score}")

# 局部比对（Smith-Waterman）
# 与NW的区别：
# 1. DP矩阵最小值为0（不允许负分）
# 2. 从最高分处开始回溯（不是右下角）
# 3. 回溯到分数为0时停止

# 面试要点：
# BLAST 不是Smith-Waterman，而是启发式的种子-延伸算法（更快但不保证最优）
# BWA-MEM 使用 BWT+后缀数组做精确匹配种子，再用Smith-Waterman做延伸

# ========== BWT变换：BWA/Bowtie2比对的核心算法 ==========

def bwt_transform(text):
    """
    Burrows-Wheeler Transform
    将文本变换为更容易压缩和搜索的形式
    """
    text = text + '$'                              # 添加终止符
    n = len(text)
    rotations = [text[i:] + text[:i] for i in range(n)]  # 所有旋转
    rotations.sort()                               # 字典序排序
    bwt = ''.join(r[-1] for r in rotations)        # 取每个旋转的最后一个字符
    return bwt

# 示例
seq = "ATCGATCG"
bwt = bwt_transform(seq)
print(f"原始序列: {seq}")
print(f"BWT变换: {bwt}")

# 面试要点：
# BWT的神奇之处：
# 1. 变换是可逆的（可以从BWT恢复原始序列）
# 2. 相同字符聚在一起，利于压缩
# 3. 配合FM-index可以O(m)时间精确匹配长度为m的模式
# 4. 参考基因组3GB，BWT索引约3-4GB，比后缀树(~45GB)小很多
# BWA建索引(bwa index)就是在构建BWT+FM-index

# ========== de Bruijn图：基因组组装的核心 ==========

from collections import defaultdict

def build_debruijn(reads, k=3):
    """
    从reads构建de Bruijn图
    节点 = (k-1)-mer
    边 = k-mer
    """
    graph = defaultdict(list)                     # 邻接表
    for read in reads:
        for i in range(len(read) - k + 1):        # 遍历所有k-mer
            kmer = read[i:i+k]
            prefix = kmer[:-1]                     # 前(k-1)个字符
            suffix = kmer[1:]                      # 后(k-1)个字符
            graph[prefix].append(suffix)           # prefix → suffix
    return dict(graph)

# 示例
reads = ["ATCGA", "TCGAT", "CGATC"]
graph = build_debruijn(reads, k=3)
print("de Bruijn图：")
for node, edges in sorted(graph.items()):
    print(f"  {node} → {edges}")

# 面试要点：
# 理想情况：de Bruijn图有欧拉路径 → 组装出完整序列
# 实际问题：重复序列导致图有分支（bubble）、错误导致tip
# 组装工具（如SPAdes、MEGAHIT）就是在解决这些图结构问题
# k值选择：k太小→图太复杂（重复多），k太大→图不连通（覆盖度不够）

# ========== 图算法：生信网络分析基础 ==========

# BFS（广度优先搜索）— 找最短路径
from collections import deque

def bfs_shortest_path(graph, start, end):
    """
    BFS找两个基因在PPI网络中的最短路径
    应用：药物靶点预测、基因功能推断
    """
    visited = {start}                             # 已访问节点
    queue = deque([(start, [start])])             # 队列：(节点, 路径)

    while queue:
        node, path = queue.popleft()              # 取出队首
        if node == end:
            return path                            # 找到目标

        for neighbor in graph.get(node, []):      # 遍历邻居
            if neighbor not in visited:
                visited.add(neighbor)
                queue.append((neighbor, path + [neighbor]))

    return None                                    # 不可达

# PPI网络示例
ppi_network = {
    "TP53": ["MDM2", "BRCA1", "ATM"],
    "BRCA1": ["TP53", "RAD51", "PALB2"],
    "MDM2": ["TP53"],
    "ATM": ["TP53", "CHK2"],
    "CHK2": ["ATM", "BRCA1"],
    "RAD51": ["BRCA1"],
    "PALB2": ["BRCA1"],
}

path = bfs_shortest_path(ppi_network, "MDM2", "RAD51")
print(f"MDM2到RAD51最短路径: {' → '.join(path)}")
# MDM2 → TP53 → BRCA1 → RAD51

# 节点度（degree）— 连接数
degrees = {gene: len(partners) for gene, partners in ppi_network.items()}
hub_genes = sorted(degrees.items(), key=lambda x: -x[1])  # 按度排序
print(f"Hub基因: {hub_genes[0][0]}（度={hub_genes[0][1]}）")
# Hub基因通常是重要的调控因子

# ========== 复杂度分析速查 ==========

# 时间复杂度（从快到慢）：
# O(1)      — 常数时间（哈希查找）
# O(log n)  — 对数时间（二分查找、BAM索引查询）
# O(n)      — 线性时间（遍历一遍、grep搜索）
# O(n log n)— 排序时间（sort、samtools sort）
# O(n²)     — 平方时间（简单比对、嵌套循环）
# O(n*m)    — Smith-Waterman比对（n和m是两条序列长度）
# O(2^n)    — 指数时间（暴力搜索，不可接受）

# 常见生信操作的复杂度：
operations = {
    "FASTA读取": "O(n) 时间, O(n) 空间",
    "K-mer计数": "O(n) 时间, O(4^k) 空间",
    "BLAST搜索": "O(n*m) 最坏, 启发式约O(n)",
    "BWA比对": "O(m) 查找种子 + O(m*w) 延伸",
    "samtools sort": "O(n log n) 时间",
    "VCF过滤(grep)": "O(n) 时间",
    "BED区间查询": "O(log n) 用tabix索引",
    "NW全局比对": "O(n*m) 时间, O(n*m) 空间",
    "SW局部比对": "O(n*m) 时间, O(n*m) 空间",
}

for op, complexity in operations.items():
    print(f"  {op}: {complexity}")

# === 生信算法速查 ===

# 数据结构选择
# 查找多 → 哈希表(dict/set)
# 有序数据 → 二分查找(bisect)
# 网络关系 → 图(邻接表)
# 序列比对 → 动态规划(DP矩阵)

# 复杂度口诀
# 一次遍历 O(n)，排序 O(n log n)
# 嵌套循环 O(n²)，二分 O(log n)
# 哈希查找 O(1)，BFS/DFS O(V+E)

# 生信算法对应
# BWA → BWT + FM-index + SW延伸
# BLAST → 种子(word) + 延伸(gapped alignment)
# SPAdes → de Bruijn图 + 欧拉路径
# GATK → Pair-HMM + 贝叶斯
# samtools sort → 外部归并排序