MIT-6-047-基因组机器学习笔记-全-

MIT 6.047 基因组机器学习笔记（全）

课程 P1：L1 - 计算生物学导论 🧬

在本节课中，我们将要学习计算生物学（也称为基因组学机器学习）的课程概述。课程将围绕剖析人类疾病的分子电路展开，涵盖从理解DNA基序、基因网络、表观基因组学，到解析疾病机制、理解跨表型疾病、解读医疗记录与数据整合，再到单细胞水平的基因组与基因分析，以及介导分析和理解遗传变异如何通过基因调控层影响疾病。我们的目标是介绍课程概览、教学团队、学生情况，并探讨课程的基础与前沿二分法、教材、作业、测验以及最终项目团队。

课程概述与目标 🎯

上一节我们介绍了课程的整体框架，本节中我们来看看课程的具体组织和内容安排。

课程由Manolis Kellis教授主讲，他是SEAS的教授，同时在Broad研究所担任副研究员。他的研究聚焦于疾病机制，运用基因组学、学习算法、统计学和人工智能技术，应用于癌症、大脑研究，以理解基因调控、进化和单细胞组学。助教Damon和Brian也分别介绍了自己的背景和研究兴趣。

课程网站提供了所有材料，包括未来所有讲座的幻灯片、往年的视频录像、复习笔记、习题集以及项目指导。课程日程表是了解课程内容安排最重要的文件。

以下是课程的核心目标：

介绍计算生物学领域的知识：理解计算生物学中的基本问题。
理解数据科学的算法与机器学习技术：掌握应用于基因组学和健康领域的通用方法。
聚焦当前研究前沿：这是一个年轻且快速发展的领域，许多技术是近年才发明的。
强调方法的原理：课程重点在于编写代码来理解和解释数据集，理解方法如何工作，从而赋能学生开发下一代工具。
实践与研究经验：通过习题集获得编程和算法实践的动手经验，通过最终项目完成一个独立的研究项目。

课程内容与模块 📚

上一节我们明确了课程目标，本节中我们来看看课程内容是如何围绕两个核心二元性组织的。

课程内容可以从两个维度理解：

计算与生物学：关注当前重要且相关的生物学问题，并运用计算机科学和机器学习的一般技术与原理。
基础与前沿：基础部分涉及定义明确的问题和通用方法（经典内容）；前沿部分涉及复杂的当前问题和开放性问题，通常需要结合多种技术。

课程按模块组织，每个模块对应一个活跃的研究领域，并同时包含基础和前沿部分。以下是主要模块：

比较基因组学：对齐和建模基因组，涉及动态规划和隐马尔可夫模型。
基因与转录组：RNA测序、聚类和结构分析，涉及聚类和分类。
基因调控与表观基因组学：转录因子、基序、网络推断。
人类遗传变异：介绍人类遗传学、人类历史、遗传力、表达数量性状位点。
进化与系统发育：系统发育学、进化特征、全基因组复制、基因组组装。
前沿主题：个人基因组学、疾病基因组学、3D基因组学等。

每个模块的前半部分学习计算技术（基础），后半部分探讨该领域当前的研究方向（前沿）。

教材、作业与测验 📖

上一节我们了解了课程模块，本节中我们来看看支持学习的教材、作业和测验安排。

推荐教材：

《生物信息学算法》
《生物序列分析》
《模式分类》
此外，课程网站提供了一本由往年学生整理的665页课堂笔记合集，涵盖了所有讲座内容。

作业：

每个作业对应一个或两个讲座，强调特定技术。
包含实践问题（编写代码、下载数据集、分析解释结果）和理论问题（用纸笔探索算法、统计或机器学习细节）。
提交格式为Jupyter Notebook。
截止时间为周一晚上11:59，有弹性的迟交政策。

测验：

测验旨在评估对前四个模块材料的掌握程度。
题型包括知识题（判断正误并证明、选择题）、深度理解题（简答题）、实践问题（解决简单算法）和设计题（设计新算法或修改现有算法）。

最终项目 🚀

上一节我们介绍了课程的基础学习部分，本节中我们来看看课程的核心实践环节——最终项目。

最终项目的目标是让学生在计算生物学领域进行原创性研究。项目过程模拟完整的科研流程：

提出并完成一个独立研究项目：学习如何构建问题、搜集相关文献和数据、使用新算法解决问题、并从生物学角度解释结果。
展示想法与研究：通过撰写研究提案、科学论文和进行学术报告来提升沟通能力。
团队合作：在具有互补技能的小组中工作。
同行评审：评审其他团队的提案，提出改进建议，并学习如何接收和整合反馈。

项目将分为多个阶段（回合），配有详细的指导：

第1周：提交自我介绍视频。
第2周：选择并介绍一篇高水平论文及其数据集。
第3周：列出最感兴趣的前三个项目想法。
第4-5周：组建团队并提交初步项目提案。
第6周：进行提案的同行评审。
第7周：展示初始的端到端分析流程演示。
第11周：提交中期报告。
第14周：提交最终报告、幻灯片和录制最终报告视频。

评分标准将围绕挑战性、相关性、原创性、完成度和展示效果。

为什么选择计算生物学？💡

上一节我们详细了解了课程要求，本节中我们回到一个根本问题：为什么计算生物学如此独特和令人兴奋？

与计算天文学、地质学或社会科学相比，计算生物学的独特之处在于：

生命本身是数字化的：DNA是一种由A、C、G、T构成的数字代码，我们正在使用计算来理解运行在每个细胞中的数字计算机。
数据量巨大且不断增长：存在海量的生物学数据集等待分析。
能够利用他人数据高效研究：可以在没有实验预算的情况下进行深入研究。
整合多组学数据的需求：结合基因组、转录组、表观基因组等多层面数据以获得更深入的见解。
提出并验证假设：可以快速生成假设并指导实验验证，极大地压缩了从假设到测试的周期。
直接影响人类健康与生活：研究成果具有重要的实际应用价值，能够真正改善人类健康。

生物学基础速成 🧪

上一节我们探讨了领域的独特性，本节中我们为初学者快速回顾本课程涉及的生物学核心概念——中心法则。

中心法则描述了遗传信息的流动方向：DNA -> RNA -> 蛋白质。

DNA（脱氧核糖核酸）：
- 承载遗传信息的双螺旋分子。
- 碱基配对原则：A与T配对，G与C配对（A-T， G-C）。这为复制和互补提供了基础。
- 人类每个细胞约有2米长的DNA，通过紧密折叠储存在细胞核中。
表观基因组学：
- 在不改变DNA序列的情况下，调节基因活性的机制。
- DNA缠绕在组蛋白上形成核小体，进而折叠成染色质。
- 组蛋白尾部和DNA本身的化学修饰（如甲基化、乙酰化）构成了“表观遗传密码”，决定了不同细胞类型中哪些基因被开启或关闭。
RNA（核糖核酸）：
- 作为信使（mRNA），将DNA的遗传信息传递到蛋白质合成场所。
- 除了信使功能，许多RNA分子本身具有重要功能（如tRNA, rRNA, miRNA），它们可以折叠成复杂结构参与基因调控。
蛋白质：
- 由20种氨基酸组成，根据mRNA的指令合成。
- 蛋白质折叠成特定三维结构以执行几乎所有的细胞功能。
基因调控网络：
- 增强子、启动子、转录因子、非编码RNA等通过复杂相互作用，精确控制基因在特定时间、特定细胞中的表达水平，从而形成多样的细胞类型和功能。
遗传与变异：
- 突变可发生在DNA、RNA或蛋白质水平，影响调控网络，导致疾病。
- 人类基因组有约32亿个碱基对，个体之间存在数百万个遗传变异。
- 全基因组关联研究致力于发现与疾病相关的遗传变异，并解读其分子机制。
进化：
- 地球上所有生命源于一个共同祖先。
- 通过比较不同物种的基因组（比较基因组学、系统发育学），可以理解基因和调控元件的进化过程，并识别功能区域。

总结 📝

本节课中我们一起学习了“计算生物学导论”课程的全貌。我们了解了课程的目标是融合计算技术与生物学问题，兼顾基础方法与前沿研究。课程通过模块化教学，涵盖基因组学、转录组学、表观基因组学、网络生物学、遗传学和进化等多个核心领域。我们明确了教材、作业、测验的要求，并详细了解了最终项目作为核心实践环节的分阶段规划。我们还探讨了计算生物学领域的独特魅力，并回顾了必要的分子生物学中心法则知识。希望本课程能带领大家进入这个激动人心的领域，通过学习和研究，为理解人类疾病和推动生命科学进步做出贡献。

课程 P10：L10 - 调控基因组学与基序发现 🧬

在本节课中，我们将要学习调控基因组学的核心概念——基序。我们将探讨什么是调控基序，如何从共调控基因集中发现它们，以及如何利用进化特征在全基因组范围内进行从头发现和单个实例识别。最后，我们会介绍如何利用高通量实验从头解析调控区域。

什么是调控基序？🧩

调控基序是基因调控的基本构建模块。例如，在酵母中代谢半乳糖的GAL1基因，其调控区域包含GAL4启动子，该启动子结合一个回文序列基序。当转录因子结合这些基序时，它们会招募大量机制，最终在基因启动子处启动转录和翻译。

基因根据环境变化开启或关闭。基因组中没有直接的寻址机制，基因内部包含被称为基序的序列标签。专门的蛋白质（即转录因子）识别这些标签。基序发现之所以困难，是因为基序通常很短，有时具有简并性（例如，允许A或T），可以包含任何核苷酸组合，并且可以在目标基因上游或下游的可变距离处发挥作用。

基序无处不在：它们存在于启动子区、增强子区、剪接信号（供体、受体、分支点）以及外显子和内含子内部。此外，3‘ UTR中也存在被microRNA识别的基序。所有这些基序的共同点是：它们是可以在不同位置出现且可能简并的模式。

如何发现基序？🔍

转录因子利用DNA结合域识别基因组中特定的DNA序列。这种识别发生在闭合的DNA双螺旋上，转录因子与特定残基接触，接触的具体情况决定了基序的形状和特异性。例如，如果一个位置标明“T或C”，则意味着转录因子在该位置喜欢识别T和C共有的原子，这导致了基序的简并性。

我们通常使用基序标识图来可视化和理解基序。基序的序列直接由转录因子的结构决定。这些调控基序在基因调控中至关重要。在全基因组关联研究中，绝大多数与复杂性状相关的遗传变异位于非编码区，它们通常通过破坏这些调控基序来发挥作用。

基序的表示：位置权重矩阵

基序通常使用位置权重矩阵（PWM）来总结转录因子的结合特异性。PWM告诉我们，在实验确定的结合位点中，每个核苷酸在每个位置出现的频率。

我们可以将基序标识图视为一个生成模型：它概括了产生基序实例的生成过程。标识图中每个位置的高度是信息量的度量，它表示在该位置有多少比特的信息。这直接源于信息论，并假设每个位置独立于其他位置（位置独立性）。

基序发现的挑战与方法 🧪

基序发现面临三个核心对象的挑战：

基序：生成模型。
靶标：基序的单个实例。
调控因子：结合并识别基序的蛋白质。

根据已知信息的不同，有不同的发现路径：

已知靶标，寻找基序：在共调控基因或结合区域中寻找富集的模式。
已知基序，寻找靶标：使用进化特征或序列组成进行全基因组扫描。
已知调控因子，寻找基序：使用SELEX或蛋白质结合微阵列等实验技术。
已知调控因子，寻找靶标：使用ChIP-chip或ChIP-seq，或通过扰动实验观察基因表达变化。

本节课我们将重点介绍两种基于区域发现基序的计算方法：期望最大化算法和吉布斯采样。

方法一：基于共调控基因集的基序发现 🧬

当我们有一组共表达或功能相关的基因时，我们假设驱动它们的是结合在每个基因启动子区的同一个转录因子。计算问题是：给定一组共调控基因，如何在它们的启动子区域找到共同的基序。

一种强大的方法是迭代优化：如果我们知道基序实例的确切起始位置，就很容易通过计数推导出PWM（基序）；反之，如果我们知道PWM，就可以通过扫描序列找到起始位置。当我们两者都不知道时，可以迭代进行估计。

期望最大化算法

期望最大化算法是一种迭代方法：

初始化：随机猜测一个基序PWM或起始位置。
E步（期望步）：给定当前PWM，估计每个序列中基序起始于每个位置的概率（后验概率）。
M步（最大化步）：利用E步计算出的加权起始位置，更新PWM（最大似然估计）。
迭代：重复E步和M步，直到PWM和起始位置收敛。

核心公式：计算序列 i 中基序起始于位置 j 的后验概率 z_ij：
P(z_ij=1 | 序列, PWM) ∝ P(序列片段 | PWM) * P(起始先验)

期望最大化算法会收敛到似然函数的一个局部极大值，但它对初始值敏感，可能陷入局部最优。

吉布斯采样

吉布斯采样是期望最大化的一种随机类比，属于马尔可夫链蒙特卡洛方法。它同样迭代，但在每一步不是使用所有起始位置的加权平均，而是随机抽样一个起始位置。

步骤：

随机选择每个序列中的一个起始位置。
随机移除一个序列，用剩余序列的起始位置计算一个临时PWM。
根据这个临时PWM，在移除的序列中抽样一个新的起始位置（概率与匹配得分成正比）。
将该序列与新起始位置放回，重复步骤2-3。

吉布斯采样通过允许“侧向移动”，更不容易陷入局部最优，有时能找到全局更优解。流行的工具如AlignACE和BioProspector就使用了吉布斯采样。

上一节我们介绍了在已知共调控基因集的情况下，如何利用期望最大化或吉布斯采样发现基序。接下来，我们将看看另一种不依赖于特定基因集的全基因组方法。

方法二：基于进化特征的从头基序发现 🌍

这种方法利用比较基因组学。其动机是：单个基序实例在基因组的一个位置可能因偶然或其它功能（如RNA结构）而保守。但真正的调控基序会在基因组多个位置反复出现，并且这些实例会表现出独特的进化特征。

核心思想：在全基因组范围内评估基序模式的保守性，而不仅仅是看单个位点。

早期研究（以酵母GAL4为例）：

发现GAL4基序在基因间区的保守性（13%）远高于具有类似属性的随机控制基序（2%），富集了6倍。
相反，GAL4基序在蛋白质编码区内的保守性反而低于控制基序，表明它倾向于不在编码区结合。

方法流程：

枚举种子：系统性地枚举所有可能的短序列模式（如三核苷酸间隔三核苷酸的模式）。
评估保守性：对于每个种子模式，检查它在全基因组多个同源区域是否保守。
扩展与聚类：将保守的种子模式向两侧扩展，纳入附近同样保守的核苷酸，形成更长的“超级基序”。然后将相似的超级基序聚类，得到最终的基序集合。
功能注释：通过与已知基序数据库比对、基因本体富集分析、基因表达谱关联等方法，推测新发现基序的功能。

这种方法能够完全从头发现一个物种中的调控基序目录，无需预先知道转录因子或共调控基因集。

识别单个基序实例：进化足迹 👣

现在我们知道了某个转录因子的结合基序（PWM），如何在全基因组中精确找到它的功能结合位点（靶标）？

实验方法（如ChIP-seq）受限于抗体、组织和发育阶段。计算方法可以利用进化信息。

挑战：基序实例并非总是完美保守，它们可能发生突变（只要在基序简并性允许范围内）、在附近移动，甚至在部分物种中丢失。

解决方案：分支长度评分法

独立扫描：在多个物种的基因组中，独立扫描给定的基序PWM，找到所有超过一定阈值的匹配位点。
构建系统发育树：对于基因组中的每个同源区域，检查每个物种中是否存在该基序的匹配。
计算分支长度：在系统发育树上，对存在基序匹配的分支长度进行求和。这个总分就是该位点的“分支长度评分”（BLS）。BLS越高，说明该基序实例在进化中保守的时间越长。

处理背景噪声：置信度评分
不同基序（长短、组成不同）的BLS直接比较没有意义。短基序可能偶然保守。因此，需要将BLS转换为统一的置信度。

步骤：

生成控制基序：通过对原始基序PWM进行多次“洗牌”（保持核苷酸频率等属性），生成数百个控制基序。
计算控制基序的BLS分布：用同样的方法计算所有控制基序在全基因组中的BLS分布，这代表了“偶然保守”的背景噪声水平。
计算置信度：对于原始基序的每个BLS值，计算达到该BLS的实例中，有多少比例超过了控制基序的预期数量。置信度 = 1 - 假发现率。

置信度评分提供了一个与基序属性无关的统一指标。利用置信度，我们可以：

发现高置信度的基序实例。
验证其生物学合理性：例如，高置信度的转录因子基序实例富集在启动子和5‘ UTR，而microRNA基序实例富集在3’ UTR。
与实验数据（如ChIP-seq峰）进行对比验证，发现置信度越高的预测位点，与实验数据的重叠富集越强。

通过这种方法，我们可以大规模地鉴定转录因子与其靶基因之间的调控关系，从而构建基因调控网络。

总结 📚

本节课我们一起学习了调控基因组学的核心——基序。

基序的定义与重要性：调控基序是短的、可能简并的DNA序列模式，是转录因子结合的位点，构成了基因调控电路的基础，并与许多疾病相关。
基序的表示：使用位置权重矩阵和基序标识图来概括基序的特异性。
基于共调控基因集的发现：
- 期望最大化算法：迭代估计基序PWM和起始位置，可能陷入局部最优。
- 吉布斯采样：一种随机采样方法，能更有效地探索解空间，避免局部最优。
基于进化特征的从头发现：通过在全基因组范围内评估短序列模式的重复性和保守性，无需先验知识即可发现新的调控基序。
单个实例的识别：利用多物种比较基因组学，通过分支长度评分和置信度计算，在全基因组中高精度地预测转录因子的功能结合位点。

这些计算方法为我们理解基因调控的序列基础、解析非编码区遗传变异的机制以及构建基因调控网络提供了强大的工具。

📚 课程 P11：L11 - 网络分析基础与应用

在本节课中，我们将学习网络分析的基础知识、核心算法及其在生物学等领域的应用。我们将从图论基础开始，探讨网络的不同类型和表示方法，然后深入到网络的全局与局部性质、线性代数在网络分析中的应用，以及社区发现等高级主题。

🌐 网络与图论基础

网络无处不在，从社交关系到生物调控，都可以用图来抽象表示。一个图由一组节点（代表实体）和一组边（代表关系）构成。

以下是图的基本类型：

无向图：边没有方向，关系是对称的。
有向图：边有方向，关系是非对称的。
加权图：边带有权重，表示关系的强度。
简单图：没有自环和多重边。

我们可以用邻接矩阵来表示图。对于一个有 n 个节点的图，邻接矩阵 A 是一个 n x n 的矩阵，其中 A[i][j] 表示从节点 i 到节点 j 的边的权重（对于无权重图，通常为1或0）。

# 一个无向无权图的邻接矩阵示例
A = [
    [0, 1, 1, 0],
    [1, 0, 1, 1],
    [1, 1, 0, 0],
    [0, 1, 0, 0]
]

节点的度是指与其相连的边的数量。在有向图中，分为入度（指向该节点的边数）和出度（从该节点指出的边数）。

🔍 基本图算法：广度优先与深度优先搜索

上一节我们介绍了图的表示，本节中我们来看看两种基础的图遍历算法。

广度优先搜索 从一个起始节点开始，逐层探索其邻居。它常用于寻找最短路径。
深度优先搜索 则沿着一条路径深入探索到底，然后回溯，常用于拓扑排序或寻找连通分量。

以下是两种算法的核心思想对比：

BFS：使用队列，保证按距离起始点的层次顺序访问节点。
DFS：使用栈（或递归），尽可能深地探索图的分支。

🌍 网络的涌现性质

现实世界中的网络，无论属于哪个领域，都展现出一些共同的宏观和微观性质。

全局性质包括：

小世界现象：网络中任意两个节点之间通常存在较短的路径。
无标度特性：网络中少数节点（枢纽）拥有大量的连接，而大多数节点连接很少。节点度的分布常遵循幂律分布。
聚类系数：衡量节点的邻居之间也相互连接的概率，反映了网络的局部紧密程度。
网络模体：在真实网络中反复出现、具有特定功能的小型连接模式。

局部性质主要指节点的中心性，用于衡量节点的重要性。常见度量包括：

度中心性：节点的连接数。
介数中心性：经过该节点的最短路径数量。
接近中心性：该节点到网络中所有其他节点平均距离的倒数。
特征向量中心性：一个节点的中心性是其邻居中心性的加权和。

📊 网络的线性代数视角

我们可以用线性代数的工具来分析和操作网络。邻接矩阵 A 本身就是一个线性算子。

将网络视为矩阵后，我们可以分析其特征值和特征向量。对于一个方阵 A，如果存在一个非零向量 v 和一个标量 λ，使得 A * v = λ * v 成立，那么 v 就是 A 的特征向量，λ 是对应的特征值。特征向量指示了网络中的主要“振动模式”或影响力扩散的方向。

对于非方阵或更一般的矩阵，我们使用奇异值分解。SVD 将任意矩阵 M 分解为三个矩阵的乘积：M = U * Σ * V^T。其中 U 和 V 是正交矩阵，Σ 是对角矩阵，其对角线上的值称为奇异值。SVD 是进行低秩近似的强大工具，可以用于降维和去噪。

🎯 降维与稀疏表示

上一节我们介绍了SVD，本节中我们来看看如何利用它进行降维，并引入稀疏性约束。

主成分分析 是常用的线性降维方法，它本质上是数据协方差矩阵的特征值分解，寻找方差最大的投影方向。然而，PCA 得到的主成分是所有原始变量的线性组合，通常不稀疏，难以解释。

稀疏PCA 通过引入 L1 正则化（Lasso）约束，迫使主成分的载荷向量变得稀疏，即只有少数变量具有非零权重。这提高了模型的可解释性。其优化问题可以表述为在重建误差和系数稀疏性之间取得平衡。

除了线性方法，还有非线性降维，如 t-SNE。t-SNE 专注于保持高维数据点之间的局部邻近关系，将其映射到低维空间（如2D或3D），特别适用于可视化复杂的高维数据集，如单细胞基因表达数据。

🧩 网络社区发现：谱聚类

网络中的一个核心任务是发现其中紧密连接的子图，即社区。谱聚类 是一种基于图拉普拉斯矩阵特征向量的有效社区发现方法。

首先，我们定义图的拉普拉斯矩阵 L = D - A，其中 D 是度矩阵（对角矩阵），A 是邻接矩阵。拉普拉斯矩阵具有一些重要性质，例如所有行和为零，并且是半正定矩阵。

谱聚类的关键步骤是：

计算拉普拉斯矩阵 L。
计算 L 的前 k 个最小特征值对应的特征向量。
将这些特征向量按列排列，形成一个 n x k 的矩阵。
将这个矩阵的每一行视为 k 维空间中的一个点，使用传统的聚类算法（如K-Means）对这些点进行聚类。

这样，图节点的聚类问题就转化为了特征向量空间中点的聚类问题。拉普拉斯矩阵的第二小特征值对应的特征向量（称为费德勒向量）通常给出了图的一个最优二分切割。

📝 总结

本节课中我们一起学习了网络分析的核心内容。我们从图的基本定义和表示出发，介绍了BFS和DFS等基础算法。然后，我们探讨了现实世界网络的小世界、无标度等涌现性质，以及度量节点重要性的各种中心性指标。接着，我们从线性代数的视角重新审视网络，学习了特征值分解、奇异值分解及其在低秩近似中的应用。在此基础上，我们介绍了带有稀疏约束的PCA和非线性降维方法t-SNE。最后，我们学习了如何使用谱聚类这一基于图拉普拉斯矩阵的强大方法来发现网络中的社区结构。这些工具为理解和分析复杂的生物网络、社交网络等提供了坚实的基础。

🧠 深度学习课程 P12：从基础到基因组学应用

在本节课中，我们将学习深度学习的基础概念、核心架构及其在基因组学领域的应用。课程内容涵盖从受监督学习到无监督学习，再到现代深度学习架构的演变。

🧬 课程概述

本节课将介绍深度学习的基本原理，包括受监督学习、无监督学习、卷积神经网络、循环神经网络、迁移学习以及生成对抗网络。我们还将探讨这些技术在调控基因组学中的应用，并提供实用的计算工具和代码示例。

🧠 受监督学习与神经网络

上一节我们介绍了课程的整体框架，本节中我们来看看受监督学习与神经网络的基础概念。

人类大脑在处理视觉信息时，能够从低级像素中抽象出高级概念，如物体和形状。这种能力启发了早期神经网络的设计。

人工神经元与感知机

人工神经元模拟生物神经元，接收多个输入，每个输入乘以一个权重，然后求和。如果总和超过某个阈值，神经元则“激活”并产生输出。

公式：
输出 = 激活函数( Σ (权重_i * 输入_i) )

激活函数

线性函数的表达能力有限。为了学习复杂的非线性关系，我们引入了激活函数。

以下是几种常见的激活函数：

Sigmoid函数：早期常用，将输入映射到0到1之间，模拟神经元的饱和状态。
Softplus函数：从0开始平滑过渡到线性区域。
修正线性单元：目前最常用，在输入小于0时输出为0，大于0时输出等于输入。

梯度下降与学习

学习的目标是调整权重，使网络的预测输出尽可能接近真实输出。

核心思想：计算预测输出与真实输出之间的误差，然后根据误差相对于每个权重的梯度来更新权重。

权重更新公式：
W_t = W_{t-1} - ε * (∂误差/∂W) + 动量项 + 权重衰减项

其中：

ε 是学习率，控制更新步长。
∂误差/∂W 是误差相对于权重的梯度。
动量项利用之前的更新方向来加速收敛并减少震荡。
权重衰减项惩罚较大的权重值，防止过拟合。

反向传播

反向传播是一种高效计算网络中所有权重梯度的方法。它首先进行前向传播计算输出和误差，然后从输出层向输入层反向传播误差，并利用链式法则计算每个权重的梯度。

避免过拟合

神经网络参数众多，容易对训练数据过拟合，导致在新数据上表现不佳。

以下是避免过拟合的几种策略：

划分数据集：将数据分为训练集、验证集和测试集。用验证集监控模型性能并决定何时停止训练，用测试集进行最终评估。
正则化：如L1、L2正则化，惩罚大的权重值。
Dropout：在训练过程中，随机“丢弃”网络中的一部分神经元。这相当于同时训练多个子网络，是一种强大的正则化方法，能防止网络对某些特定神经元过度依赖。

🔍 无监督学习与深度信念网络

上一节我们介绍了受监督学习，本节中我们来看看如何在没有标签的情况下学习数据的内在表示。

受限玻尔兹曼机

受限玻尔兹曼机是一种无监督学习模型，包含一层可见单元和一层隐藏单元。它学习可见单元和隐藏单元之间的联合概率分布，从而发现数据中的特征。

深度信念网络

通过堆叠多个RBM，可以构建深度信念网络。较低层的RBM学习低级特征，较高层的RBM基于这些低级特征学习更抽象的特征。

学习技巧

模拟退火：在训练初期引入“噪声”，帮助模型跳出局部最优解；随着训练进行逐渐“降温”，使模型稳定在更好的解上。
吉布斯采样：一种从RBM中采样的方法，用于学习和推理。
醒眠算法：交替进行“醒”阶段和“睡”阶段来训练网络。

🎨 自编码器

自编码器是一种通过“自我重建”来学习数据表示的神经网络。

核心思想：将输入数据压缩到一个低维的“瓶颈”层，然后再重建回原始维度。通过最小化重建误差，网络被迫学习输入数据最重要、最具代表性的特征。

自编码器可以用于数据降维、去噪或作为其他任务的预训练步骤。

🖼️ 现代深度学习架构：卷积神经网络

上一节我们探讨了无监督表示学习，本节中我们来看看专门为图像数据设计的强大架构——卷积神经网络。

CNN的灵感来源于视觉皮层，它通过卷积核在图像上滑动来检测局部特征。

核心操作

卷积：使用一组可学习的滤波器在输入图像上滑动。每个滤波器负责检测一种特定的局部模式。
激活函数：对卷积结果应用非线性激活函数。
池化：对局部区域进行下采样，常用最大池化或平均池化。这能提供平移不变性并减少参数数量。

一个典型的CNN由多个“卷积-激活-池化”块堆叠而成，后面连接一个或多个全连接层用于最终分类。

关键优势

参数共享：同一个滤波器在整个图像上使用，极大地减少了参数量。
局部连接：每个神经元只与输入图像的局部区域连接，这符合图像的空间局部性。
层次化特征：底层学习边缘、颜色等低级特征，高层组合这些低级特征形成物体部件等高级特征。

经典网络演进

从AlexNet、VGGNet到GoogleNet、ResNet，网络层数不断加深。ResNet引入了“残差连接”，允许信息跨层直接传播，有效缓解了深度网络中的梯度消失问题，使得训练成百上千层的网络成为可能。

⏳ 现代深度学习架构：循环神经网络

上一节我们学习了处理空间数据的CNN，本节中我们来看看处理序列数据的循环神经网络。

RNN专为处理序列数据设计，如文本、语音、时间序列等。其核心思想是网络不仅接收当前输入，还接收上一个时间步的“隐藏状态”，从而拥有记忆能力。

长期依赖问题与LSTM

基本RNN难以学习长序列中的长期依赖关系，因为梯度在反向传播时会消失或爆炸。

长短期记忆网络通过引入“门”机制来解决这个问题：

遗忘门：决定从细胞状态中丢弃哪些信息。
输入门：决定哪些新信息存入细胞状态。
输出门：基于细胞状态决定输出什么。

LSTM单元能够有选择地保留和传递信息，从而有效地学习长期依赖。

双向RNN与架构组合

双向RNN同时从序列的前向和后向处理信息，能获得更丰富的上下文。RNN/LSTM常与CNN组合使用，例如用CNN处理视频的每一帧，再用RNN学习帧之间的时序关系。

🔄 迁移学习

迁移学习旨在将一个任务上训练好的模型知识，迁移到另一个相关任务上。

典型做法：

在一个大型通用数据集上预训练一个深度网络。
移除网络的最后几层。
将预训练网络的其余部分作为新任务的固定特征提取器，或者进行微调。
为新任务添加并训练一个新的输出层。

这种方法在目标领域数据量较少时特别有效。

🎭 生成对抗网络

GAN由两个相互对抗的网络组成：

生成器：学习从随机噪声生成逼真的数据。
判别器：学习区分真实数据和生成器产生的“假”数据。

两个网络在对抗中共同进步：生成器努力生成更逼真的数据以骗过判别器，判别器则努力提高鉴别能力。GAN能生成极其逼真的图像、音频等，也可用于数据增强。

🧬 深度学习在调控基因组学中的应用

现在，我们将上述工具应用于基因组学领域。

序列表示与卷积神经网络

DNA序列可以通过独热编码表示为一个4行（A、C、G、T）的矩阵。将这个矩阵视为一个“图像”，CNN的卷积核就相当于在自动学习DNA上的“基序”。

低级卷积层：学习类似于转录因子结合位点的简单序列模式。
高级卷积层：组合简单模式，识别更复杂的调控元件。
全连接层：基于学习到的特征预测特定生物学事件，如蛋白质结合、染色质可及性等。

多任务学习与蛋白质结构预测

可以训练一个网络同时预测多个细胞类型中转录因子的结合，共享的低层特征有助于提升各任务的性能。类似地，可以从蛋白质氨基酸序列预测其三维结构域和功能。

💻 深度学习计算工具与实践

现在有许多库简化了深度学习的实现：

底层框架：PyTorch, TensorFlow。它们提供自动微分和GPU加速。
高级API：Keras。它构建在TensorFlow等之上，提供了更简洁的接口。

以下是一个使用Keras构建简单CNN来预测DNA序列功能的示例代码框架：

import keras
from keras.layers import Conv2D, MaxPooling2D
from keras.models import Sequential
from keras.layers import Activation, Dropout, Flatten, Dense

# 定义模型参数
input_shape = (4, 150, 1) # 4个通道（A,C,G,T），长度150
num_filters = 10
conv_width = 4
pool_width = 5

# 构建模型
model = Sequential()
model.add(Conv2D(num_filters, (conv_width, 4), padding='valid', input_shape=input_shape))
model.add(Activation('relu'))
model.add(MaxPooling2D(pool_size=(pool_width, 1)))
model.add(Flatten())
model.add(Dense(1))
model.add(Activation('sigmoid'))

# 编译模型
model.compile(loss='binary_crossentropy', optimizer='adam')

# 训练模型
# model.fit(X_train, y_train, batch_size=1000, epochs=20, validation_data=(X_val, y_val))

# 进行预测
# predictions = model.predict(X_test)

📚 课程总结

本节课我们一起学习了深度学习的广阔领域：

神经网络基础：感知机、激活函数、梯度下降、反向传播以及应对过拟合的技术。
无监督学习：深度信念网络和自编码器如何学习数据的内在表示。
现代架构：
- 卷积神经网络：用于处理图像和序列等具有空间/局部结构的数据。
- 循环神经网络与LSTM：用于处理时间序列和文本等序列数据。
- 迁移学习：利用预训练模型解决数据有限的新任务。
- 生成对抗网络：用于生成新数据和数据增强。
基因组学应用：如何将CNN等模型应用于DNA和蛋白质序列分析，以预测调控功能。
实践工具：介绍了PyTorch、TensorFlow和Keras等库，并给出了一个简单的代码示例。

深度学习是一个强大且快速发展的工具集，正在推动从计算机视觉到基因组学等多个领域的进步。

群体遗传学课程 P13：人类遗传变异与群体差异 🧬

在本节课中，我们将学习群体遗传学的基础知识，重点是人类遗传变异的测量、理解与量化。我们将回顾人类遗传学的历史，探讨不同类型的遗传变异，并学习如何检测和解读这些变异。此外，我们还将了解单倍型、连锁不平衡、基因型填充和定相等核心概念，以及它们如何帮助我们理解人类群体历史和关联研究。

遗传学历史与基础概念

上一节我们概述了课程内容，本节中我们来看看遗传学思想的发展历程。

遗传的概念早已被理解。人类长期对动植物进行选择性育种，例如将具有攻击性的动物驯化为人类最好的朋友，或将不可食用的类蜀黍培育成今天的主要粮食作物玉米。

古希腊时期已有关于物种起源和变化的早期理论。
亚里士多德是第一个对物种进行分类的学者。
拉马克提出了“用进废退”和获得性遗传的观点，但这缺乏现实基础。

达尔文的理论具有革命性。他提出了物种连续性的概念，认为随机突变是多样性的来源，自然选择在此基础上发挥作用。然而，达尔文的理论并不完整，他仍持有一些融合遗传和泛生论的观点。

与达尔文同时代，孟德尔首次认识到非融合的、离散的遗传单位（基因），以及显性与隐性等位基因的概念。他的独立分配定律揭示了遗传的“数字化”本质。然而，孟德尔的理论被忽视了约50年，因为当时观察到的大多是连续的表型性状，难以与离散遗传的概念调和。

直到20世纪，费希尔等人将孟德尔遗传与连续性状变异结合起来，认识到连续的表型变异可以由多个具有微小效应的孟德尔位点共同解释。

在DNA双螺旋结构被发现之前，理论上的重大飞跃已经发生。摩尔根和他的学生斯特蒂文特通过连锁作图，首次构建了遗传图谱，将基因定位在染色体上。

基于连锁作图的方法，20世纪80年代发生了孟德尔性状作图的革命，可以定位强效应的单基因孟德尔突变。随后，人类基因组测序的完成，以及通过大规模、非亲缘个体的全基因组关联研究，使我们能够发现具有微小效应的常见遗传变异。我们现在认识到，基因对我们的生理和认知有着巨大影响。

遗传变异的类型与影响

上一节我们回顾了遗传学思想的演进，本节中我们来看看遗传变异的具体形式及其生物学意义。

人类个体间99.9%的DNA是相同的。个体间的所有差异和疾病易感性都存在于剩下的不到0.1%的基因组中。这些变异有多种形式：

单核苷酸多态性：最常见，约每1000个碱基出现一次。公式表示为：SNP位置：参考碱基 -> 替代碱基。
插入与缺失：约每10,000个核苷酸出现一次。
短串联重复序列：约每1,000个碱基出现一次。
结构变异/拷贝数变异：通常大于5000个碱基对，缺失约每百万个碱基出现一次。

SNP的影响取决于其位置。在蛋白质编码区，约三分之一的变异会改变氨基酸，可能严重影响功能。例如，导致镰状细胞贫血的变异（HBB基因：Glu6Val）虽然有害，但在疟疾流行地区却被正向选择，因为杂合子个体对疟疾有抵抗力。

其他类型的变异也能导致严重疾病，例如亨廷顿病（CAG重复扩展）和囊性纤维化（CFTR基因缺失）。

等位基因的分类与频率谱

上一节我们介绍了变异的类型，本节中我们来看看如何描述和分类这些变异。

描述等位基因有多种方式：

参考等位基因 vs. 替代等位基因：是否与参考基因组序列匹配。
主要等位基因 vs. 次要等位基因：在群体中的频率高低。
祖先等位基因 vs. 衍生等位基因：与人类-黑猩猩最近共同祖先的序列比较。
风险等位基因 vs. 保护等位基因：与疾病关联的方向（依赖于环境）。

根据等位基因频率，可将其分类：

常见等位基因：频率 > 5%
低频等位基因：频率 0.5% - 5%
罕见等位基因：频率 < 0.5%
私有或新生变异：仅存在于单一个体或部分细胞（体细胞突变）

等位基因的频率与其效应大小通常存在权衡。常见等位基因通常效应微弱，而强效应等位基因通常罕见。这是因为自然选择会清除有害的常见变异。这是一个连续谱，大多数罕见变异效应也微弱，只是我们更容易发现强效应的罕见变异。一个罕见的例外是载脂蛋白E ε4等位基因，它常见且对阿尔茨海默病有强风险效应。

遗传变异的表示、检测与项目

上一节我们讨论了等位基因的频谱，本节中我们来看看如何表示、存储和检测这些变异。

人类是二倍体生物，每个个体在每个变异位点有两个等位基因（分别来自父母）。这些变异共同存在于作为单位遗传的单倍型中。基因型告诉我们在一个位点上拥有0、1或2个替代等位基因。从基因型推断单倍型的过程称为定相。

检测变异的第一步是对大量个体进行测序以发现变异。将测序读数比对到参考基因组后，需要区分真正的变异与测序错误，这个过程称为变异检测。GATK等工具使用生成模型和隐马尔可夫模型来推断最可能的单倍型和基因型。

发现变异后，可以针对已知变异设计更经济的基因分型芯片进行大规模分型。国际人类基因组单体型图计划（HapMap）和千人基因组计划等重大项目系统性地绘制了人类遗传变异图谱。千人基因组计划对全球26个亚群体的2500个低深度全基因组进行了测序，建立了庞大的变异目录和用于定相、填充的统计工具。

数据显示，非洲群体保留了最多的遗传多样性，而非洲以外群体的变异只是非洲多样性的一部分，这支持了“走出非洲”的迁徙模型。

单倍型、连锁不平衡与重组热点

上一节我们了解了如何检测变异，本节中我们来看看变异在基因组中是如何被共同继承的。

在减数分裂过程中，同源染色体之间会发生重组，这是染色体正确分离所必需的。重组并非随机发生，而是集中在特定的重组热点。这些热点由PRDM9蛋白结合特定DNA基序所决定，但PRDM9的结合会导致该基序被破坏，从而驱动了PRDM9蛋白及其识别基序的快速协同进化。

由于重组热点稀少，基因组被分割成单倍型区块。区块内的变异倾向于共同遗传，这种非随机关联称为连锁不平衡。

衡量LD的常用指标包括：

D值：衡量观察到的单倍型频率与随机组合预期频率的偏差。公式：D = P(AB) - P(A)P(B)
D‘值：D值标准化后的结果，范围在0到1之间。公式：D' = D / Dmax
r²值：两个位点等位基因之间的相关系数平方。公式：r² = D² / (P(A)P(a)P(B)P(b))

单倍型区块的存在意味着人类群体中实际存在的单倍型数量远少于理论上所有SNP组合的可能性。因此，我们只需检测每个区块中的少数几个标签SNP，就能推断出整个区块的单倍型信息。这既是祝福（降低成本），也是诅咒（难以精确定位致病变异），因为关联信号往往属于整个单倍型区块。

基因型填充与定相

上一节我们认识到只需检测部分SNP即可推断单倍型，本节中我们来看看实现这一推断的具体计算方法。

基因型填充的目标是根据已测量的部分SNP基因型，推断未测量的SNP基因型。这通常需要先解决定相问题，即确定每个等位基因分别来自哪条亲本染色体。

定相可以利用家系信息（如父母基因型）无歧义地解决部分位点。对于无亲缘关系的个体，则需要利用参考单倍型面板和连锁不平衡信息。核心思想是：个体的单倍型可以被视为从参考面板中拷贝片段拼接而成。重组事件对应于拷贝来源的切换。

基于隐马尔可夫模型的方法可以沿基因组滑动窗口，计算从不同参考单倍型拷贝的概率，从而推断出最可能的单倍型路径。一旦单倍型被确定，填充缺失的基因型就变得简单，只需从推断出的单倍型中拷贝相应等位基因即可。

群体亲缘关系、历史与祖先推断

上一节我们学习了如何推断个体的单倍型，本节中我们来看看如何利用遗传变异理解群体之间的关系和历史。

不同群体的遗传多样性不同。非洲群体拥有最高的遗传多样性，而其他群体由于“走出非洲”过程中的奠基者效应，多样性较低。通过分析基因组中变异位点的分布，可以估算历史上不同时期的有效群体大小，揭示群体扩张和瓶颈事件。

对于具有混合祖先的个体，我们可以将其基因组“描绘”成来自不同祖先群体的片段。这可以通过基于条件随机场等模型的程序（如RFMix）实现，该程序利用参考群体的等位基因频率来预测每个局部片段的来源群体。

对基因型矩阵进行主成分分析，可以揭示群体的遗传结构。例如，在欧洲样本中，前两个主成分往往对应地理上的南北和东西梯度，反映了历史上的迁徙模式。在关联研究中，这些主成分常作为协变量以控制群体分层造成的假阳性。

群体间分化指数用于量化群体间的遗传分化。通过分析线粒体DNA（母系遗传）和Y染色体DNA（父系遗传），可以追溯人类迁徙历史。研究还发现，现代非非洲人群的基因组中含有约2%的尼安德特人DNA，而某些亚洲人群则含有丹尼索瓦人的遗传成分，这证明了古代人类群体间的基因交流。

总结与展望

本节课中，我们一起学习了群体遗传学的核心内容。我们回顾了人类遗传学的历史，从孟德尔定律到现代全基因组关联研究。我们探讨了不同类型的遗传变异（SNP、Indel、STR、CNV）及其生物学影响，理解了等位基因频率与效应大小的权衡关系。

我们深入研究了单倍型、连锁不平衡和重组热点的生物学基础与计算方法，学习了如何通过基因型填充和定相从部分数据推断完整信息。最后，我们了解了如何利用遗传数据推断群体亲缘关系、历史迁徙事件和个体祖先成分，并认识到主成分分析在控制群体分层中的重要性。

这些关于遗传变异的知识为我们下一节课学习表型变异和疾病关联研究打下了坚实的基础。我们将学习如何将遗传变异与疾病、性状联系起来，开展全基因组关联研究，并理解其背后的生物学机制。

🧬 P14：L14 - GWAS与疾病解析

在本节课中，我们将学习如何理解遗传变异与疾病之间的关系。我们将从孟德尔遗传学开始，探讨连锁分析如何用于定位疾病基因，然后转向现代遗传学，研究复杂性状和全基因组关联研究。最后，我们将深入探讨如何解释GWAS结果，将遗传位点与疾病机制联系起来，并通过具体案例展示如何从遗传变异走向功能解析。

🧬 第一部分：孟德尔性状与连锁分析

上一节我们介绍了遗传变异的基础。本节中，我们来看看如何利用孟德尔遗传学的原理，通过家族研究来定位与疾病相关的基因变异。

孟德尔遗传学的核心是认识到存在离散的遗传单位，这些单位的变异是可遗传的，并导致表型差异。它区分了基因型（个体生殖细胞内的遗传信息）和表型（基因型编码的可见特征）。关键概念包括显性与隐性（例如，大写字母代表显性等位基因，小写字母代表隐性等位基因），以及不同性状的独立分配。

然而，孟德尔遇到了“连锁”的概念，即某些性状对并非独立遗传，这违反了他的第二定律（独立分配定律）。原因在于这些基因在染色体上的物理距离不够远，在减数分裂中不易发生交叉互换。连锁强度反映了基因在染色体上的空间距离，这成为了构建遗传图谱的基础。

以下是连锁分析在孟德尔疾病研究中的应用要点：

家族连锁研究：利用已知的染色体标记（多态性位点），追踪这些标记与疾病状态的共分离现象，从而定位包含致病基因和突变的染色体区域。
定位克隆：通过连锁分析确定致病基因在染色体上的大致位置后，进一步克隆并鉴定该基因。
适用性：这种方法在定位高外显率的单基因遗传病（如亨廷顿病、囊性纤维化）方面取得了巨大成功。

🧬 第二部分：复杂性状与全基因组关联研究

上一节我们介绍了用于单基因病的连锁分析。本节中，我们来看看面对大多数常见疾病和复杂性状时，我们需要不同的方法。

大多数复杂性状并非由单个高外显率变异决定，而是具有极高的多基因性，每个基因的效应很小，并受到环境因素的影响，表现为连续分布。这意味着许多基因各自施加微小的影响，与环境因素结合共同决定个体结果。

进入21世纪初，人类基因组计划的完成、高通量基因分型芯片技术的发展以及大规模国际合作，共同推动了全基因组关联研究时代的到来。

进行GWAS的基本步骤包括严格的质量控制，然后对基因组中的每一个SNP进行统计检验。

核心的统计检验是卡方检验。其基本思想是，比较病例组和对照组中某个等位基因的频率是否存在显著差异。公式可以简化为观察值与期望值偏差的平方和除以期望值：

χ² = Σ (Oᵢ - Eᵢ)² / Eᵢ

其中，Oᵢ 是观察到的计数，Eᵢ 是在无效假设（即无关联）下的期望计数。

然而，GWAS面临一个关键挑战：多重假设检验。当对数百万个SNP进行检验时，即使没有真实关联，仅凭偶然性也会产生许多看似显著的结果。为此，领域内设定了严格的“全基因组显著性”阈值（通常为 p < 5 × 10⁻⁸），以控制假阳性率。

结果通常以“曼哈顿图”展示，其中每个点代表一个SNP，X轴是染色体位置，Y轴是关联显著性的负对数。显著超越阈值的峰值代表可能的疾病关联位点。

另一个重要工具是QQ图，用于评估观察到的检验统计量分布与期望分布（在无效假设下）的一致性，帮助判断是否存在系统性误差或真实的遗传信号。

🧬 第三部分：精细定位与多研究整合

上一节我们介绍了如何发现GWAS信号。本节中，我们来看看如何从关联区域 pinpoint 到可能的因果变异，以及如何通过整合多个研究来增强发现能力。

GWAS发现的关联信号通常指向一个染色体区域，其中可能包含多个高度连锁的SNP。精细定位的目标是区分哪些SNP仅仅是“搭便车”（与真正的因果变异连锁），哪些更可能是功能性的。

以下是几种主要的精细定位方法：

基于连锁不平衡：计算区域内每个SNP与最强关联SNP（lead SNP）之间的连锁不平衡程度（如r²）。通常，与lead SNP高度连锁（r² > 0.8）的SNP被归入同一个可信区间。
惩罚回归模型：在模型中校正lead SNP的效应后，检查其他SNP是否仍有独立的关联信号（非零的β值）。
后验包含概率：利用贝叶斯方法，考虑所有SNP之间的连锁结构，计算每个SNP是因果变异的概率，并汇总成可信集合。
跨族群研究：利用不同人群（如亚洲人、非洲人）中不同的连锁不平衡结构，帮助缩小因果变异的范围。
功能注释优先：如果发现多个独立关联位点的可信区间都富集某种特定的基因组注释（如胚胎干细胞增强子），那么在新位点中，重叠该注释的SNP更可能被优先视为因果变异。

那么，GWAS发现的位点与之前连锁分析发现的位点有何关系？这主要取决于等位基因的频率和效应大小（比值比）。连锁分析对低频、强效应的风险变异有较好的检测能力，而GWAS则擅长发现常见、效应较弱（无论是风险还是保护作用）的变异。

由于大多数复杂性状涉及大量小效应位点，增加样本量是发现更多关联的关键。通过整合多个独立研究的“荟萃分析”，可以极大地提高统计效能。这在许多疾病（如精神分裂症、2型糖尿病）的研究中都观察到了一个“拐点”：当样本量累积到一定程度后，发现的位点数量急剧增加。

🧬 第四部分：从遗传位点到疾病机制

上一节我们讨论了如何找到并精细定位关联信号。本节中，我们来看看最大的挑战：如何解释这些非编码区的遗传位点，并阐明其导致疾病的生物学机制。

GWAS的主要目标是揭示疾病的生物学基础，并指导治疗。然而，超过90%的疾病关联位点位于蛋白质编码区之外，这带来了解释上的挑战：目标基因未知、因果变异未知、发挥作用的细胞类型未知、相关通路未知、作用机制未知。

为了克服这些挑战，我们需要综合利用多组学数据：

表观基因组注释：利用ENCODE、Roadmap等项目的数据，确定SNP是否落在特定细胞类型的增强子、启动子等调控元件内，从而推断其作用的细胞环境。
序列 motif 分析：检查SNP是否破坏或创建了转录因子结合位点的DNA序列模式。
进化保守性：分析SNP是否位于跨物种保守的基因组区域，这通常暗示功能重要性。
系统水平分析：
- 基因集富集分析：查看关联位点附近的基因是否富集于特定的生物学过程或通路。
- 染色质状态富集：检查关联位点是否在特定细胞类型的活跃染色质区域（如增强子）中富集。
- 整合数量性状基因座：结合表达数量性状基因座数据，寻找遗传变异与基因表达水平的关联，直接提示可能的靶基因。

通过系统性地应用这些工具，我们可以构建从遗传变异到表型的因果链。例如，利用Epigenomics Roadmap数据，可以将数百种性状的GWAS位点与特定细胞类型的增强子活性进行关联，绘制出“疾病-组织”关联图谱。研究发现，身高相关变异富集于胚胎干细胞增强子，免疫疾病富集于免疫细胞增强子，而阿尔茨海默病则意外地富集于免疫细胞（如小胶质细胞）而非神经元的增强子，这提示了该疾病的免疫基础。

🧬 第五部分：案例研究：解析FTO肥胖位点

上一节我们介绍了系统解析GWAS位点的通用策略。本节中，我们通过一个经典案例——FTO基因座与肥胖的关联，来具体展示如何完成从遗传位点到分子、细胞和生物体表型的完整解析。

FTO基因座是已知与肥胖关联最强的遗传区域，但其中包含89个常见变异，且都不改变蛋白质序列。早期研究曾认为FTO基因本身是靶标，但功能研究未能明确其与肥胖的机制联系。

我们的研究遵循了以下步骤：

确定相关细胞类型：利用表观基因组数据，发现风险单倍型位于一个在间充质干细胞（脂肪细胞的前体细胞）中活跃的“超级增强子”内。
确定靶基因：结合染色质构象捕获数据（显示基因组三维空间相互作用）和eQTL分析，发现该增强子通过长距离相互作用调控远端（相距60万和120万个碱基）的IRX3和IRX5基因，而非其邻近的FTO基因。
确定因果核苷酸：通过进化保守性和转录因子motif分析，锁定了一个将富含AT的序列变为C的特定SNP。该变异破坏了一个抑制性转录因子（ARID5B）的结合位点。
确定上游调控因子：证实ARID5B在脂肪前体细胞中高表达，并与该motif结合，共同抑制增强子活性。风险等位基因破坏了这种结合，导致增强子去抑制。
阐明细胞与生物体表型：风险等位基因导致IRX3和IRX5表达升高，这抑制了脂肪细胞的“米色化”（产热过程），使其更倾向于储存能量的“白色”脂肪细胞表型。在小鼠中，降低Irx3表达可显著抵抗高脂饮食诱导的肥胖。最后，使用CRISPR基因编辑技术将风险等位基因（C）纠正为保护性等位基因（T），成功恢复了脂肪细胞的产热能力。

这个案例完整地展示了如何整合多种基因组学工具，将一个非编码区的遗传关联，解析为具体的细胞类型、靶基因、因果变异、调控因子和生理表型，为理解和干预复杂疾病提供了范例。

🎯 总结

在本节课中，我们一起学习了：

孟德尔遗传与连锁分析：如何利用家族研究和连锁图谱定位高外显率单基因病的致病基因。
复杂性状与GWAS基础：面对多基因、小效应的复杂性状，如何设计并进行全基因组关联研究，包括质量控制、统计检验和多重检验校正。
结果深化：如何通过精细定位区分连锁与因果变异，以及如何通过荟萃分析整合多个研究来发现更多位点。
机制解析：如何利用表观基因组学、motif分析、eQTL和系统生物学工具，将非编码的GWAS信号转化为对疾病生物学机制的理解，包括确定作用的细胞类型、靶基因和因果变异。
完整案例：通过对FTO肥胖位点的深入剖析，展示了从遗传关联到功能机制的全链条研究范式。

通过这门课，我们看到了现代遗传学如何将海量的遗传变异数据与多层次的功能基因组学信息相结合，逐步揭开复杂疾病的神秘面纱，并最终指向新的治疗靶点和干预策略。

课程 P15：L15 - eQTLs 表达数量性状位点 🧬

在本节课中，我们将学习表达数量性状位点的概念、动机、方法论及其在疾病研究中的应用。我们将从回顾上一讲内容开始，逐步深入到eQTLs的生物学意义、发现方法以及相关的补充分析技术。

上一讲我们介绍了群体遗传学与疾病作图的基础，包括全基因组关联研究的挑战、连锁不平衡对精细定位的影响，以及利用表观基因组学注释来解析疾病位点功能的方法。本节中，我们将重点探讨分子层面的数量性状，特别是表达数量性状位点。

概述：为何研究eQTLs？

数千个遗传位点各自对疾病表型产生微小影响。直接研究这些位点与复杂疾病的关系非常困难，因为单个位点的效应量通常极小。然而，这些遗传变异对中间分子表型（如基因表达水平）的影响可能更强、更易检测。通过研究遗传变异如何影响基因表达，我们可以搭建起连接遗传变异与最终疾病表型的桥梁，从而更有效地发现疾病相关基因和通路。

核心概念：eQTL（Expression Quantitative Trait Locus）即表达数量性状位点，指的是一个与基因表达水平这一数量性状相关联的遗传区域。

1. eQTLs的动机与概念框架

理解非编码位点的功能回路及其中间表型至关重要。eQTL研究的目标是促进对疾病易感基因的定位，即找到非编码变异的下游靶基因。

基本逻辑是：一个遗传变异影响某个基因的表达水平，而这种表达水平的改变最终导致疾病表型。例如，在之前解析的FTO肥胖相关位点中，风险基因型导致了IRX3和IRX5基因表达的升高，而FTO基因本身的表达并未改变，这提示IRX3和IRX5才是该位点的真正靶基因。

因果关系方向：从遗传变异到分子表型（如基因表达、DNA甲基化）的箭头是单向的，因为基因型在个体出生时即已确定。然而，分子表型与疾病表型之间的关联可能是双向的，需要借助特定统计方法（如孟德尔随机化）来推断因果。

2. eQTLs的研究案例：阿尔茨海默病中的甲基化QTLs

本节我们通过一个具体案例来了解如何研究分子QTLs。该研究在阿尔茨海默病患者和对照中，检测了大脑多个区域数十万个DNA甲基化位点的水平。

以下是数据处理和发现QTLs的关键步骤：

数据预处理与协变量校正：首先需要消除对数据有较大影响的全局协变量。这些可能包括实验批次、个体年龄、性别、细胞类型比例等。通过主成分分析等方法识别并校正这些协变量，可以凸显遗传变异的局部效应。
不同染色质状态的甲基化特性：增强子区域通常呈现中等程度的甲基化水平，并且是变异度最高的染色质状态。相比之下，启动子区域通常低甲基化且变异度很低。
发现甲基化QTLs：mQTL是指特定CpG位点的甲基化水平与附近基因型之间的关联。研究发现，大多数mQTLs都位于增强子区域，而非启动子区域。
样本量与统计效能：随着样本量增加，能够发现的QTLs数量增多，可检测到的效应量阈值降低，并且QTLs与靶标之间的距离分布也更广。
等位基因频率的影响：样本量较小时，只能发现较高频率的变异位点；随着样本量增大，才能发现更多低频变异位点。
与疾病的关联：研究发现，位于增强子区域的、受遗传控制的甲基化变异，比启动子区域的变异更能预测阿尔茨海默病。这提示增强子可能是遗传变异影响疾病的重要媒介。

3. eQTLs的方法学基础

本节我们来看看发现eQTLs的具体方法。基本概念是将基因表达视为数量性状，寻找附近遗传变异与表达水平之间的关联。

核心分析框架：通常采用线性回归模型。表达式为：
表达水平 = β0 + β1 * 基因型 + β2 * 协变量1 + ... + βn * 协变量n + 误差
我们通过检验回归系数β1是否显著不为零，来判断该遗传位点是否为eQTL。

以下是进行eQTL分析时需要考虑的关键因素：

表达测量与预处理：使用微阵列或RNA测序技术。需进行标准化，并注意探针可能被SNP干扰的问题。
遗传变异搜索空间：通常搜索基因转录起始位点上下游一定范围内的变异（如100 kb或1 Mb）。搜索半径影响多重检验负担和统计效能。
协变量校正：需要校正已知协变量（如年龄、性别、种群结构）以及通过PCA等方法从数据中推断出的未知混杂因素（如批次效应）。
统计效能考量：基因表达的总体方差、次要等位基因频率阈值等因素都会影响发现eQTL的能力。
参数优化：可以通过网格搜索等方法，优化MAF阈值、搜索半径、使用的PC数量等参数，以最大化发现的eQTL基因数量。

4. eQTLs的生物学洞察

随着研究深入和样本量扩大（如GTEx项目），我们对eQTLs有了更多认识：

普遍性：eQTLs广泛存在，随着统计效能提升，几乎每个基因都能在某些组织中发现受遗传调控的证据。
位置：eQTLs通常非常靠近转录起始位点。
组织特异性与共享性：eQTLs具有一定组织特异性，但大脑不同区域之间、或某些功能相关的组织之间也存在大量共享的eQTLs。
反式eQTLs：影响转录因子编码基因或其结合亲和力的变异，可以以反式作用方式影响多个靶基因的表达。
助力GWAS发现：将搜索范围限制在显著的eQTLs区域内，可以大大提高发现疾病相关GWAS位点的能力。

5. 等位基因特异性分析：一种补充方法

传统的eQTL分析比较不同个体间的表达差异。等位基因特异性分析则着眼于杂合子个体内部，比较来自父本和母本染色体的等位基因表达比例是否失衡。

这种方法优势在于：

它是在同一个体、同一细胞环境下进行比较，消除了个体间混杂因素的影响。
可以与传统的eQTL信息结合，进一步提高发现调控变异的统计效能。
有助于研究表观等位基因效应，即附近的表观遗传状态如何影响遗传变异的效应。

6. 中介分析与贝叶斯方法

最后，我们探讨如何利用中介分析来理解遗传变异影响疾病的路径。核心问题是：一个遗传变异影响疾病，是通过哪些中间表型介导的？

方法一：基于基因型预测表型：首先在部分样本中建立遗传变异到中间表型的预测模型，然后利用该模型在全基因组队列中估算每个人的中间表型水平，最后检验估算值与疾病的关联。由于估算值完全由遗传决定，此关联可提示因果方向。

方法二：基于汇总数据的转录组关联研究：这是一种更高效的方法。它利用GWAS和eQTL研究的汇总统计数据，通过线性模型的组合，直接推断基因表达对疾病的潜在因果效应，而无需个体层面的基因型数据。

此外，贝叶斯方法在线性回归框架中引入先验分布，有助于：

选择相关变量，实现稀疏建模。
利用连锁不平衡信息进行精细定位。
通过线性混合模型更好地校正种群结构等复杂效应。

总结

本节课我们一起学习了表达数量性状位点的核心概念。我们从研究eQTLs的动机出发，回顾了其在解析疾病位点中的应用案例。接着，我们深入探讨了发现eQTLs的方法学框架，包括数据处理、回归模型和参数优化。然后，我们了解了eQTLs的普遍性和生物学特性。最后，我们介绍了两种强大的补充分析技术：等位基因特异性分析可以从个体内部提供调控证据，而中介分析与贝叶斯方法则能帮助我们推断遗传变异影响疾病的因果路径和中间媒介。这些工具共同构成了解析复杂性状遗传架构的强大武器库。

课程 P16：L16- 系统遗传学与遗传力 🧬

在本节课中，我们将要学习系统遗传学的核心概念，特别是遗传力。我们将探讨如何定义和计算遗传力，如何将表型变异分解为遗传和环境成分，以及如何利用全基因组关联研究（GWAS）数据来估计和划分遗传力。此外，我们还将介绍多基因风险评分、线性混合模型、LD分数回归等关键技术，并探讨复杂性状的遗传结构，包括多基因性和“全基因”模型。

什么是遗传力？🧬

上一节我们介绍了表型变异与基因型变异的关系。本节中，我们来看看如何量化遗传因素对表型变异的贡献，即遗传力。

我们假设表型（P）是基因型（G）和环境（E）的函数。表型变异（Vp）可以分解为遗传变异（Vg）、环境变异（Ve）以及基因与环境的协方差。为简化模型，我们通常假设基因与环境之间没有交互作用，因此表型变异可以表示为：

Vp ≈ Vg + Ve

遗传力（h²）被定义为遗传变异占表型总变异的比例：

h² = Vg / Vp

这个比例取决于我们测量表型的精确度。如果测量误差很大，那么归因于环境的变异（Ve）就会增加，从而导致估计的遗传力降低。

遗传力的划分：狭义与广义

我们讨论了表型变异的总体分解。现在，我们可以进一步将遗传变异（Vg）本身划分为不同的成分。

遗传变异可以进一步分解为：

加性效应（Va）：多个基因位点效应简单线性相加的贡献。
显性效应（Vd）：等位基因间非加性互作的贡献。
上位效应（Vi）：不同基因位点间互作的贡献。

因此，广义遗传力（H²） 指的是遗传变异（Vg）占表型总变异的比例，它包含了所有遗传效应：

H² = Vg / Vp

而狭义遗传力（h²） 则特指加性遗传效应（Va）占表型总变异的比例：

h² = Va / Vp

在人类遗传学中，关于加性、显性和互作效应的相对重要性一直存在讨论。

为什么研究遗传力？🔍

理解了遗传力的定义和划分后，一个自然的问题是：研究遗传力有什么实际意义？

研究遗传力非常重要，原因如下：

它允许我们量化遗传与环境对特定性状的相对重要性。
它可以揭示性状的遗传结构，例如有多少个因果变异、它们的效应大小和等位基因频率分布如何。
遗传力是表型预测的基础参数，它定义了线性模型理论上最佳的预测性能。
在GWAS研究趋于成熟的今天，遗传力帮助我们判断是否已经找到了大部分因果位点，或者是否还存在大量未被发现的微小效应位点。

从GWAS估计遗传力

我们知道了遗传力的理论定义。那么，如何利用实际的GWAS数据来估计它呢？

在GWAS中，我们估计每个SNP的效应大小（β）。每个SNP所能解释的表型方差取决于其效应大小（β）和次要等位基因频率（MAF, f）：

方差贡献 ≈ 2 * f * (1 - f) * β²

将所有SNP的贡献相加，理论上可以得到加性遗传力（h²）的估计值。然而，一个著名的现象是，仅使用达到全基因组显著性水平的SNP所解释的方差，远低于通过双胞胎或家系研究估计的总遗传力。这被称为“遗传力缺失”问题。

可能的解释包括：

存在成千上万个效应极微小的SNP尚未被GWAS发现。
存在未被GWAS标签SNP捕获的稀有变异。
存在非加性遗传效应（如显性、互作效应）。

目前的证据更倾向于第一种解释，即存在大量效应微小的SNP共同贡献了“缺失”的遗传力。

多基因风险评分

既然GWAS发现了许多与性状相关的SNP，我们能否利用它们来预测个体的表型或疾病风险？这就是多基因风险评分的目的。

多基因风险评分 是通过累加个体在所有相关SNP位点上的风险等位基因数量（加权以该SNP的效应大小）来计算的一个综合分数：

PRS = Σ (β_i * G_i)

其中，β_i 是SNP i的效应大小估计值，G_i 是个体在该SNP位点的基因型剂量（0, 1, 2）。

实际操作中面临挑战：

连锁不平衡：SNP之间并非独立，直接求和会导致重复计算。
效应估计误差：对于未达到显著性的SNP，其β估计值噪音很大。

解决方法包括：

SNP剪除：在高度相关的SNP中只保留一个。
LD调整：使用主成分分析等方法将SNP空间转换为正交的“特征基因座”空间，以获取独立的效应成分。

研究表明，即使纳入远低于全基因组显著性阈值的SNP，PRS的预测能力仍会持续提升，这进一步支持了复杂性状受大量微小效应SNP影响的模型。

疾病遗传结构与遗传力划分

多基因风险评分表明许多SNP都有贡献。我们能否更精细地探究遗传力在基因组中的分布？这就需要对遗传力进行划分。

我们可以不计算整个基因组的遗传相关性，而是只针对特定的SNP子集（例如特定染色体、特定功能注释区域）来计算。这被称为遗传力划分。

例如，按染色体划分遗传力时发现，较长的染色体倾向于解释更多的遗传力，这暗示因果变异在基因组中几乎是均匀分布的，而非集中在少数热点。

更精细的划分可以基于功能注释（如编码区、启动子、增强子）。通过分层LD分数回归等方法，我们可以评估不同功能类别对遗传力的贡献。研究发现：

与功能相关的区域（如DNA可及性区域、增强子）的遗传力富集程度远高于其基因组占比。
与疾病相关的细胞类型中活跃的增强子，其遗传力富集尤为明显。

这种方法仅需GWAS汇总统计数据和一个参考群体的LD矩阵，无需个体基因型数据。

线性混合模型

为了更准确地估计遗传力并校正样本间的复杂关系，我们需要引入更强大的统计模型：线性混合模型。

标准的线性模型假设样本间独立同分布。然而，在遗传学研究中，样本间可能存在亲缘关系或群体分层，导致假关联。线性混合模型 通过引入一个随机效应来捕获这些由样本间遗传相似性引起的协方差。

模型可以表示为：

y = Xβ + u + ε

其中：

y 是表型向量。
X 是基因型矩阵。
β 是固定效应（SNP效应）。
u 是随机效应，服从分布 N(0, σ_g² * K)，K是亲缘关系矩阵，衡量个体间的遗传相似性。
ε 是残差。

利用限制性最大似然估计等方法，我们可以直接估计加性遗传方差σ_g²，而无需估计每个SNP的具体效应β。LMMs能够解释比仅使用显著SNP更多的遗传力，是遗传力估计和校正混淆因素的强大工具。

多基因性与全基因模型

通过划分遗传力，我们发现许多功能类别在很宽松的P值阈值下仍然富集。这引出了关于性状多基因性程度的深刻问题。

对多种疾病的研究表明，与疾病生物学相关的功能注释（如免疫疾病中的免疫增强子，精神疾病中的神经增强子）的富集信号会延伸到成千上万个SNP位点，形成很长的“尾巴”。这意味着有数千个独立的功能性位点集中在相同的通路上。

更令人惊讶的是，几乎所有基因本体论类别在多种疾病的关联基因中都有一定程度的富集。这促使了“全基因模型”的提出。该模型认为：

存在一个“核心基因”集合，直接与疾病病理相关。
存在一个更庞大的“外围基因”集合，它们通过调控网络与核心基因相互作用，间接影响性状。
虽然核心基因的富集程度最高，但大部分遗传力实际上由外围基因贡献。GWAS早期发现的大多是核心基因，而更大规模的GWAS则不断揭示出外围基因的贡献。

这个模型解释了为什么遗传力如此广泛地分布在基因组中，以及为什么看似不相关的功能类别也会与疾病关联。

总结 📚

本节课中我们一起学习了系统遗传学的核心——遗传力。我们从定义和计算遗传力开始，探讨了如何将表型变异分解为遗传和环境成分，以及狭义与广义遗传力的区别。我们深入分析了如何从GWAS数据估计遗传力，并探讨了“遗传力缺失”问题。

接着，我们介绍了利用GWAS结果进行表型预测的工具——多基因风险评分，以及用于估计和划分遗传力的关键统计模型——线性混合模型和LD分数回归。这些技术帮助我们理解复杂性状的遗传结构，揭示了遗传力在基因组不同功能区域间的分布。

最后，我们探讨了复杂性状的多基因性本质，并介绍了最新的“全基因模型”。该模型认为，除了直接相关的核心基因外，大量通过调控网络间接作用的外围基因共同贡献了性状的大部分遗传力。这为我们从系统层面理解遗传对复杂性状的影响提供了全新的视角。

🧬 课程 P17：L17 - 比较基因组学

在本节课中，我们将要学习比较基因组学，特别是用于基因组注释的进化特征。这是本课程关于比较基因组学和进化的最后一个模块。今天的讲座将聚焦于比较基因组学和进化特征。

下一次课，我们将讨论基因组尺度的进化和基因组复制，然后是系统发育学和系统发育基因组学，接着是测验，最后我们将转向未来方向。

之前我们讨论了很多关于基因组注释、基因表达、表观基因组学、网络和调控基因组学，以及疾病遗传学和疾病基因组学的内容，这些主要关注人类谱系内部的变异。今天，我们将讨论跨物种的变异，特别是研究如何利用全基因组研究来推断进化特征，然后利用这些特征来注释元件。我们将把这些特征应用于注释蛋白质编码基因、非编码基因、microRNA基因和调控基序。

下一次课，我们将讨论基因组重排和基因组复制，如何基于同线性构建全局比对，以及如何检测结构变化。然后，我们将有两节关于系统发育学和系统发育基因组学的讲座，即如何利用系统发育来理解进化速率和进化模型，回到树上的动态规划和贝叶斯进化模型。系统发育基因组学与今天的讲座非常相似，都是利用全局特征来理解局部的个体事件，我们将研究基因树、物种树和协调。

概述：利用比对揭示功能元件

之前我们展示过这张幻灯片，用以激发构建基因组比对的动机。本次讲座基本上从那里开始，即我们可以利用近缘物种的比对来研究不同区域的保守模式。

这里，星号表示这四个物种间的完美保守。基于连续的星号，可以看到保守的“岛屿”。事实上，这些保守岛屿与这些区域内已知的调控基序实例非常吻合。因此，核心思想是：我们现在能否系统地使用比较基因组学来研究这些特征并发现这些元件？

具体来说，今天的讲座将侧重于“解读进化以揭示功能元件”。下一次课，我们将研究如何利用基因组来理解进化本身。基本上，我们越理解进化，就越理解基因组；反之亦然。

本节课目标

首先，观察核苷酸保守性和进化约束。
其次，关注进化特征，重点是变化的模式，而非保守的数量。
然后，研究蛋白质编码基因的特征。
接着，研究RNA、基序、转录因子结合位点和microRNA的特征。
最后，了解如何测量人类谱系内部的选择。

比较基因组学的力量与核心概念

比较基因组学极其强大。我们可以通过比较相关物种来发现功能元件。其核心概念自达尔文以来就已理解：一方面是随机突变，另一方面是自然选择。

随机突变通过随机突变元件来探索空间，这并非有意为之，只是犯错误。如果由工程师负责选择，我们可能仍然是完美复制的细菌，没有错误。因此，突变是一把双刃剑：一方面它引入问题，另一方面它引入新事物。

我们今天主要关注的是纯化选择如何运作，即那些破坏元件的突变是如何被维持的。在系统发育学讲座中，我们会更多地讨论正选择。但现在，我们主要关注随机突变，它大多会破坏功能，但偶尔也会使功能变得更好。

自然选择则是残酷的决策，在环境背景下选择那些维持功能的序列。因此，一方面是完全盲目的随机突变，另一方面是维持那些仍保留功能的序列的自然选择。

这两种力量在进化时间尺度上作用的结果是：一方面，非功能区域会积累突变并被保留；另一方面，功能区域也会积累突变，但这些突变会降低适应性，因此随着时间的推移，那些具有破坏性元件的生物体适应性会降低，其基因最终会减少。

尽管功能和非功能区域发生相同数量的突变，但随着时间的推移，功能区域的突变会从基因库中被排除，而非功能区域的突变则会积累。

检测能力与物种选择策略

我们关心的一个问题是检测能力：我们有多大能力来发现这些功能元件？一个非常简单的概念是：系统发育树上的分支长度越长，非功能区域相对于功能区域积累的事件就越多，你区分非功能区域突变积累和功能区域突变排除的能力就越强。

因此，你的检测能力应大致与所比对物种树中捕获的进化距离成线性比例。如果我们的目标是根据突变数量区分功能与非功能区域，那么在非常近的距离下，两个区域都没有什么突变；在足够的距离下，我们可以区分功能与非功能区域之间替换的差异积累；而在非常遥远的比较中，功能区域不再保守，因此你将无法真正捕捉到功能与非功能之间的区别。

比较基因组学中常用的解决方案是比较许多近缘物种，这比比较少数远缘物种要强大得多。也就是说，对于相同的总分支长度，我们更倾向于许多近缘物种，因为功能区域在每对物种间都是保守的，任何一对物种间的进化比较都没有足够的时间积累足够的新环境适应。在较短的进化距离内，环境变化不大，因此成对保守性仍然很高，但非功能区域会独立地积累“噪音”。

我喜欢用这个类比：当你用一个麦克风录制交响乐时，如果麦克风周围有噪音，你录制的噪音和信号一样多。如果有多个麦克风，信号会线性传输到所有麦克风，但每个麦克风的噪音是独立的，因此我们能够真正捕捉到进化信号，因为不同谱系中噪音的独立积累。

检测约束元件：从全局比对着手

现在我们已经研究了中性分支能力和纯化选择如何通过排除那些积累了破坏性突变的生物体来维持功能元件不变。接下来，让我们看看如何检测约束元件，特别是观察可以检测这些核苷酸的单个核苷酸或模型。

现在，你可以构建这些全基因组比对。对于任何一个物种，你可以看到不同颜色的基因。然后，在进化时间尺度上，你可以比对这些物种，并看到顺序大致是保守的。

这很重要，因为通过选择染色体片段顺序大致保守的物种，意味着我们不仅可以构建蛋白质编码区域的比对，还可以构建这些基因间区域的比对。这意味着我们实际上可以研究整个基因组中编码和非编码区域的进化模式。

这可以同时应用于数百个基因。我们现在可以开始使用这些比较来识别其中的功能元件。

在第一讲或第二讲中，我们也展示过这张幻灯片，通过显示例如外显子（用浅蓝色高亮显示的蛋白质编码外显子）实际上是非常深度保守的。当然，挑战在于许多其他元件也高度保守，问题是这些到底是外显子还是调控区域。

如果我们试图理解这个特定区域是编码还是非编码，我们可以观察该序列的变化模式。我们可以看到，它在哺乳动物之外就不再保守了，在鸡、河豚或斑马鱼中不保守。这表明这可能实际上是一个非编码区域，因为蛋白质编码外显子往往在哺乳动物树之外也非常深度保守。

一般来说，蛋白质编码基因进化得更慢。现在我们可以首先开发检测约束的方法，然后开发区分不同类型约束的方法。

检测约束的方法

对于检测约束，我们讨论过如何计算编辑操作的总数或替换和空位的数量，这是我们动态规划讲座的内容。或者，我们可以估计突变的数量，包括对回复突变的估计，这将是我们系统发育学讲座的重点。

我们还可以结合邻域信息，寻找保守窗口。例如，当我们观察不同的核苷酸片段时，我们可能会说，这里有一个保守的核苷酸，但它是单独保守的；那里有另一个保守的核苷酸，但也是单独保守的。也许连续多个核苷酸保守的事实实际上很重要。

因此，我们可以使用保守窗口的概念。你可以编写一个隐马尔可夫模型来寻找高度保守区域和低度保守区域之间的最佳转换。我们在HMM讲座中讨论过隐马尔可夫模型在理解人类基因组中蛋白质编码约束岛屿或非编码约束岛屿位置的应用。

具体来说，这里Y轴上的信息显示了处于与保守与非保守元件对应的隐马尔可夫模型状态的后验概率。同样，你可以估计约束隐藏状态的后验概率，可以使用维特比解码来找到与保守/非保守最佳部分对应的路径，或者使用后验解码来定义该区域每个核苷酸最可能的隐藏状态。

同样，你可以使用系统发育来估计树突变率或拒绝的替换。这基本上是询问：沿着我的树，是否存在我预期会发生但并未发生的单个SNP突变？我们也可以允许不同的速率，这将在系统发育学部分讨论。

但今天我要关注的主要概念是选择模式与选择速率的区别。

选择模式 vs. 选择速率

这里有三个区域进化中性或非中性的例子。在中性区域，你可以看到有很多替换。每个保守位置显示为一个点，每个变化显示为该序列已变为的核苷酸。

因此，一个中性序列基本上有多个独立的变化，在每个核苷酸之间切换。这是一个中性序列。约束序列可以显示变化速率的降低，这基本上是说：我可以发生任何类型的突变，但总体突变更少。所以这里有更多的点，有更多约束序列。

因此，我们可以建立一个替换率的概率模型，对该进化速率ω进行最大似然估计，并报告速率ω，以及该速率非中性的对数优势得分。然后，我们可以使用基于窗口的或位点方式的应用ω。

到目前为止，大家都理解了吗？

分支是如何发生的？

分支基本上是物种形成。例如，当一个物种的祖先染色体发生重排，使得它们的染色体不再对齐，即使果蝇可能有后代，这些后代也可能无法存活，因为它们的染色体无法正确对齐。

另一个可能性是地理隔离，例如夏威夷的果蝇物种，熔岩流可以分隔物种，使它们无法再相互交配，从而导致物种形成事件。

所以，分支的发生要么是由于地理隔离，要么是由于合子后不兼容。

约束强度与特定模式

我们讨论了约束的强度，但我们也可以研究这个区域。在这个区域，我们发现有很多独立事件，其中C变成了G，又在不同的独立谱系中变成了G。这告诉我们什么？这是否意味着这个核苷酸在随机进化？是中性进化吗？不，不完全是中性。它被约束为要么是C，要么是G。一半的物种是C，一半是物种是G。

因此，与其研究一个速率（这里速率会很高，表明该核苷酸在中性进化），我们可以检测不寻常的替换模式。我们可以建立一个四个字符A、C、G、T之间马尔可夫链的平稳分布的概率模型。

然后，随着字母在该马尔可夫链中切换，我们可以问：它们是否以非随机的方式相互切换？我们可以构建该区域被约束的核苷酸向量的最大似然估计量，用于基因组中的每个k-mer，以及该k-mer非中性的对数优势得分。

到目前为止，大家都理解整体变化量与特定模式之间的区别了吗？

应用：检测特定约束与基序

利用这一点，我们现在可以遍历基因组，不仅判断什么是高度保守或低度保守，还可以判断什么是保守为T或G，保守为A或C或T或G，等等。

因此，我们基本上可以匹配全基因组的位置权重矩阵。我们已经用它来揭示单个转录因子结合位点，以及在基序实例内的位置特异性偏好。有了更多物种，你实际上可以直接从序列推导出基序共识，你可以看到多个核苷酸实际上在特定位置受到约束，并可以用它来识别破坏进化保守基序的遗传变异。

我们在讨论GWAS功能解释时谈到过这一点。

估计基因组中受约束的比例

这是关于通过序列特异性（位置特异性突变量）或位置特异性变化模式来检测约束元件。现在，我们也可以用这个来估计基因组中受约束的比例。

也就是说，即使我们无法检测到单个核苷酸，我们是否可以检测到整个基因组似乎比偶然预期更受约束的事实？这里的想法是，我们将观察进化约束的分布。

我们可以通过ω（速率）或π（模式）来测量约束。然后，基于这两条曲线之间的区别，识别实际上有多少保守是可检测的。这是关于在这里设置一个阈值，然后询问高于该阈值的是什么。

在每个阈值处，我们都有一些假阳性率。因此，我们可以基于中性元件选择这个假阳性率。我们可以问：我预期中性进化区域有多少进化？这就是这里的蓝色部分。

因此，我们基本上可以观察保守分数在红色（实际）与蓝色（预期背景，若无约束）中的分布。然后，在任何截止值处，我们基本上都有一组真阳性（红色分布中高于阈值的所有部分）和假阳性预测（蓝色分布中高于该阈值的部分）。

在每个阈值处，我们基本上都有某个假阳性率。假阳性率基本上告诉你，我可以基于蓝色分布选择阈值，以获得5%的假阳性率之类。基于这两个分布之间的差异，我可以选择一个阈值和一个假阳性率。

但问题是我们无法检测到所有的约束元件，因为曲线重叠。基本上，真实的信号相对于蓝色的中性信号略有偏移，但分布仍然高度重叠。因此，我们将永远无法在不引入一堆假阳性的情况下检测到这里的东西。

然而，我们实际上可以做的是，通过对两条曲线之间的整个区域进行积分来估计约束的总过量。利用这一点，我们基本上可以说，基因组中有多少比例不仅可检测为受约束，而且总体上估计有多少比例的基因组受约束，尽管我根本没有能力检测到那些。

在哺乳动物比对的背景下，这种分析表明，蛋白质编码转录本内部，以及5‘ UTR、3’ UTR和启动子中，约束有非常强的富集。在其他基因组区域也有，但新的元件主要落在内含子和基因间区域，即富集较少，但那里有大量的约束。

利用这一点，我们基本上可以通过比较蓝色和红色曲线来估计存在这种过量的约束。但如果你放大这条曲线，比较全基因组（蓝色）和祖先元件（红色），这些是祖先重复元件，我们预计它们不受很强的约束。因此，我们可以使用这些红色元件来基本估计过量约束。

然后，这种过量约束实际上有两种模式：一种是过量纯化约束，另一种实际上是过量正选择，那里基本上有新元件被发现、进化出来。

到目前为止，大家都理解这个不仅使用特定阈值检测元件，还观察过量元件的整体概念了吗？

覆盖深度作为约束的额外度量

我们在这里讨论的约束是基于对不同元件约束的实际评估，但你也可以观察不同区域的整体覆盖深度。你可以看到，例如，蛋白质编码外显子在29个物种中约有20个物种的覆盖，而对于祖先重复序列，覆盖要少得多。

这意味着什么？这意味着不仅比对物种在那里高度保守，而且在基因组的功能区域有更多比对物种。因此，你可以使用整体覆盖深度作为约束的额外度量。

如果你观察在四种哺乳动物中保守的元件，你可以看到它们在29种哺乳动物中对齐得更多。然后，我们发现的这些新元件也对齐度很高，且高度保守。

目前，我们只使用了约束信号，但在比对的存在或不存在中，还有额外的残留信号。你可以看到，从4个物种（人、小鼠、大鼠、狗）到29个物种，检测能力如何增加。我们基本上可以看到能检测到更多的约束，但也可以估计更多的约束。

这里在50个核苷酸分辨率下的约束大约是5%。通过变化模式，我们发现了更多的约束估计。

进化特征：关注变化模式

以上都是关于在核苷酸水平上测量约束。但现在，当我们观察更大的元件时，我们实际上可以寻找涉及多个核苷酸的模式。现在，我们将讨论进化特征的概念，然后特别关注变化模式。

这里的想法是，我们可以开发特征来检测，也许这个区域实际上与那些其他区域的进化方式不同。这里的核心概念是，一个区域的具体功能——无论是蛋白质编码基因、非编码RNA还是调控基序——决定了作用于该区域的选择压力，而这些选择压力实际上决定了我们将要发现的突变、插入和缺失的模式。

因此，我们可以开发对每种功能类型具有特征性的进化特征。例如，如果你观察蛋白质编码基因，并比较蛋白质编码与非编码区域，密码子替换频率本身告诉你，核苷酸三联体以与氨基酸功能选择一致的速率进行交换，而在非编码区域，三联体以不同的模式交换。

蛋白质编码进化的一个特征是密码子替换频率。第二个特征是阅读框保守，即保持翻译阅读框的压力，因此空位总是3的倍数（例如6和3），而不是像这里的14、7和17。这些是与蛋白质编码功能相关的进化特征。

我们也可以讨论与RNA结构相关的进化特征。这里的功能不是在氨基酸水平，而是在碱基配对水平。因此，那些改变核苷酸序列但保留配对的变化将被容忍。所以，如果你有一个将GC碱基对变为AT碱基对的替换，这完全没问题，折叠仍然发生。

因此，保留折叠的补偿性变化实际上是RNA结构的特征。还有沉默的GU替换，其中G在RNA结构中既可以与C配对，也可以与U配对，这些GU/GC替换是另一个进化特征。

注意，我们不是在问有多少保守，而是在问变化是否专门保留了特定功能选择的特征。因此，对于蛋白质编码区域，白色核苷酸是没有变化的地方。我们甚至不使用白色核苷酸，只使用彩色核苷酸（即发生变化的地方）来寻找这些进化特征。

这对于检测特定类别的功能元件非常具有特征性。对于调控基序，特征是突变实际上保留了共识，并允许在任何位置增加分支长度得分，以及全基因组保守。这些特征共同使我们能够从进化特征中检测调控基序，并检测单个基序实例。

我展示这张幻灯片的原因，即使我们将深入探讨每一个细节，是希望你们认识到，对于每种类型的功能元件，存在非常不同类别的进化特征。其背后的关键思想是，突变将随机探索这个进化空间，而这个空间的形状对于蛋白质编码基因、RNA、microRNA等将非常不同。

我们正在观察编码相同功能的所有序列的空间，而这个空间使我们能够检测这些进化特征。在那个空间中，你基本上必须考虑更多分支长度的影响。如果你试图推断形状是正方形还是圆形，在纸上只扔了几个点，你可能无法分辨。但随着更多事件的发生，可容忍突变的形状实际上变得更加清晰。

这就是进化特征的概念：不同类别的元件以不同的模式、不同的方式进化。我们不会关注变化量或保守量，而是关注特定的变化模式。

到目前为止，大家都理解进化特征的概念以及它们与单纯的整体变化量有何不同了吗？

应用进化特征

利用这一点，我们现在可以深入研究蛋白质编码基因，发现新的蛋白质编码基因或修订现有的蛋白质编码基因，或检测RNA内部不寻常的基因结构。我们可以开发方法来检测新的结构家族或识别靶向、编辑或稳定性信号的方式，以及核糖开关。

在microRNA进化特征中，我们可以发现新的microRNA家族，扩展现有的microRNA家族，寻找两个臂协同作用或正义链和反义链实际上都有功能的特定类别的microRNA。我们也可以在基序领域发现新的调控基序、新的调控基序实例，或我们在调控基序和调控基因组学讲座中讨论过的转录因子和microRNA网络，并且我们还可以推断其中许多的单结合位点分辨率。

深入探讨蛋白质编码进化的特征

现在让我们深入探讨蛋白质编码进化的特征。我们将讨论阅读框保守和密码子替换频率。

这里有一个测试：这是一个12种果蝇物种的比对，它们跨越的距离与29种哺乳动物大约6000万年的分歧相似。我们现在可以问：这里的蛋白质编码区域在哪里？这个区域的一部分编码蛋白质，一部分是非编码。你们可以猜猜哪个区域是蛋白质编码，哪个是非编码。

如果我们在这里画一条垂直线，我会声称这条线左侧的序列与右侧的序列进化方式不同。我们讨论过的一些进化特征是什么？

第一个特征是空位。如果你看这里的空位，长度是5，长度是1，这不是你在蛋白质编码区域预期的那种空位。相比之下，这是一个长度为6的空位，

🧬 课程 P18：L18- 基因组进化

在本节课中，我们将要学习基因组进化。我们将从如何组装和比对基因组开始，探讨基因组快速进化的机制，特别是全基因组复制事件，最后研究近期人类进化中的自然选择。

🧩 基因组组装

上一节我们介绍了基因组进化研究的背景，本节中我们来看看如何从测序片段中组装出完整的基因组。

测序技术每次只能读取几百个核苷酸。因此，我们无法直接从一端测序到另一端。常见的策略是使用双端测序。这种方法将DNA打断成固定长度的片段，然后对片段的两端进行测序。通过已知的片段长度，我们可以推断两端序列在基因组中的大致距离。

以下是基因组组装的传统方法步骤：

寻找重叠序列：通过比对或哈希等方法，找到不同读段之间重叠的部分。
构建重叠群：将重叠的读段逐步连接起来，形成更长的连续序列，称为重叠群。
搭建支架：利用双端测序提供的距离和方向信息，将重叠群排列并连接成更大的支架。支架之间的未知序列用字母“N”填充。
生成一致序列：对于每个位置，根据覆盖该位置的所有读段的碱基和质量分数，通过加权投票等方式，确定最终的共识碱基。

组装过程中的主要挑战是基因组中的重复序列。如果重复序列的长度超过读段长度且序列高度相似，组装算法可能错误地将来自基因组不同位置的重复序列“折叠”到一起。

为了处理重复序列，我们可以构建一个重叠图。图中的节点是读段，边表示读段之间的重叠关系。通过分析图的拓扑结构和测序深度（即某个区域被测序覆盖的次数），我们可以推断并“展开”那些被错误折叠的重复区域。

另一种现代组装方法是字符串图法。它直接基于读段构建图结构，并通过寻找图中的最优路径来推导基因组序列。这种方法可以更有效地利用长读段数据来解决重复序列问题。

🔗 全基因组比对

现在我们已经有了组装好的基因组，接下来看看如何将不同物种的基因组进行比对。

对于两个序列的全局比对，我们可以使用动态规划算法，但其时间复杂度是 O(n²)。对于多个物种的全基因组比对，如果使用多维动态规划（例如，3个物种就是三维立方体），计算将变得不可行。

因此，我们采用渐进式比对策略。首先，根据物种间的进化关系构建一棵系统发育树。然后，从树的叶子节点开始，两两比对关系最近的序列，产生一个共识序列。接着，将这个共识序列与下一个最近的序列进行比对，如此逐步向上，直到树的根部。这样就将一个复杂的多重比对问题，分解为一系列两两比对问题。

对于很长的基因组序列，直接进行全局动态规划仍然开销巨大。因此，实际流程通常是：

首先寻找高得分的局部匹配。
将这些局部匹配连接成链。
只在链所定义的区域内进行受限的动态规划，从而得到最终的全局比对。

这种比对方法可以检测到序列之间的大规模进化事件，如倒位、易位和复制。

除了基于序列相似性的方法，我们还可以进行基于基因的全基因组比对。其核心思想是利用基因作为“锚点”，先确定不同基因组间基因的对应关系，再填充基因间的序列比对。

BOS算法（最佳明确子群算法）是解决基因对应关系的一种方法。它构建一个加权二分图，节点是不同基因组中的基因，边的权重是蛋白质序列相似性。算法通过迭代筛选，优先保留那些位于同线区块（基因顺序保守的区域）且相似性高的边，从而解析出一对一、一对多或多对多的同源关系。

通过比对，我们可以发现基因组中快速进化的区域。例如，在酵母基因组中，染色体末端（端粒附近）存在大量基因家族的扩增，这些区域结构不稳定，但有利于快速适应环境变化。

⚡ 快速进化机制与全基因组复制

通过比较亲缘关系较近的物种，我们可以观察到基因组快速进化的具体机制。

研究发现，许多基因组重排（如易位、倒位）的断点附近常常有转座元件或tRNA基因。这些重复元件可能作为“进化热点”，促进了基因组结构的改变。

此外，还存在一些“进化捷径”，例如：

蛋白质内含子：一种可以自我剪切并移动的蛋白质片段，能插入到其他蛋白质中。
通读终止密码子：在特定条件下，核糖体可以越过终止密码子继续翻译。

接下来，我们探讨一个重大的进化跳跃事件：全基因组复制。

通过对酿酒酵母和其近缘物种克莱沃毕赤酵母的基因组比较，研究者发现了一个关键模式。克莱沃毕赤酵母的许多基因组区域，会交错地映射到酿酒酵母的两条不同染色体上。这种基因的“交织”模式，是远古发生的一次全基因组复制的有力证据。

在复制发生后，祖先基因组中的每个基因都变成了两个拷贝。随后，绝大多数冗余的基因拷贝被丢失，只有约8%的基因保留了两个拷贝。通过分析基因顺序，我们可以“描绘”出复制后基因丢失的历史，从而重建这一重大进化事件。

全基因组复制为新功能的产生提供了原材料。对酵母的研究发现，在保留的两个拷贝基因中，约有20%的基因对表现出不对称进化速率，其中一个拷贝进化加速，这可能是新功能化的迹象。这些加速进化的拷贝往往在特定条件（如胁迫）下表达，而进化速率较慢的“祖先功能”拷贝则通常执行必需的基础功能。

👨‍👩‍👧‍👦 近期人类进化与自然选择

前面我们看了跨越百万年的进化，现在把目光聚焦到近期，看看如何在人类群体中检测自然选择的信号。

首先，我们需要一个基线。在理想条件下（无突变、无迁移、无选择、大群体、随机交配），群体中的基因型频率会达到哈迪-温伯格平衡。设等位基因A的频率为 p，等位基因a的频率为 q (且 p + q = 1)，则基因型频率为：

AA: p²
Aa: 2pq
aa: q²

当实际观测值偏离这个预期时，就可能存在进化力量（如选择、遗传漂变）的作用。

我们有多种统计方法，可以在不同时间尺度上检测自然选择：

物种间比较 (长尺度)：使用 dN/dS 检验。比较非同义替换率(dN)和同义替换率(dS)。如果 dN/dS > 1，表明氨基酸改变受到正向选择。
群体内多态性与物种间分化比较 (中尺度)：如 McDonald-Kreitman 检验。比较物种内多态性位点和物种间固定位点中同义与非同义替换的比例。如果固定位点中非同义替换的比例异常高，提示正向选择。
高频率衍生等位基因 (较短尺度)：通过与黑猩猩基因组比较，确定每个位点的“祖先型”和“衍生型”碱基。如果在人类群体中，某个衍生等位基因的频率异常高，可能是近期正向选择的结果。
群体间分化 (群体尺度)：计算 群体固定指数 Fst。Fst 衡量群体间的遗传分化程度。某个位点在群体间Fst值极高，表明该位点可能受到局域适应选择。
长单倍型 (近期选择)：这是检测近期强正向选择的有力工具。当一个有益突变在群体中快速上升时，它周围紧密连锁的染色体区域还来不及被重组事件打断，因此会形成一个长而高频的单倍型。

一个经典案例是欧洲人群中的乳糖耐受性。负责消化乳糖的LCT基因上游的一个突变，使成人在摄入乳制品后仍能表达乳糖酶。这个突变在畜牧业发展后提供了巨大优势，因此受到了强烈的正向选择。在现代欧洲人群中，携带此突变的染色体单倍型不仅频率高（~77%），而且长度非常长（超过1Mb），这正是近期强选择的典型特征。

通过整合这些检测方法，科学家们已经在人类基因组中发现了许多近期受到选择的基因，涉及消化代谢、疾病抵抗、肤色、身高适应等多个方面。

📚 总结

本节课中我们一起学习了基因组进化的多个层面：

我们了解了如何从短测序读段组装出完整基因组，并处理重复序列的挑战。
我们探讨了如何进行全基因组比对，包括渐进式比对和基于基因的比对，以识别保守区域和进化事件。
我们深入研究了快速进化机制，如转座元件介导的重排，并重点剖析了全基因组复制这一重大进化事件及其后续的功能分化。
最后，我们将尺度拉近到人类群体内部，学习了多种检测自然选择信号的统计学方法，并看到了这些方法如何揭示人类近期适应性的进化历史。

从基因组组装到跨物种比对，再到群体遗传分析，这些工具共同帮助我们利用基因组数据来解读生命演化的壮丽篇章。

课程 P19：L19- 系统发育学 🧬

在本节课中，我们将学习分子进化与系统发育学的基础知识。我们将探讨如何从序列数据推断进化历史，包括构建距离矩阵、使用不同算法构建系统发育树，以及如何直接基于序列比对来评估和搜索最优树形结构。

1. 系统发育学基础 📚

上一节我们介绍了课程的整体目标，本节中我们来看看系统发育学的基本概念和定义。

系统发育学研究的是物种、基因或其他实体之间的进化关系。我们利用可观察的特征（性状）来推断它们的历史谱系关系。

性状：可以是形态特征（如牙齿形状），也可以是分子特征（如DNA序列中的特定核苷酸）。
节点：代表分类单元（如物种、基因）。
- 根节点/祖先节点：代表所有分类单元的共同祖先。
- 内部节点：代表假想的祖先，表示分化事件。
- 叶节点/终端节点：代表我们实际观察到的现存分类单元。
分支/谱系：连接节点的线段，代表进化路径和时间。

系统发育树主要有三种类型：

分支图：只关注拓扑结构（谁和谁关系更近），不包含分支长度信息。
时序图：假设所有谱系从共同祖先开始以相同速率进化，分支长度代表实际时间。
系统发育图：结合了拓扑结构和分支长度，分支长度可以代表进化速率或遗传距离。

我们的核心挑战是：如何从一组已比对的核苷酸或蛋白质序列中，推断出最能代表其进化历史的树形结构。

2. 从序列比对到距离矩阵 📏

上一节我们定义了系统发育树的基本要素，本节中我们来看看如何量化序列之间的进化差异，即构建距离矩阵。

我们首先需要测量两条已比对序列之间的进化距离。最简单的方法是计算总体核苷酸差异百分比。例如，如果两条序列90%相同，则距离为10%。

然而，简单的计数会低估实际的突变数量，因为存在回复突变的可能性：一个位点可能从A变为C，后来又变回A，这样在最终比对中我们就观察不到变化。因此，我们需要进化模型来校正这种效应。

以下是校正模型的核心思想：我们将核苷酸的替换过程建模为一个马尔可夫链。在极短的时间 Δt 内，一个核苷酸以特定概率变为其他核苷酸。

Jukes-Cantor模型是一个单参数模型，它假设所有核苷酸之间的相互替换速率相同。其替换概率矩阵可以表示为：

给定时间 t 后，一个位点保持不变的概率 P(相同) 和发生改变的概率 P(不同) 的公式为：

P(相同) = 1/4 + (3/4) * exp(-4αt)
P(不同) = 3/4 - (3/4) * exp(-4αt)

其中 α 是替换速率参数。

通过观测到的差异比例 p，我们可以反推校正后的遗传距离 d：

d = -(3/4) * ln(1 - (4/3)p)

Kimura双参数模型则更进一步，它考虑了转换（嘌呤间或嘧啶间的替换，如A↔G）和颠换（嘌呤与嘧啶间的替换，如A↔C）通常具有不同速率的事实，因此使用两个参数进行建模。

实际上，存在一整套模型层次结构，从简单到复杂（如考虑不同碱基频率的GTR模型），用于更精确地估计序列间的进化距离。

3. 从距离矩阵到系统发育树 🌳

上一节我们学会了如何计算序列间的校正距离，本节中我们来看看如何利用这些距离来构建系统发育树。

距离矩阵包含了所有序列对之间的差异信息。构建树的目标是找到一个树形结构，使得树上任意两叶节点之间的路径距离（即沿途分支长度之和）尽可能接近距离矩阵中对应的观测值。

首先，我们需要了解距离矩阵的类型：

超度量距离：满足三点条件，d(i,k) = d(j,k) > d(i,j)。这对应于一个所有叶节点到根节点距离相等的有根树（即存在分子钟）。
可加距离：满足四点条件，d(i,j) + d(k,l) ≤ d(i,k) + d(j,l) = d(i,l) + d(j,k)。这保证存在一个树能精确反映这些距离。
一般距离：由于噪声、模型不准确或进化速率不均等原因，实际数据往往既非超度量也非严格可加。

以下是几种基于距离的建树算法：

UPGMA（非加权组平均法）
这是一种层次聚类方法。

找到距离最近的两个类群（初始为单个序列），合并它们，形成一个新节点。
新节点的高度设为两个类群距离的一半。
重新计算新类群与其他所有类群的平均距离。
重复上述步骤，直到所有类群合并完毕。

UPGMA假设分子钟存在，如果数据是超度量的，它能构建出正确的树；如果进化速率不均，它可能会得出错误拓扑。

邻接法
邻接法通过引入一个“净分歧度”来修正UPGMA的缺陷，该度量反映了每个节点到所有其他节点的总距离。
其修正后的距离 D(i,j) 计算公式为：

D(i,j) = d(i,j) - (r_i + r_j) / (N - 2)

其中 r_i 是序列 i 到所有其他序列的距离之和，N 是序列总数。

邻接法选择使 D(i,j) 最小的两个节点进行合并。它对于可加距离矩阵能保证得到正确树，对于一般距离矩阵也通常表现良好。

除了这些循序渐进的算法，我们还可以使用最优性标准来搜索最佳树，例如最小二乘法（最小化观测距离与树距之差的平方和）或最小进化法（最小化树的总分支长度）。但这些方法需要结合树空间搜索策略。

4. 从序列比对直接评估系统发育树 ⚖️

上一节我们介绍了基于距离的建树方法，本节中我们绕过距离矩阵，直接学习如何根据序列比对数据来评估一个给定树形结构的优劣。

这类方法不直接构建树，而是为不同的树拓扑结构进行评分，需要与树搜索算法（下一节介绍）结合来寻找最优树。

最大简约法
核心思想是选择所需进化事件（如核苷酸替换）最少的树。对于比对中的每一列（位点），我们尝试为所有内部节点分配核苷酸状态，使得整棵树上该位点所需的状态变化总数最小。

可以使用动态规划算法高效计算：

从叶节点开始，每个叶节点的状态集合已知（即观测到的核苷酸）。
向根节点方向计算。对于一个内部节点，查看其两个子节点的状态集合。
- 如果两个集合有交集，则父节点状态设为该交集中的一个，不增加代价。
- 如果无交集，则父节点状态设为两个集合的并集，代价增加1。
遍历所有位点，将代价相加，总代价最小的树被认为是最优的。

最大似然法与最大后验概率法
这些是概率模型方法，不仅考虑拓扑，还考虑分支长度 t。

最大似然：寻找使观测数据 D（序列比对）出现概率最大的树拓扑 T 和分支长度 t。
```
(T, t)_ML = argmax P(D | T, t)
```
最大后验概率：在最大似然的基础上，结合了关于树和分支长度的先验知识 P(T, t)，寻找后验概率最大的树。
```
(T, t)_MAP = argmax P(T, t | D) = argmax P(D | T, t) * P(T, t)
```

计算 P(D | T, t) 需要使用进化模型（如Jukes-Cantor）和Felsenstein的剪枝算法。该算法利用动态规划，从叶节点到根节点递归计算每个节点处于每种核苷酸状态的可能性，最终得到在给定树和分支长度下，产生观测数据的整体概率。

这些概率方法的优点在于其坚实的统计学和进化论基础，但计算量通常非常庞大。

5. 系统发育树空间的搜索与评估 🔍

上一节我们学会了如何给一棵树打分，本节中我们来看如何在所有可能的树中寻找得分最高的那一棵。

对于包含 n 个物种的树，可能的无根二叉树拓扑数量是 (2n-5)!!，这是一个天文数字，无法进行穷举搜索。因此，我们需要启发式搜索策略。

马尔可夫链蒙特卡洛（MCMC）搜索
这是一种常用的在树空间中进行随机采样的方法。

初始化：通常从一棵邻接法构建的树开始。
提议新树：通过特定的树形调整操作在当前树附近产生一棵新树。
- 最近邻交换：交换一个内部节点连接的两个子树。
- 子树剪枝与重接：剪下一棵子树，将其连接到另一条分支上。
评估与接受：计算新树的得分（似然值或后验概率）。根据Metropolis-Hastings准则，以一定概率接受新树作为当前状态（即使它比旧树差，但接受概率较低）。这保证了长期采样会收敛到后验分布。
迭代：重复步骤2和3数百万次，在丢弃初始的“老化”样本后，剩余的样本可以用于估计最可信的树或计算分支支持度。

自举法评估
为了评估所建树的可靠性，特别是各个分支的可信度，我们使用自举法。

从原始比对数据中有放回地随机抽取列，生成许多个（如1000个）相同大小的新比对数据集。
为每个自举数据集构建一棵系统发育树。
观察原始树中的每个分支在这些自举树中出现的频率（即自举支持率）。支持率高的分支（如>95%）被认为是可靠的，而支持率低的分支则可信度低。

总结 📝

本节课中我们一起学习了系统发育学的核心内容：

基础概念：了解了系统发育树的基本组成部分和类型。
距离矩阵构建：学习了如何通过进化模型（如Jukes-Cantor模型）从序列比对中计算校正后的遗传距离。
距离法建树：掌握了UPGMA和邻接法等算法，以及它们各自的假设和适用条件。
特征法评分：理解了最大简约法和最大似然/后验概率法如何直接基于序列数据评估树的优劣。
树空间搜索：认识了如何使用MCMC等启发式方法在巨大的树空间中搜索最优树，并用自举法评估结果的可靠性。

通过这些方法，我们可以从分子序列数据中推断出描绘进化历史的系统发育树，这是比较基因组学和进化生物学中至关重要的工具。

P2：L2- 动态规划 - ShowMeAI - BV1RM4y1g76r

在本节课中，我们将要学习序列比对和动态规划。这是课程第一个模块的开始。这个模块将全部关于动态规划，以及代理推断和如何使用隐马尔可夫模型和快速字符串搜索。模块一关注基因组。这些模块的计算基础将是动态规划。具体来说，我们将探讨如何探索指数级空间，但只使用多项式时间，这听起来有些神奇。我们希望你能真正欣赏它的精妙之处。我们将强调这些空间有多大、多复杂。下周我们将介绍隐马尔可夫模型，这是计算机科学的核心工具。我们将研究隐马尔可夫模型算法，如编码、评估、解析和评分。我们还将研究哈希，这在计算机科学中同样无处不在，用于实现基于内容的索引。本周我们将专注于序列比对和比较基因组学。我们将研究局部比对、全局比对，以及如何推断核苷酸水平的进化事件。下节课我们将研究数据库搜索，如何扫描可能具有快速匹配的区域。现在，我们将研究它为何是一个指数级问题，以及我们如何用多项式时间解决它。然后下周我们将使其实际达到线性时间，这很惊人。接着，我们将研究如何用隐马尔可夫模型对基因组建模。这就是模块一。

首先，环顾四周，整个世界充满了运动的事物。它们要么在攻击我们并试图杀死我们，比如冠状病毒；要么是食物，比如栖息在你体内的大约30万亿个细菌细胞；要么是你的一部分，比如另外30万亿个真正的人类细胞，它们与你的细菌一起携带。就像我们用来酿酒（如果你放太久）的酵母；就像第一个拥有脊髓和中央处理单元的有机体，如果蝇；以及比新冠病毒致命得多的蚊子，例如疟疾每年导致数百万人死亡；就像第一个脊椎动物基因组，以及我们所有的脊椎动物亲属，包括我们可爱的毛茸茸的朋友，它们有助于生物医学研究；还有我们的表亲，它们实际上有很多与我们相似的特征，只是它们的额头在眉毛处停止，而我们的额头则继续向上，拥有这个巨大的大脑。所有的生命，你从窗外、海底（如果你去潜水和浮潜）和空中看到的这种令人难以置信的多样性，都从一个共同祖先进化而来。这些共同祖先产生了许多许多谱系，其中一些如你所知已经灭绝。在尤卡坦半岛奇琴伊察以北的希克苏鲁伯曾有一次巨大的小行星撞击，基本上消灭了绝大多数恐龙，除了一些进化成鳄鱼，一些进化成鸟类。灭绝是生命的一部分，很遗憾。这使我们能够开始研究所有这些不同物种之间的进化关系。我们可以做的是，实际上获取密切相关的物种的基因组。这里有一系列不同的果蝇物种。我的叔叔迈克尔·坎比塞利斯实际上是夏威夷的研究人员之一，他研究这些果蝇物种，它们非常有趣和复杂。但每个人都说，为什么人类的基因组比苍蝇的优越那么多？事实证明，我们不能飞，我们没有一百只眼睛，我们不能以毫秒级的速度移动。所以苍蝇有很多优势。但它们的基因组与人类非常相似，基因数量也非常可比。如果你观察越来越多密切相关物种的基因组，你会发现这些物种的基因大致顺序相同，方向也大致相同。在这里，你可以看到基因编码在沃森链上，指向这个方向，显示为黄色；基因编码在克里克链上，指向另一个方向，显示为某种粉紫色。然后你可以看到，它们的顺序大致保留，基于这些显示不同物种基因组坐标匹配的线条。这类比对一次跨越数百个基因，这些基因基本上是共线的。所以你意识到，所有这些物种基本上都是彼此的不同版本。就像一个是运行Windows 98，另一个是运行Windows 2000。它们都只是略有不同，但或多或少保留了相同的功能。这非常重要，因为我们可以使用比较基因组学来识别功能元件。我们基本上可以比较Windows 98和Windows 2000，寻找被保留的指令，这些指令更可能是功能性的。而许多随机的非编码位点是非功能性的。

为了更具体，如果你观察人类基因组，绝大多数人类基因组不编码蛋白质。只有1.5%的人类基因组实际编码构成人类约20，000个基因的氨基酸，并最终构成人类基因组的所有蛋白质功能。这令人难以置信。另外98.5%不编码蛋白质。如果你观察密切相关的物种，如狗、小鼠、大鼠，你会发现这另外98.5%在很大程度上是不保守的。所以人类基因组的绝大部分，大约5%到15%，可以明确识别为保守的。这基本上可以帮助我们精确定位这些基因组中的功能元件。因此，我们基本上可以观察保守区域，并说也许保守区域可以帮助我们找到基因，也许它们还可以帮助我们找到其他可能具有基因调控或控制这些基因功能的东西。所以我们可以使用比较基因组学，这基本上意味着排列不同物种的基因组，以揭示功能元件。例如，蛋白质编码外显子，这些是在剪接发生时保留的用于编码外显子的部分，在小鼠、鸡、鱼等物种中深度保守。你可以看到，每次这里有一个蛋白质编码外显子（用深蓝色表示），在我们所有的脊椎动物亲属中。相比之下，还有许多其他元件也高度保守。问题是它们是外显子还是调控元件等等。当我们讲到课程的进化部分时，我们将研究这些，这实际上稍晚一些。我们将开发估计约束水平的方法。我们将计算编辑操作的数量、替换的数量和空位的数量，作为约束的度量。我们将使用概率模型估计包括回复突变在内的突变数量。我们将结合邻域信息，使用隐马尔可夫模型寻找保守窗口。我们将估计约束隐藏状态的概率。事实上，这就是我在这里展示的内容。对于那些听说过隐马尔可夫模型的人（我们将在两节课中介绍），隐马尔可夫模型基本上允许你基于观察变量估计隐藏状态。这里的观察变量是你看到的保守程度。然后隐藏状态具有是保守状态或非保守状态的后验概率。然后，对于真正保守的区域，该概率更高。这是另一种估计约束的方法。或者，你实际上可以使用将这些物种联系在树中的系统发育关系，来基本理解沿着树发生了多少次替换。我将读出数字，因为我每次投票后都会擦掉它们。所以如果你听到我说一串五个数字，基本上是跟随人数：80以上、60到80、40到60、20到40、0到20。这里的数字是：28， 10， 0， 0， 0。这很棒。然后我给你们的下一个问题是：节奏如何？很好。关于节奏，28人认为正好，5人认为太快，4人认为太慢。聊天中似乎有一个问题：如果你有完整的遗传数据，你会有什么隐藏状态？这是一个很好的问题。问题是为什么你需要一个隐藏状态来估计？答案是，你拥有的是某物是否保守的观察结果。你推断的是该区域是否受到进化约束，因为它可能只是偶然出现保守。基于你拥有的证据量，你能够权衡那个概率。这回答了你的问题吗？这是一个很好的问题。基本上，隐藏状态是它是否真正具有功能，或者是否不具有功能。但观察结果是保守的核苷酸数量。莉莉问什么是功能基因？如果你破坏它，生物体可能会有劣势。它可能不会立即死亡，那个基因可能只在每17年当有巨大的暴风雪时才有所帮助。在进化时间尺度上，前16年可能没问题，第17年所有拥有该基因缺陷版本的生物体可能会死亡。这就是进化的工作方式，它基本上在巨大的进化跨度上求和，以那种方式测量功能。曼努埃尔问，许多不编码蛋白质的基因做什么？人类基因组中大约有20，000个基因，它们编码蛋白质。我们可能知道其中大约80%的功能。另外20%信不信由你，我们实际上不知道它们做什么，这相当引人注目。除了20，000个蛋白质编码基因外，还有大约2，000个长非编码RNA。其中许多作为染色质折叠的支架，因此它们并不真正具有功能，但被转录。也许转录行为本身具有功能，也许它们具有我们尚未理解的功能。它们通常以非常低的丰度表达，因此很难找到它们的功能。除了长非编码RNA，还有无数短非编码RNA。这些可以是帮助将mRNA翻译成蛋白质的tRNA，这些可以是剪接RNA，这些可以是帮助翻译的核糖体RNA。我希望这回答了你的问题。好了，还有一个问题：谁觉得这些问题有帮助？如果你觉得问题有帮助，请给我一个5；如果你觉得一般，给4；3， 2， 1，否则。对于那些提问的人，这应该是一些积极的强化，因为17人说对问题超级兴奋，13人说好，7人在中间，只有1人在底部。所以提问很重要，请继续提问。好了，这就是为什么我们需要比较基因组学的基础。但我告诉你的一切都假设我们实际上可以比对基因组。今天，我们将学习如何实际比对基因组，如何基于这些比对来研究不同类别元件的进化特征。我们可以识别蛋白质编码基因与非编码区域、RNA结构、microRNA和调控基序以及单个基序实例的特定进化模式。这些是我们可以区分的特定类别的元件，不仅基于它们的保守水平，还基于它们的进化模式。我们将在第13或17讲左右更多地学习这些。但今天，我们将专注于如何甚至生成这些比对，如何获取一个物种的一个序列，并与另一个物种的序列对齐。这实际上非常贴近我的心，因为这是我博士论文中的一张幻灯片。所以如果你搜索“kellis nature03。pdf”，你会找到它。是的，很久以前，在一个遥远的星系，我也曾是一名学生。这是我的第一篇论文。我的名字是凯利斯。

如果你看这篇论文的图六，我现在很自豪能在讲座中教授它，这很酷。总之，这是酵母酿酒酵母（我是希腊人，所以我会给你很多词源：糖，然后是霉菌，即真菌，所以糖真菌，然后酿酒酵母与啤酒相同，所以制造啤酒的糖真菌，也称为面包酵母，更政治正确）与奇异酵母、尼克松酵母和比安努斯酵母的比对。我们在这里展示的是一个多序列比对，它跨越多行。然后星号表示比对中该核苷酸完全保守的所有位置。真正令人兴奋的是，如果你看TATA基序，它是完全保守的。基因GAL10，如果你看GAL4基序，它是完全保守的。所以3， 11， 3，这是双螺旋的一圈。如果你看MIG-1基序，它相当保守。然后有第二个BIG1基序，不是特别保守。然后有第二个保守岛，它不匹配。为什么这令人兴奋？因为它意味着我们实际上可以通过观察比对来解读进化。这基本上也意味着我们可以使用这些比对来解码功能元件。我们基本上可以说，保守岛在哪里？那里一定有什么重要的东西，因为所有其他区域都不保守，但那些岛是保守的，所以这可能意味着它们很重要。因此，我们可以解读进化以揭示功能元件。我们今天的目标是使这成为可能，实际比对两个基因。我们将通过首先形式化问题来实现这一点。我在这部分讲得慢一点，因为这是你们所有人将在项目中做的练习，所以在实践中看到它非常重要。因此，我们将提出一个生物学问题，将其转化为计算术语，以不同的方式，我们将研究不同的形式化。然后我们将看到所有这些形式化实际上相当困难，实际上存在指数级数量的可能比对。然后我们将介绍我们将用于解决这些指数级数量问题的计算技术，然后将该计算技术直接应用于序列比对。首先，我们将介绍一个问题，然后我们将介绍该技术，然后我们将两者结合起来。最后，我们将研究一些关于该锚点变体的很酷的高级主题。到目前为止，我们已经看到了比较基因组学和分子进化，以及为什么存在进化及其如何运作的简要介绍，以及序列比对。现在我们将研究如何将这个生物学问题转化为计算机科学术语。首先，基因组随时间变化，存在突变、缺失和插入。这些发生是因为当DNA聚合酶愉快地沿着DNA行走并将DNA复制成DNA时，它基本上打开双螺旋，然后在这里制作第二个副本，在这里制作第二个副本。这是DNA的半保守复制，你在一侧有旧链和新链，在另一侧有旧链和新链。每次，偶尔会出错，它基本上有时插入错误的核苷酸，有时插入一个空位等等。这些事情就这样发生了。听起来像是错误，但突变也是进化的动力。我总是喜欢说，如果我们没有突变，如果工程师设计了进化和复制机制，那么我们就会是完美复制的细菌，我们基本上永远不会发生突变。突变既有害，但也是进化的动力。它们也允许多样性，然后选择保留或拒绝。这些是进化的基本操作。问题是，如果你只有起始序列（即酿酒酵母）和结束序列（即奇异酵母），那么你如何将一个转化为另一个？你如何推断导致从一个到另一个的一系列突变、缺失和插入操作？这就是形式化问题的部分。实际上有一个问题：这是否意味着基因长度在短时间内不保守？不，基因长度不是问题，只要它保持功能，你就不需要担心基因长度。我们要做的第一件事是定义一组进化操作：插入、缺失和突变。使这些操作对称允许时间可逆性，这是我们设计选择的一部分。原因是我们并不真正比较酿酒酵母和奇异酵母，我们并不真正比较人类和老鼠，因为人类是从老鼠进化而来的吗？不，不，我们都从一个共同祖先进化而来。基本上，时间可逆性的原因是，当我们将人类与小鼠比对时，我们希望能够讨论方向性的操作，无论它们是朝一个方向还是另一个方向移动，这都不应该重要。这就是为什么我们有这种时间可逆性。有一些例外，比如甲基化的CpG二核苷酸，绝对是非甲基化的，失去甲基化的胞嘧啶变成胸腺嘧啶，但你不必担心那个。接下来我们需要问的是，我们的最优性标准是什么？我们怎么知道我们做得好？这是一个好的解决方案吗？也许我们想要的是最少数量的操作，即有多少插入、缺失和突变。或者，我们可能想要最小成本，这些操作实际上可能有不同的成本。根据成本，我们可能更喜欢三个廉价操作，而不是两个非常昂贵的操作。第三，可能实际上不可能推断出确切的操作序列，因为当人类转化为小鼠，反之亦然，在这个时间可逆模型中，你基本上可能走一条非常迂回的路径，你将A变成T，然后变回C，然后变回A。看起来从未改变，但实际上我们进行了几次不可见的操作。但我们将选择一条最小路径，一个最小成本转换。这就是奥卡姆剃刀原理。奥卡姆剃刀基本上说，在多个同样好地解释数据的假设中，我将用我的剃刀选择最尖锐的假设，最小的模型。所以从所有我可以拥有的、同样好地解释数据的超级复杂模型中，我将选择这些模型中最小的一个。一旦我们有了最优性标准，我们还必须设计一个实现该最优性或近似该最优性的算法。该解决方案的可处理性将取决于我们在形式化中所做的具体假设。例如，我们可能做出一些捕捉生物学每一个方面的决策，例如CpG二核苷酸不可逆性，或者这种随时间来回的奇怪情况等等。这可能使其与生物学更相关、更正确，可能处理更多特殊情况，但会使我们远离可处理性。所以基本上存在权衡。我们可以做出一些简化假设，使计算更容易，或者处理一些使这更困难的特殊情况。在你们的所有最终项目中，当你们设计这个时，都会有权衡。当然，并非所有决策都是冲突的。一些决策可能既与生物学更相关，又更易处理。有一个著名的例子，佩夫兹纳与桑科夫关于染色体倒位方向性的争论。事实证明，将方向性作为标准并使其更具生物学意义，也使其更易处理，并将其从NP中移除。现在，让我们沿着这个形式化连续体看看，我们可以为进化问题推断使用哪些不同的形式化。我们将看的第一个形式化是最长公共子串。给定两个可能相关的字符串（没有空位的字符串），子串是一组连续的字符。那么你在这里看到的最长公共子串的长度是多少？我看到长度为3，例如TCA匹配DCA相当好，然后有一个长度为4的，即G T C A。所以如果我只偏移一个位置，我得到TCA，但如果我再偏移一个位置，我实际上在这里得到GTCA。这很酷，对吧？我们基本上能够搜索两个序列，将它们彼此扫描，这是你们在脑海中解决的计算挑战，能够找到这三个，然后找到这四个。你使用的那个算法的运行时间是多少？在每个偏移处，你可能比较了一定数量的字符，然后这是一个n平方操作，因为你必须比较序列的每个部分。但你们中的许多人可能更聪明地解决了它，通过对齐单个字符，然后从该对齐扩展，这就是我们将在周二研究的BLAST算法的基础。总之，这是寻找最长公共子串的基本算法的基础。但这不处理空位，插入和删除无法用此建模。所以我们需要稍微增加复杂性。下一个形式化将是最长公共子序列。子序列是最长的不一定连续的字符集。注意，在这里我可以捕获我之前有的TCA，如果我在这里只插入一个空位，我有G和T，如果我在这里插入另一个空位，我也有A。所以现在我可以捕获更多具有共同祖先的字符，因为我容忍插入和删除事件。这在序列二中推断一个空位，一个T的插入，在序列一中推断C和A的删除，然后C到G的突变。如果你注意到，这正是我们之前的例子。所以我们基本上创建了这些之间的最长公共子序列，并发现了一个插入、两个删除和一个突变。这就是我们想要做的。我们想通过一个可以做到这一点的很酷的算法来解决我们提出的原始问题。同样，每个这样的比对立即推断、暗示和关联一系列事件：两个删除、一个插入和一个突变。所以最长公共子序列实际上与编辑距离相同，其中使用统一评分函数在s1和s2之间所需的更改次数，实际上就是最长公共子序列所解决的。但我们可以更进一步。我们可以开始建模空位可能有固定惩罚，例如每个插入和删除对每个字符有唯一成本。或者你可能实际上对不同的空位有不同的惩罚。你可能还对不同类型的操作有不同的惩罚。例如，A和G之间的转换（两者都是相当大的碱基）比A与较小碱基（T或C）之间的转换发生得容易得多，仅仅因为它们在DNA聚合酶内的化学计量学。这些类型的错误，大碱基被大碱基替换更容易，小碱基被小碱基替换更容易，而不是大碱基被小碱基替换，后者更罕见。所以我们可能实际上想要编码一个评分矩阵，而不是对每个突变惩罚完全相同的量。我们可能想要包括一个评分矩阵，基本上对E和G之间的错配只惩罚一半，然后对T和C之间的错配惩罚全额，然后对匹配奖励加一。很好。所以31， 5， 0， 0， 0相当不错。这些是越来越复杂的形式化，因为聚合酶更常将A与G、G与T混淆。所以让我们编码那个。然后另一个形式化可能是改变空位成本。基本上，线性空位成本可能与以前相同。现在，仿射空位成本可能对开始或结束一个空位有很大的初始成本，并对每个额外字符有小的增量成本。因为聚合酶有时会滑动，它可能滑过多个字符，多个核苷酸。或者一般空位惩罚基本上可以允许任何类型的成本，但那不再能用相同的模型计算。你还可以构建一个阅读框感知的空位惩罚，其中三的倍数破坏编码区域，保留编码区域的没问题，但破坏编码区域的非三的倍数（

系统基因组学课程 P20：L20 🧬🌳

在本节课中，我们将学习系统基因组学。我们将探讨如何将之前学习的关于复制、丢失、分化等事件的基因树，整合到物种树和溯祖树中，并利用这些信息来推断基因复制、丢失以及深层溯祖事件。本节课的核心是理解“树中之树”的概念。

概述：基因树与物种树

上一节我们介绍了系统发育分析。本节中，我们来看看如何区分基因树和物种树，讨论两者的“调和”问题，以及种群差异的溯祖过程。

首先，我们将讨论“调和”：如何将基因树映射到物种树上，如何利用这种映射推断直系同源、旁系同源、基因复制和丢失事件，以及调和问题的具体案例和最大简约法调和。

基因树与物种树的调和 🧩

基因树是在物种谱系内部进化的。我们之前研究的系统发育树，其分支通常是无标签的。但现在，我们需要考虑序列所属的物种。例如，一段序列来自狗，另一段来自人。通过了解比对序列的物种归属，我们可以理解它们在物种树内部的进化过程。

核心概念与定义

物种树：描述了物种之间的进化关系。每个物种可以被看作一个“基因袋”。
基因树：描述了特定基因在物种间的进化历史。
调和：将基因树的节点映射到物种树节点上的过程，以推断进化事件。

以下是调和过程中的关键事件定义：

复制事件：发生在基因树节点上，当该节点映射到的物种树位置与其任一子节点相同时。公式表示为：reconciliation(v) == reconciliation(left_child(v)) 或 reconciliation(v) == reconciliation(right_child(v))。
丢失事件：当基因树中一个节点的父节点映射到的物种树位置，与该节点映射的位置不同，且中间跨越了物种树分支时被推断出来。
直系同源：两个基因通过物种形成事件找到最近共同祖先。它们在不同物种中很可能执行相同功能。
旁系同源：两个基因通过基因复制事件找到最近共同祖先。它们可能获得新的功能或产生基因剂量效应。

调和算法简介

最简单的算法是最大简约法调和。其递归过程如下：

对于基因树的叶节点（物种），直接映射到物种树对应的物种。
对于基因树的内部节点 v，其映射位置是其左右子节点映射位置的最近共同祖先。
根据映射结果判断事件：如果节点 v 的映射位置与其任一子节点相同，则为复制事件；如果节点 v 的父节点映射位置与 v 不同，且跨越了物种树分支，则推断存在丢失事件。

这个算法从叶节点开始，自底向上遍历基因树，完成映射和事件推断。

系统基因组学：利用全基因组数据 🌍

上一节我们介绍了单个基因树的调和。本节中，我们来看看如何利用全基因组范围内的多个基因树来构建更可靠的物种树，并反过来利用物种树提高基因树重建的准确性。

构建物种树的方法

给定一组基因树，目标是找到一个能最好概括这些基因树信息的物种树。主要方法有：

超级矩阵法：将所有基因的序列比对拼接成一个大的超级矩阵，然后基于此矩阵构建系统发育树。
超级树法：分析每个基因树支持的“分裂”（即物种分组），然后整合所有基因树的分裂信息，构建一个共识树。
最小化进化事件法：选择一个物种树假设，计算将所有基因树调和到该物种树所需的总的复制和丢失事件数（或深层溯祖事件数），选择使总成本最小的物种树。

利用物种树优化基因树重建

传统的流程是：重建基因树 -> 将其调和到物种树。但基因树重建可能出错，导致错误推断大量复制和丢失事件。

更优的方法是：将物种树信息作为先验知识，整合到基因树重建过程中。我们可以建立一个生成式概率模型，该模型考虑：

物种树拓扑结构和分歧时间。
基因特异性进化速率和物种特异性进化速率。分支长度可以建模为：分支长度 = 基因速率 * 物种速率 * 分歧时间。
复制与丢失模型。
序列进化模型（如HKY、Jukes-Cantor等）。

通过贝叶斯推理，我们可以计算在给定物种树和序列数据的情况下，某个基因树拓扑结构的后验概率，从而找到最可能的基因树。这种方法能显著提高直系同源基因推断的准确性和减少错误推断的进化事件。

种群水平的进化：溯祖理论 👥

之前我们关注的是物种间的进化。本节中，我们来看看种群内部、个体之间的进化过程，这需要适应更短的时间尺度。

赖特-费希尔模型（前进时间）

这是一个前进时间模型，用于研究遗传漂变等效应。它假设：

种群大小固定为 N。
随机交配。
中性进化（无选择）。
非重叠世代。

在每一代，每个个体会产生大量配子，下一代个体从这些配子中随机抽取形成。因此，一个等位基因可能没有后代，也可能有多个后代。

溯祖模型（后退时间）

与赖特-费希尔模型等价但更高效的是溯祖模型，它向后追溯时间。其核心思想是：随机抽取当代的两个等位基因，追溯它们的祖先，直到找到最近共同祖先。

溯祖时间：k 个谱系追溯到其最近共同祖先所需的时间。
概率分布：在种群大小 N 很大时，k 个谱系在 t 代内未发生溯祖的概率近似服从指数分布：P(T > t) ≈ exp( - (k choose 2) * t / (2N) )。

这个模型可以高效地模拟种群中基因谱系的合并历史。

多物种溯祖与不完全谱系分选

将溯祖模型扩展到物种树，就得到了多物种溯祖模型。物种树的每个分支代表一个种群，有其特定的种群大小和存在时间。

不完全谱系分选：当种群很大或物种形成时间很短时，来自不同物种的两个等位基因可能没有足够时间在祖先种群中溯祖，它们的溯祖事件可能发生在更古老的祖先种群中。
结果：这会导致基因树的拓扑结构与物种树的拓扑结构不一致，而这种不一致并非由于基因复制或丢失引起。

ILS是导致基因树与物种树差异的重要原因之一，可用于推断物种分化时间和祖先种群大小。

统一模型：等位基因、基因座与物种 🧬

我们已看到两种解释基因树/物种树冲突的视角：复制/丢失 vs. 溯祖/ILS。然而，传统模型无法同时处理这两种情况。

一个统一的框架是等位基因在基因座中进化，基因座在物种树中进化的模型。在这个模型中：

基因座树在物种树中向前进化，经历复制和丢失事件。
在每个基因座内部，等位基因通过溯祖过程向后进化。
复制事件产生的新基因座，其等位基因可以在种群中固定、丢失或保持多态性。

这个DL-溯祖模型结合了前进时间的复制/丢失过程和后退时间的溯祖过程，能够更全面地建模基因组进化历史。

重组与祖先重组图 🔀

之前的模型都假设无重组。但真实情况中，重组会打断谱系的历史。

重组的影响：基因组上不同区域可能具有不同的进化历史（拓扑结构）。
祖先重组图：一种表示包含重组事件的进化历史的图结构，它不是一棵树，而是一个有向无环图。
建模方法：可以使用序列马尔可夫溯祖模型等，将基因组视为从左到右的序列，在不同位置之间，基因树的拓扑结构以一定的概率（重组率）发生切换。这类似于一个隐马尔可夫模型，其中隐藏状态是局部的树拓扑。

总结 📚

本节课中我们一起学习了：

基因树与物种树的调和：如何映射并推断复制、丢失事件，定义直系/旁系同源。
系统基因组学：如何利用全基因组数据构建更稳健的物种树和基因树，通过整合物种树信息的生成式模型显著提升推断准确性。
种群遗传学模型：包括前进时间的赖特-费希尔模型和后退时间、更高效的溯祖模型。
不完全谱系分选：多物种溯祖模型下的重要现象，能导致基因树与物种树冲突，并用于推断历史参数。
统一进化模型：整合基因复制、丢失和溯祖过程的DL-溯祖模型。
重组：重组如何使进化历史非树状，以及祖先重组图的概念。

通过系统基因组学，我们将微观的基因进化与宏观的物种形成、种群历史联系起来，获得了对生命进化历史更完整、更深入的理解。

课程 P21：L21- 癌症基因组学 🧬

在本节课中，我们将要学习癌症基因组学的基础知识。我们将从理解癌症的核心特征开始，探讨如何通过外显子组和全基因组测序来发现驱动癌症的基因突变，并进一步了解超越基因序列的癌症驱动因素，包括表观遗传学改变和肿瘤与免疫系统的相互作用。

癌症的基石：癌基因、抑癌基因与癌症特征

上一节我们介绍了课程的整体框架，本节中我们来看看理解癌症的基础概念。癌症并非单一疾病，而是一类具有共同特征的细胞异常增殖性疾病。

2000年，Bob Weinberg 和 Doug Hanahan 在《细胞》杂志上发表了一篇名为“癌症的特征”的综述，将癌细胞区别于正常细胞所获得的能力总结为六大类，并在2011年补充了四个新特征。这些能力是肿瘤形成所必需的“技能”。

以下是癌症的六大核心特征：

避免细胞凋亡：癌细胞能够逃避程序性细胞死亡信号。
实现生长信号自给自足：癌细胞不再依赖外部生长信号，其“油门”被卡住。
对生长抑制信号不敏感：癌细胞无视周围环境发出的停止生长信号，其“刹车”失灵。
诱导持续血管生成：肿瘤能够发出信号，诱导周围环境形成新的血管，为其提供氧气和营养。
无限复制潜能：癌细胞能够突破正常细胞分裂次数的限制（如端粒缩短），实现无限增殖。
组织浸润和转移：癌细胞能够侵入周围组织并扩散到身体其他部位，这是被正向选择的进化优势。

此外，2011年提出的四个新兴特征包括：

细胞能量代谢异常：癌细胞线粒体功能失调，产生大量能量。
避免免疫摧毁：癌细胞能够与免疫系统相互作用，逃避免疫系统的攻击。
基因组不稳定性和突变：癌细胞倾向于制造基因组不稳定性和高突变率，通过“数量游戏”增加获得有益突变的机会。
肿瘤促进性炎症：肿瘤可以利用炎症反应，招募免疫细胞并加以利用。

肿瘤获得这些特征的顺序并不重要，关键在于最终全部获得。由于存在多种基因和调控路径，每个肿瘤都通过独特的突变组合来实现这些特征，这导致了肿瘤的异质性。

驱动突变与乘客突变：识别癌症的“元凶”

理解了癌症的特征后，我们需要找到导致这些特征的分子变化。在肿瘤发生过程中，细胞会积累大量突变，但并非所有突变都同等重要。

核心概念在于区分驱动突变和乘客突变：

驱动突变：直接赋予癌细胞生长优势，促进肿瘤发展的突变。它们通常影响癌基因、抑癌基因等关键通路。
乘客突变：在肿瘤发展过程中随机积累，不直接贡献于肿瘤生长优势的突变。它们是进化过程的“搭便车者”。

识别驱动突变的主要挑战在于，它们在所有突变中只占极少数（通常每个肿瘤仅有5-20个），而乘客突变数量庞大。常用的策略是寻找复发性突变，即在许多不同患者的相同基因或通路上反复出现的突变。

驱动突变主要影响四类基因：

癌基因：正常时促进细胞生长（原癌基因），突变后过度活跃，成为“卡住的油门”。例如 RAS 基因家族。
- 原癌基因 ->（突变）-> 癌基因 -> 过度生长
抑癌基因：正常时抑制细胞分裂，充当“刹车”。突变导致功能丧失，细胞生长失控。例如 TP53 基因。
- 抑癌基因 ->（功能丧失性突变）-> 生长失控
突变基因：参与DNA修复或维持基因组稳定。突变会导致基因组不稳定性和突变率升高。
表观突变基因：调控全局表观基因组。其功能失调会导致广泛的基因表达改变。

此外，基因组结构变异（如染色体易位产生的融合基因）也是常见的驱动事件。例如，9号与22号染色体易位产生的 BCR-ABL 融合基因是慢性髓系白血病的驱动因素。

从外显子组测序中发现复发突变与异质性

上一节我们介绍了驱动突变的概念，本节中我们来看看如何利用外显子组测序来发现它们。识别癌症驱动基因的主要方法包括全基因组关联研究（GWAS）、遗传连锁分析和体细胞突变分析。

体细胞突变分析是发现肿瘤特异性驱动突变的关键。其核心方法是配对测序：同时测序患者的肿瘤组织和正常组织（如血液），通过比较找出只在肿瘤中出现的体细胞突变。

以下是分析体细胞突变的主要步骤：

突变检测：使用工具（如 Mutect2）比较肿瘤和正常样本的测序数据，统计富集地识别体细胞突变。
过滤与注释：应用一系列过滤器排除测序错误、常见多态性位点等，然后对剩余突变进行功能注释。
寻找复发信号：在大量患者中寻找在基因水平（同一基因多次突变）、氨基酸水平（同一蛋白位点反复突变）或通路水平（同一通路多个基因被破坏）上复发的突变。例如，BRCA1 和 BRCA2 基因在遗传性乳腺癌中就有数百种不同的复发突变。

然而，癌症也是一个进化过程，这导致了肿瘤内异质性。一个肿瘤并非由完全相同的细胞组成，而是包含多个具有不同突变谱的细胞亚群（克隆）。这种异质性使得驱动突变可能只存在于一部分肿瘤细胞中，在测序数据中表现为较低的等位基因频率。

克隆进化模型描述了肿瘤的发展：起始驱动事件导致一个优势克隆扩增，随后该克隆内细胞继续积累新的突变（包括新的驱动事件），形成不同的亚克隆。治疗就像一个选择压力，可能消灭大部分克隆，但残留的、具有耐药性的亚克隆会再次扩增，导致复发。

通过分析多个肿瘤区域或不同时间点的样本，我们可以构建肿瘤的系统发育树，追踪不同克隆的进化历史，并识别在克隆扩增中起关键作用的驱动事件。

全基因组测序：探索非编码区驱动突变

外显子组测序主要关注蛋白质编码区，但人类基因组中超过98%的区域是非编码的。本节中我们来看看如何利用全基因组测序探索非编码区的癌症驱动突变。

绝大多数体细胞突变位于非编码区，且主要是乘客突变。在非编码区寻找驱动突变面临巨大挑战：

功能单元不明确：在编码区，基因是天然的功能单元。而在非编码区，功能单元可能是增强子、启动子或更复杂的调控网络，其边界难以界定。
背景突变率差异大：不同基因组区域的突变率受染色质状态、复制时间、基因表达水平等因素影响，存在巨大差异。例如，不活跃的染色质区域（如异染色质）通常比活跃的开放染色质区域突变率更高。
识别收敛性更困难：在编码区，不同患者在同一基因的不同位置突变可视为收敛于该基因。在非编码区，我们需要将散布在不同位置的突变关联到它们共同调控的靶基因上，这需要借助基因调控网络的知识。

为了应对这些挑战，我们需要：

校正背景突变率：建立统计模型，考虑染色质可及性、复制时间、序列背景等多种因素，预测每个基因组区域的“预期”突变率，然后找出突变显著高于预期的区域。
定义调控单元：利用表观基因组学数据（如组蛋白修饰、染色质开放区域）定义增强子、启动子等调控元件。
寻找调控收敛：即使突变位于基因组不同位置，如果它们都影响同一个靶基因或同一个信号通路，也可以被认为是驱动事件。这需要整合染色质相互作用（如Hi-C）和基因调控网络数据。

非编码驱动突变可以通过多种方式发挥作用，例如改变转录因子结合位点、破坏或创造microRNA结合位点、通过结构变异重写基因调控环境等。

超越突变：表观遗传改变与肿瘤微环境

癌症驱动因素不仅限于DNA序列的改变。本节中我们来看看表观遗传改变和肿瘤微环境在癌症中的作用。

表观遗传改变，如DNA甲基化、组蛋白修饰和染色质结构的全局性变化，也可以驱动癌症。这些改变通常由表观突变基因（如DNA甲基转移酶、组蛋白去乙酰化酶）的功能失调引起。它们能导致全基因组范围的基因表达重编程，影响细胞分化状态，从而促进肿瘤发生。针对这些表观遗传改变的药物（如表观遗传疗法）正在成为新的治疗策略。

肿瘤微环境是指肿瘤细胞周围的复杂生态系统，包括成纤维细胞、免疫细胞、血管内皮细胞等。肿瘤与微环境的相互作用对其生长、侵袭和转移至关重要。例如，肿瘤可以通过释放信号分子来抑制免疫细胞功能，营造一个免疫抑制性的微环境，从而逃避免疫系统的攻击。

单细胞基因组学技术使我们能够以前所未有的分辨率解析肿瘤内和微环境中的细胞异质性。通过单细胞RNA测序，我们可以：

识别不同的肿瘤细胞亚型和它们所处的细胞周期状态。
刻画肿瘤浸润免疫细胞（如T细胞、B细胞、巨噬细胞）的组成和功能状态。
推断肿瘤细胞的拷贝数变异。

肿瘤免疫学与免疫治疗 🛡️

上一节我们提到了肿瘤微环境中的免疫细胞，本节中我们深入探讨肿瘤与免疫系统的相互作用，以及如何利用这种相互作用进行治疗。

肿瘤与免疫系统的关系是一个动态的“编辑”过程。免疫系统会识别并清除具有高免疫原性（即能被免疫系统识别为“非己”）的肿瘤细胞。然而，在进化压力下，肿瘤细胞会通过多种机制逃避免疫监视：

降低免疫原性：选择抗原性较低的突变。
下调抗原呈递：减少主要组织相容性复合体（MHC）分子的表达。
营造免疫抑制微环境：上调抑制性免疫检查点分子（如PD-L1），招募抑制性免疫细胞（如调节性T细胞）。

免疫治疗旨在打破这种免疫抑制，重新激活患者自身的免疫系统来攻击肿瘤。主要策略包括：

免疫检查点抑制剂：使用抗体阻断PD-1/PD-L1或CTLA-4等抑制性通路，相当于“松开T细胞的刹车”。这类药物的疗效与肿瘤的突变负荷和新抗原数量正相关，因为更多的突变意味着免疫系统有更多可以识别的靶点。
过继性细胞疗法（如CAR-T）：从患者体内提取T细胞，在体外进行基因工程改造，使其表达能够特异性识别肿瘤抗原的嵌合抗原受体（CAR），然后扩增并回输到患者体内，相当于给免疫细胞装上“GPS导航”去攻击肿瘤。
癌症疫苗：根据患者肿瘤特有的新抗原设计个性化疫苗，刺激机体产生针对这些新抗原的特异性免疫反应。

预测免疫治疗疗效是一个活跃的研究领域，涉及计算预测新抗原、分析T细胞受体（TCR）库的多样性以及评估肿瘤免疫微环境的状态。

总结

本节课中我们一起学习了癌症基因组学的核心内容。我们从癌症的六大（及新兴）特征出发，理解了癌细胞获得的关键能力。我们探讨了如何通过外显子组和全基因组测序来区分驱动突变与乘客突变，并利用复发性和克隆进化分析来识别关键的癌症基因。我们认识到，驱动因素不仅限于编码区突变，非编码区调控元件和全局性表观遗传改变也扮演着重要角色。最后，我们深入了解了肿瘤与免疫系统复杂的相互作用，以及如何通过免疫检查点抑制剂、CAR-T疗法和癌症疫苗等免疫治疗策略来对抗癌症。癌症基因组学是一个快速发展的领域，其研究成果正在不断转化为更精准、更有效的癌症诊断和治疗方法。

📚 课程 P22：单细胞基因组学与深度学习方法 🧬

在本节课中，我们将学习单细胞基因组学的基础知识，特别是用于单细胞基因组学的深度学习方法。我们的目标是能够在单细胞水平上推断信息。我们将讨论单细胞RNA测序、单细胞ATAC测序以及单细胞多组学技术，探讨如何扩展这些技术，以及如何处理噪声、双联体和其他单细胞测序问题。课程将涵盖单细胞RNA测序数据分析的计算挑战、已应用于此的深度学习方法，并展望超越单细胞RNA测序的未来方向。

🧫 为何要进行单细胞分析？

上一节我们介绍了课程概述，本节中我们来看看为何需要分析单细胞。当我们观察批量数据时，例如脑组织活检或肝脏活检中测量的RNA，我们得到的是许多维度的平均值。这包括不同细胞类型、不同细胞状态、不同发育阶段、不同响应以及细胞周期中不同时间点的平均值。在显微镜下，细胞是极其异质性的。因此，平均这些细胞的表达量会得到一个无法代表其中任何一个细胞的结果。

第二个原因是细胞会分化。例如，它们从多能干细胞开始，然后分化成所有不同的血统谱系。即使在同一种细胞类型内，个体细胞之间也存在巨大差异，这是批量分析无法捕捉的。此外，当细胞对环境做出反应时，它们拥有的特定受体只能部分地捕捉细胞环境。在组织或器官水平上，你得到的是这些反应的平均值，而实际上，每个细胞可能根据其检测到的病毒蛋白量而有非常不同的反应。

最后，进行单细胞分析的另一个原因是，如果你观察循环肿瘤细胞或斑马鱼早期胚胎，可能只有八个细胞。如果我们使用批量RNA捕获方法，当样本中只有百分之一的细胞是肿瘤细胞时，我们可能会完全错过它。同样，如果某个区域细胞非常少，也可能被遗漏。

批量数据的另一个主要问题是，当你观察一个基因的表达分布时，真实数据可能呈现非常双峰性的分布。一些细胞可能大量表达该基因，而另一些细胞可能表达极少。但当你观察批量数据时，你得到的是这些的平均值，完全错过了这种双峰性。批量数据的另一个问题是稀有事件可能丢失。如果你观察平均值的分布，最终分布会非常接近平均值，而实际上可能有一些个体细胞对该基因有极强的信号，这些信号可能在批量数据中完全丢失。

🔬 单细胞分析的传统与现代方法

上一节我们探讨了单细胞分析的必要性，本节中我们来看看实现单细胞分析的方法。单细胞分析并非新概念，人们长期以来一直试图研究单个细胞。传统方法使用显微镜观察带有大约五个探针的单个标记RNA，这可以捕获丰富的空间信息，但只能捕获很少的基因。

随着电信号、FISH技术和原位测序技术的发展，这种方法得到了一定程度的扩展，可以观察大量探针。通过单细胞RNA测序，我们现在可以开始观察成千上万个基因在成千上万个细胞中的表达。

许多单细胞RNA测序技术的底层基础技术是RT-PCR。概念是从活细胞开始，提取RNA，进行逆转录，添加引物，然后测量群体的RNA表达。早期研究通过流式细胞术或微孔板捕获等方法费力地分离单个细胞，然后使用这些单个细胞测量尽可能多的基因。当然，这需要在孔内进行大量扩增，以通过单个孔中的单个细胞的RNA测序捕获大量基因。

单细胞RNA测序的关键优势在于，你现在不仅可以对五个或十个探针进行此操作，还可以对成千上万个基因在成千上万个细胞中进行此操作。基本思路是：获取组织，可以进行整个组织的RNA测序，也可以创建包含单细胞和RNA制备的单细胞悬浮液，然后关联单个细胞，分离这些单细胞。你可以使用移液器吸取一个细胞，也可以使用细胞分选仪逐个发射细胞，可视化其特性并捕获它们，或者让它们通过。最近，微流控技术通过一系列腔室引导细胞移动，实现了对单个细胞的操控。

🧪 单细胞RNA测序技术流程

以下是单细胞RNA测序的基本步骤：

组织获取与处理：从生物样本中获取组织，并制备成单细胞悬浮液。
单细胞分离：使用微流控、液滴或孔板等技术将单个细胞物理分离。
细胞裂解与RNA捕获：裂解细胞释放RNA，并通过带有条形码的磁珠或溶液捕获mRNA。
逆转录与扩增：将RNA逆转录为cDNA，并通过PCR进行扩增，为测序准备足够的材料。
文库构建与测序：添加测序接头，构建测序文库，然后进行高通量测序。
数据分析：对测序产生的数据进行生物信息学分析，包括基因表达定量、细胞聚类和轨迹推断等。

🚀 高通量单细胞技术的演进

上一节我们介绍了基本流程，本节中我们来看看技术如何实现高通量化。早期技术可能一次只能分析十个细胞。但近年来，技术已显著改进，现在可以一次分析数万甚至数十万个细胞。让我们看看哪些底层技术实现了单细胞RNA测序技术的这种指数级扩展。

我们经历了从使用移液器捕获单个细胞并将其放入孔中，到使用显微镜和毛细管移液器，最终到使用流式细胞分选，通过激光决定关心哪些细胞。最近，我们能够进行激光捕获显微切割，选择特定的细胞类型。更近期的技术包括基于微流控的单个细胞分离，有效地捕获这些细胞并将其与特定的条形码一起包裹在凝胶液滴中，这些液滴在油相介质中彼此分离。最近还有基于液体的收集和基底辅助的组合条形码技术。

无论使用何种具体技术，关键思想是相同的：将细胞物理分离到某种“孔”中。这可以是一个大管，也可以是96孔板或384孔板中的单个孔。传统方法是通过将这些细胞物理分离到微流控芯片中，然后每次对整孔进行测序，进行全长RNA测序，从而能够检测每个孔中单个细胞的基因表达、剪接和其他附加信息。

最近，通过像Drop-seq和inDrop这样的方法，将细胞捕获在水凝胶液滴内部，而不是进行分选，从而实现了更高的通量。经过物理分离的细胞也使用简化方法（如CEL-seq、MARS-seq等）进行了分析。在这些方法中，通常只测序基因的3‘端UTR区域，逆转录从polyA尾开始，使用 oligo-dT 引物启动逆转录，从而测量基因表达，但不一定能捕获全长转录本、剪接模式等信息。

💡 主流单细胞技术平台比较

以下是几种主流单细胞技术的简要比较：

Smart-seq：将单个细胞分离到独立的孔或微流控腔室中，进行全长RNA测序。通量较低（通常一次几十到几百个细胞），成本较高，但能获得更完整的转录本信息。
10x Genomics (Drop-seq原理)：基于液滴的高通量技术。将单个细胞与带有独特条形码的凝胶磁珠共同包裹在油包水液滴中，在液滴内进行逆转录和条形码标记。通量极高（一次可处理数千至数万个细胞），成本低，但通常只捕获3‘端序列。
Split-seq/sci-RNA-seq (组合索引)：一种基于组合条形码索引的技术。不进行物理分离，而是通过多轮细胞池化、条形码标记和再混合，为每个细胞的RNA分子分配一个独特的组合条形码序列。通过有限的条形码池（如几百个）的多次组合，理论上可以标记海量细胞，成本极低，且无需特殊设备。

⚠️ 单细胞数据分析的挑战与噪声处理

上一节我们介绍了高通量技术，本节中我们来看看数据分析中的挑战，特别是噪声处理。单细胞RNA测序面临几个主要挑战。

第一个挑战是rRNA污染。每个细胞内都有大量核糖体RNA（rRNA），其丰度远高于mRNA。如果不处理，测序数据中绝大部分 reads 将来自rRNA。常用的去除方法包括 polyA 选择（捕获带有polyA尾的mRNA）或使用探针特异性去除rRNA。

第二个挑战是PCR偏好性。当扩增特定的mRNA分子时，第一个被扩增的分子可能会在“富者愈富”的模式中占据优势，导致该分子的 reads 数被严重高估。为了缓解这个问题，可以在建库时加入唯一分子标识符（UMI）。UMI是连接在每条原始RNA分子上的短随机序列，通过统计具有相同UMI的 reads，可以区分来自同一原始分子的PCR重复副本和来自不同原始分子的 reads，从而更准确地定量基因表达。

第三个挑战是质量控制。我们需要判断一个细胞是否是高质量的。常用的质控指标包括：比对到基因组的 reads 比例、检测到的基因数量、线粒体基因表达比例（高比例可能表示细胞状态不佳或死亡）、每个细胞的总 reads 数等。不符合质控标准的细胞（如基因数过少、线粒体基因比例过高）通常会被过滤掉。

一个特别的挑战是双联体（Doublets）。在液滴技术中，一个液滴可能偶然包裹了两个细胞。这两个细胞的所有RNA分子都会被打上相同的细胞条形码，在数据分析中会被误认为是一个“细胞”，但其基因表达谱是两个细胞类型的混合，会干扰下游的细胞分群和差异表达分析。识别双联体的方法包括：利用已知的细胞类型标记基因（混合表达谱）、使用样本多重标记（如果两个细胞来自不同样本的混合，其条形码组合会异常），或使用专门的算法进行预测。

🧮 单细胞RNA-seq数据分析流程

处理完原始数据后，我们得到一个基因（行）x 细胞（列）的计数矩阵。从这个矩阵出发，标准的分析流程包括：

质量控制与过滤：根据上述指标过滤低质量细胞和低表达基因。
标准化：消除由于测序深度不同导致的细胞间总 reads 数的差异。常用方法包括将每个细胞的计数除以一个大小因子，然后进行对数转换（如 log1p）。
高变基因选择：并非所有基因都提供有用的信息。选择在细胞间表达变异程度高的基因进行下游分析，通常能捕获更多的生物学信号。
降维：单细胞数据维度极高（数万个基因）。为了可视化和计算，需要使用主成分分析（PCA）、t-分布随机邻域嵌入（t-SNE）或均匀流形近似与投影（UMAP）等方法将数据投影到低维空间（如2维或3维）。
细胞聚类：在低维空间中，根据细胞之间的相似性（如欧氏距离）将细胞分组。这些簇通常对应于不同的细胞类型或状态。
差异表达分析：比较不同簇之间或不同实验条件（如疾病 vs 对照）下细胞的基因表达，找出显著差异表达的基因。
轨迹推断：对于发育或分化过程，细胞状态是连续的。轨迹推断算法（如Monocle, PAGA）试图根据表达相似性，将细胞排列成一条或多条假时间轨迹，以模拟分化或激活过程。
细胞类型注释：利用已知的细胞类型标记基因，为每个聚类分配生物学上有意义的细胞类型名称。

🤖 用于单细胞分析的深度学习方法

上一节我们概述了标准分析流程，本节中我们重点看看深度学习方法如何应用于单细胞数据分析。深度学习模型，特别是自编码器，在解决单细胞数据分析的特定挑战方面显示出巨大潜力。

批次效应校正：当数据来自不同实验批次时，批次间的技术差异可能掩盖真实的生物学差异。MMD-ResNet 等方法使用自编码器学习细胞的一个低维表示（潜空间），并在此过程中通过最大均值差异（MMD）等损失函数，强制不同批次细胞在潜空间中的分布对齐，从而消除批次效应。
细胞聚类：DESC 是一种深度嵌入聚类方法。它使用自编码器学习细胞表示，并同时优化一个聚类目标函数。该方法以分层和迭代的方式进行，能够获得细胞类型特异性的批次校正和清晰的聚类边界。
数据插补（Imputation）：单细胞数据含有大量“漏失（dropout）”的零值（基因在某个细胞中未被检测到，可能是技术原因，也可能是该细胞确实不表达）。AutoImpute 等基于自编码器的方法，通过学习数据的低维表示并重建完整的表达矩阵，可以智能地填充这些零值，恢复基因-基因之间的相关性，并改善下游分析。
变分自编码器（VAE）用于建模：scVI 是一个强大的概率框架，它使用变分自编码器对单细胞RNA-seq数据进行建模。scVI 将观察到的计数数据建模为零膨胀负二项分布等，并同时推断细胞潜变量、批次效应和文库大小因子。它可以用于批次校正、降维、差异表达分析和插补等多种任务，提供了一个统一的分析框架。

🌌 超越RNA：单细胞多组学与未来展望

单细胞分析不仅限于RNA。现在可以对同一细胞同时或关联分析多种分子层面，即单细胞多组学。

单细胞ATAC-seq：测量单个细胞中染色质的可及性（开放程度），从而推断其活跃的调控元件和潜在的转录因子活性。
单细胞DNA测序：分析单个细胞的基因组变异。
单细胞蛋白组学：例如通过质谱流式细胞术（CyTOF）或抗体条形码（如CITE-seq）测量细胞表面或内部的蛋白质丰度。
空间转录组学：在保留组织空间位置信息的前提下，测量局部区域的基因表达，将细胞类型映射回其原始组织结构。

整合这些多组学数据（如同时分析同一细胞的基因表达和染色质可及性）是当前的前沿。计算方法如 Multi-Omics Factor Analysis (MOFA) 和 Seurat v4 的整合工具，旨在将不同模态的数据对齐，以获得对细胞状态更全面、更深入的理解。

📝 总结

在本节课中，我们一起学习了单细胞基因组学的核心内容。我们从为何需要单细胞分析开始，探讨了批量分析的局限性。接着，我们回顾了单细胞RNA测序技术的发展历程，从低通量方法到如今主流的基于液滴（如10x Genomics）和组合索引的高通量平台。我们深入了解了单细胞实验和数据面临的主要挑战，包括rRNA污染、PCR偏好性、质量控制以及双联体问题。

然后，我们梳理了标准单细胞RNA-seq数据分析流程，包括质控、标准化、降维、聚类、注释和轨迹推断。重点介绍了深度学习方法在这一领域的应用，例如利用自编码器进行批次效应校正（MMD-ResNet）、深度聚类（DESC）、数据插补（AutoImpute）以及使用变分自编码器进行概率建模（scVI）。最后，我们展望了超越RNA的单细胞多组学世界，包括ATAC-seq、蛋白组学和空间转录组学，这些技术正在共同推动我们对细胞异质性和复杂生物系统的理解达到前所未有的深度。

单细胞基因组学是一个飞速发展的领域，它正在彻底改变我们对生物学和疾病的认识。掌握其基本原理和计算方法，是进入这一前沿领域的钥匙。

🧬 P3：L3 - 局部对齐、哈希与BLAST对齐分数

在本节课中，我们将深入学习序列比对和数据库搜索。我们将探讨全局比对与局部比对的区别，学习如何通过哈希技术加速字符串匹配，并深入了解BLAST算法的核心思想。最后，我们将探讨比对分数背后的概率学基础，理解这些分数如何反映序列间的进化关系。

🔄 全局比对与局部比对

上一节我们介绍了使用动态规划进行全局序列比对的方法。本节中，我们来看看如何将其调整为局部比对。

局部比对的目标是寻找两个序列中相似度最高的子串进行比对，而不是强制比对整个序列。这在生物学上很有意义，例如寻找基因中的保守结构域或在长染色体中定位短基因。

以下是全局比对（Needleman-Wunsch算法）与局部比对（Smith-Waterman算法）的核心区别：

初始化：全局比对从矩阵左上角（0,0）开始。局部比对可以从矩阵任何位置开始。
递推规则：全局比对的递推公式为 M[i,j] = max(M[i-1,j-1] + s(i,j), M[i-1,j] + gap, M[i,j-1] + gap)。局部比对在此基础上增加了一个选项：M[i,j] = max(0, M[i-1,j-1] + s(i,j), M[i-1,j] + gap, M[i,j-1] + gap)。这个“0”选项允许比对在任何位置重新开始。
终止：全局比对在矩阵右下角结束。局部比对可以在矩阵中任何得分最高的位置结束。

通过以上调整，算法就能找出两个序列内部最相似的片段，而不是强迫进行端到端的比对。

⚡ 半数值字符串匹配与哈希

我们了解了如何高效地进行比对，但当需要在超长文本（如染色体）中快速搜索一个短模式（如基因）时，直接使用动态规划（O(nm)复杂度）仍然太慢。本节中我们来看看如何利用哈希技术实现接近线性的搜索速度。

核心思想是将字符串视为数字，从而将字符间的比较转化为数字间的比较。Karp-Rabin算法正是基于这一思想。

以下是该算法的基本步骤和关键优化：

将模式串转换为数字：例如，模式串“31415”可以解释为数字31415。
滑动窗口计算：在文本串上滑动一个与模式等长的窗口，计算每个窗口对应的数字。
快速更新：通过数学推导，可以从上一个窗口的数字Y(i-1)快速计算出当前窗口的数字Y(i)，公式为：Y(i) = 10 * (Y(i-1) - 10^(m-1) * digit_out) + digit_in。这避免了每次重新计算整个数字，使窗口更新成为常数时间操作。
哈希处理：直接处理长数字（如1000位）在计算机中效率很低。我们通过取模运算（哈希）将大数字映射到一个小范围内，例如 H(x) = x mod 13。这样，我们只需要在小整数上进行运算。
处理碰撞：哈希可能导致不同的原始字符串映射到相同的哈希值（碰撞）。因此，当哈希值匹配时，我们必须回退到原始字符串进行逐字符验证，以确保是真正的匹配，而非伪命中。

通过哈希，我们将每次比较的复杂度从O(m)降到了O(1)（期望情况），使得整个搜索过程的期望时间复杂度接近O(n+m)。

🚀 BLAST算法：快速数据库搜索

基于哈希的思想，我们可以构建更强大的工具来搜索生物序列数据库。本节中我们来看看生物信息学中最著名的工具之一——BLAST算法。

BLAST用于解决数据库搜索问题：给定一个查询序列，快速在庞大的序列数据库中找出与之相关的序列。其核心策略是“先筛选，后验证”。

以下是BLAST算法的工作流程：

分词：将查询序列切分成一系列重叠的短单词（例如，长度为3个氨基酸的单词）。
构建邻域：对于每个单词，根据替换矩阵（如BLOSUM62）找出所有得分超过一定阈值的相似单词，构成一个“邻域”。这允许在匹配时容忍一定程度的错配。
哈希与搜索：预先为数据库中的所有序列建立哈希表，记录每个短单词出现的位置。然后，用查询序列所有单词的“邻域”去扫描这个哈希表，寻找命中点。
延伸与评估：对找到的命中点向两侧延伸，进行更细致的局部比对，并计算一个最终得分。只报告那些得分超过显著性阈值的匹配。

BLAST之所以快速，是因为它利用哈希实现了对数据库的常数时间查找，并且通过邻域搜索和双命中策略等技巧，在保持高灵敏度的同时最大限度地减少了不必要的详细比对计算。

📊 比对分数的概率学基础

我们一直在使用匹配、错配和空位罚分，但这些分数是如何制定的？它们代表什么？本节中我们将揭示比对分数背后的概率学解释。

比对分数本质上衡量的是“观察到的比对是由于进化相关（同源）而非随机巧合”的可能性。我们可以用对数似然比来形式化这一概念。

给定两个序列X和Y的一个比对，我们比较两个假设：

假设R（相关）：序列源于共同祖先。
假设U（不相关）：序列是独立随机产生的。

比对的概率在假设R下是 P(X,Y | R) = Π_i P(x_i, y_i)，即每个位置氨基酸对联合出现的概率的乘积。
在假设U下是 P(X,Y | U) = Π_i P(x_i) P(y_i)，即每个位置氨基酸独立出现概率的乘积。

对数似然比为：
S = log( P(X,Y | R) / P(X,Y | U) ) = Σ_i log( P(x_i, y_i) / (P(x_i)P(y_i)) )

这个公式的妙处在于：我们常用的加和式比对得分，恰好就是这个对数似然比的和。因此，替换矩阵（如BLOSUM62）中的分数 s(a, b)，实际上就是 log( P(a,b) / (P(a)P(b)) )，即观察到氨基酸a和b配对相对于随机期望的对数几率。

所以，当我们用动态规划计算出一个最优比对的总分时，这个总分直接对应了这个比对由进化关系产生的（对数）可能性相对于随机产生的可能性有多大。这使得我们可以为比对结果计算统计显著性（p值）。

✨ 总结

本节课中我们一起深入学习了序列比对和数据库搜索的多个核心主题。

我们首先比较了全局比对和局部比对的算法差异，理解了如何修改动态规划规则来寻找序列内部最相似的区域。接着，我们探索了半数值字符串匹配和哈希技术，看到了如何将字符串比较转化为数字运算，并利用哈希实现快速的近似搜索。在此基础上，我们剖析了BLAST算法，了解了它如何通过分词、邻域搜索和哈希表等技术，在庞大的生物序列数据库中实现高速且灵敏的搜索。最后，我们探讨了比对分数的概率学基础，明白了常用的替换矩阵分数实际上源于对数似然比，这使得我们的比对得分具有明确的统计意义。

这些概念将算法效率、生物学意义和统计推断紧密结合，是生物信息学分析的基础。

📚 课程 P4：L4 - 隐马尔可夫模型 (HMM) 第1部分

在本节课中，我们将学习隐马尔可夫模型及其在生物序列分析中的应用。我们将首先介绍贝叶斯推断和贝叶斯规则的基本概念，然后将这些概念与马尔可夫链的时间动态特性相结合。接着，我们会从一个简单的背景序列与CpG富集启动子序列的例子开始，逐步探索隐马尔可夫模型。我们将学习一系列用于处理隐马尔可夫模型的算法，包括如何对序列进行评分、如何推断最佳路径、如何对所有路径进行评分。之后，我们会通过增加状态空间来扩展模型，并转向后验解码和学习。本节课是两部分中的第一部分，旨在为初学者提供清晰、易懂的讲解。

🧠 贝叶斯推断简介

上一节我们介绍了课程的整体目标。本节中，我们来看看如何对世界进行推断。

生成模型允许你表达在给定世界隐藏状态下，事件发生的正向概率。我们通过感官收集观察结果，但永远无法直接了解世界的真实状态。例如，你坐在房间里，只能看到走廊外的景象，但看不到天空。你观察到有阳光、下雨或下雪，并试图推断出是什么季节。你是在根据一组观察结果，推断世界的真实状态。

贝叶斯推断允许我们逆转这个箭头。生成模型提供了在给定季节（假设）下观察到某种天气（数据）的概率。贝叶斯推断则允许我们计算在观察到特定数据（如下雨）后，某个假设（如存在大型风暴系统）为真的概率。

为了做到这一点，我们将计算假设在给定数据下的后验概率。这依赖于似然（数据在给定假设下的概率）、先验（假设在收集数据前的概率）以及数据在所有可能假设下的总概率（即数据的边际似然）。

以下是贝叶斯规则的公式：

P(H|D) = [P(D|H) * P(H)] / P(D)

其中：

P(H|D) 是后验概率（观察到数据D后，假设H为真的概率）。
P(D|H) 是似然（假设H为真时，观察到数据D的概率）。
P(H) 是先验概率（在观察数据前，假设H为真的概率）。
P(D) 是证据或边际似然（在所有可能假设下，观察到数据D的总概率）。

我们关心概率序列建模，因为生物数据是嘈杂的，概率为操作这些模型提供了一套演算方法。它不局限于“是”或“否”的答案，可以提供信念的程度。许多常见的计算工具都基于这些概率模型，而今天我们将重点讨论马尔可夫链和隐马尔可夫模型。

⛓️ 从马尔可夫链到隐马尔可夫模型

上一节我们介绍了贝叶斯推断。本节中，我们来看看如何将时间动态引入模型。

隐马尔可夫模型与简单贝叶斯推断的关键区别在于，我们处理的不是孤立的观察，而是一系列观察。例如，如果每天早上我通过长走廊的窗户观察，我可以利用从一天到另一天的信息来改进我对当前季节的后验概率判断。

隐马尔可夫模型允许你将贝叶斯推断的概念（世界的隐藏状态 vs. 观察状态）与状态间的依赖关系和转移耦合起来。这些状态间的转移由一个马尔可夫链控制。

马尔可夫链 是一个概率模型，它根据状态间的转移概率在不同状态间移动（例如，从夏季到秋季，从秋季到冬季）。它没有记忆性，下一个状态只取决于当前状态。
隐马尔可夫模型 与马尔可夫链的区别在于，在HMM中，链（隐藏状态序列）与观察序列是解耦的。隐藏的世界状态（例如，季节或风暴系统）决定了观察结果的发射概率。你看到的（如下雪）和世界的真实状态（如冬季）之间存在脱节。

🔬 隐马尔可夫模型的数学描述

现在，让我们为所有这些概念赋予数学形式。

一个隐马尔可夫模型包含以下部分：

观察序列 X：这是你可以观察到的序列（例如，阳光、多云、雨、雪）。
隐藏路径 π：这是生成每个观察结果的隐藏状态序列（例如，夏季、冬季）。
发射概率 e_k(b)：这是在隐藏状态 k 下发射出符号 b（观察结果）的概率。例如，P(x_i | π_i = k)，即在位置 i 的状态为 k 时，观察到字符 x_i 的概率。
转移概率 a_{kl}：这是从状态 k 转移到状态 l 的概率。例如，P(π_i = l | π_{i-1} = k)，即从位置 i-1 的状态 k 转移到位置 i 的状态 l 的概率。转移矩阵 A 是一个 K x K 的矩阵（K 是状态数）。发射概率可以组织成一个矩阵，其中每行对应一个状态，每列对应一个可能的观察符号。

我们将使用贝叶斯规则来估计在给定观察结果下，隐藏状态的概率。但我们的生成模型将只包含这些正向概率：在给定隐藏状态下发射每个字符的概率。

🧬 生物学示例：GC富集区域检测

让我们回到生物学领域，看一个具体例子。

如果你观察人类基因组中转录起始位点周围的GC含量（G和C核苷酸的比例），你会发现靠近转录起始位点的地方，GC含量和CpG二核苷酸频率显著富集。根据物种不同，这个信号可能不同。

我们想要建模的是：给定一段DNA序列的GC含量，它属于启动子区域的概率。我们将基因组建模为两种状态：

状态 P：启动子（Promoter）
状态 B：背景（Background）

我们将对背景区域和启动子区域的不同核苷酸组成进行建模。

在背景区域，我们可能期望四个核苷酸（A, C, G, T）的频率各为25%。
在GC富集的启动子区域，我们可能期望G和C的频率高达80%，而A和T各占10%。

此外，从对基因组的研究中，我们可能知道GC富集区域平均持续约20个核苷酸，而非启动子（背景）区域平均持续约100个核苷酸。

现在，我们将这些生物学知识编码到一个隐马尔可夫模型中：

发射概率：编码了背景区各字符25%的均等频率，以及GC富集区80%为G/C、20%为A/T的特性。
转移概率：编码了状态的平均持续时间。例如，从背景状态转移出去的概率为1%（意味着平均停留100个时间步），从启动子状态转移出去的概率为5%（意味着平均停留20个时间步）。

这个简单的两状态HMM可以用于检测GC富集区域。类似地，HMM可以用于检测保守区域（状态：保守 vs. 不保守；发射：保守水平）、蛋白质编码外显子（状态：编码 vs. 非编码；发射：三核苷酸组成或三联体）以及更复杂的基因结构或染色质状态。

📊 评估给定路径的联合概率

上一节我们构建了一个HMM模型。本节中，我们来看看如何评估一个特定的序列解析（路径）的好坏。

给定一个观察序列 X 和一个特定的隐藏路径 π，它们的联合概率 P(X, π) 是序列和该路径同时出现的概率。这可以通过生成模型计算：

P(X, π) = P(X | π) * P(π)

计算过程是：从起始状态概率开始，然后沿着路径，依次乘以从当前状态发射当前观察字符的概率，以及转移到下一个状态的概率。

让我们用之前的启动子/背景模型和一个示例序列 “ATTAGGCT” 来测试。我们比较三种解析：

全部为背景状态（BBBBBBBB）。
全部为启动子状态（PPPPPPPP）。
中间部分是启动子，两边是背景（例如，BBBBPPBB）。

通过计算每种解析的联合概率，我们发现：

全部背景的解析具有最高的联合概率（尽管绝对值很小，例如 5.2e-9），因为序列中A/T较多，更符合背景的发射概率，且背景状态本身更常见（先验高）。
全部启动子的解析概率较低（例如 2.0e-9），因为序列中的A/T在启动子模型下发射概率低。
混合解析的概率甚至更低（例如 1.6e-9），因为状态转移本身概率较低，引入两次转移带来了较大惩罚。

这个比较实际上是在计算不同假设（解析）下的似然与先验的乘积。由于观察序列 X 对于所有假设是相同的，在比较时，边际似然 P(X) 会被抵消掉。因此，我们通常直接比较 P(X|π) * P(π)。

🧭 维特比算法：寻找最佳路径

上一节我们学会了如何评估一个特定路径。本节中，我们来看看如何从指数级数量的可能路径中找到最佳的那一条。

对于一个长序列，可能路径的数量是状态数的指数级（对于长度为 n 的序列和 K 个状态，有 K^n 条路径）。逐一评估是不可行的。

解决方案是使用动态规划，具体算法称为维特比算法。其核心思想是存储中间计算结果，并利用最优子结构性质：到达当前位置 i 的状态 k 的最佳路径，必然由到达前一个位置 i-1 的某个状态 j 的最佳路径，加上从 j 到 k 的转移构成。

我们定义一个维特比变量 v_k(i)，它表示在位置 i 以状态 k 结束的所有可能路径中的最大联合概率（或对数概率）。

算法步骤如下：

初始化：对于每个状态 k，计算 v_k(1) = initial_prob(k) * emission_k(x_1)。
递归：对于位置 i = 2 到 n，对于每个状态 k，计算：
v_k(i) = emission_k(x_i) * max_{j} [ v_j(i-1) * transition_{j->k} ]
同时，记录下是哪个前驱状态 j 带来了这个最大值（用于后续回溯）。
终止：在序列末尾，最佳路径的概率是 max_{k} [ v_k(n) ]。
回溯：从获得最大概率的最终状态开始，根据记录的前驱指针，反向追踪到序列开头，即可得到最佳隐藏路径 π*。

该算法的时间复杂度为 O(K^2 * n)，空间复杂度为 O(K * n)，远优于穷举搜索。

∑ 前向算法：计算所有路径的总概率

上一节我们找到了最可能的单一路径。本节中，我们来看看如何计算生成整个观察序列的总概率，即对所有可能路径进行求和。

总概率 P(X) 表示在模型参数下，观察到序列 X 的可能性，它考虑了所有可能的隐藏路径。这对于模型比较和后续的学习算法非常重要。

计算 P(X) 同样面临路径数量指数级的问题。解决方案是另一个动态规划算法——前向算法。

我们定义一个前向变量 f_k(i)，它表示在位置 i 以状态 k 结束，并且生成了前 i 个观察字符 x_1...x_i 的所有可能路径的概率之和。

算法步骤如下：

初始化：对于每个状态 k，计算 f_k(1) = initial_prob(k) * emission_k(x_1)。
递归：对于位置 i = 2 到 n，对于每个状态 k，计算：
f_k(i) = emission_k(x_i) * sum_{j} [ f_j(i-1) * transition_{j->k} ]
注意，这里用的是求和 sum，而不是维特比算法中的取最大值 max。
终止：序列 X 的总概率是 P(X) = sum_{k} [ f_k(n) ]。

前向算法的时间复杂度也是 O(K^2 * n)。它高效地计算了所有可能路径对总概率的贡献，而不需要显式枚举每一条路径。

🎯 总结

在本节课中，我们一起学习了隐马尔可夫模型的第一部分内容。

我们首先从贝叶斯推断入手，理解了如何通过观察数据来推断世界的隐藏状态，并掌握了贝叶斯规则的核心公式。接着，我们引入了时间动态，将贝叶斯推断与马尔可夫链结合，从而得到了隐马尔可夫模型的框架。我们明确了HMM的组成部分：观察序列、隐藏路径、发射概率矩阵和转移概率矩阵。

然后，我们探讨了HMM在生物序列分析中的一个经典应用：区分GC富集的启动子序列和背景序列。我们学习了如何计算一个特定序列与隐藏路径的联合概率。

面对指数级数量的可能路径，我们介绍了两个核心的动态规划算法：

维特比算法：用于高效地找到最有可能的隐藏状态路径（解码问题）。
前向算法：用于高效地计算观察到整个序列的总概率，即对所有可能路径求和（评估问题）。

这些算法是理解和使用HMM的基石。在下节课中，我们将继续探讨HMM的更多内容，包括增加状态空间的复杂性、后验解码以及如何从数据中学习HMM的参数。

📚 课程 P5：L5 - 隐马尔可夫模型（HMM）第二部分教程

在本节课中，我们将深入学习隐马尔可夫模型（HMM）的核心算法与应用。我们将回顾HMM的基础知识，探讨如何扩展模型状态空间以处理更复杂的生物序列，并详细介绍参数学习的方法，包括监督学习和无监督学习。课程内容旨在让初学者能够理解并掌握HMM在生物信息学中的关键应用。

🔄 回顾：HMM基础与核心算法

上一节我们介绍了HMM的基本概念，包括评估、解析和后验解码。我们探讨了观测模型、贝叶斯规则以及贝叶斯推断。核心思想是，我们现在能够对序列进行建模，而不仅仅是简单的比对。

建模通过三个主要任务实现，对应我们矩阵的行：

发射特定类型的DNA序列（即从不同类别的元素生成发射）。
识别特定类型的DNA序列（即根据观测序列推断最可能产生它的隐藏状态）。
推断这些状态的概率分布（本节课将重点学习）。

我们分离了隐藏状态和观测序列。隐藏状态代表我们对世界的推断模型，观测序列则是我们实际看到的数据。我们利用贝叶斯规则来逆转这些条件概率的方向性。

贝叶斯规则公式：
P(H|D) = [P(D|H) * P(H)] / P(D)
其中，P(H|D) 是后验概率，P(D|H) 是似然，P(H) 是先验，P(D) 是证据或边际概率。

HMM将这种贝叶斯推断与顶层的马尔可夫链（控制状态转移）结合起来。我们讨论了观测向量 X、路径或解析向量 π（即隐藏状态序列），以及发射概率和转移概率。

发射概率：对于每个状态，发射每个字符的概率，即 P(字符 | 状态)。
转移概率：给定前一个位置 i-1 的状态是 k，下一个位置状态是 l 的概率，即 P(状态_l | 状态_k)。

我们看到了一个简单的HMM示例，用于基于核苷酸频率和状态间转移概率的差异来检测背景序列与GC富集启动子区域。我们计算了给定路径和观测序列的总概率，并比较了不同解析方式的可能性。

由于可能路径的数量是指数级的，直接计算不可行。因此我们定义了维特比算法，通过动态规划递归地计算结束于每个状态的最大可能概率路径。

维特比变量递归公式：
v_k(i) = e_k(x_i) * max_j [ v_j(i-1) * a_{j,k} ]
其中，v_k(i) 是在位置 i 结束于状态 k 的所有路径中的最大概率，e_k(x_i) 是状态 k 发射字符 x_i 的概率，a_{j,k} 是从状态 j 转移到状态 k 的概率。

我们还定义了前向算法，用于计算所有可能路径的总概率（即边际化所有路径）。

前向变量递归公式：
f_k(i) = e_k(x_i) * sum_j [ f_j(i-1) * a_{j,k} ]
其中，f_k(i) 是在位置 i 结束于状态 k 的所有路径的概率总和。

前向算法可以按列从左到右计算，最后对所有状态的概率求和，得到观测序列的总概率 P(X)。在实际计算中，通常使用对数空间以避免数值下溢。

🧠 扩展状态空间：增加模型“记忆”

第一个HMM示例只是一个玩具模型，只有两个状态，发射四种字符，状态间转移简单。现在，让我们讨论如何检测二核苷酸频率。既然HMM基于的马尔可夫链是无记忆的，我们如何增加系统的“记忆”呢？

解决方案是增加状态的数量。通过将更多信息编码到当前状态中，我们可以创建更多的“记忆”。例如，为了记住前一个核苷酸，我们可以将原来的“启动子”状态扩展为四个子状态：P_A, P_C, P_G, P_T，分别表示前一个字符是A、C、G、T的启动子状态。对“背景”状态也进行同样的扩展。这样，我们就得到了一个八状态的HMM。

这种扩展允许我们编码二核苷酸频率。转移概率现在可以捕获二核苷酸模式。例如，在CpG岛（启动子）模型中，C后接G的频率很高，而在非CpG岛（背景）区域，这种频率则低得多。

为什么CpG岛重要？
CpG二核苷酸中的胞嘧啶（C）可以被甲基化，这是一种常见的基因表达抑制标记。在进化过程中，甲基化的C容易脱氨基变成T，因此基因组中大部分CpG位点会逐渐消失。然而，启动子等区域通常保持低甲基化状态，从而保留了CpG岛。因此，寻找CpG岛有助于定位可能的启动子区域。

以下是扩展状态空间的其他应用示例：

蛋白质编码基因检测：需要更复杂的结构，包括起始密码子（ATG）、外显子、内含子、剪接位点（GT-AG）、终止密码子（TGA, TAG, TAA）以及5‘和3’非翻译区（UTR）。关键在于必须记住阅读框，因为翻译是按三个核苷酸（一个密码子）进行的。从一个外显子过渡到下一个外显子时，必须保持阅读框的连续性。这通常通过创建多个外显子状态（如E0, E1, E2）来编码阅读框的偏移量来实现。
染色质状态检测（如ChromHMM）：使用HMM根据组蛋白修饰的组合模式来推断不同的染色质状态（如增强子、启动子、转录区域、抑制区域）。每个隐藏状态发射一个染色质标记的向量，转移矩阵捕获这些状态的空间关系。
蛋白质编码保守性检测：使用两状态HMM，基于密码子替换频率的差异来区分编码和非编码外显子。

所有这些应用都展示了通过增加状态数量，HMM能够捕获生物序列中复杂而多样的模式。

🎯 后验解码

之前我们讨论了维特比解码，它找到最可能的单一路径（最大似然路径）。后验解码则提供了另一种方式：它找到在每个位置上最可能的状态，这个判断考虑了所有可能路径的概率。

后验解码定义：
对于序列的每个位置 i，选择状态 k，使得在所有路径中，位置 i 处于状态 k 的总概率最大。即：
π_i = argmax_k P(π_i = k | X)

如何计算后验概率 P(π_i = k | X)？
这需要计算通过位置 i 状态 k 的所有路径的总概率。我们可以利用前向-后向算法：

前向算法：计算从序列开头到位置 i 并结束于状态 k 的所有路径的概率 f_k(i)。
后向算法：计算从位置 i 的状态 k 开始，到序列末尾的所有路径的概率 b_k(i)。
那么，通过位置 i 状态 k 的所有路径的概率为 f_k(i) * b_k(i)。
后验概率为：P(π_i = k | X) = [f_k(i) * b_k(i)] / P(X)，其中 P(X) 可由前向或后向算法在序列末端/始端得到。

后向算法递归公式（从右向左计算）：
b_k(i) = sum_l [ a_{k,l} * e_l(x_{i+1}) * b_l(i+1) ]

后验解码的优点：

比维特比路径提供更精细的、每个位置的信度度量。
对于分类任务非常有用，能更好地反映状态的不确定性。

后验解码的缺点：

得到的状态序列可能不是一个有效的路径（即某些相邻状态之间的转移概率可能为零或极低，不符合模型约束）。而维特比路径始终是一个连贯的、符合模型转移规则的有效路径。

📈 参数学习：监督与无监督

现在我们来探讨如何学习HMM的参数（发射概率和转移概率）。这对应我们最初的第三个任务：学习每个状态的区分性特征。

监督学习

如果我们有已标注的序列（即已知每个位置对应的隐藏状态），那么参数学习非常简单，只需计数即可。

最大似然估计方法：

发射概率：e_k(b) = （在状态 k 下观察到字符 b 的次数） / （处于状态 k 的总次数）。
转移概率：a_{k,l} = （从状态 k 转移到状态 l 的次数） / （离开状态 k 的总次数）。

为了避免零概率（在取对数时会导致无穷大的惩罚），通常加入伪计数（加性平滑），这相当于引入了一个弱的先验分布。

无监督学习（鲍姆-韦尔奇算法）

更常见的情况是，我们只有观测序列，没有标注。我们需要同时推断隐藏状态序列和模型参数。这可以通过期望最大化（EM）算法来实现，具体到HMM就是鲍姆-韦尔奇算法。

基本思想（迭代优化）：

初始化：随机猜测一组模型参数（或基于先验知识初始化）。
E步（期望步）：使用当前参数，计算所有位置处于每个状态的后验概率 P(π_i = k | X)（前向-后向算法），以及所有相邻位置状态对的后验概率 P(π_i = k, π_{i+1} = l | X)。这相当于对隐藏状态进行“软”赋值或概率估计。
M步（最大化步）：利用E步计算出的后验概率（作为“软计数”），重新估计模型参数，类似于监督学习中的计数，但现在是加权计数。
- 发射概率：e_k(b) ∝ 对所有 x_i = b 的位置 i 求和 P(π_i = k | X)。
- 转移概率：a_{k,l} ∝ 对所有位置 i 求和 P(π_i = k, π_{i+1} = l | X)。
重复：用新的参数替换旧参数，回到E步。迭代直到参数收敛（变化很小）。

为什么EM算法有效？
即使从随机参数开始，观测序列本身的结构也会引导E步中的后验概率计算。在M步中，这些后验概率会更新参数，使其更符合序列中实际存在的模式（如不同的核苷酸组成区域）。每次迭代，模型的对数似然值都会增加（或保持不变），最终收敛到一个局部最优解。

简化版本：维特比训练
一种更简单的无监督学习方法是维特比训练，它是EM算法的一个近似：

E步：使用当前参数运行维特比算法，找到最可能的单一路径（硬赋值）。
M步：像监督学习一样，在这条路径上进行计数，更新参数。
重复。
这种方法计算更快，但不如使用所有路径信息的完整鲍姆-韦尔奇算法精确。

✅ 总结

在本节课中，我们一起深入学习了隐马尔可夫模型（HMM）的核心内容：

回顾了HMM基础：包括贝叶斯推断、维特比算法（用于寻找最可能路径）和前向算法（用于计算序列总概率）。
探讨了状态空间扩展：通过增加状态数量，HMM可以编码更复杂的“记忆”，例如二核苷酸频率、蛋白质编码基因的阅读框信息等，从而应用于CpG岛检测、基因预测、染色质状态分析等多个生物信息学领域。
学习了后验解码：这是一种不同于维特比解码的方法，它找出每个位置上最可能的状态，考虑了所有路径的信息，提供了更细致的概率解释。
掌握了参数学习：
- 监督学习：在已知状态序列时，直接通过计数获得参数。
- 无监督学习（鲍姆-韦尔奇算法）：在只有观测序列时，使用期望最大化（EM）算法迭代地估计状态后验概率并更新模型参数，直至收敛。

通过掌握这些概念和算法，你已经具备了使用HMM对复杂生物序列进行建模、解析和学习的基本能力。

🧬 P6：L6- 基因表达分析：聚类与分类

在本节课中，我们将学习基因表达分析的基础知识。我们将探讨两种核心技术——微阵列和RNA测序，并深入理解两种核心的计算方法：无监督学习（聚类）和有监督学习（分类）。课程将涵盖K均值聚类、层次聚类以及朴素贝叶斯分类等核心算法。

🧬 基因表达分析基础

基因表达分析是生物学和计算机科学的基础。在完成基因组注释模块后，我们现在进入模块二，该模块关注基因表达的动态变化。模块一讨论了静态基因组，模块二将关注转录本动态和表观基因组修饰。后续模块将涉及基因调控网络、疾病基因组学和进化基因组学。

模块二的计算基础将包括有监督和无监督学习、快速读段比对以及RNA的二维结构。生物学前沿将涵盖基因表达分析、转录本结构、表观基因组学和三维基因组学。

本节课的目标是：

介绍基因表达分析的基础。
讲解两大技术支柱：微阵列和RNA测序。
探讨有监督和无监督分析的计算方法。

在介绍之后，我们将讨论K均值聚类和层次聚类，然后是分类的生成式方法。如果时间允许，我们还将探讨分类的判别式方法。

🔬 基因表达分析技术

基因表达谱分析的两大技术支柱是微阵列技术和下一代测序技术。

微阵列技术 的工作原理是：针对基因组中的一组固定基因，为每个基因制作探针。通过逆转录将RNA转化为互补DNA（cDNA），然后将这些cDNA与微阵列上的探针进行杂交。通过DNA-DNA杂交，标记了荧光染料的cDNA会与互补探针结合，通过测量每个DNA点的荧光强度，可以推断每个基因的表达水平。

探针可以是每个基因的一个长探针，也可以是覆盖基因的多个短探针（平铺探针）。微阵列技术的优势在于，即使每个细胞的分子数非常少，它也能让你专注于感兴趣的特定小区域。

RNA测序技术 是过去十年的革命性进展。它直接对成熟的mRNA分子进行测序，然后使用我们之前讨论过的技术（以及即将讨论的更新技术）将读段映射回基因组。通过这种映射方法，可以计算映射到基因组中每个已知基因的读段数量，从而推断基因的表达量。

与微阵列的模拟信号（荧光强度）不同，RNA测序提供的是数字化的分子计数。需要进行各种标准化，例如根据基因长度和总测序读段数进行校正（如RPKM）。这种方法的优势在于无需预先知道感兴趣的特定基因，可以进行无偏的mRNA测序实验，从中推断所有已知基因的计数，甚至可能发现新的高表达基因。

📊 数据分析的两个维度

完成测量后，你通常拥有成千上万个基因在数百个条件下的表达数据。这形成了一个数据矩阵。

你可以从两个维度分析数据：

基因相似性：将每个基因在所有条件下的表达水平视为一个向量。通过比较这些基因向量，可以询问哪些基因在表达模式上行为相似。
条件相似性：将每个条件下所有基因的表达水平视为一个向量。通过比较这些条件向量，可以询问哪些实验条件最为相似。例如，比较精神分裂症与某种医学干预后的表达谱。

这两个维度为研究基因表达数据提供了不同的视角。

🧩 无监督学习：聚类

聚类意味着将相似的项目分组，这些项目可能来自同一类别，从而揭示数据的隐藏结构。我们可以对基因进行聚类，也可以对实验条件进行聚类。

聚类不依赖于任何先验注释，它是一种无监督的方法，旨在寻找相似性模式。聚类结果可以通过后续分析进行验证，例如检查聚类出的基因是否富集了特定的生物学功能。

实验条件可以是不同的组织（如肝脏、大脑）、不同的生理状态（如餐后、住院）、疾病状态或年龄等。聚类完成后，我们可以检查每个簇中的基因是否在特定功能类别上显著富集。

🏷️ 有监督学习：分类

与聚类不同，分类属于有监督学习。我们预先知道一些类别（例如，已知某些基因参与胰腺β细胞功能）。分类的目标是从数据中提取特征（即特定条件下的表达水平），这些特征能最好地将新元素分配到一个或多个明确定义的类别中。

在有监督学习中，类别是预先定义好的，算法学习基因的表达模式如何映射到这些类别。而在无监督学习中，我们只在聚类完成后才用生物学数据验证结果。

分类的一个子问题是特征选择，即从高维向量中识别出最重要的特征。评估分类效果的指标包括准确率、敏感性、特异性等。

📐 聚类与分类的形式化描述

聚类和分类是机器学习和计算机科学中普遍使用的方法，用于通过一个或多个特征来表征对象。

在我们的例子中，对象是基因，特征可以是基因在不同组织或数百万个细胞中的表达水平。每个基因点都存在于一个超高维空间中。

分类是指我们的数据点带有标签（如红色或绿色），我们希望找到一个“规则”，能够准确地将新数据点分配标签。特征选择是分类的一个子问题。

聚类是无监督的，没有标签。它根据数据点之间的邻近程度将其分组，目标是识别数据中的结构。评估指标可以是某种独立的验证特征，例如聚类后检查基因的功能富集情况。

🔢 聚类方法一：K均值聚类（划分式）

我们将探讨两种聚类方法：划分式聚类和凝聚式聚类。

K均值聚类是一种划分式方法，它将对象划分为K个不重叠的簇，每个数据对象只属于一个簇。其基本思想是：假设存在固定数量（K）的簇，目标是将点划分到这些簇中，并使簇内紧凑。

算法步骤如下：

随机初始化：随机选择K个点作为初始簇中心（质心）。
分配点：将每个数据点分配给距离它最近的簇中心。
更新质心：将每个簇的中心移动到属于该簇的所有点的中心（平均值）。
迭代：重复步骤2和3，直到分配不再发生变化（收敛）。

K均值算法试图最小化一个成本标准，即所有点到其所属簇质心的距离平方和。这实质上是在使每个簇尽可能紧凑。

模糊K均值 是K均值的一个变体。它不进行硬分配（一个点100%属于一个簇），而是进行软分配。每个点以一定的概率属于所有簇，概率大小取决于该点到各簇中心的距离。然后，质心更新为所有点的加权平均值，权重就是各点属于该簇的概率。标准的K均值是模糊K均值在概率为100%或0%时的特例。

🤖 K均值的机器学习公式：生成式模型

我们可以从生成式模型的角度来看待K均值聚类。这被称为高斯混合模型。

生成过程如下：首先抛一枚硬币来选择使用哪个高斯分布（例如，以30%的概率选择蓝色分布，以70%的概率选择红色分布），然后从选定的高斯分布中采样生成数据点。

在实际中，我们观察到数据点（灰色），但不知道每个点是由哪个高斯分布生成的（这是隐变量）。我们的目标是：仅给定样本，估计模型参数（两个高斯分布的均值和方差）以及点的分配标签。

这可以通过期望最大化算法来解决：

E步（期望步）：给定当前估计的簇中心，计算每个数据点属于每个簇的概率（即“期望”的分配）。
M步（最大化步）：根据E步计算出的概率权重，重新计算每个簇的中心（即“最大化”似然性的参数）。

迭代进行E步和M步，直到收敛。可以证明，在高斯分布方差相等且为先验知识的情况下，EM算法的解与K均值算法的解是等价的。

模糊K均值对应于使用所有可能路径（全路径）的EM算法，而标准K均值对应于只选择最可能路径（维特比路径）的EM算法。此外，还可以使用吉布斯采样从后验分布中采样单个标签。

🌳 聚类方法二：层次聚类（凝聚式）

K均值需要预先指定簇的数量K，而层次聚类则不需要。层次聚类会产生一组嵌套的簇，组织成树状结构（树状图）。

算法步骤如下（凝聚式）：

开始：将每个点视为一个单独的簇。
合并最接近的簇：在每一步，找到并合并两个最接近的簇。
重复：重复步骤2，直到所有点都合并成一个簇。

结果是一个树状图。要获得离散的簇，只需在树的特定高度“切割”即可。切割的位置决定了簇的数量。

关键问题是如何定义簇与簇之间的距离，常见方法有：

单链法：两个簇中最近的两个点之间的距离。
全链法：两个簇中最远的两个点之间的距离。
平均链法：两个簇中所有点对之间的平均距离。
质心法：两个簇的质心之间的距离。

不同的距离度量会影响聚类结果和运行时间。同样，点与点之间的距离（如欧氏距离、曼哈顿距离、皮尔逊相关系数等）也会影响最终结果。

📈 聚类结果评估

如何评估聚类结果的好坏？

紧凑性：可以测量簇内的紧凑程度。
鲁棒性：通过从数据中随机抽样并重新聚类，检查簇是否重复出现。
富集分析：使用独立的注释信息（如基因功能类别）来评估。对于一个聚类结果，可以计算某个功能类别在特定簇中是否显著富集。这通常使用超几何分布来计算p值，以评估观察到的富集程度是否随机。

🎯 分类方法：朴素贝叶斯分类器

分类方法可分为生成式方法和判别式方法。朴素贝叶斯分类器是一种常见的生成式方法。

生成式方法对每个类别的特征分布进行建模。例如，假设我们有两个类别（红色和绿色），每个类别下，某个特征（如基因表达水平）服从一个高斯分布（类条件分布）。

分类过程：对于一个新数据点，我们使用贝叶斯规则计算它属于每个类别的后验概率。后验概率正比于 （类条件概率 × 类先验概率）。选择后验概率较高的类别作为预测结果。

决策边界可以表示为一个判别函数。当判别函数值大于0时，选择类别1；否则选择类别2。

训练过程：使用已标记的数据（训练集）来估计类条件分布的参数（如均值和方差）以及每个类别的先验概率。

“朴素”的假设：为了简化计算，特别是在特征维度很高时，朴素贝叶斯假设在给定类别的条件下，各个特征之间是相互独立的。因此，联合的类条件概率可以分解为每个特征单独的概率的乘积。这使得计算变得非常高效。

尽管这个假设很强，但在许多实际应用中，朴素贝叶斯分类器表现得出奇地好。评估分类性能可以使用准确率、敏感性、特异性等指标。

🎓 总结

在本节课中，我们一起学习了基因表达分析的核心内容。

我们首先介绍了基因表达分析的技术基础：微阵列和RNA测序。然后，我们探讨了分析基因表达数据的两个计算范式：无监督学习（聚类） 和 有监督学习（分类）。

在聚类部分，我们深入讲解了：

K均值聚类：一种划分式方法，通过迭代优化使簇内紧凑。我们看到了它的算法步骤、最优性准则，以及其模糊变体和生成式模型（高斯混合模型）解释。
层次聚类：一种凝聚式方法，无需预先指定簇数，通过构建树状图来揭示数据的层次结构。我们讨论了不同的簇间距离度量方法。

在分类部分，我们重点介绍了：

朴素贝叶斯分类器：一种生成式方法，基于贝叶斯定理和特征条件独立性假设，通过建模类条件概率分布来进行分类。

通过本节课的学习，你应该对如何利用计算工具从基因表达数据中提取生物学见解有了基本的认识。这些聚类和分类方法是生物信息学数据分析的基石。

🧬 P7：L7 - RNA折叠、RNA世界与RNA结构

在本节课中，我们将学习RNA生物学、RNA转录本重建以及RNA结构预测。我们将探讨RNA在生命起源中的核心作用，了解如何从测序数据中重建RNA转录本，并学习预测RNA二级结构的计算方法。

📚 RNA生物学：为什么RNA如此重要

上一节我们介绍了基因表达分析，本节中我们来看看RNA在分子生物学中的核心作用。传统观点认为，RNA仅仅是DNA和蛋白质之间的信使，但现代研究发现，RNA的功能远不止于此。

以下是RNA扮演的多种角色：

信使RNA (mRNA)：携带遗传信息，用于蛋白质合成。
核糖体RNA (rRNA)：构成核糖体的催化核心，负责翻译过程。
转运RNA (tRNA)：作为适配器分子，在翻译过程中携带特定氨基酸。
参与剪接的RNA：如核小RNA (snRNA)，在mRNA前体的剪接中起关键作用。
调控RNA：如微RNA (miRNA) 和长链非编码RNA (lncRNA)，参与转录后基因表达调控。
结构RNA：作为支架，参与染色质组织等过程。
核酶：某些RNA分子具有酶活性，可以催化化学反应，例如自我剪接的内含子。

RNA甚至可以作为某些病毒（如SARS-CoV-2）的遗传信息存储载体。这表明，RNA能够扮演现代生物学“中心法则”中DNA（信息存储）、RNA（信息传递）和蛋白质（功能执行）的全部三种角色。

🌍 RNA世界假说

基于RNA的多功能性，科学家提出了“RNA世界”假说。该假说认为，在生命进化早期，存在一个以RNA为主导的世界。

在这个假想的世界中：

RNA是核心：RNA同时承担信息存储、自我复制和催化化学反应的功能。
自我复制与进化：能够自我复制的RNA分子通过自然选择不断进化，复杂性增加。
向现代世界过渡：RNA后来“发明”了蛋白质来分担催化工作，并最终将信息存储的任务交给了更稳定的DNA，形成了我们今天所知的“DNA → RNA → 蛋白质”的中心法则。

现代细胞中保留的多种RNA功能（如rRNA的催化作用）被认为是“RNA世界”留下的遗迹。

🧬 从测序数据重建RNA转录本

接下来，我们探讨如何利用现代测序技术来识别和量化RNA转录本。核心挑战在于，高通量测序（如RNA-seq）产生的读段很短（约30-70个核苷酸），而转录本可能很长且存在可变剪接。

方法一：基于基因组的比对

该方法需要一个参考基因组作为“地图”。

基本流程如下：

短读段比对：首先使用超快速比对工具（如Bowtie），将所有能连续比对到基因组上的读段进行定位。
识别剪接位点：对于无法连续比对的读段，其两端可能分别比对到基因组上相距很远的位置，这提示存在剪接事件。使用TopHat等工具来识别这些潜在的剪接位点。
转录本组装与定量：将比对上的读段连接起来，重建出完整的转录本（Cufflinks）。由于一个读段可能来自多个不同的剪接变体，需要使用期望最大化（EM）等算法，概率性地分配读段，并估算每个转录本的表达量。表达量通常用 FPKM（每千碱基每百万片段数）进行标准化。
差异表达分析：在量化基础上，使用统计模型（如Cuffdiff）比较不同样本间的基因表达差异，考虑技术变异和生物学变异。

方法二：从头组装

当没有参考基因组时（例如研究非模式生物），可以采用从头组装方法。

基本思路如下：

构建重叠图：直接将所有测序读段相互比较，构建一个德布鲁因图。该图能展示读段之间的重叠关系，并将序列差异（如测序错误、单核苷酸多态性、可变剪接）表示为图中的不同路径。
遍历路径推断转录本：通过遍历图中所有可能的路径，可以推导出不同的转录本异构体。
定量与比对：估算各转录本的丰度，并可将组装出的转录本序列比对到近缘物种的基因组上进行分析。

🧮 RNA二级结构预测

了解了RNA的序列和表达后，我们进一步探索其功能的关键——RNA结构。RNA分子可以通过碱基配对形成复杂的二级和三级结构，这些结构决定了其功能。

我们主要关注二级结构的预测，因为它决定了大部分折叠自由能，且计算上更易处理。二级结构不允许假结（即碱基配对交叉）。

纽辛格算法：最小自由能折叠

最经典的RNA二级结构预测算法是纽辛格算法，其目标是找到自由能最低的二级结构。该算法基于动态规划思想。

核心步骤如下：

能量模型：为不同的碱基对（如A-U， C-G， G-U）以及结构单元（如茎环、凸环、内环、多分支环）分配自由能值。这些值来自实验测量。
动态规划矩阵：定义一个二维矩阵 E(i, j)，表示从序列位置 i 到 j 的子序列所能形成的最佳二级结构的最小自由能。
递归计算：从短序列片段开始（矩阵对角线附近），逐步向外计算更长的片段。对于子序列 (i, j)，其最小自由能 E(i, j) 可以通过考虑以下几种情况递归计算：
- i 和 j 配对，则能量为 E(i+1, j-1) 加上配对 (i, j) 的能量。
- i 不参与配对，则能量为 E(i+1, j)。
- j 不参与配对，则能量为 E(i, j-1)。
- 在 i 和 j 之间某个位置 k 将结构分成两部分，则能量为 E(i, k) + E(k+1, j)。
- 计算 (i, j) 形成各种环结构（如发夹环、内环、凸环）的能量。
回溯：填充完整个矩阵后，从 E(1, n) 开始回溯，即可得到自由能最低的二级结构。

从单一结构到结构集合：分区函数

纽辛格算法只给出一个最优结构。然而，RNA在溶液中可能以多种构象存在。我们可以计算分区函数 Z，即所有可能结构的玻尔兹曼权重之和。通过计算 Z，我们可以得到每个碱基对的配对概率，从而了解结构的稳定性和动态性。

基于语法的预测：随机上下文无关文法

另一种强大的方法是使用随机上下文无关文法来建模和解析RNA二级结构。

基本概念如下：

产生式规则：定义一套文法规则来描述如何“生成”一个RNA结构。例如，一条规则可以是：S → a S u，表示在结构S的两端添加一对配对的碱基a和u。
随机性：为每条产生式规则赋予一个概率。
解析与推断：给定一个RNA序列，可以使用CYK算法等动态规划方法，找到概率最高的结构解析（即最可能生成该序列的结构），或计算所有可能解析的总概率（分区函数）。
系统发育SCFG：可以进一步扩展，结合多个同源RNA序列的进化信息，提高结构预测的准确性。

📝 总结

本节课中我们一起学习了：

RNA的多功能性：RNA远不止是信使，它在催化、调控、结构等方面扮演着核心角色。
RNA世界假说：生命可能起源于一个以RNA为核心的自我复制和催化系统。
RNA转录本重建：我们学习了如何从短读段测序数据中，通过“比对到基因组”或“从头组装”两种策略来识别和量化RNA转录本。
RNA二级结构预测：我们探讨了使用纽辛格算法基于最小自由能模型预测RNA结构，以及使用随机上下文无关文法进行基于概率模型的结构解析。这些方法是理解RNA功能机制的基础。

通过本课的学习，我们看到了RNA从生命起源到现代细胞复杂调控中的核心地位，以及计算生物学在解析其序列、表达和结构中的强大作用。

课程P8：L8- 表观基因组学L1部分 🧬

在本节课中，我们将要学习表观基因组学的基础知识。表观基因组学是研究在DNA序列之上发生的、可遗传的化学修饰如何调控基因表达的领域。我们将从介绍该领域开始，逐步深入到数据生成、处理和分析的核心计算方法。

引言：什么是表观基因组学？🧬

大家好，今天我们将讨论表观基因组学。基因组学领域存在许多挑战，而表观基因组学是我最喜欢的领域之一。

这是模块二关于基因表达和表观基因组学的第三讲，内容介于基因组注释和即将开始的调控基因组学、网络与疾病以及比较基因组学之间。

今天，我们将聚焦于表观基因组学的计算方面。首先，我将介绍这个领域：表观基因组学是关于什么的？染色质修饰的多样性是怎样的？然后，我们将探讨如何找到这些修饰的位置，即如何使用抗体和染色质免疫沉淀结合下一代测序技术。

接着，我们会介绍数据生成项目、原始数据以及初级数据处理：如何比对测序读段、如何识别峰以及如何进行质量控制。之后，我们将重点介绍整合多种表观基因组修饰以推断染色质状态的方法，并定义一个用于此目的的多变量隐马尔可夫模型来注释这些状态。

很多人问过关于模型复杂性的问题：如何选择染色质状态的数量以及标记的数量？如何捕捉状态间的依赖性和条件标记依赖性？如何跨多种细胞类型进行联合学习？以及如何利用这些染色质状态的活动模式来进行我们下周将要讨论的调控基因组学研究，特别是将增强子与其上游调控因子和下游靶基因联系起来。

最后，我们将讨论如何通过利用染色质标记之间的相关性来进行表观基因组推算。

让我们开始吧。

表观基因组学简介 🧬

表观基因组是为什么你可以在几乎所有细胞中共享一个几乎相同的、静态的基因组，却拥有如此惊人的细胞类型多样性的原因。

如果你观察血细胞类型，有大量不同类型的免疫细胞、携氧细胞、凝血细胞等等。我们最近通过新冠疫情新闻学到了很多免疫学知识。大脑拥有令人难以置信的细胞类型多样性，有主要的神经元、星形胶质细胞、少突胶质细胞、小胶质细胞等类型，还有早期大脑和人类新皮层中大量的组织多样性。

如果你观察我们所有的内部器官，如心脏、肝脏、脾脏等，同样存在令人难以置信的细胞类型和腔室多样性。即使只看你的皮肤，也有令人难以置信的神经支配，神经元直接与你的皮肤、毛发、血液供应、神经支配以及其他环境相互作用。如果你有伤口，它会愈合，所有相同的细胞类型会重新构建。这是一个相互连接的系统，不仅仅是细胞类型，还包括通讯、供给、反馈和控制的组织网络，这些在你身体的任何部分都有所体现。

所有这些都是通过表观基因组实现的。那么，这是如何实现的呢？

表观基因组本质上是将你的DNA维系在一起的结构，它首先扮演结构角色，但也扮演功能角色，并使得这种细胞类型的多样性成为可能。为了让你了解每个细胞中DNA的包装程度，你大约有两米长的DNA被包装在一个尺寸小六个数量级的细胞里。而DNA链本身的直径又比其长度小九个数量级。这是一个惊人的包装壮举，将两米长的DNA包装进你万亿个细胞中的每一个。如果你将你30万亿个细胞首尾相连，其DNA总长度足以到达木星10次。地球上的每个人体内都包装了如此多的DNA。

这就是表观基因组的第一个功能：它使得所有这些DNA能够组合在一起。这是如何可能的呢？通过许多不同层次的包装。

首先，DNA双螺旋缠绕在这些核小体“碗”周围。每个核小体包含147个核苷酸，加上大约50个核苷酸的连接区。你可以将基因组中200个核苷酸的片段想象成包装在这些“碗”里。这些“碗”本身又被包装在染色质纤维中，而染色质纤维本身又在更高的组织层次上被包装成A区和B区，这些区域更活跃或更受抑制。我们将在下一讲详细讨论这些。

所有这些最终凝聚成我们所熟知的染色体。顺便说一下，你通常看到的染色体图片只出现在细胞周期的特定阶段，即染色体排列准备在子细胞间分离时。但在大多数时候，染色体的包装程度要低得多，并且是漂浮的，但仍然保持它们自己的结构域。

再次强调，DNA非常长，细胞非常小，压缩程度比伸展的DNA小几个数量级。为了使用DNA，这种紧凑的结构必须在局部打开，以便转录因子结合。染色质状态的可及性和三维相互作用在实现这一切中扮演着巨大角色。

当我们谈论表观基因组学时，“epi”意味着“之上”，“genomics”意味着“基因组之上”，所以表观基因组学是研究在遗传信息之上发生的修饰。因此，表观基因组修饰是可以施加在人类基因组主要信息之上的修饰。

主要DNA信息可以被调节和修饰。我喜欢给出的类比是：就像乐谱上的注释可以让你演奏得更响亮或更柔和、增加或减少强度等。这些修饰可以通过三种类型的修饰来实现，这些修饰涉及DNA本身以及DNA缠绕的蛋白质。

第一种修饰是在DNA本身上，我们之前简要提到过，即CpG甲基化。“CpG”是因为它是C后面跟着磷酸基团，再跟着G。在连续的CG核苷酸中，你基本上有一个甲基化的C，这通常在调控区域指示抑制。如果你有这个“第五个碱基”——甲基化的C碱基，该区域更可能被抑制。事实上，人类基因组的大部分是甲基化的，只有一小部分在我们之前讨论过的CpG岛中是非甲基化的。

所以，这些甲基化是表观基因组修饰的第一种形式，它直接发生在DNA上。

第二种形式是DNA可及性。这基本上是DNA在哪些区域实际上是可以接近的，即没有紧密包装在核小体周围。这种可及性基本上告诉你调控因子在哪里结合。如果你有一个转录因子，它通常会结合DNA的可及区域。如果你进行可及性测定（例如，切割DNA的可及区域，然后读取哪些区域被切割，或者以某种方式标记这些开放的可及区域，然后读取哪些区域被标记），你会看到这些区域被核小体大小的不可及区域隔开。在这些可及区域内，你实际上会以超高分辨率看到单个调控因子的足迹，这些足迹有时一次覆盖10个核苷酸。围绕这些足迹，你会看到可及性，然后是缺乏可及性，等等。

因此，可及性有两个层次：第一个是在峰的水平，但当你放大到峰内部时，你会注意到峰实际上具有这种双峰形状，其中转录因子结合在那里。在某些情况下，你会看到可及区域的多个凹陷，这表明这是一个普遍可及的区域，但在这里和那里结合了两个不同的调控因子。

所以，可及性可以在多个分辨率层次上考虑。这是第二种类型的修饰：第一种是DNA甲基化，第二种是DNA可及性。

第三种是组蛋白修饰。这些表观遗传因子基本上可以修饰DNA本身的解读方式。这是迄今为止最多样化和最丰富的表观基因组量化形式。它们是如何发生的呢？它们基本上发生在构成核小体的组蛋白的尾部。

我介绍了很多术语：DNA缠绕两圈，这就是我之前谈到的147个碱基对。每一个核小体由八种蛋白质组成，这些蛋白质被称为组蛋白。这些组蛋白通常是：H2A两份、H2B两份、H3两份和H4两份。但是这些组蛋白有很多变体，它们可以在基因组特定位置交换，形成修饰的核小体和核小体变体。

就表观基因组信息而言，在基因组信息之上，你有DNA甲基化、DNA可及性和组蛋白修饰。这些修饰是什么？

这些组蛋白尾部有长的氨基酸尾巴，可以伸出并被甲基化、泛素化、乙酰化、磷酸化等修饰。这些修饰可以被专门的蛋白质识别，然后通过标记在细胞类型中重要的区域来影响该DNA区域的解读。

如果你回到人体不同的细胞类型，也许神经元会有一组特定的区域被表观基因组标记为活跃，而相同的区域在例如小胶质细胞或星形胶质细胞中被标记为非活跃。

这就是通过这种表观基因组记忆实现的。当我说它们调节时，它们基本上可以直接物理地使其不可接近，或者改变“序列”，使得通常识别非甲基化C的转录因子无法再识别甲基化的C，或者通过招募其他因子来竞争掉通常结合的调控因子，或者通过自身使染色质更凝聚或更松散，等等。

因此，你可以将基因调控想象为转录因子（直接结合DNA的调控因子，这些是序列特异性因子，结合我们将在两讲后讨论的序列基序）、组蛋白（对DNA有固有亲和力）、它们构成的核小体、染色质调控因子等之间不断的竞争。

我们将基于特定的蛋白质（如H2A、H2B、H3或H4）来命名这些修饰。例如，这是H3，然后是从末端计数的特定氨基酸残基，所以赖氨酸4是这个，我们称之为K4，或者赖氨酸36更靠下。然后是特定的化学修饰，例如是甲基化、磷酸化、乙酰化等。有时你可以在同一个位点有多个甲基基团，我们称之为三甲基化。

所以，我们将讨论H3，即组蛋白H3，赖氨酸4三甲基化，简称为H3K4me3。或者H2B赖氨酸5乙酰化，称为H2BK5ac。除了DNA修饰（CpG上的甲基化C）、核小体定位和DNA可及性之外，这些都是组蛋白修饰。

总的来说，表观基因组的标记方式是通过DNA甲基化、组蛋白修饰、DNA可及性共同协作，用不同的“颜色”标记DNA。这些“颜色”对应于在感兴趣的细胞类型中活跃的增强子、活跃的启动子、转录区域、抑制区域、重复区域等。

已知有数百种这样的修饰，还有许多新的修饰正在被发现。我们可以使用染色质免疫沉淀、亚硫酸氢盐测序和DNA可及性测定来系统地绘制它们。

DNA可及性我之前已经提到过，基本上是找出哪些区域实际上在物理上可以被特定的酶（如DNA酶I，一种切割可及DNA的切割酶）接近。所以，如果DNA被切割，那么你就知道它是可及的。这是其中一种方法。

第二种方法是，你通过化学方式处理DNA，以不同的方式修饰甲基化C和非甲基化C核苷酸，查看差异，并基于差异推断哪些区域实际上是甲基化的与非甲基化的。为了获取DNA甲基化信息，你基本上进行亚硫酸氢盐测序，这是一种化学修饰然后测序的方法。

然后是染色质免疫沉淀，我们将会讨论，但这是使用抗体来下拉具有特定修饰的基因组区域。然后你切割DNA，下拉这些区域，对你得到的东西进行测序，然后查看它们落在基因组的哪个位置。

这三种方法（以及许多其他方法）被用来系统地绘制数百个样本，涵盖成人和人类胚胎身体的不同区域、大脑的不同区域、不同类型的原代细胞、体外分化细胞等。

其中一个例子是ENCODE项目（DNA元件百科全书），另一个是Roadmap表观基因组学项目（这是NIH路线图计划的一部分），还有modENCODE、mouseENCODE等。我们的研究小组在这个领域非常活跃，我们帮助领导了其中一些联盟，当然是在计算方面，整合所有这些数据集以推断人类表观基因组图谱。

今天的许多内容实际上与我的研究密切相关，因为许多方法是由我小组的学生和博士后开发的，然后广泛应用于表观基因组学领域。

Roadmap项目绘制了超过100种不同的原代组织和细胞，并绘制了多种组蛋白修饰、DNA可及性、DNA甲基化和基因表达数据，然后整合所有这些数据以产生人类表观基因组图谱。

你会注意到我在这里使用这些颜色。这些颜色是我们发明并在许多不同项目中重复使用的，所以它们有点标准化。这些颜色基本上代表：远端基因调控区域（极其细胞类型特异性，它们环化到启动子附近以帮助招募多种因子类型，当需要时）——这些是增强子，因为它们增强表达；启动子（近端调控区域，这是RNA聚合酶最终结合以启动转录的地方）——它们同样由组蛋白修饰、DNA甲基化和可及性的组合标记；转录区域我们将用绿色表示，启动子用红色，增强子用橙色，抑制区域用灰色，等等。

这里我展示了每种类型区域中一些普遍的标记，但再次强调，有数十种这样的修饰。

这个项目的目标是绘制多种修饰、多种细胞类型、多个个体、多个物种、多种条件、使用多种抗体、跨越整个基因组的数据。已经发布了多波这样的数据，我们最近刚刚发布了ENCODE的最新一波数据，这在Nature杂志上有多篇论文报道。

你们很幸运，因为你们拥有大量在我还是学生时根本不可用的数据。

这是表观基因组学的基本概述：为什么表观基因组有用？染色质修饰的多样性。现在我们将更多地讨论抗体、ChIP-seq、数据生成项目和原始数据。

数据处理：比对、质量控制和峰识别 🔬

现在让我们深入探讨我们实际上如何处理数据：如何推断峰、读段去向以及这些修饰的位置。

这就是染色质免疫沉淀结合测序的用武之地。染色质免疫沉淀结合测序是什么意思？首先，染色质基本上是DNA和蛋白质的组合，它们维系着你的遗传信息。染色质通常指构成组蛋白的蛋白质、缠绕在它们周围的DNA以及将其维系在一起的辅助因子。

“免疫”基本上意味着我将使用抗体，即利用免疫系统来做这件事。“沉淀”意味着下拉某些东西，使其沉淀，有点像雨水沉淀。

我们如何找出基因组中特定修饰所在的所有位置？方法是：你基本上构建一种抗体（我们用这种可爱的抗体结构表示），然后你基于该修饰的许多副本来训练这种抗体。你生成大量该修饰，将其喂给兔子、山羊、骆驼或其他具有良好免疫系统的动物，这些动物会产生大量这些抗体，然后你生成这种抗体库存，世界各地的人订购它。

如果你想绘制例如启动子标记H3K4me3或增强子标记H3K27ac，你基本上只需订购大量那种抗体。当你得到抗体后，你使用该抗体下拉你的染色质。你首先片段化DNA，这样每次下拉时，你只下拉实际被该修饰结合的区域。然后你使用抗体将其下拉，然后纯化DNA。你首先进行交联以将蛋白质和DNA捆绑在一起，然后片段化，然后下拉，用你的抗体选择，然后逆转交联并纯化DNA，然后使用我们上次讨论的不同类型的测序技术对DNA进行测序。

你通常还会做一个输入DNA对照，即不使用抗体或使用像IgG这样的不结合任何东西的抗体。然后你基本上从你的对照中得到一堆读段，从例如抗体到Pol II或STAT2或CSE4或任何转录因子中得到一堆读段，或者你可以专门针对组蛋白修饰构建你的抗体。完成后，你进行测序，然后将它们比对到基因组。

所以，再次强调：你有你的细胞核，你进行交联，添加你的蛋白质或修饰特异性抗体，用这种抗体下拉，逆转交联，得到一堆DNA，对DNA进行测序，并将DNA比对回基因组。

你最终得到的是：这基本上是对应于x轴上每个基因组位置，y轴上我得到的测序读段数量的表示，对应于该修饰。

这里真正酷的是，在这种下拉数据中有大量的数据。你首先会注意到这种凹凸不平的性质，这基本上是个别核小体的定位。所以你基本上在核小体所在的位置下拉了很多，它们具有大约200个碱基对的周期性。

你注意到的第二件事是，存在这些无核小体区域。你在这里看到很多核小体，你下拉了组蛋白修饰标记H3K4me1，然后有一个低谷，你不再看到它们，然后又是更多。这个低谷也存在于这个H3K4me3数据中。那个低谷是什么？那个低谷基本上是启动子区域，通常是一个无核小体区域，即核小体通常不占据的区域。所以你基本上可以在那里看到那个特征。

你注意到的另一件事是，H3K4me1（这个增强子标记）可以看到它侧翼于转录起始点，所有增强子都在转录的上游和下游，在这个特定例子中，位于第一个和第二个内含子内。你看到这个启动子标记H3K4me3，它相当精确地定位在该基因的启动子和起始点。然后你看到H3K36三甲基化，这是我之前用绿色标记转录区域的标记。

所以，H3K36me3是转录区域的标记，启动子用H3K4me3，增强子用H3K4me1和H3K27ac。还有与延伸、Polycomb抑制、异染色质抑制等相关的额外标记。

这基本上是你对所有那些组蛋白修饰染色质免疫沉淀实验或ChIP-seq实验得到的主要数据。每个测序标签大约有30个碱基长，就像我的小片段一样。然后我将这些标签比对到30亿个碱基的参考基因组中的唯一位置。我得到的读段数量取决于测序深度、读段的密度和定位、浓度以及这些读段的分布，告诉我该修饰在基因组中的位置。

当然，我们不仅用一种标记做这个，而是用许多不同的标记。所以，这里是你之前看到的H3K4me3，这是你之前看到的H3K4me1，这是你之前看到的H3K36me3，但然后我们在许多额外标记的背景下看到这些。这是一个Polycomb抑制标记H3K27me3。相比之下，H3K27ac是增强子区域的最佳标记之一。注意它们都在组蛋白H3的第27位赖氨酸上。这基本上意味着该位置只能有一种标记，因为每个基因组有两条染色体，每个细胞有每条染色体的两个拷贝。当然，不仅仅是两条染色体，而是每条染色体的两个拷贝。两个拷贝实际上可能有不同的标记，所以一个拷贝可能有H3K27ac，另一个拷贝实际上可能有H3K27me3。所以这些被称为二价区域，其中两个姐妹染色体或混合细胞中的不同细胞基本上具有活跃或抑制标记的组合。

这就是数据的样子。这是原始数据。我们今天将要学习的一件事是如何将这些数据转化为染色质状态，使用我们之前看到的隐马尔可夫模型框架，你可以基于观察到的不同组蛋白修饰标记的定位，推断每个位置表观基因组的隐藏状态。

我们讨论了很多探测表观基因组的不同技术。我们特别讨论了染色质免疫沉淀，以及我们如何使用这些抗体在下拉染色质后找出其位置并相应地进行比对。

当然，问题是围绕这些的一些计算挑战是什么？其中一个计算挑战是，我们如何通过这种表观基因组下拉（这种抗体染色质免疫沉淀下拉）的视角，超快速地对人类基因组进行重测序比对？我们基本上只是对基因组片段进行重测序，我们需要将数十万甚至数百万的测序读段比对到相同的参考基因组。

你最不想做的事情就是动态规划。动态规划基本上会进行匹配，但速度不够快。BLAST是一个非常好的解决方案，你可能想要使用它，因为它使用了哈希方法。但BLAST的挑战在于它具有巨大的内存占用。因此，当每个人都想从他们的所有实验中BLAST 1亿个读段时，最终会导致服务器崩溃，人们无法在自己的机器上完成，因为哈希表的内存占用太大了。

所以我们要看的第一件事是如何将数百万个短读段比对到基因组。首先，我们将讨论传统的哈希方案，然后我们将讨论一种非常酷的计算技术，称为Burrows-Wheeler变换。

问题是：大家

🧬 P9：L9- 表观基因组学L2 部分和 3D 基因组

在本节课中，我们将继续学习表观基因组学。我们将探讨如何确定染色质状态模型的复杂度，如何跨多种细胞类型联合学习染色质状态，以及如何利用多细胞类型数据来解析增强子调控网络。最后，我们将转向三维基因组学，学习如何研究细胞核内基因组的空间组织。

上一节我们介绍了表观基因组学的基础知识、数据生成和处理流程。本节中，我们来看看更高级的分析方法。

📊 模型复杂度选择

在上一讲中，我们学习了如何使用多元隐马尔可夫模型来发现不同的染色质状态。一个关键问题是：我们如何确定模型中染色质状态的数量和所需的表观遗传标记？

一个简单的答案是，这取决于输入数据的丰富程度以及我们解释生物学差异的能力。输入数据越多，区分不同状态的能力就越强。然而，传统的机器学习模型选择标准（如贝叶斯信息准则）在基因组这种极其复杂的数据上往往倾向于选择越来越多的状态，难以确定一个合理的阈值。

因此，我们需要其他指标。一个有用的技巧是评估模型是否捕获了不同标记之间的依赖关系。

以下是评估模型捕获标记间依赖关系的方法：

计算每个染色质状态下，单个标记的预期出现频率（即发射概率）。
基于这些单个频率的乘积，预测任意两个标记在该状态下共同出现的频率。
将预测的共同出现频率与在该状态注释的基因组位置上实际观察到的共同出现频率进行比较。

如果模型能很好地捕获标记间的依赖关系，预测值和观测值将落在对角线上。如果存在大量偏离对角线的点，则说明当前模型复杂度不足以捕获这些标记组合关系，可能需要增加状态数量。

为了解决随机初始化对模型结果的影响，并更稳健地选择模型复杂度，可以采用一种称为“嵌套初始化”的计算技巧。

以下是嵌套初始化的步骤：

首先，学习一个包含大量状态（例如80个）的模型，以捕获所有可能的相关状态。
然后，采用贪心策略，逐步剔除信息量最少的冗余状态。
最终，基于生物学的可解释性，选择一个主观的截止点来确定最终的状态数量。

这种方法比随机初始化更稳健，能更一致地捕获关键的染色质状态。

🔗 跨细胞类型的联合学习

之前的学习都是针对单一细胞类型。现在，我们希望利用多种细胞类型的数据进行联合学习，以确保染色质状态定义的一致性，并研究其动态变化。

主要有三种策略来实现跨细胞类型的染色质状态联合学习。

以下是三种主要的联合学习策略：

独立学习后聚类：在每个细胞类型中独立学习一个染色质状态模型（例如15个状态），然后基于发射概率矩阵或基因组注释的后验概率对所有模型进行聚类，将相似的州合并。
堆叠法：将所有细胞类型的所有标记数据堆叠成一个很长的向量，让每个隐藏状态发射这个完整向量。这种方法会组合性地增加状态数量（同时捕获元件类别和细胞类型特异性），在细胞类型很多时难以扩展。
串联法：将不同细胞类型的基因组数据像串联不同染色体一样串联起来，形成一个“巨型基因组”，然后在此之上学习一套共同的染色质状态定义。这种方法要求分析的标记相同，或能将缺失数据作为缺失值处理。这是目前最广泛采用的方法。

串联法确保了所有细胞类型使用共同的状态定义，便于直接比较。通过这种方法，我们可以将每个细胞类型总结为一条染色质状态轨道，然后比较不同细胞类型间相同基因组区域的状态变化（如从“预备启动子”变为“活跃启动子”或“抑制状态”）。

多细胞类型分析的强大之处在于，它允许我们寻找基因组不同位置之间相关的活性模式。例如，一个基因间增强子与某个基因在不同细胞类型中表现出同步的开启/关闭模式，这提示该增强子可能调控那个基因，而不是其物理位置上更近的其他基因。这种相关性分析有助于我们解析基因调控回路。

🧩 利用活性谱解析调控网络

通过跨细胞类型联合学习得到的染色质状态活性谱（即每个元件在不同细胞类型中是活跃还是沉默的向量），是解析增强子调控网络的强大工具。

我们可以利用这些活性谱做两件事：1）将增强子与其靶基因连接；2）将转录因子与其调控的增强子连接。

连接增强子与靶基因：通过计算增强子活性谱与附近基因表达谱在不同细胞类型间的相关性，可以推断增强子的靶基因。当增强子活性与某个基因的表达高度正相关时，该基因很可能是其靶标。

连接转录因子与增强子：这需要一个三向相关性分析。我们同时考察：1) 增强子的活性谱；2) 特定转录因子结合基序在该增强子区域的富集程度；3) 该转录因子自身的表达谱。如果三者呈正相关，则预测该因子是激活子；如果增强子活性与因子表达或基序富集呈负相关，则预测该因子是抑制子。

通过系统性地进行上述分析，可以预测调控特定细胞类型增强子模块的转录因子，并构建组织特异性的调控网络。

🔮 表观基因组插补

在实际项目中，我们可能无法在所有细胞类型中检测所有表观遗传标记。表观基因组插补技术旨在利用已测序标记在不同细胞类型间的相关性，来预测缺失的标记数据。

核心方法是使用随机森林回归等计算技术。预测一个目标细胞类型中的目标标记时，输入特征包括：

同一细胞类型，不同标记：其他标记在该细胞类型中的信号。
不同细胞类型，同一标记：目标标记在其他细胞类型中的信号。
不同细胞类型，不同标记：其他标记在其他细胞类型中的信号，它们与目标标记的相关性可作为信息。

通过这种方法，即使某些标记在特定细胞类型中未被实验检测，也能高精度地预测其信号。有趣的是，插补得到的数据有时甚至比低质量的实验观测数据更可靠，并且能更好地揭示细胞类型间的生物学关系。

🧵 三维基因组学简介

现在，我们将话题转向三维基因组学，研究细胞核内基因组的空间组织。DNA在细胞核内并非随机分布，而是高度有序地折叠和包装。

研究三维基因组的主要实验方法分为两类：基于标记的方法和基于相互作用的方法。

以下是研究三维基因组的主要实验方法：

基于标记的方法（如DamID）：以核纤层等细胞核结构为“地标”，通过类似于染色质免疫沉淀的技术，找出与这些地标相互作用的基因组区域，从而定义如核纤层关联域。
基于相互作用的方法（如Hi-C）：通过甲醛交联将空间上临近的DNA片段“粘”在一起，随后将DNA片段化并重新连接，形成嵌合片段，最后通过高通量测序鉴定这些相互作用。这可以系统性地绘制全基因组范围的染色质互作图谱。

通过Hi-C数据分析，可以揭示三维基因组的几个关键特征：

以下是三维基因组的关键特征：

区室化：基因组被分为两个主要区室。A区室（活跃）基因丰富、染色质开放，通常位于核中央；B区室（抑制）基因贫乏、染色质紧凑，常位于核周边。
拓扑关联域：染色体被进一步划分为一些亚区域，域内的DNA片段彼此频繁互作，而与域外的片段互作较少。这些域被称为拓扑关联域。
染色质环：由CTCF蛋白和黏连蛋白复合物介导，形成更局部的、增强子-启动子之间的特异性环状结构，这对基因调控至关重要。

目前主流的“环挤压”模型认为，CTCF蛋白像锚点一样结合在DNA特定方向位点上，黏连蛋白复合物则像“转轴”一样将中间的染色质丝环化挤出，直至遇到CTCF锚点为止，从而形成TAD和染色质环。

研究表明，尽管基因组序列在进化中重排，但TAD边界和核纤层关联在人类和小鼠之间具有相当程度的保守性。有趣的是，预测一个区域是否与核纤层关联的最佳特征之一，竟然是其简单的AT碱基含量。

本节课中我们一起学习了如何选择染色质状态模型的复杂度，以及通过跨细胞类型联合学习来解析染色质状态动态和基因调控网络。我们还介绍了利用标记相关性进行表观基因组插补的强大技术。最后，我们探讨了三维基因组学的基本概念、研究方法和主要特征，包括区室、拓扑关联域和染色质环，以及它们对基因调控的重要意义。

posted @ 2026-02-03 19:59 绝不原创的飞龙阅读(4) 评论(0) 收藏举报

刷新页面返回顶部

龙哥盟

掠夺·扩张·投机·博弈

MIT-6-047-基因组机器学习笔记-全-

MIT 6.047 基因组机器学习笔记（全）

课程 P1：L1 - 计算生物学导论 🧬

课程概述与目标 🎯

课程内容与模块 📚

教材、作业与测验 📖

最终项目 🚀

为什么选择计算生物学？💡

生物学基础速成 🧪

总结 📝

课程 P10：L10 - 调控基因组学与基序发现 🧬

什么是调控基序？🧩

如何发现基序？🔍

基序的表示：位置权重矩阵

基序发现的挑战与方法 🧪

方法一：基于共调控基因集的基序发现 🧬

期望最大化算法

吉布斯采样

方法二：基于进化特征的从头基序发现 🌍

识别单个基序实例：进化足迹 👣

总结 📚

📚 课程 P11：L11 - 网络分析基础与应用

🌐 网络与图论基础

🔍 基本图算法：广度优先与深度优先搜索

🌍 网络的涌现性质

📊 网络的线性代数视角

🎯 降维与稀疏表示

🧩 网络社区发现：谱聚类

📝 总结

🧠 深度学习课程 P12：从基础到基因组学应用

🧬 课程概述

🧠 受监督学习与神经网络

人工神经元与感知机

激活函数

梯度下降与学习

反向传播

避免过拟合

🔍 无监督学习与深度信念网络

受限玻尔兹曼机

深度信念网络

学习技巧

🎨 自编码器

🖼️ 现代深度学习架构：卷积神经网络

核心操作

关键优势

经典网络演进

⏳ 现代深度学习架构：循环神经网络

长期依赖问题与LSTM

双向RNN与架构组合

🔄 迁移学习

🎭 生成对抗网络

🧬 深度学习在调控基因组学中的应用

序列表示与卷积神经网络

多任务学习与蛋白质结构预测

💻 深度学习计算工具与实践

📚 课程总结

群体遗传学课程 P13：人类遗传变异与群体差异 🧬

遗传学历史与基础概念

遗传变异的类型与影响

等位基因的分类与频率谱

遗传变异的表示、检测与项目

单倍型、连锁不平衡与重组热点

基因型填充与定相

群体亲缘关系、历史与祖先推断

总结与展望

🧬 P14：L14 - GWAS与疾病解析

🧬 第一部分：孟德尔性状与连锁分析

🧬 第二部分：复杂性状与全基因组关联研究

🧬 第三部分：精细定位与多研究整合

🧬 第四部分：从遗传位点到疾病机制

🧬 第五部分：案例研究：解析FTO肥胖位点

🎯 总结

课程 P15：L15 - eQTLs 表达数量性状位点 🧬

概述：为何研究eQTLs？

1. eQTLs的动机与概念框架

2. eQTLs的研究案例：阿尔茨海默病中的甲基化QTLs

3. eQTLs的方法学基础

4. eQTLs的生物学洞察