生物信息学概述 (第6章)

6.1.1 基因组学的概念

基因组学研究的范畴可以从三个层次理解：
（1）结构基因组学（Structural Genomics）
结构基因组学关注基因组的物理结构，包括DNA序列的测定、基因组图谱的构建、基因的定位和结构分析等。其核心目标是获得基因组的完整序列，并理解基因在基因组中的排列方式。
（2）功能基因组学（Functional Genomics）
功能基因组学研究基因组中基因的功能、基因之间的相互作用、基因表达调控等。与结构基因组学回答"是什么"不同，功能基因组学回答"做什么"和"怎么做"。主要技术包括RNA-Seq、ChIP-Seq、ATAC-Seq等。
（3）比较基因组学（Comparative Genomics）
比较基因组学通过比较不同物种或不同个体基因组的异同，揭示基因组的进化规律、基因的起源和功能分化、物种之间的亲缘关系等。

基因组学中所研究的问题主要分为三类：
第一类：如何获得基因组序列
这是基因组学最基础的问题。随着DNA测序技术的发展，从一代Sanger测序到二代高通量测序（NGS），再到三代单分子长读长测序，获得基因组序列的成本急剧下降、速度大幅提升。一个物种的基因组项目通常包括：基因组大小估计、测序策略设计、序列组装、质量评估等步骤。
第二类：如何解读/解码基因组
获得序列后，需要解读这本"生命之书"。这包括：基因预测与结构注释（识别编码区、内含子、外显子、调控元件等）、功能注释（确定基因的功能、参与的通路）、比较分析（与近缘物种比较，理解进化关系）、变异分析（识别SNP、Indel、SV等）。
第三类：如何重写/编写新的基因组
这是合成生物学的前沿方向。通过基因编辑技术（如CRISPR-Cas9）对现有基因组进行精确修饰，甚至从头设计和合成全新的基因组。2010年，Craig Venter团队首次合成了支原体的人工基因组，标志着人类开始具备"编写生命"的能力。

人类基因组计划采用了层次测序策略：通过构建DNA大片段克隆文库、重叠克隆群和物理图谱，再进行逐个克隆鸟枪法测序。中国科学家承担了1%的测序任务（3号染色体端部约3000万个碱基对），成功跟进世界先进测序技术，推动了国内测序仪研发和生物信息学软件开发。
HGP的完成不仅提供了人类基因组的参考序列，更重要的是：
1. 技术推动：测序技术从Sanger法发展到高通量测序
2. 数据共享：建立了基因组数据开放共享的文化
3. 疾病研究：为疾病基因定位、药物基因组学奠定基础
4. 进化研究：为理解人类起源和进化提供基础数据

Sanger测序的核心是利用ddNTP终止DNA链延伸：
1. 反应体系：模板DNA、引物、DNA聚合酶、四种dNTP（其中一种带有放射性或荧光标记）和少量ddNTP
2. 链延伸与终止：DNA聚合酶沿模板合成新链，遇到dNTP时正常延伸，遇到ddNTP时因缺少3'-OH而终止
3. 四个独立反应：分别加入ddATP、ddCTP、ddGTP、ddTTP，产生四组不同长度的片段
4. 电泳分离：通过聚丙烯酰胺凝胶电泳（PAGE）或毛细管电泳按大小分离片段
5. 读取序列：从短到长读取末端碱基，获得互补链序列
自动化Sanger测序使用四种不同荧光标记的ddNTP，四个反应在同一管中进行，通过激光检测荧光颜色确定碱基类型。
局限性：通量低（每次反应只能读取一条序列），成本高，不适合大规模基因组项目。

Illumina测序流程包括四个关键步骤：
（1）文库制备
- 将基因组DNA随机打断成200-500 bp片段
- 末端修复、加A尾、连接测序接头（adapter）
- PCR扩增文库（可选）
（2）Flow Cell杂交与桥式PCR
- Flow Cell表面固定有两种与接头互补的引物
- 单链DNA文库与引物杂交
- 通过桥式PCR扩增：DNA片段弯曲成桥状，以另一端为模板进行扩增
- 经过多轮扩增，每个原始分子形成一个约1000-2000个克隆的"簇"（cluster）
- 变性后得到单链簇，作为测序模板
（3）边合成边测序（SBS）
- 加入DNA聚合酶和四种带有不同荧光标记的dNTP（每种dNTP的3'-OH被可切除的化学基团阻断）
- 每个循环：
a. 四种dNTP竞争掺入，每次只有一个碱基被加入
b. 激光扫描，根据荧光颜色确定掺入的碱基
c. 切除荧光基团和3'阻断基团，使链可以继续延伸
- 重复上述循环，逐碱基读取序列
（4）双端测序（Paired-End Sequencing）
- 从簇的两端分别测序，获得两段各约150 bp的序列
- 两段序列之间距离已知（约200-500 bp）
- 双端信息对于基因组组装和结构变异检测至关重要
数据格式与质量控制
测序数据以FASTQ格式存储，每条记录包含四行：序列标识符、碱基序列、"+"分隔符、碱基质量值（Phred分数）。Phred质量值Q与碱基错误率ε的关系为：Q = -10log₁₀(ε)。Q30表示碱基错误率约为0.1%，即准确性为99.9%。

（1）PacBio SMRT测序
PacBio测序的核心是零模波导孔（ZMW）技术：
- 在芯片上蚀刻数十万个直径约70 nm、深度约100 nm的纳米孔
- 每个孔底部固定一个DNA聚合酶分子
- 待测DNA与引物结合，被聚合酶捕获
- 加入四种带有不同荧光标记的dNTP（标记在磷酸基团上，而非碱基上）
- 当dNTP被掺入新链时，荧光标记被释放到溶液中，产生闪光信号
- 记录每个ZMW中的闪光序列和时间，即可获得DNA序列
PacBio Sequel II可产生约15-20 kb的长读长，最新的Revio系统进一步提高了通量。通过环化一致性测序（CCS，Circular Consensus Sequencing），同一个分子多次循环测序，可获得准确性>99.9%的HiFi reads。
（2）Nanopore测序
Nanopore测序的核心原理：
- 蛋白质纳米孔嵌入脂质双分子层膜中
- 在膜两侧施加电压，产生离子电流
- DNA/RNA分子在动力蛋白（helicase）驱动下穿过纳米孔
- 不同碱基（或k-mer）穿过纳米孔时引起特征性的电流变化
- 通过深度学习模型（如Bonito、Dorado）将电流信号解码为碱基序列
Nanopore的优势：
- 超长读长：目前最长记录超过4 Mb
- 便携性：MinION设备仅U盘大小，可在野外、太空等极端环境使用
- 直接测序RNA：无需逆转录，可检测RNA修饰
- 实时测序：数据边测边出，可即时分析
局限性：原始读取准确性较低（约92-97%），但R10.4.1 flow cell配合Dorado basecaller的准确性已大幅提升。

由于二代测序准确性高但读长短，三代测序读长长但原始准确性低，混合测序（Hybrid-Seq）成为当前基因组项目的标准配置：
- 用三代长读长进行基因组骨架组装
- 用二代短读长对组装结果进行纠错和填补缺口
- Hi-C数据辅助染色体级别组装

传统的基因组组装在高度重复序列区域（如着丝粒、端粒、rDNA阵列）存在大量缺口，主要原因是：
- 二代短读长无法跨越长重复单元
- 重复序列使组装算法难以确定正确的路径
T2T联盟实现人类完整基因组的关键技术组合：
（1）PacBio HiFi长读长
- 读长约15-20 kb，准确性>99.9%
- 足以跨越大多数重复单元（如着丝粒卫星重复单元约171 bp）
- 高准确性使得区分高度相似的重复序列成为可能
（2）Oxford Nanopore超长读长
- 读长可达100 kb至数Mb
- 直接跨越最大的重复阵列（如5S rDNA阵列可长达数Mb）
- 提供长距离的序列连续性信息
（3）Hi-C辅助scaffold连接
- 将HiFi组装产生的contig提升到染色体级别
- 确定contig的顺序、方向和间隙大小
（4）手动审校与验证
- 对复杂区域（如着丝粒）进行手动检查
- 使用荧光原位杂交（FISH）和光学图谱进行验证

Hi-C实验流程：
1. 甲醛交联：将空间上接近的DNA片段共价交联
2. 酶切：用限制性内切酶（如HindIII、DpnII）切割DNA
3. 末端修复和生物素标记：补平粘性末端并掺入生物素标记的dNTP
4. 连接：稀释后连接，使空间上接近的片段连接在一起
5. 去交联和测序：打断DNA，富集生物素标记的嵌合片段，进行双端测序
Hi-C数据分析用于基因组组装：
- 构建接触矩阵（contact matrix）：染色体各区域之间的相互作用频率
- 同一染色体上的区域相互作用频率高（呈现对角线富集）
- 相邻区域相互作用频率高于远端区域
- 使用算法（如3D-DNA、SALSA、ALLHiC）将contig聚类到染色体组、确定顺序和方向

泛基因组分为两部分：
- 核心基因组（Core Genome）：物种所有个体共有的基因，通常与基本生命功能相关
- 可变基因组/附属基因组（Accessory/Dispensable Genome）：仅在部分个体中存在的基因，通常与环境适应、致病性、代谢多样性等相关
泛基因组构建方法：
- 迭代式：将新个体的基因组与已有泛基因组比较，添加新序列
- 图基因组（Graph Genome）：将所有序列变异表示为图结构，而非线性序列
泛基因组的应用：
- 捕获参考基因组中缺失的遗传多样性
- 改进变异检测（减少比对偏差）
- 识别与重要性状相关的基因（如作物抗病基因）

一个典型的T2T基因组项目测序策略：
| 数据类型 | 测序深度 | 用途 |
|----------|---------|------|
| PacBio HiFi | 30-50× | 主要组装，产生高质量contig |
| Nanopore Ultra-long | 10-20× | 跨越超长重复区域，解决复杂区域 |
| Illumina PE150 | 50-100× | 纠错、质量验证 |
| Hi-C | 60-100× | 染色体级别scaffold |
| 光学图谱 | 5-10× | 验证组装结构 |

（1）线性序列展示
- 以染色体为坐标轴，从左到右展示DNA序列
- 支持不同尺度的缩放（从单碱基到全染色体）
- 显示正反链的序列和注释
（2）多轨道（Track）整合
- 每个track代表一类数据（基因注释、测序reads、变异、保守性等）
- 用户可自定义加载、排列和配置track
- 支持track的叠加和比较
（3）交互式查询
- 通过基因名、染色体坐标、rsID等定位基因组区域
- 支持序列搜索（如查找特定基序）
- 支持不同物种间的基因组比对（synteny view）
（4）数据导出
- 导出特定区域的序列（FASTA）
- 导出track数据（BED、GTF等）
- 生成高质量图片用于发表

（1）UCSC Genome Browser
- 最经典的基因组浏览器之一
- 整合了大量注释数据和工具
- 支持自定义track加载（通过URL或本地文件）
- 强大的Table Browser用于数据查询和下载
（2）Ensembl Genome Browser
- 欧洲生物信息研究所（EBI）维护
- 注释质量高，更新及时
- 强大的比较基因组学工具（如Ensembl Compara）
- 提供API和BioMart数据查询工具
（3）IGV（Integrative Genomics Viewer）
- 本地运行的桌面软件
- 支持大规模数据（如BAM、VCF、Hi-C矩阵）
- 特别适合查看测序reads比对情况和变异位点
- 支持基因组重排和结构变异的可视化
（4）JBrowse/GBrowse
- 可部署在本地服务器的网页浏览器
- 适合构建物种特异性的基因组数据库
- 支持插件扩展

BED格式
BED文件至少包含3列，最多12列：
- chrom：染色体名称
- chromStart：起始位置（0-based）
- chromEnd：终止位置（1-based，不包含）
- name：名称（可选）
- score：分数（0-1000，可选）
- strand：方向（+/-/.）
- thickStart：编码区起始（可选，用于显示）
- thickEnd：编码区终止（可选）
- itemRgb：颜色（可选）
- blockCount：外显子数量（BED12）
- blockSizes：外显子大小列表（逗号分隔）
- blockStarts：外显子相对起始位置列表
WIG格式
variableStep格式示例：

variableStep chrom=chr1 span=50
1001 12.5
1051 13.2
1101 11.8

表示chr1上1001-1050位置的值为12.5，1051-1100位置的值为13.2。
GTF格式
GTF是GFF2的变体，每行9列：
1. seqname：染色体或scaffold
2. source：注释来源
3. feature：特征类型（gene, transcript, exon, CDS等）
4. start：起始位置（1-based）
5. end：终止位置（1-based，包含）
6. score：分数（通常为".")
7. strand：方向（+/-/.)
8. frame：阅读框（0/1/2/.）
9. attributes：属性（gene_id, transcript_id等）

现代基因组研究往往涉及多组学数据的整合：
- 基因组层：基因结构、变异位点
- 转录组层：RNA-Seq信号、剪接异构体
- 表观组层：ChIP-Seq峰（组蛋白修饰、转录因子结合）、DNA甲基化、ATAC-Seq开放染色质
- 三维结构层：Hi-C相互作用矩阵、TAD边界
- 进化层：跨物种保守性分数（phyloP, phastCons）
整合可视化可以帮助发现跨层面的关联，例如：
- 某个增强子区域（H3K27ac峰+ATAC-Seq开放区）与远处基因的启动子有Hi-C相互作用
- 一个SNP位于保守非编码区，且与附近基因的表达量显著相关

6.2.1 基因组的组装

（1）重复序列
基因组中普遍存在重复序列，包括串联重复（如微卫星）和散在重复（如转座子、segmental duplication）。重复序列导致组装歧义——相同的reads可能来源于基因组的不同位置，组装算法无法确定其真实来源。
（2）测序错误
Illumina测序错误率约0.1-1%（Q30-Q20），Nanopore原始错误率可达5-15%。测序错误会干扰k-mer的正确连接，在DBG中产生"气泡"（bubble）和"死胡同"（tip）。
（3）覆盖度不均
基因组不同区域的测序深度差异很大。高GC区域、AT富集区域、重复序列区域的覆盖度往往偏低，导致这些区域组装困难。
（4）杂合度
二倍体生物的两条同源染色体存在差异。组装软件需要区分同源序列并可能分别组装（phased assembly），否则杂合位点会导致组装碎片化。

DBG算法是二代短读长组装的核心方法，其工作流程如下：
（1）k-mer分解
将所有reads拆分为长度为k的短片段（k-mer）。例如，序列"ATGCGGAT"在k=5时产生：ATGCG、TGCGG、GCGGA、CGGAT。通常k取奇数（如21、31、51、71），且不超过read长度的2/3。
（2）构建DBG
- 每个k-mer作为图中的一个顶点（node）
- 若两个k-mer之间有k-1个碱基重叠，则在它们之间建立一条有向边
- 例如，ATGCG和TGCGG共享TGCG（4个碱基），建立从ATGCG到TGCGG的边
- 将不同reads产生的相同k-mer合并为一个顶点，记录覆盖度
（3）图简化
测序错误会在DBG中产生两类错误结构：
- Tips（死胡同）：由于测序错误产生的低频k-mer连接到主路径后又很快终止。处理：去除长度和覆盖度低于阈值的tips。
- Bubbles（气泡）：由于杂合位点或测序错误，主路径在某处分叉为两条路径，在下游又汇合。处理：保留覆盖度较高的路径，或在杂合组装中保留两条路径。
（4）contig提取
在简化后的DBG中，一条连续的路径（从没有入度的顶点开始到没有出度的顶点结束）对应一个contig。
k-mer大小的选择：
- k值较小：对重复序列的分辨率低（不同重复单元可能共享k-mer），但覆盖度高（每个read产生更多k-mer）
- k值较大：对重复序列的分辨率高，但覆盖度低（要求更大的测序深度才能保证每个k-mer被覆盖）
- 多k-mer策略（如SPAdes）：从较小的k开始组装，逐步增大k值，利用小k的连续性和大k的准确性优势

OLC算法更适合长读长数据（PacBio/Nanopore），其工作流程：
（1）Overlap（重叠检测）
计算所有reads两两之间的重叠关系。对于N条reads，需要O(N²)次比对，计算量巨大。对于长读长数据，通常使用minimap2等快速比对工具进行all-vs-all比对。
（2）Layout（布局构建）
根据重叠关系构建布局图（string graph）：
- 每个read是图中的一个节点
- 如果两个reads有显著重叠，则在它们之间建立边
- 边的方向和长度由重叠的位置和长度确定
- 简化图结构（去除 transitive edges，压缩 unambiguous paths）
（3）Consensus（一致性序列生成）
对布局图中的每个unitig（unambiguous path），将对应的reads进行多序列比对，生成一致性序列作为contig。

利用二代和三代数据的互补性：
（1）二代纠错三代
- 用高准确性的二代reads对三代reads进行纠错
- 代表工具：PBcR、LoRDEC、Jabba
- 纠错后的三代reads准确性可提升至>99%
（2）三代搭建骨架，二代填充
- 用三代长读长构建初步的contig骨架
- 用二代reads对三代contig进行 polish（填补gap、纠正碱基错误）
- 代表工具：Pilon、Racon
（3）graph-based混合组装
- 将二代和三代数据同时输入组装算法
- 在组装图中同时利用短读长的准确性和长读长的连续性
- 代表工具：HybridSPAdes、MaSuRCA

连续性指标：
- N50/L50：衡量contig/scaffold长度的分布
- 最大contig长度：最长的连续序列
- gap数量：scaffold中N碱基的数量和总长度
完整性指标：
- BUSCO（Benchmarking Universal Single-Copy Orthologs）：检测基因组中保守单拷贝基因的存在情况。高BUSCO完整性（如>90%）表示基因组组装较完整。
- CEGMA（Core Eukaryotic Genes Mapping Approach）：类似BUSCO，基于核心真核基因集。
准确性指标：
- QUAST：计算N50、错误组装（misassembly）数量、基因组覆盖度等
- Pilon polish后的二代比对一致性：检测碱基错误
- k-mer spectrum分析：与原始测序数据的k-mer谱比较，评估完整性

原核生物基因特征：
- 基因密度高（大肠杆菌基因组约85%为编码区）
- 无内含子（少数古菌有内含子）
- 无5'端帽子和3'端polyA尾
n- 操纵子结构：多个基因共用一个启动子，转录为多顺反子mRNA
- 起始密码子通常为ATG（少数GTG、TTG）
- 启动子元件：-10区（Pribnow框：TATAAT）和-35区（TTGACA）
真核生物基因特征：
- 基因密度低（人类基因组仅约1-2%为编码区）
- 含内含子（人类基因平均含8个外显子、7个内含子）
- 5'端有帽子结构，3'端有polyA尾
- 单顺反子mRNA（一个转录本对应一个蛋白质）
- 复杂的剪接机制：内含子5'端为GT（供体位点），3'端为AG（受体位点）
- 核心启动子元件：TATA框（约-30区）、CAAT框（约-80区）
- 存在大量的可变剪接（人类约95%的多外显子基因发生可变剪接）

从头预测的核心是利用基因组的统计特征和信号特征：
（1）开放阅读框（ORF）检测
- 寻找起始于ATG、终止于TAA/TAG/TGA的连续密码子序列
- 计算ORF长度：真核生物外显子通常>300 bp（100个氨基酸），短ORF可能是假阳性
- 计算密码子使用频率（Codon Usage Frequency, CUF）：真实基因通常偏好使用特定密码子
（2）隐马尔可夫模型（HMM）
HMM是基因预测中最常用的统计框架：
- 隐藏状态：基因间区、外显子（起始/中间/终止）、内含子（起始/中间/终止）、UTR等
- 观测符号：A、C、G、T
- 状态转移概率：从一个状态转移到另一个状态的概率（如外显子→内含子的概率）
- 发射概率：在某个状态下观察到特定碱基的概率
- Viterbi算法：寻找最可能的隐藏状态序列（即基因结构）
经典HMM基因预测软件：
- Genscan：针对人类和脊椎动物，基于广义HMM
- Augustus：支持真核和原核生物，可整合外部证据
- Glimmer/GlimmerHMM：Glimmer用于原核，GlimmerHMM用于真核
- GeneMark/GeneMark-ES：自训练HMM，适合新物种
- SNAP：可训练参数，适合非模式生物
（3）支持向量机（SVM）和深度学习方法
近年来，机器学习方法也被用于基因预测：
- 将序列特征（k-mer频率、GC含量、密码子偏好等）编码为特征向量
- 用已知基因训练分类器（SVM、随机森林、神经网络）
- 对新序列进行预测

（1）蛋白质同源比对
- 将已知物种的蛋白质序列与目标基因组进行比对（如tblastn）
- 根据蛋白质比对结果推断基因的编码区
- 考虑剪接位点信号，将基因组DNA与外显子对应
- 代表工具：GeneWise、Exonerate（protein2genome模式）
（2）cDNA/EST比对
- 将已知的cDNA或EST序列与基因组比对（如BLAT、GMAP）
- 直接确定外显子-内含子边界
- EST数据质量参差不齐，需要过滤和聚类
（3）跨物种比较基因组学
- 利用进化保守性：编码区通常比非编码区更保守
- 利用密码子替换模式：同义替换率（Ks）和非同义替换率（Ka）
- 代表方法：Twinscan（整合人类-小鼠比较信息）、CONTRAST

（1）RNA-Seq数据比对与组装
- 将RNA-Seq reads比对到基因组（如HISAT2、STAR）
- 使用StringTie、Cufflinks等工具组装转录本
- 获得GTF格式的转录本结构注释
（2）全长转录本测序
- PacBio Iso-Seq或Nanopore direct RNA测序获得全长cDNA
- 无需组装，直接获得完整的外显子-内含子结构
- 可发现新的转录本异构体
（3）整合策略（Evidence Modeler, EVM）
- 将从不同来源（从头预测、同源比对、转录本组装）获得的基因模型进行整合
- 根据证据的可靠性和一致性进行加权评分
- 输出一致性最高的基因模型集合
- 代表工具：EVM（Evidence Modeler）、PASA（Program to Assemble Spliced Alignments）、BRAKER

非编码RNA（ncRNA）是不编码蛋白质的RNA分子，包括：
- miRNA（微小RNA）：~22 nt，调控基因表达
- lncRNA（长链非编码RNA）：>200 nt，多种调控功能
- circRNA（环状RNA）：共价闭合环，miRNA海绵等
- tRNA、rRNA、snRNA、snoRNA等
预测方法：
- miRNA：基于发夹结构、种子区域保守性、热力学稳定性（RNAfold、miRDeep2）
- lncRNA：基于开放阅读框长度、编码潜能（CPAT、CPC2）、序列保守性
- circRNA：基于反向剪接位点（backsplice junction）的reads支持（CIRI、find_circ）

GO数据库的核心是一个有向无环图（DAG）结构，其中：
- 每个节点（GO term）代表一个功能概念
- 边代表"is_a"（是一种）或"part_of"（是...的一部分）关系
- 子节点比父节点更具体（如"蛋白质激酶活性" is_a "转移酶活性"）
GO注释的常用方法：
（1）基于序列相似性的GO注释
- 用BLAST将查询蛋白与GO注释数据库（如UniProt-GOA）比对
- 将相似蛋白的GO term转移给查询蛋白
- 常用工具：Blast2GO、eggNOG-mapper
（2）基于结构域的GO注释
- 通过InterProScan识别蛋白质结构域
- 利用结构域与GO term的映射关系进行注释
- InterPro数据库提供了Pfam、SUPERFAMILY等结构与GO的映射
（3）基于机器学习的方法
- 使用蛋白序列特征（氨基酸组成、物理化学性质、亚细胞定位信号等）训练分类器
- 预测GO term
- 代表方法：DeepGO（基于深度学习）

KEGG数据库的核心是通路图（Pathway Map），包括：
- 代谢通路（Metabolism）：碳水化合物代谢、能量代谢、脂质代谢、核苷酸代谢、氨基酸代谢等
- 遗传信息处理（Genetic Information Processing）：转录、翻译、折叠/分选/降解、复制与修复
- 环境信息处理（Environmental Information Processing）：膜转运、信号转导、配体-受体相互作用
- 细胞过程（Cellular Processes）：细胞生长与死亡、细胞群落、运输与分解代谢
- 生物体系统（Organismal Systems）：免疫系统、神经系统、内分泌系统等
- 人类疾病（Human Diseases）：癌症、神经退行性疾病、心血管疾病等
KEGG注释方法：
- 序列比对法：用BLAST/KAAS（KEGG Automatic Annotation Server）将基因与KEGG GENES数据库比对
- KO（KEGG Orthology）注释：将基因映射到KO编号（功能单位），再根据KO编号映射到通路
- 工具：KAAS、KOFAM、eggNOG-mapper

InterPro整合了多个蛋白质特征数据库：
- Pfam：基于HMM的蛋白质家族和结构域数据库
- PRINTS：基于指纹（保守基序组）的数据库
- PANTHER：基于系统发育树的蛋白质家族分类
- ProSite：基于模式（pattern）和谱（profile）的功能位点数据库
- SUPERFAMILY：基于SCOP结构分类的远缘同源识别
- SMART：信号肽、跨膜区、结构域识别
InterProScan工作流程：
1. 用多种方法（HMM、模式匹配、profile等）扫描蛋白质序列
2. 整合不同数据库的预测结果
3. 为每个预测的domain/feature提供GO term映射
4. 输出GFF3或TSV格式的注释结果

（1）COG/KOG/EggNOG
- COG（Clusters of Orthologous Groups）：原核生物蛋白质直系同源簇
- KOG：真核生物版本
- EggNOG：扩展版本，覆盖更多物种，提供GO、KEGG、CAZy等注释
- 将蛋白质分为功能类别（如J: Translation, K: Transcription, S: Unknown）
（2）NR/UniProt/Swiss-Prot
- NR（Non-Redundant）：NCBI的非冗余蛋白数据库，用于序列相似性搜索
- UniProt：最全面的蛋白质数据库，包含Swiss-Prot（人工审编）和TrEMBL（自动注释）
- 注释流程：BLASTp比对NR/UniProt → 取top hit → 继承功能描述
（3）专用数据库
- CAZy：糖类活性酶分类（糖苷水解酶、糖基转移酶等）
- TCDB：膜转运蛋白分类
- VFDB：细菌毒力因子
- CARD：抗生素耐药基因
- antiSMASH：生物合成基因簇（次级代谢产物）

现代基因组项目通常使用自动化流程进行功能注释：

预测蛋白质序列
    ↓
1. InterProScan（结构域识别 + GO注释）
    ↓
2. KAAS/KOFAM（KEGG通路注释）
    ↓
3. eggNOG-mapper（COG/GO/KEGG综合注释）
    ↓
4. BLASTp against NR/Swiss-Prot（功能描述）
    ↓
5. 整合所有注释结果

6.3.1 基本原理

基因组变异按尺度可分为：
点突变尺度：
- SNP（Single Nucleotide Polymorphism）：单个碱基的替换，是最常见的变异类型
- Indel（Insertion/Deletion）：小片段（通常<50 bp）的插入或缺失
结构变异尺度：
- 缺失（Deletion）：>50 bp的序列缺失
- 插入（Insertion）：>50 bp的序列插入
- 倒位（Inversion）：序列方向的翻转
- 易位（Translocation）：序列从一个染色体位置移动到另一个位置
- 重复（Duplication）：序列拷贝数的增加
- CNV（Copy Number Variation）：>1 kb的拷贝数变化

原始测序数据 (FASTQ)
    ↓
质量控制 (FastQC, Trimmomatic)
    ↓
比对到参考基因组 (BWA, Bowtie2)
    ↓
比对后处理 (排序, 标记重复, 碱基质量校正)
    ↓
变异检测 (GATK, Samtools, FreeBayes)
    ↓
变异过滤和质量控制
    ↓
变异注释 (功能影响预测)
    ↓
最终变异集合 (VCF格式)

（1）排序（Sorting）
比对后的SAM/BAM文件需要按染色体坐标排序，以便后续处理。
（2）标记PCR重复（Mark Duplicates）
PCR扩增会引入系统性偏差——来自同一模板分子的多个reads并非独立采样。使用Picard MarkDuplicates等工具识别并标记PCR重复reads，避免重复reads对变异检测的干扰。
（3）Indel区域重比对（Indel Realignment）
比对软件在Indel附近可能产生比对偏差（如将Indel附近的SNP错误比对）。GATK的IndelRealigner通过识别Indel附近不一致的比对，进行局部重比对，提高Indel和周围SNP的检测准确性。
（4）碱基质量校正（Base Quality Score Recalibration, BQSR）
测序仪给出的原始碱基质量值可能存在系统性偏差（如某些测序循环质量系统性地偏低）。BQSR通过将碱基质量与已知变异数据库（如dbSNP）进行比较，建立校正模型，重新校准碱基质量值。

变异检测的核心是判断某个位置的碱基差异是真实的遗传变异还是测序/比对错误。这本质上是一个统计假设检验问题。
零假设（H₀）：该位置没有变异，所有观察到的差异都是测序错误
备择假设（H₁）：该位置存在变异
测序错误通常被建模为二项分布或贝叶斯模型：
- 设测序错误率为ε（通常ε≈0.001，对应Q30）
- 在覆盖度为D的位置，观察到k个非参考碱基
- 如果H₀成立，k ~ Binomial(D, ε)
- 计算P(k | H₀)，如果该P值小于阈值（如0.05），则拒绝H₀，判定为变异
基因型似然模型：
更精确的方法是对三种可能的基因型（参考纯合RR、杂合RA、变异纯合AA）分别计算似然值：
- L(RR) = ∏ P(read_i | RR)
- L(RA) = ∏ P(read_i | RA)
- L(AA) = ∏ P(read_i | AA)
选择似然值最大的基因型作为该位点的基因型。
GATK HaplotypeCaller使用更先进的单倍型模型，在局部区域（通常数百bp）同时考虑所有可能的变异，通过局部de novo组装生成候选单倍型，然后计算每个单倍型的似然值。

即使使用了严格的统计模型，变异检测结果中仍包含大量假阳性。需要通过多重过滤来提高准确性：
基于测序质量的过滤：
- QD（Quality by Depth）：质量值除以覆盖度，过低提示质量不可靠
- FS（Fisher Strand）：检验正负链reads是否均衡支持变异，不均衡提示比对偏差
- MQ（Mapping Quality）：比对质量值，过低提示比对不可靠
- ReadPosRankSum：检验变异是否主要出现在reads末端（末端碱基质量通常较低）
基于覆盖度的过滤：
- DP（Depth）：覆盖度过低（<5-10×）的变异不可靠
- 覆盖度过高（>平均深度的2-3倍）可能提示重复区域
基于群体频率的过滤：
- 使用已知变异数据库（如gnomAD）过滤常见变异
- 在罕见病研究中，过滤在人群中频率>1%的变异

SNP主要来源于DNA复制过程中的自发突变：
（1）脱氨基作用
- 胞嘧啶（C）自发脱氨基转变为尿嘧啶（U），DNA修复系统将U识别为T，导致C→T突变
- 这是转换突变的主要来源，解释了为什么转换/颠换比（Ti/Tv ratio）通常大于2
- 5-甲基胞嘧啶（5mC）的脱氨基率更高，导致CpG岛中的C→T突变频率显著升高
（2）复制错误
- DNA聚合酶在复制过程中偶尔掺入错误的碱基
- 尽管有校对功能，错误率仍约为10⁻⁹/碱基/复制
（3）氧化损伤
- 活性氧（ROS）可导致碱基氧化修饰，如8-oxo-guanine可导致G→T颠换

按功能位置分类：
| 位置 | 比例 | 功能影响 |
|------|------|---------|
| 基因间区 | ~90% | 通常无直接功能影响 |
| 内含子 | ~5% | 可能影响剪接调控 |
| 外显子同义 | ~1% | 不改变氨基酸（密码子简并性） |
| 外显子错义 | ~1.5% | 改变氨基酸序列 |
| 外显子无义 | ~0.1% | 引入提前终止密码子 |
| 调控区 | ~2% | 可能影响基因表达 |
| 剪接位点 | <0.1% | 可能导致剪接异常 |
按群体频率分类：
- 常见变异（Common Variant）：群体频率>5%，通常与复杂疾病的易感风险相关
- 低频变异（Low-frequency Variant）：群体频率1-5%
- 稀有变异（Rare Variant）：群体频率<1%，可能与孟德尔遗传病相关

（1）基于pileup的方法（Samtools/bcftools）
- 统计每个基因组位置上比对reads的碱基组成
- 使用贝叶斯模型计算每个可能基因型的后验概率
- 优点：速度快，计算资源需求低
- 缺点：不处理Indel附近的比对偏差，对Indel不敏感
（2）基于单倍型的方法（GATK HaplotypeCaller）
- 在局部窗口（数百bp）内进行de novo组装，生成候选单倍型
- 将reads与候选单倍型进行比对，计算每个单倍型的似然值
- 选择最可能的单倍型组合作为基因型
- 优点：对Indel和复杂区域的SNP检测更准确
- 缺点：计算量大，耗时较长
（3）基于深度学习的方法（DeepVariant）
- 将每个候选变异位点转化为"图像"（read pileup的视觉表示）
- 使用卷积神经网络（CNN）进行变异识别
- 优点：准确性极高，在多个基准测试中优于传统方法
- 缺点：计算资源需求高，训练数据依赖

QD（Quality by Depth）：
QD = QUAL / DP
- QD过低（<2）提示变异可能被低质量reads或重复区域支持
- QD过高（>40）可能提示比对到重复区域
MQ（Mapping Quality）：
- 反映reads在该区域的比对可靠性
- MQ < 40提示比对可能存在问题
FS（Fisher Strand Bias）：
- 检验变异是否被正负链均衡支持
- FS > 60提示严重的链偏倚，变异可能不可靠
SOR（Strand Odds Ratio）：
- FS的替代指标，对覆盖度不均更稳健
- SOR > 3提示链偏倚
MQRankSum：
- 检验支持变异的reads与支持参考的reads的比对质量是否存在显著差异
- 负值过大提示支持变异的reads比对质量较差
ReadPosRankSum：
- 检验变异是否主要出现在reads末端
- 末端碱基质量通常较低，可能导致假阳性

转换/颠换比率（Ti/Tv ratio）是评估SNP集合质量的重要指标：
- 全基因组范围的期望Ti/Tv：约2.0-2.2
- 外显子区域的期望Ti/Tv：约2.8-3.0（因为同义SNP只发生在特定密码子位置）
如果检测到的SNP集合Ti/Tv显著低于期望值，可能提示：
- 存在大量假阳性（特别是颠换类型的假阳性较多）
- 过滤条件不够严格
- 样本存在污染

（1）DNA聚合酶滑移（Replication Slippage/Strand Slippage）
- 在重复序列区域（尤其是微卫星），模板链和新生链可能发生相对滑动
- 导致重复单元的插入或缺失
- 这是Indel最主要的产生机制，解释了为什么Indel在重复序列区域富集
（2）不等交换（Unequal Crossing Over）
- 减数分裂过程中的同源重组发生错位
- 导致一个染色体上重复序列增加，另一个上减少
（3）DNA损伤修复错误
- 双链断裂修复过程中可能发生小片段的插入或缺失

Indel检测比SNP检测困难的主要原因：
（1）比对困难
- 含有Indel的reads与参考基因组比对时，需要在Indel位置引入"空位"（gap）
- 比对算法对空位有罚分，倾向于将Indel区域的碱基错配解释为SNP
- 例如，参考序列"ATCGATCG"，reads为"ATC---GATCG"（3bp缺失），比对软件可能错误地比对为"ATCGATCG"（多个错配）
（2）Indel周围的碱基质量下降
- Illumina测序在Indel附近容易产生碱基识别错误
- 需要专门的Indel重比对步骤来纠正
（3）Indel大小的不确定性
- 对于较大的Indel（>10 bp），单个read可能无法完全跨越，导致Indel边界不确定
- 需要多个reads的信息来确定Indel的精确边界

（1）基于局部重比对的Indel检测（GATK HaplotypeCaller）
- HaplotypeCaller在局部窗口内生成候选单倍型（包含可能的SNP和Indel组合）
- 将reads与所有候选单倍型比对，选择最优解释
- 优点：同时考虑SNP和Indel，避免"SNP vs Indel"的比对歧义
（2）基于split-read的Indel检测（Pindel）
- 识别"split reads"——一条read的两部分分别比对到Indel两侧
- 通过分析未比对上的read片段来推断Indel序列
- 优点：对较大Indel（>20 bp）的边界定位准确
- 缺点：只能检测被reads跨越的Indel
（3）基于de novo局部组装的Indel检测（Platypus, Dindel）
- 在候选Indel区域提取相关reads
- 进行局部de novo组装
- 将组装结果与参考基因组比对，识别Indel

Indel的功能影响通常比SNP更严重：
| Indel类型 | 功能影响 | 示例 |
|----------|---------|------|
| 编码区非3倍数Indel | 移码突变，通常导致功能丧失 | BRCA1移码突变导致乳腺癌风险 |
| 编码区3倍数Indel | 氨基酸插入/缺失，可能影响功能 | CFTR ΔF508导致囊性纤维化 |
| 剪接位点Indel | 影响pre-mRNA剪接 | 深度内含子Indel激活隐蔽剪接位点 |
| 调控区Indel | 影响转录因子结合 | 胰岛素基因启动子Indel影响表达 |
| UTR区Indel | 可能影响mRNA稳定性或翻译 | 5'UTR的uORF相关Indel |
经典病例：
- 囊性纤维化：CFTR基因最常见的突变是ΔF508（一个3bp缺失），导致苯丙氨酸缺失，约占所有囊性纤维化病例的70%
- 亨廷顿舞蹈症：HTT基因中CAG重复序列的扩增（可视为一种特殊Indel），当CAG重复数>36时导致疾病

（1）基于Read深度（Read Depth, RD）的方法
原理：测序深度与基因组拷贝数成正比。
- 缺失区域：reads深度显著低于平均值（约50%对于杂合缺失）
- 重复区域：reads深度显著高于平均值
代表软件：CNVnator、BIC-seq、CNV-seq、DELLY（结合RD和PE）
优点：
- 可以检测大尺度的CNV
- 计算相对简单
局限：
- 无法精确定位断点（只能确定区域）
- 无法检测平衡性SV（如倒位、易位不改变拷贝数）
- 受GC含量、比对偏好等因素影响
（2）基于Split Read（SR）的方法
原理：识别一条read的两部分分别比对到不同位置的"split reads"，从而精确定位断点。
代表软件：Pindel、DELLY、LUMPY、Manta
工作流程：
1. 提取未完全比对或部分比对的reads（soft-clipped reads）
2. 将未比对部分重新比对到基因组其他位置
3. 如果找到匹配，确定断点坐标
优点：
- 断点定位精确（可达bp级别）
- 可以检测小至50 bp的SV
局限：
- 只能检测被reads跨越的SV
- 需要足够长的reads（对Illumina短读长有限制）
（3）基于Paired-End（PE）的方法
原理：利用paired-end测序中insert size和方向的异常来检测SV。
正常配对特征：
- 两个read在同一染色体上
- 方向相对（read1正向，read2反向）
- 距离在预期insert size范围内
异常配对类型及对应的SV：
| 异常类型 | 提示的SV |
|---------|---------|
| 距离远大于预期 | 缺失（reads跨越了被删除的区域） |
| 距离远小于预期 | 插入（reads之间的区域被插入序列替代） |
| 方向相同 | 倒位 |
| 位于不同染色体 | 易位 |
代表软件：DELLY、LUMPY、BreakDancer、Manta
优点：
- 可以检测大尺度SV
- 对read长度要求不高
局限：
- 断点定位不精确（通常在数百bp范围内）
- 受insert size分布的影响
（4）基于de novo组装的方法
原理：对SV区域进行局部de novo组装，将组装结果与参考基因组比对来识别SV。
代表软件：TIGRA、MindTheGap、SMRT-SV（PacBio专用）、Sniffles（Nanopore专用）
优点：
- 可以检测最复杂的SV，包括串联重复扩增、倒位等
- 可以获得SV的完整序列
局限：
- 计算资源需求大
- 组装质量依赖于数据质量和覆盖度

现代SV检测工具通常整合多种信号：
DELLY：整合PE、SR和RD信号
LUMPY：整合PE和SR信号，使用概率框架
Manta：结合PE、SR和局部组装，速度快，准确性高
GRIDSS：基于基因组图（assembly graph）的SV检测

长读长测序（PacBio/Nanopore）对于SV检测具有革命性优势：
1. 跨越复杂SV：长读长可以直接跨越整个SV区域，获得完整的SV序列
2. 精确定位断点：无需依赖PE或SR推断，直接比对即可定位
3. 检测重复区域SV：短读长难以比对的重复区域，长读长可以跨越
专用工具：
- Sniffles：Nanopore SV检测
- PBSV：PacBio SV检测
- SVIM：同时支持PacBio和Nanopore

SV检测的假阳性率通常高于SNP/Indel，验证至关重要：
实验验证方法：
- PCR+Sanger测序：对断点区域进行PCR扩增和测序，是SV验证的金标准
- qPCR：验证拷贝数变化
- Southern blot：验证大尺度重复/缺失
- 光学图谱（Bionano）：验证大尺度SV
计算验证方法：
- 使用不同检测工具的结果取交集
- 检查SV区域的reads比对情况（IGV可视化）
- 与已知SV数据库（如gnomAD-SV、DGV）比对

6.4.1 微生物组概述

人体微生物组的规模：
- 总重量：约1-2 kg（主要在肠道）
- 细胞数量：约3×10¹³个（与人体自身细胞数量相当或略多）
- 基因数量：约300-500万个基因（是人类基因组基因数的150-250倍）
人体各部位的微生物组特征：
| 部位 | 主要门类 | 细菌密度 | 主要功能 |
|------|---------|---------|---------|
| 结肠 | Firmicutes, Bacteroidetes | 10¹¹-10¹²/g | 发酵膳食纤维、合成维生素 |
| 口腔 | Firmicutes, Proteobacteria | 10⁹/mL唾液 | 参与氮循环、免疫调节 |
| 皮肤 | Actinobacteria, Firmicutes | 10⁶-10⁷/cm² | 屏障保护、免疫训练 |
| 阴道 | Lactobacillus | 10⁸-10⁹/mL | 维持酸性环境、防御病原 |
| 鼻腔 | Firmicutes, Actinobacteria | 10⁶-10⁷/cm² | 免疫防御 |

传统微生物学估计，环境中可培养的微生物仅占1%左右（"Great Plate Count Anomaly"）。主要原因包括：
1. 营养需求未知：很多微生物需要特殊的生长因子或复杂的共生关系
2. 环境条件难以模拟：如极端pH、温度、压力、厌氧条件等
3. 生长缓慢：某些微生物代时可达数周甚至数月
4. 群体感应依赖：需要达到一定密度才能生长
5. 严格的共生关系：依赖其他微生物产生的代谢产物
宏基因组学通过直接提取环境DNA进行测序，彻底 bypass 了培养瓶颈。

16S rRNA（原核生物）和18S rRNA/ITS（真核微生物）是微生物分类鉴定的"条形码"。
16S rRNA的优势：
1. 普遍存在：所有细菌和古菌都含有16S rRNA基因
2. 功能保守：核糖体功能对生命至关重要，序列演化相对缓慢
3. 大小适中：约1,500 bp，既包含保守区（设计通用引物），又包含可变区（区分物种）
16S rRNA的可变区：
- V1-V9共9个可变区
- 常用测序区域：V3-V4（约465 bp，Illumina MiSeq 2×300可覆盖）
- 不同可变区的分辨能力不同：V1-V2适合属水平，V3-V4可区分到种/株水平
分类阈值：
- 16S rRNA序列相似性 >97%：通常认为是同一个种
- 16S rRNA序列相似性 >95%：通常认为是同一个属

原始测序数据 (FASTQ)
    ↓
质量控制 (FastQC, cutadapt)
    ↓
去接头、去低质量、去嵌合体 (DADA2/Deblur/QIIME2)
    ↓
生成ASV/OTU表
    ↓
物种注释 (SILVA/Greengenes/RDP数据库)
    ↓
多样性分析 (α, β)
    ↓
统计分析 (LEfSe, ANCOM, DESeq2)
    ↓
可视化 (phyloseq, ggplot2)

（1）DADA2去噪算法
DADA2（Divisive Amplicon Denoising Algorithm 2）通过以下步骤生成ASV：
1. 质控过滤：去除低质量和含N的reads
2. 去重复：将相同序列合并，记录丰度
3. 错误学习：从数据中估计测序错误模型
4. 去噪：将含错误的序列纠正回其真实序列
5. 去嵌合体：识别并去除PCR嵌合体
6. 合并双端reads
（2）OTU聚类 vs ASV
- OTU聚类以97%相似性为阈值，将序列分组。缺点：阈值是人为设定的；近缘种可能因高于阈值而无法区分；信息在聚类过程中丢失
- ASV保留了每个独特的真实序列，分辨率更高，且可重复性更好

原始测序数据 (FASTQ)
    ↓
质量控制 (FastQC, Trimmomatic)
    ↓
去宿主DNA (Bowtie2比对宿主参考基因组)
    ↓
宏基因组组装 (MEGAHIT, metaSPAdes)
    ↓
基因预测 (Prodigal, MetaGeneMark)
    ↓
基因 catalog构建 (CD-HIT聚类)
    ↓
功能注释 (eggNOG-mapper, KEGG)
    ↓
分箱 (MetaBAT2, MaxBin2, CONCOCT)
    ↓
MAG质量评估 (CheckM)
    ↓
物种分类和系统发育分析

分箱的目标是将来自同一基因组的contig聚类到一起。
（1）基于组成的分箱（Composition-based Binning）
- 原理：不同物种的基因组具有不同的碱基组成特征（如GC含量、四核苷酸频率）
- 方法：将contig表示为高维特征向量（如136维四核苷酸频率向量），使用聚类算法（如K-means、高斯混合模型）进行聚类
- 代表工具：MetaBAT2、MaxBin2
（2）基于覆盖度的分箱（Coverage-based Binning）
- 原理：来自同一物种的contig在不同样本中的覆盖度变化模式一致（因为它们来自同一个基因组）
- 方法：在多个样本中计算每个contig的覆盖度，构建覆盖度矩阵，进行聚类
- 代表工具：CONCOCT、MaxBin2
（3）混合分箱
- 同时使用组成和覆盖度信息进行分箱
- 代表工具：MetaBAT2（使用TNF和覆盖度）、VAMB（使用变分自编码器）
MAG质量评估：
- CheckM：基于保守单拷贝基因（SCG）评估MAG的完整性和污染度
- 高质量MAG：完整性>90%，污染度<5%
- 中等质量MAG：完整性≥50%，污染度<10%

α多样性指数：
| 指数 | 衡量 | 公式 | 特点 |
|------|------|------|------|
| 物种数（Richness） | 丰富度 | S = 物种数量 | 不考虑丰度 |
| Shannon | 丰富度+均匀度 | H' = -Σ(pi × ln pi) | 对丰富度敏感 |
| Simpson | 丰富度+均匀度 | D = 1 - Σ(pi²) | 对优势种敏感 |
| Chao1 | 估计总丰富度 | Chao1 = S_obs + n₁²/(2n₂) | 考虑稀有物种 |
| Pielou's J | 均匀度 | J' = H'/ln(S) | 纯度的衡量 |
β多样性距离度量：
| 距离 | 公式 | 特点 |
|------|------|------|
| Bray-Curtis | 1 - 2Σmin(xi, yi) / Σ(xi + yi) | 考虑丰度，最常用 |
| Jaccard | 1 - |A∩B| / |A∪B| | 仅考虑有无 |
| Unweighted UniFrac | 考虑系统发育距离 | 考虑进化关系 |
| Weighted UniFrac | 考虑丰度和系统发育 | 综合考量 |
β多样性可视化：
- PCoA（主坐标分析）
- NMDS（非度量多维尺度分析）
- UPGMA聚类树

识别在不同组间丰度显著差异的微生物或功能基因：
- LEfSe（Linear Discriminant Analysis Effect Size）：结合统计检验和LDA评分，筛选组间差异标志物
- ANCOM（Analysis of Composition of Microbiomes）：考虑微生物组数据的组成性特点
- DESeq2/edgeR：基于负二项分布的差异分析（适用于OTU/基因计数数据）
- ALDEx2：考虑组成性和稀疏性的差异分析

微生物组数据具有以下"大数据"特征：
1. 高维度：数千个OTU/ASV × 数百个样本
2. 稀疏性：很多物种只在少数样本中出现（零膨胀分布）
3. 组成性：相对丰度总和为100%，变量间存在约束
4. 高变异性：个体差异大，受饮食、环境等多种因素影响
5. 多层次：物种、功能、代谢物等多个层面的数据

（1）松散的整合（Loose Integration）
- 分别分析各组学数据
- 比较不同组学层面的结果，寻找一致性
- 例如：比较差异物种和差异代谢物，寻找相关性
（2）统计学整合（Statistical Integration）
- 使用多变量统计方法整合多个数据矩阵
- CCA（Canonical Correlation Analysis）：寻找两个数据集之间的最大相关
- sPLS（sparse Partial Least Squares）：寻找预测关系的同时进行变量选择
- MCIA（Multiple Co-Inertia Analysis）：整合多个组学数据集
- MOFA（Multi-Omics Factor Analysis）：分解共享和特异的变异源
（3）网络整合（Network Integration）
- 构建跨组学的关联网络
- 例如：物种-代谢物关联网络、基因-代谢物关联网络
- 识别网络中的关键节点（hub）
（4）机器学习整合
- 将多组学特征作为输入，训练预测模型
- 使用特征选择识别最重要的预测因子

（1）分类任务：疾病预测
- 输入：微生物组组成特征（物种丰度、功能通路丰度）
- 输出：疾病状态（如IBD vs 健康）
- 常用算法：
- 随机森林（Random Forest）：对高维稀疏数据表现好，可输出特征重要性
- 支持向量机（SVM）：适合小样本高维数据
- 深度学习（神经网络）：可学习复杂的非线性关系
案例：使用随机森林预测CRC（结直肠癌）
- 训练集：数百个粪便样本的16S/宏基因组数据
- 特征：物种丰度（或MetaCyc通路丰度）
- 模型：随机森林分类器
- 结果：AUC可达0.8-0.9，某些研究报道准确率>90%
- 关键标志物：Fusobacterium nucleatum、Peptostreptococcus anaerobius等
（2）回归任务：表型预测
- 输入：微生物组特征 + 宿主因素
- 输出：连续变量（如BMI、药物响应、血糖水平）
- 常用算法：LASSO回归、弹性网络、随机森林回归
（3）聚类任务：亚型鉴定
- 输入：微生物组组成
- 输出：患者亚型（enterotypes，肠型）
- 肠型概念：人类肠道微生物组可分为几种主要的组成模式
- 肠型1：以Bacteroides为主
- 肠型2：以Prevotella为主
- 肠型3：以Ruminococcus为主

微生物组研究面临的一个核心挑战是因果关系的确定：
- 相关性不等于因果性：微生物组变化可能是疾病的原因，也可能是疾病的结果
- 混杂因素：饮食、药物、生活方式等同时影响微生物组和疾病状态
- 策略：
- 纵向研究：追踪疾病发生前后的微生物组变化
- 动物模型：无菌小鼠或抗生素处理后的定植实验
- 孟德尔随机化：利用宿主基因型作为工具变量
- 干预研究：益生菌、FMT等干预后的效果评估

（1）疾病诊断
- 结直肠癌（CRC）筛查：F. nucleatum等细菌标志物+粪便免疫化学检测（FIT）联合筛查
- 艰难梭菌感染（CDI）诊断：结合临床表现和微生物组特征
- 非酒精性脂肪性肝病（NAFLD）评估：肠道菌群组成与肝纤维化程度相关
（2）疾病治疗
- 复发性CDI：FMT治愈率>90%，已被FDA批准
- 溃疡性结肠炎（UC）：多项RCT显示FMT对UC有一定疗效
- 自闭症谱系障碍（ASD）：小规模临床试验显示FMT可改善胃肠道症状和行为
- 黑色素瘤免疫治疗：肠道菌群组成影响PD-1抑制剂的疗效
（3）药物代谢
- 肠道微生物组参与多种药物的代谢，影响药效和毒性
- 例如：伊立替康（irinotecan）的毒性由肠道细菌的β-葡萄糖醛酸酶介导
- 通过抑制特定细菌酶，可以降低药物毒性

（1）植物微生物组与农业
- 根际微生物组：影响植物养分吸收、抗病性和抗逆性
- 生物肥料：固氮菌、解磷菌等促进植物生长
- 生物防治：利用拮抗微生物防治植物病原菌
- 案例：印度矮生小麦与根际促生菌（PGPR）联合使用，可减少化肥用量30%
（2）动物微生物组与畜牧业
- 反刍动物瘤胃微生物组：决定饲料转化效率
- 益生菌：改善家禽和家畜的生长性能和免疫力
- 替代抗生素：欧盟禁用抗生素促生长剂后，微生物组产品成为替代方案

（1）生物修复
- 利用环境微生物组降解污染物
- 石油污染修复：海洋微生物组中的烃降解菌可降解石油烃
- 重金属污染修复：某些微生物可将有毒重金属转化为低毒形态
（2）功能基因挖掘
- 从环境宏基因组中发现新的酶和生物活性物质
- 纤维素酶：从瘤胃、白蚁肠道等环境中发现高效纤维素降解酶
- 新型抗生素：从土壤宏基因组中发现新结构抗生素（如Teixobactin）
（3）碳循环和气候变化
- 土壤微生物组在全球碳循环中发挥关键作用
- 冻土融化释放古老微生物，可能加速温室气体排放
- 海洋微生物组（特别是浮游植物共生菌）影响全球碳固定

（1）益生菌设计
- 传统益生菌：Lactobacillus, Bifidobacterium等
- 下一代益生菌（Next-Generation Probiotics）：Akkermansia muciniphila, Faecalibacterium prausnitzii等
- 工程化益生菌：通过基因工程赋予益生菌特定功能（如产丁酸盐、分泌治疗性蛋白）
（2）合成微生物组
- 目标：构建具有稳定功能的最小微生物群落
- 方法：
- 自上而下：从复杂群落中简化出核心功能群
- 自下而上：从单一菌株开始，逐步构建稳定群落
- 挑战：微生物间互作复杂，群落稳定性难以预测

（1）样本采集和保存
- 微生物组在采样后迅速变化（如氧气暴露导致厌氧菌死亡）
- 保存方法影响结果：冷冻、RNAlater、FTA卡等各有优缺点
- 建议：采样后尽快冷冻（-80°C），或使用稳定化试剂
（2）DNA提取偏差
- 不同试剂盒对不同类型细菌的提取效率差异很大
- 革兰氏阳性菌（厚壁菌门）细胞壁含大量肽聚糖，比革兰氏阴性菌更难裂解
- 策略：在同一研究中统一使用同一种提取方法；使用裂解对照（spike-in controls）
（3）PCR扩增偏差
- 通用引物对不同物种的扩增效率不同
- 嵌合体（chimera）和异源双链（heteroduplex）的形成
- 策略：使用高保真酶、优化PCR循环数、使用去嵌合体算法
（4）测序深度不足
- 复杂样本（如土壤）可能需要>10 Gb数据才能充分覆盖稀有物种
- 深度不足会导致稀有物种被遗漏，影响多样性估计
（5）生物信息学分析的异质性
- 不同分析流程（QIIME vs mothur）、不同数据库（SILVA vs Greengenes）会产生不同结果
- 缺乏统一的金标准

（1）样本量计算
- 微生物组研究通常需要比传统临床研究更大的样本量
- 因为个体差异大，效应量通常较小
- 建议：初步探索性研究每组至少20-30个样本；验证性研究每组50-100个
（2）对照设置
- 阴性对照（空白对照）：检测试剂污染
- 阳性对照（mock community）：评估技术流程的准确性
- 裂解对照（spike-in）：评估DNA提取效率
（3）批次控制
- 随机化：将不同组的样本随机分配到不同批次
- 平衡设计：确保每批次包含所有组的样本
- 记录所有批次信息，以便后续统计校正
- 统计校正：使用ComBat等方法校正批次效应
（4）纵向设计
- 对于疾病研究，纵向采样（疾病前-中-后）比横断面设计更能揭示因果关系
- 建议：基线（T0）、急性期（T1）、恢复期（T2）

标准化倡议：
- MIxS标准（Minimum Information about any (x) Sequence）：包括MIMARKS、MIMS等，规定了宏基因组研究需要报告的元数据
- Earth Microbiome Project（EMP）：制定了标准的采样和测序方案
数据共享平台：
- NCBI SRA/EBI ENA：原始测序数据
- Qiita/MG-RAST：处理后的微生物组数据
- 数据共享是提高研究可重复性和促进元分析的关键

（1）方法学进展
- 长读长宏基因组学：Nanopore/PacBio长读长有助于解决重复序列和分箱问题
- 单细胞宏基因组学：在单细胞分辨率上解析微生物群落
- 空间宏基因组学：保留微生物的空间分布信息
（2）人工智能与微生物组
- 深度学习用于物种分类和功能预测
- 生成式AI用于设计合成微生物组
- 大语言模型整合文献知识和微生物组数据
（3）精准微生物组医学
- 基于个体微生物组特征的个性化益生菌
- 微生物组作为疾病早期诊断的生物标志物
- 微生物组靶向治疗（如噬菌体疗法、CRISPR编辑肠道菌）
（4）全球微生物组计划
- 类似于人类基因组计划，全球微生物组计划旨在系统性地绘制地球微生物多样性图谱
- 应用：生物勘探、环境监测、气候变化预测

基因组学作为生物信息学的核心领域，在过去二十年中经历了前所未有的技术革命和方法学进步。从人类基因组计划的完成到个人基因组测序进入千元时代，从短读长测序主导到长读长技术崛起，基因组学的发展深刻改变了生命科学研究的面貌。

本章核心要点回顾

测序技术层面：Illumina短读长测序以其高准确性和高通量继续占据主流地位，而PacBio HiFi和Nanopore长读长测序则在复杂基因组组装和结构变异检测中展现出独特优势。混合测序策略正在成为高质量基因组项目的标准配置。
组装与注释层面：de Bruijn图算法是二代测序组装的核心方法，而OLC算法在长读长时代重新焕发活力。Hi-C技术使染色体级别组装成为可能。基因预测已从单纯的从头预测发展为整合多种证据的共识预测策略。功能注释依赖于同源性比对、结构域识别和通路映射等多种方法。
变异检测层面：SNP检测已达到极高的准确性，Indel检测需要特殊的比对处理，而SV检测仍是挑战性最大的领域。长读长测序为SV检测带来了突破性进展。
宏基因组学层面：宏基因组学通过bypass培养瓶颈，揭示了微生物世界的巨大多样性。从16S扩增子测序到宏基因组shotgun测序，从多样性分析到功能挖掘，宏基因组学正在从描述性研究向机制性研究和应用转化迈进。

未来展望

完整基因组时代：端粒到端粒（T2T）组装技术的成熟将使每个物种都能获得真正完整的基因组参考序列，包括着丝粒、端粒和核糖体DNA等 previously intractable 区域。
泛基因组学：单一参考基因组已不足以代表物种内的遗传多样性。泛基因组（pan-genome）研究将揭示核心基因组和可变基因组的动态演化，为作物改良和疾病研究提供新的视角。
单细胞基因组学：单细胞测序技术的进步将使研究者能够在单细胞分辨率上解析基因组变异，揭示肿瘤异质性、发育谱系和微生物群落结构。
人工智能与基因组学：深度学习已经在变异检测（DeepVariant）、蛋白质结构预测（AlphaFold）等领域取得突破。未来，AI将进一步整合基因组、转录组、蛋白质组和表观遗传组数据，实现从序列到功能的端到端预测。
精准基因组医学：随着测序成本的进一步下降和数据分析方法的成熟，全基因组测序有望成为临床常规检测手段，为罕见病诊断、肿瘤精准治疗和药物基因组学提供个体化的遗传信息。
基因组学的发展远未止步。正如人类基因组计划的领导者之一Francis Collins所言："基因组测序只是开始，真正的挑战和机遇在于解读这本生命之书。"作为生物信息学工作者，我们正处于这一伟大事业的核心，肩负着从海量基因组数据中提取生物学洞见、推动精准医学发展的历史使命。

本章思考题

1. 回顾基因组学发展的历史，你认为哪三项技术突破对领域的影响最大？为什么？
2. 如果你要启动一个全新物种的基因组项目，请设计一个完整的技术路线，包括测序策略、组装方案和注释流程。
3. 比较SNP、Indel和SV在检测难度、功能影响和进化意义方面的异同。
4. 宏基因组学如何改变了我们对微生物世界和宿主-微生物互作的理解？
5. 展望基因组学的未来，你认为下一个重大突破可能出现在哪个方向？

推荐阅读
1. Goodwin S, McPherson J D, McCombie W R. Coming of age: ten years of next-generation sequencing technologies [J]. Nature reviews genetics, 2016, 17: 333-351.
2. Nurk S, Koren S, Rhie A, et al. The complete sequence of a human genome [J]. Science, 2022, 376: 44-53. （T2T-CHM13）
3. Quince C, Walker A W, Simpson J T, et al. Shotgun metagenomics, from sampling to analysis [J]. Nature biotechnology, 2017, 35: 833-844.
4. Turner T N, Hormozdiari F, Duyzend M H, et al. Genome sequencing of autism-affected families reveals disruption of putative noncoding regulatory DNA [J]. American journal of human genetics, 2016, 98: 58-74.

（1）测序技术的三次革命
第一代测序（Sanger, 1977）：
- 原理：双脱氧链终止法
- 特点：准确性极高（>99.99%），读长可达1000 bp，但通量极低
- 贡献：完成首批模式生物基因组（噬菌体φX174、流感嗜血杆菌、酿酒酵母）
- 局限：人类基因组计划耗时13年、耗资27亿美元
第二代测序（Illumina, 2005）：
- 原理：桥式PCR + 可逆终止子边合成边测序
- 特点：通量提升百万倍，读长100-300 bp，准确性>99.9%，成本降至$1000/人基因组
- 贡献：千人基因组计划、肿瘤基因组图谱（TCGA）、宏基因组学兴起
- 局限：短读长无法跨越重复序列和复杂区域
第三代测序（PacBio/Nanopore, 2011-2015）：
- 原理：单分子实时测序（ZMW荧光检测 / 纳米孔电流检测）
- 特点：读长可达数Mb，无需PCR扩增，可检测碱基修饰
- 贡献：T2T基因组、复杂结构变异解析、直接RNA测序
- 局限：原始错误率较高（PacBio HiFi已解决），成本仍高于二代
（2）组装算法的演进
- 2001-2008：基于OLC的Sanger组装（Celera Assembler, Arachne）
- 2008-2015：基于DBG的二代短读长组装（Velvet, SOAPdenovo, SPAdes）
- 2015-2020：长读长OLC组装（Canu, FALCON）
- 2020至今：HiFi精准长读长组装（Hifiasm, Verkko）+ T2T整合策略
（3）变异检测的精度跃升
- 早期：基于pileup的统计方法（Samtools mpileup）
- 中期：单倍型模型（GATK HaplotypeCaller）
- 近期：深度学习（DeepVariant, Google Brain）——在多个基准测试中达到或超越人类专家水平

（1）数据层面的挑战
- 数据量爆炸：全球测序数据已超100 EB（Exabyte），存储和传输成本急剧上升
- 数据异质性：不同平台（Illumina/PacBio/Nanopore）、不同版本（hg19/hg38/T2T-CHM13）的数据整合困难
- 数据质量不均：公共数据库中大量低质量组装和错误注释存在
（2）分析层面的挑战
- 计算资源需求：T2T组装、泛基因组构建、单细胞分析需要PB级内存和万核级计算
- 算法瓶颈：高度重复序列、复杂结构变异、多倍体基因组仍是组装的难点
- 标准化缺失：不同分析流程的结果可比性差，缺乏统一的金标准
（3）解读层面的挑战
- 功能未知基因：即使在人类基因组中，仍有数千个基因功能不明
- 非编码区解读：增强子、沉默子、绝缘子等调控元件的精确预测和功能验证困难
- 变异临床意义：大量罕见变异的致病性难以判定（ACMG指南仍依赖专家判断）
- 多基因复杂疾病：大多数常见疾病（糖尿病、高血压、精神疾病）受数千个微效变异影响，单一基因解释力有限

（1）完整基因组时代（T2T for All）
- 目标：为地球上所有物种（>1000万种真核生物）提供T2T参考基因组
- 驱动力：地球生物基因组计划（EBP）、Zoonomia项目（哺乳类）、达尔文树计划（中国）
- 技术：HiFi + ONT + Hi-C的标准化流程，自动化组装和审校
- 意义：理解生物多样性、保护濒危物种、挖掘新资源
（2）泛基因组与图基因组
- 单一参考基因组已不能满足研究和临床应用的需求
- 图基因组将线性参考升级为包含群体变异的图结构，减少比对偏差
- 人类泛基因组参考联盟（HPRC）正在构建包含>40个多样化单倍型的泛基因组
- 应用：改进变异检测、发现参考基因组中缺失的序列、群体遗传学研究
（3）单细胞基因组学
- 从"bulk平均信号"到"单细胞分辨率"
- 单细胞基因组测序揭示肿瘤异质性、发育谱系、微生物群落结构
- 技术挑战：全基因组扩增的偏倚、低覆盖度下的变异检测、海量数据的分析
- 前沿方向：单细胞多组学（同时测基因组+转录组+表观组+蛋白质组）
（4）人工智能与基因组学
- 变异检测：DeepVariant、Clair3等深度学习工具已达到金标准水平
- 蛋白质结构预测：AlphaFold2革命性地解决了蛋白质折叠问题
- 调控元件预测：Enformer、Basenji等模型从序列预测基因表达和调控
- 端到端预测：从DNA序列直接预测功能（如SpliceAI预测剪接位点）
- 大语言模型：基于Transformer的基因组语言模型（如DNABERT、Nucleotide Transformer）正在学习DNA序列的"语义"
（5）精准基因组医学
- 罕见病诊断：全基因组测序（WGS）已成为未确诊遗传病的首选诊断工具，诊断率约25-40%
- 肿瘤精准治疗：基于ctDNA的液体活检、MRD（微小残留病灶）监测、个性化疫苗
- 药物基因组学：基于CYP2D6、HLA-B57:01等基因型指导用药，避免严重不良反应
- 新生儿筛查：扩展性基因组筛查（NICUSeq）在重症新生儿中快速诊断遗传病
- 预防医学：基于多基因风险评分（PRS）的疾病风险预测和早期干预
（6）合成基因组学*
- 从"读取"到"编写"：基因组合成成本从2003年的$4/bp降至2020年的$0.001/bp以下
- 最小基因组：Mycoplasma laboratorium（2010年，首个合成细胞）
- 酵母基因组合成计划（Sc2.0）：重新设计并合成16条酵母染色体
- 人类基因组的合成规划（GP-write）正在讨论中
- 应用：设计新型生物制造底盘、创造抗病毒/抗逆作物、开发活体疗法

中国在基因组学领域已从"跟跑"进入"并跑"甚至部分"领跑"阶段：
测序技术自主化：
- 华大智造（MGI）的DNBSEQ技术成为全球第二大测序平台
- 真迈生物、齐碳科技等在纳米孔测序领域取得突破
大科学计划：
- 中国十万人基因组计划
- 百万微生态计划
- 万种鸟类基因组计划（B10K中国部分）
- 地球生物基因组计划中国节点
数据基础设施建设：
- 国家基因组科学数据中心（NGDC）
- 中国国家生物信息中心（CNCB）
- 基因组数据安全存储和共享的国家标准

6.1.1 基因组学的概念

基因组学研究的范畴可以从三个层次理解：
（1）结构基因组学（Structural Genomics）
结构基因组学关注基因组的物理结构，包括DNA序列的测定、基因组图谱的构建、基因的定位和结构分析等。其核心目标是获得基因组的完整序列，并理解基因在基因组中的排列方式。
（2）功能基因组学（Functional Genomics）
功能基因组学研究基因组中基因的功能、基因之间的相互作用、基因表达调控等。与结构基因组学回答"是什么"不同，功能基因组学回答"做什么"和"怎么做"。主要技术包括RNA-Seq、ChIP-Seq、ATAC-Seq等。
（3）比较基因组学（Comparative Genomics）
比较基因组学通过比较不同物种或不同个体基因组的异同，揭示基因组的进化规律、基因的起源和功能分化、物种之间的亲缘关系等。

基因组学中所研究的问题主要分为三类：
第一类：如何获得基因组序列
这是基因组学最基础的问题。随着DNA测序技术的发展，从一代Sanger测序到二代高通量测序（NGS），再到三代单分子长读长测序，获得基因组序列的成本急剧下降、速度大幅提升。一个物种的基因组项目通常包括：基因组大小估计、测序策略设计、序列组装、质量评估等步骤。
第二类：如何解读/解码基因组
获得序列后，需要解读这本"生命之书"。这包括：基因预测与结构注释（识别编码区、内含子、外显子、调控元件等）、功能注释（确定基因的功能、参与的通路）、比较分析（与近缘物种比较，理解进化关系）、变异分析（识别SNP、Indel、SV等）。
第三类：如何重写/编写新的基因组
这是合成生物学的前沿方向。通过基因编辑技术（如CRISPR-Cas9）对现有基因组进行精确修饰，甚至从头设计和合成全新的基因组。2010年，Craig Venter团队首次合成了支原体的人工基因组，标志着人类开始具备"编写生命"的能力。

人类基因组计划采用了层次测序策略：通过构建DNA大片段克隆文库、重叠克隆群和物理图谱，再进行逐个克隆鸟枪法测序。中国科学家承担了1%的测序任务（3号染色体端部约3000万个碱基对），成功跟进世界先进测序技术，推动了国内测序仪研发和生物信息学软件开发。
HGP的完成不仅提供了人类基因组的参考序列，更重要的是：
1. 技术推动：测序技术从Sanger法发展到高通量测序
2. 数据共享：建立了基因组数据开放共享的文化
3. 疾病研究：为疾病基因定位、药物基因组学奠定基础
4. 进化研究：为理解人类起源和进化提供基础数据

Sanger测序的核心是利用ddNTP终止DNA链延伸：
1. 反应体系：模板DNA、引物、DNA聚合酶、四种dNTP（其中一种带有放射性或荧光标记）和少量ddNTP
2. 链延伸与终止：DNA聚合酶沿模板合成新链，遇到dNTP时正常延伸，遇到ddNTP时因缺少3'-OH而终止
3. 四个独立反应：分别加入ddATP、ddCTP、ddGTP、ddTTP，产生四组不同长度的片段
4. 电泳分离：通过聚丙烯酰胺凝胶电泳（PAGE）或毛细管电泳按大小分离片段
5. 读取序列：从短到长读取末端碱基，获得互补链序列
自动化Sanger测序使用四种不同荧光标记的ddNTP，四个反应在同一管中进行，通过激光检测荧光颜色确定碱基类型。
局限性：通量低（每次反应只能读取一条序列），成本高，不适合大规模基因组项目。

Illumina测序流程包括四个关键步骤：
（1）文库制备
- 将基因组DNA随机打断成200-500 bp片段
- 末端修复、加A尾、连接测序接头（adapter）
- PCR扩增文库（可选）
（2）Flow Cell杂交与桥式PCR
- Flow Cell表面固定有两种与接头互补的引物
- 单链DNA文库与引物杂交
- 通过桥式PCR扩增：DNA片段弯曲成桥状，以另一端为模板进行扩增
- 经过多轮扩增，每个原始分子形成一个约1000-2000个克隆的"簇"（cluster）
- 变性后得到单链簇，作为测序模板
（3）边合成边测序（SBS）
- 加入DNA聚合酶和四种带有不同荧光标记的dNTP（每种dNTP的3'-OH被可切除的化学基团阻断）
- 每个循环：
a. 四种dNTP竞争掺入，每次只有一个碱基被加入
b. 激光扫描，根据荧光颜色确定掺入的碱基
c. 切除荧光基团和3'阻断基团，使链可以继续延伸
- 重复上述循环，逐碱基读取序列
（4）双端测序（Paired-End Sequencing）
- 从簇的两端分别测序，获得两段各约150 bp的序列
- 两段序列之间距离已知（约200-500 bp）
- 双端信息对于基因组组装和结构变异检测至关重要
数据格式与质量控制
测序数据以FASTQ格式存储，每条记录包含四行：序列标识符、碱基序列、"+"分隔符、碱基质量值（Phred分数）。Phred质量值Q与碱基错误率ε的关系为：Q = -10log₁₀(ε)。Q30表示碱基错误率约为0.1%，即准确性为99.9%。

（1）PacBio SMRT测序
PacBio测序的核心是零模波导孔（ZMW）技术：
- 在芯片上蚀刻数十万个直径约70 nm、深度约100 nm的纳米孔
- 每个孔底部固定一个DNA聚合酶分子
- 待测DNA与引物结合，被聚合酶捕获
- 加入四种带有不同荧光标记的dNTP（标记在磷酸基团上，而非碱基上）
- 当dNTP被掺入新链时，荧光标记被释放到溶液中，产生闪光信号
- 记录每个ZMW中的闪光序列和时间，即可获得DNA序列
PacBio Sequel II可产生约15-20 kb的长读长，最新的Revio系统进一步提高了通量。通过环化一致性测序（CCS，Circular Consensus Sequencing），同一个分子多次循环测序，可获得准确性>99.9%的HiFi reads。
（2）Nanopore测序
Nanopore测序的核心原理：
- 蛋白质纳米孔嵌入脂质双分子层膜中
- 在膜两侧施加电压，产生离子电流
- DNA/RNA分子在动力蛋白（helicase）驱动下穿过纳米孔
- 不同碱基（或k-mer）穿过纳米孔时引起特征性的电流变化
- 通过深度学习模型（如Bonito、Dorado）将电流信号解码为碱基序列
Nanopore的优势：
- 超长读长：目前最长记录超过4 Mb
- 便携性：MinION设备仅U盘大小，可在野外、太空等极端环境使用
- 直接测序RNA：无需逆转录，可检测RNA修饰
- 实时测序：数据边测边出，可即时分析
局限性：原始读取准确性较低（约92-97%），但R10.4.1 flow cell配合Dorado basecaller的准确性已大幅提升。

由于二代测序准确性高但读长短，三代测序读长长但原始准确性低，混合测序（Hybrid-Seq）成为当前基因组项目的标准配置：
- 用三代长读长进行基因组骨架组装
- 用二代短读长对组装结果进行纠错和填补缺口
- Hi-C数据辅助染色体级别组装

传统的基因组组装在高度重复序列区域（如着丝粒、端粒、rDNA阵列）存在大量缺口，主要原因是：
- 二代短读长无法跨越长重复单元
- 重复序列使组装算法难以确定正确的路径
T2T联盟实现人类完整基因组的关键技术组合：
（1）PacBio HiFi长读长
- 读长约15-20 kb，准确性>99.9%
- 足以跨越大多数重复单元（如着丝粒卫星重复单元约171 bp）
- 高准确性使得区分高度相似的重复序列成为可能
（2）Oxford Nanopore超长读长
- 读长可达100 kb至数Mb
- 直接跨越最大的重复阵列（如5S rDNA阵列可长达数Mb）
- 提供长距离的序列连续性信息
（3）Hi-C辅助scaffold连接
- 将HiFi组装产生的contig提升到染色体级别
- 确定contig的顺序、方向和间隙大小
（4）手动审校与验证
- 对复杂区域（如着丝粒）进行手动检查
- 使用荧光原位杂交（FISH）和光学图谱进行验证

Hi-C实验流程：
1. 甲醛交联：将空间上接近的DNA片段共价交联
2. 酶切：用限制性内切酶（如HindIII、DpnII）切割DNA
3. 末端修复和生物素标记：补平粘性末端并掺入生物素标记的dNTP
4. 连接：稀释后连接，使空间上接近的片段连接在一起
5. 去交联和测序：打断DNA，富集生物素标记的嵌合片段，进行双端测序
Hi-C数据分析用于基因组组装：
- 构建接触矩阵（contact matrix）：染色体各区域之间的相互作用频率
- 同一染色体上的区域相互作用频率高（呈现对角线富集）
- 相邻区域相互作用频率高于远端区域
- 使用算法（如3D-DNA、SALSA、ALLHiC）将contig聚类到染色体组、确定顺序和方向

泛基因组分为两部分：
- 核心基因组（Core Genome）：物种所有个体共有的基因，通常与基本生命功能相关
- 可变基因组/附属基因组（Accessory/Dispensable Genome）：仅在部分个体中存在的基因，通常与环境适应、致病性、代谢多样性等相关
泛基因组构建方法：
- 迭代式：将新个体的基因组与已有泛基因组比较，添加新序列
- 图基因组（Graph Genome）：将所有序列变异表示为图结构，而非线性序列
泛基因组的应用：
- 捕获参考基因组中缺失的遗传多样性
- 改进变异检测（减少比对偏差）
- 识别与重要性状相关的基因（如作物抗病基因）

一个典型的T2T基因组项目测序策略：
| 数据类型 | 测序深度 | 用途 |
|----------|---------|------|
| PacBio HiFi | 30-50× | 主要组装，产生高质量contig |
| Nanopore Ultra-long | 10-20× | 跨越超长重复区域，解决复杂区域 |
| Illumina PE150 | 50-100× | 纠错、质量验证 |
| Hi-C | 60-100× | 染色体级别scaffold |
| 光学图谱 | 5-10× | 验证组装结构 |

（1）线性序列展示
- 以染色体为坐标轴，从左到右展示DNA序列
- 支持不同尺度的缩放（从单碱基到全染色体）
- 显示正反链的序列和注释
（2）多轨道（Track）整合
- 每个track代表一类数据（基因注释、测序reads、变异、保守性等）
- 用户可自定义加载、排列和配置track
- 支持track的叠加和比较
（3）交互式查询
- 通过基因名、染色体坐标、rsID等定位基因组区域
- 支持序列搜索（如查找特定基序）
- 支持不同物种间的基因组比对（synteny view）
（4）数据导出
- 导出特定区域的序列（FASTA）
- 导出track数据（BED、GTF等）
- 生成高质量图片用于发表

（1）UCSC Genome Browser
- 最经典的基因组浏览器之一
- 整合了大量注释数据和工具
- 支持自定义track加载（通过URL或本地文件）
- 强大的Table Browser用于数据查询和下载
（2）Ensembl Genome Browser
- 欧洲生物信息研究所（EBI）维护
- 注释质量高，更新及时
- 强大的比较基因组学工具（如Ensembl Compara）
- 提供API和BioMart数据查询工具
（3）IGV（Integrative Genomics Viewer）
- 本地运行的桌面软件
- 支持大规模数据（如BAM、VCF、Hi-C矩阵）
- 特别适合查看测序reads比对情况和变异位点
- 支持基因组重排和结构变异的可视化
（4）JBrowse/GBrowse
- 可部署在本地服务器的网页浏览器
- 适合构建物种特异性的基因组数据库
- 支持插件扩展

BED格式
BED文件至少包含3列，最多12列：
- chrom：染色体名称
- chromStart：起始位置（0-based）
- chromEnd：终止位置（1-based，不包含）
- name：名称（可选）
- score：分数（0-1000，可选）
- strand：方向（+/-/.）
- thickStart：编码区起始（可选，用于显示）
- thickEnd：编码区终止（可选）
- itemRgb：颜色（可选）
- blockCount：外显子数量（BED12）
- blockSizes：外显子大小列表（逗号分隔）
- blockStarts：外显子相对起始位置列表
WIG格式
variableStep格式示例：

variableStep chrom=chr1 span=50
1001 12.5
1051 13.2
1101 11.8

表示chr1上1001-1050位置的值为12.5，1051-1100位置的值为13.2。
GTF格式
GTF是GFF2的变体，每行9列：
1. seqname：染色体或scaffold
2. source：注释来源
3. feature：特征类型（gene, transcript, exon, CDS等）
4. start：起始位置（1-based）
5. end：终止位置（1-based，包含）
6. score：分数（通常为".")
7. strand：方向（+/-/.)
8. frame：阅读框（0/1/2/.）
9. attributes：属性（gene_id, transcript_id等）

现代基因组研究往往涉及多组学数据的整合：
- 基因组层：基因结构、变异位点
- 转录组层：RNA-Seq信号、剪接异构体
- 表观组层：ChIP-Seq峰（组蛋白修饰、转录因子结合）、DNA甲基化、ATAC-Seq开放染色质
- 三维结构层：Hi-C相互作用矩阵、TAD边界
- 进化层：跨物种保守性分数（phyloP, phastCons）
整合可视化可以帮助发现跨层面的关联，例如：
- 某个增强子区域（H3K27ac峰+ATAC-Seq开放区）与远处基因的启动子有Hi-C相互作用
- 一个SNP位于保守非编码区，且与附近基因的表达量显著相关

6.2.1 基因组的组装

（1）重复序列
基因组中普遍存在重复序列，包括串联重复（如微卫星）和散在重复（如转座子、segmental duplication）。重复序列导致组装歧义——相同的reads可能来源于基因组的不同位置，组装算法无法确定其真实来源。
（2）测序错误
Illumina测序错误率约0.1-1%（Q30-Q20），Nanopore原始错误率可达5-15%。测序错误会干扰k-mer的正确连接，在DBG中产生"气泡"（bubble）和"死胡同"（tip）。
（3）覆盖度不均
基因组不同区域的测序深度差异很大。高GC区域、AT富集区域、重复序列区域的覆盖度往往偏低，导致这些区域组装困难。
（4）杂合度
二倍体生物的两条同源染色体存在差异。组装软件需要区分同源序列并可能分别组装（phased assembly），否则杂合位点会导致组装碎片化。

DBG算法是二代短读长组装的核心方法，其工作流程如下：
（1）k-mer分解
将所有reads拆分为长度为k的短片段（k-mer）。例如，序列"ATGCGGAT"在k=5时产生：ATGCG、TGCGG、GCGGA、CGGAT。通常k取奇数（如21、31、51、71），且不超过read长度的2/3。
（2）构建DBG
- 每个k-mer作为图中的一个顶点（node）
- 若两个k-mer之间有k-1个碱基重叠，则在它们之间建立一条有向边
- 例如，ATGCG和TGCGG共享TGCG（4个碱基），建立从ATGCG到TGCGG的边
- 将不同reads产生的相同k-mer合并为一个顶点，记录覆盖度
（3）图简化
测序错误会在DBG中产生两类错误结构：
- Tips（死胡同）：由于测序错误产生的低频k-mer连接到主路径后又很快终止。处理：去除长度和覆盖度低于阈值的tips。
- Bubbles（气泡）：由于杂合位点或测序错误，主路径在某处分叉为两条路径，在下游又汇合。处理：保留覆盖度较高的路径，或在杂合组装中保留两条路径。
（4）contig提取
在简化后的DBG中，一条连续的路径（从没有入度的顶点开始到没有出度的顶点结束）对应一个contig。
k-mer大小的选择：
- k值较小：对重复序列的分辨率低（不同重复单元可能共享k-mer），但覆盖度高（每个read产生更多k-mer）
- k值较大：对重复序列的分辨率高，但覆盖度低（要求更大的测序深度才能保证每个k-mer被覆盖）
- 多k-mer策略（如SPAdes）：从较小的k开始组装，逐步增大k值，利用小k的连续性和大k的准确性优势

OLC算法更适合长读长数据（PacBio/Nanopore），其工作流程：
（1）Overlap（重叠检测）
计算所有reads两两之间的重叠关系。对于N条reads，需要O(N²)次比对，计算量巨大。对于长读长数据，通常使用minimap2等快速比对工具进行all-vs-all比对。
（2）Layout（布局构建）
根据重叠关系构建布局图（string graph）：
- 每个read是图中的一个节点
- 如果两个reads有显著重叠，则在它们之间建立边
- 边的方向和长度由重叠的位置和长度确定
- 简化图结构（去除 transitive edges，压缩 unambiguous paths）
（3）Consensus（一致性序列生成）
对布局图中的每个unitig（unambiguous path），将对应的reads进行多序列比对，生成一致性序列作为contig。

利用二代和三代数据的互补性：
（1）二代纠错三代
- 用高准确性的二代reads对三代reads进行纠错
- 代表工具：PBcR、LoRDEC、Jabba
- 纠错后的三代reads准确性可提升至>99%
（2）三代搭建骨架，二代填充
- 用三代长读长构建初步的contig骨架
- 用二代reads对三代contig进行 polish（填补gap、纠正碱基错误）
- 代表工具：Pilon、Racon
（3）graph-based混合组装
- 将二代和三代数据同时输入组装算法
- 在组装图中同时利用短读长的准确性和长读长的连续性
- 代表工具：HybridSPAdes、MaSuRCA

连续性指标：
- N50/L50：衡量contig/scaffold长度的分布
- 最大contig长度：最长的连续序列
- gap数量：scaffold中N碱基的数量和总长度
完整性指标：
- BUSCO（Benchmarking Universal Single-Copy Orthologs）：检测基因组中保守单拷贝基因的存在情况。高BUSCO完整性（如>90%）表示基因组组装较完整。
- CEGMA（Core Eukaryotic Genes Mapping Approach）：类似BUSCO，基于核心真核基因集。
准确性指标：
- QUAST：计算N50、错误组装（misassembly）数量、基因组覆盖度等
- Pilon polish后的二代比对一致性：检测碱基错误
- k-mer spectrum分析：与原始测序数据的k-mer谱比较，评估完整性

原核生物基因特征：
- 基因密度高（大肠杆菌基因组约85%为编码区）
- 无内含子（少数古菌有内含子）
- 无5'端帽子和3'端polyA尾
n- 操纵子结构：多个基因共用一个启动子，转录为多顺反子mRNA
- 起始密码子通常为ATG（少数GTG、TTG）
- 启动子元件：-10区（Pribnow框：TATAAT）和-35区（TTGACA）
真核生物基因特征：
- 基因密度低（人类基因组仅约1-2%为编码区）
- 含内含子（人类基因平均含8个外显子、7个内含子）
- 5'端有帽子结构，3'端有polyA尾
- 单顺反子mRNA（一个转录本对应一个蛋白质）
- 复杂的剪接机制：内含子5'端为GT（供体位点），3'端为AG（受体位点）
- 核心启动子元件：TATA框（约-30区）、CAAT框（约-80区）
- 存在大量的可变剪接（人类约95%的多外显子基因发生可变剪接）

从头预测的核心是利用基因组的统计特征和信号特征：
（1）开放阅读框（ORF）检测
- 寻找起始于ATG、终止于TAA/TAG/TGA的连续密码子序列
- 计算ORF长度：真核生物外显子通常>300 bp（100个氨基酸），短ORF可能是假阳性
- 计算密码子使用频率（Codon Usage Frequency, CUF）：真实基因通常偏好使用特定密码子
（2）隐马尔可夫模型（HMM）
HMM是基因预测中最常用的统计框架：
- 隐藏状态：基因间区、外显子（起始/中间/终止）、内含子（起始/中间/终止）、UTR等
- 观测符号：A、C、G、T
- 状态转移概率：从一个状态转移到另一个状态的概率（如外显子→内含子的概率）
- 发射概率：在某个状态下观察到特定碱基的概率
- Viterbi算法：寻找最可能的隐藏状态序列（即基因结构）
经典HMM基因预测软件：
- Genscan：针对人类和脊椎动物，基于广义HMM
- Augustus：支持真核和原核生物，可整合外部证据
- Glimmer/GlimmerHMM：Glimmer用于原核，GlimmerHMM用于真核
- GeneMark/GeneMark-ES：自训练HMM，适合新物种
- SNAP：可训练参数，适合非模式生物
（3）支持向量机（SVM）和深度学习方法
近年来，机器学习方法也被用于基因预测：
- 将序列特征（k-mer频率、GC含量、密码子偏好等）编码为特征向量
- 用已知基因训练分类器（SVM、随机森林、神经网络）
- 对新序列进行预测

（1）蛋白质同源比对
- 将已知物种的蛋白质序列与目标基因组进行比对（如tblastn）
- 根据蛋白质比对结果推断基因的编码区
- 考虑剪接位点信号，将基因组DNA与外显子对应
- 代表工具：GeneWise、Exonerate（protein2genome模式）
（2）cDNA/EST比对
- 将已知的cDNA或EST序列与基因组比对（如BLAT、GMAP）
- 直接确定外显子-内含子边界
- EST数据质量参差不齐，需要过滤和聚类
（3）跨物种比较基因组学
- 利用进化保守性：编码区通常比非编码区更保守
- 利用密码子替换模式：同义替换率（Ks）和非同义替换率（Ka）
- 代表方法：Twinscan（整合人类-小鼠比较信息）、CONTRAST

（1）RNA-Seq数据比对与组装
- 将RNA-Seq reads比对到基因组（如HISAT2、STAR）
- 使用StringTie、Cufflinks等工具组装转录本
- 获得GTF格式的转录本结构注释
（2）全长转录本测序
- PacBio Iso-Seq或Nanopore direct RNA测序获得全长cDNA
- 无需组装，直接获得完整的外显子-内含子结构
- 可发现新的转录本异构体
（3）整合策略（Evidence Modeler, EVM）
- 将从不同来源（从头预测、同源比对、转录本组装）获得的基因模型进行整合
- 根据证据的可靠性和一致性进行加权评分
- 输出一致性最高的基因模型集合
- 代表工具：EVM（Evidence Modeler）、PASA（Program to Assemble Spliced Alignments）、BRAKER

非编码RNA（ncRNA）是不编码蛋白质的RNA分子，包括：
- miRNA（微小RNA）：~22 nt，调控基因表达
- lncRNA（长链非编码RNA）：>200 nt，多种调控功能
- circRNA（环状RNA）：共价闭合环，miRNA海绵等
- tRNA、rRNA、snRNA、snoRNA等
预测方法：
- miRNA：基于发夹结构、种子区域保守性、热力学稳定性（RNAfold、miRDeep2）
- lncRNA：基于开放阅读框长度、编码潜能（CPAT、CPC2）、序列保守性
- circRNA：基于反向剪接位点（backsplice junction）的reads支持（CIRI、find_circ）

GO数据库的核心是一个有向无环图（DAG）结构，其中：
- 每个节点（GO term）代表一个功能概念
- 边代表"is_a"（是一种）或"part_of"（是...的一部分）关系
- 子节点比父节点更具体（如"蛋白质激酶活性" is_a "转移酶活性"）
GO注释的常用方法：
（1）基于序列相似性的GO注释
- 用BLAST将查询蛋白与GO注释数据库（如UniProt-GOA）比对
- 将相似蛋白的GO term转移给查询蛋白
- 常用工具：Blast2GO、eggNOG-mapper
（2）基于结构域的GO注释
- 通过InterProScan识别蛋白质结构域
- 利用结构域与GO term的映射关系进行注释
- InterPro数据库提供了Pfam、SUPERFAMILY等结构与GO的映射
（3）基于机器学习的方法
- 使用蛋白序列特征（氨基酸组成、物理化学性质、亚细胞定位信号等）训练分类器
- 预测GO term
- 代表方法：DeepGO（基于深度学习）

KEGG数据库的核心是通路图（Pathway Map），包括：
- 代谢通路（Metabolism）：碳水化合物代谢、能量代谢、脂质代谢、核苷酸代谢、氨基酸代谢等
- 遗传信息处理（Genetic Information Processing）：转录、翻译、折叠/分选/降解、复制与修复
- 环境信息处理（Environmental Information Processing）：膜转运、信号转导、配体-受体相互作用
- 细胞过程（Cellular Processes）：细胞生长与死亡、细胞群落、运输与分解代谢
- 生物体系统（Organismal Systems）：免疫系统、神经系统、内分泌系统等
- 人类疾病（Human Diseases）：癌症、神经退行性疾病、心血管疾病等
KEGG注释方法：
- 序列比对法：用BLAST/KAAS（KEGG Automatic Annotation Server）将基因与KEGG GENES数据库比对
- KO（KEGG Orthology）注释：将基因映射到KO编号（功能单位），再根据KO编号映射到通路
- 工具：KAAS、KOFAM、eggNOG-mapper

InterPro整合了多个蛋白质特征数据库：
- Pfam：基于HMM的蛋白质家族和结构域数据库
- PRINTS：基于指纹（保守基序组）的数据库
- PANTHER：基于系统发育树的蛋白质家族分类
- ProSite：基于模式（pattern）和谱（profile）的功能位点数据库
- SUPERFAMILY：基于SCOP结构分类的远缘同源识别
- SMART：信号肽、跨膜区、结构域识别
InterProScan工作流程：
1. 用多种方法（HMM、模式匹配、profile等）扫描蛋白质序列
2. 整合不同数据库的预测结果
3. 为每个预测的domain/feature提供GO term映射
4. 输出GFF3或TSV格式的注释结果

（1）COG/KOG/EggNOG
- COG（Clusters of Orthologous Groups）：原核生物蛋白质直系同源簇
- KOG：真核生物版本
- EggNOG：扩展版本，覆盖更多物种，提供GO、KEGG、CAZy等注释
- 将蛋白质分为功能类别（如J: Translation, K: Transcription, S: Unknown）
（2）NR/UniProt/Swiss-Prot
- NR（Non-Redundant）：NCBI的非冗余蛋白数据库，用于序列相似性搜索
- UniProt：最全面的蛋白质数据库，包含Swiss-Prot（人工审编）和TrEMBL（自动注释）
- 注释流程：BLASTp比对NR/UniProt → 取top hit → 继承功能描述
（3）专用数据库
- CAZy：糖类活性酶分类（糖苷水解酶、糖基转移酶等）
- TCDB：膜转运蛋白分类
- VFDB：细菌毒力因子
- CARD：抗生素耐药基因
- antiSMASH：生物合成基因簇（次级代谢产物）

现代基因组项目通常使用自动化流程进行功能注释：

预测蛋白质序列
    ↓
1. InterProScan（结构域识别 + GO注释）
    ↓
2. KAAS/KOFAM（KEGG通路注释）
    ↓
3. eggNOG-mapper（COG/GO/KEGG综合注释）
    ↓
4. BLASTp against NR/Swiss-Prot（功能描述）
    ↓
5. 整合所有注释结果

6.3.1 基本原理

基因组变异按尺度可分为：
点突变尺度：
- SNP（Single Nucleotide Polymorphism）：单个碱基的替换，是最常见的变异类型
- Indel（Insertion/Deletion）：小片段（通常<50 bp）的插入或缺失
结构变异尺度：
- 缺失（Deletion）：>50 bp的序列缺失
- 插入（Insertion）：>50 bp的序列插入
- 倒位（Inversion）：序列方向的翻转
- 易位（Translocation）：序列从一个染色体位置移动到另一个位置
- 重复（Duplication）：序列拷贝数的增加
- CNV（Copy Number Variation）：>1 kb的拷贝数变化

原始测序数据 (FASTQ)
    ↓
质量控制 (FastQC, Trimmomatic)
    ↓
比对到参考基因组 (BWA, Bowtie2)
    ↓
比对后处理 (排序, 标记重复, 碱基质量校正)
    ↓
变异检测 (GATK, Samtools, FreeBayes)
    ↓
变异过滤和质量控制
    ↓
变异注释 (功能影响预测)
    ↓
最终变异集合 (VCF格式)

（1）排序（Sorting）
比对后的SAM/BAM文件需要按染色体坐标排序，以便后续处理。
（2）标记PCR重复（Mark Duplicates）
PCR扩增会引入系统性偏差——来自同一模板分子的多个reads并非独立采样。使用Picard MarkDuplicates等工具识别并标记PCR重复reads，避免重复reads对变异检测的干扰。
（3）Indel区域重比对（Indel Realignment）
比对软件在Indel附近可能产生比对偏差（如将Indel附近的SNP错误比对）。GATK的IndelRealigner通过识别Indel附近不一致的比对，进行局部重比对，提高Indel和周围SNP的检测准确性。
（4）碱基质量校正（Base Quality Score Recalibration, BQSR）
测序仪给出的原始碱基质量值可能存在系统性偏差（如某些测序循环质量系统性地偏低）。BQSR通过将碱基质量与已知变异数据库（如dbSNP）进行比较，建立校正模型，重新校准碱基质量值。

变异检测的核心是判断某个位置的碱基差异是真实的遗传变异还是测序/比对错误。这本质上是一个统计假设检验问题。
零假设（H₀）：该位置没有变异，所有观察到的差异都是测序错误
备择假设（H₁）：该位置存在变异
测序错误通常被建模为二项分布或贝叶斯模型：
- 设测序错误率为ε（通常ε≈0.001，对应Q30）
- 在覆盖度为D的位置，观察到k个非参考碱基
- 如果H₀成立，k ~ Binomial(D, ε)
- 计算P(k | H₀)，如果该P值小于阈值（如0.05），则拒绝H₀，判定为变异
基因型似然模型：
更精确的方法是对三种可能的基因型（参考纯合RR、杂合RA、变异纯合AA）分别计算似然值：
- L(RR) = ∏ P(read_i | RR)
- L(RA) = ∏ P(read_i | RA)
- L(AA) = ∏ P(read_i | AA)
选择似然值最大的基因型作为该位点的基因型。
GATK HaplotypeCaller使用更先进的单倍型模型，在局部区域（通常数百bp）同时考虑所有可能的变异，通过局部de novo组装生成候选单倍型，然后计算每个单倍型的似然值。

即使使用了严格的统计模型，变异检测结果中仍包含大量假阳性。需要通过多重过滤来提高准确性：
基于测序质量的过滤：
- QD（Quality by Depth）：质量值除以覆盖度，过低提示质量不可靠
- FS（Fisher Strand）：检验正负链reads是否均衡支持变异，不均衡提示比对偏差
- MQ（Mapping Quality）：比对质量值，过低提示比对不可靠
- ReadPosRankSum：检验变异是否主要出现在reads末端（末端碱基质量通常较低）
基于覆盖度的过滤：
- DP（Depth）：覆盖度过低（<5-10×）的变异不可靠
- 覆盖度过高（>平均深度的2-3倍）可能提示重复区域
基于群体频率的过滤：
- 使用已知变异数据库（如gnomAD）过滤常见变异
- 在罕见病研究中，过滤在人群中频率>1%的变异

SNP主要来源于DNA复制过程中的自发突变：
（1）脱氨基作用
- 胞嘧啶（C）自发脱氨基转变为尿嘧啶（U），DNA修复系统将U识别为T，导致C→T突变
- 这是转换突变的主要来源，解释了为什么转换/颠换比（Ti/Tv ratio）通常大于2
- 5-甲基胞嘧啶（5mC）的脱氨基率更高，导致CpG岛中的C→T突变频率显著升高
（2）复制错误
- DNA聚合酶在复制过程中偶尔掺入错误的碱基
- 尽管有校对功能，错误率仍约为10⁻⁹/碱基/复制
（3）氧化损伤
- 活性氧（ROS）可导致碱基氧化修饰，如8-oxo-guanine可导致G→T颠换

按功能位置分类：
| 位置 | 比例 | 功能影响 |
|------|------|---------|
| 基因间区 | ~90% | 通常无直接功能影响 |
| 内含子 | ~5% | 可能影响剪接调控 |
| 外显子同义 | ~1% | 不改变氨基酸（密码子简并性） |
| 外显子错义 | ~1.5% | 改变氨基酸序列 |
| 外显子无义 | ~0.1% | 引入提前终止密码子 |
| 调控区 | ~2% | 可能影响基因表达 |
| 剪接位点 | <0.1% | 可能导致剪接异常 |
按群体频率分类：
- 常见变异（Common Variant）：群体频率>5%，通常与复杂疾病的易感风险相关
- 低频变异（Low-frequency Variant）：群体频率1-5%
- 稀有变异（Rare Variant）：群体频率<1%，可能与孟德尔遗传病相关

（1）基于pileup的方法（Samtools/bcftools）
- 统计每个基因组位置上比对reads的碱基组成
- 使用贝叶斯模型计算每个可能基因型的后验概率
- 优点：速度快，计算资源需求低
- 缺点：不处理Indel附近的比对偏差，对Indel不敏感
（2）基于单倍型的方法（GATK HaplotypeCaller）
- 在局部窗口（数百bp）内进行de novo组装，生成候选单倍型
- 将reads与候选单倍型进行比对，计算每个单倍型的似然值
- 选择最可能的单倍型组合作为基因型
- 优点：对Indel和复杂区域的SNP检测更准确
- 缺点：计算量大，耗时较长
（3）基于深度学习的方法（DeepVariant）
- 将每个候选变异位点转化为"图像"（read pileup的视觉表示）
- 使用卷积神经网络（CNN）进行变异识别
- 优点：准确性极高，在多个基准测试中优于传统方法
- 缺点：计算资源需求高，训练数据依赖

QD（Quality by Depth）：
QD = QUAL / DP
- QD过低（<2）提示变异可能被低质量reads或重复区域支持
- QD过高（>40）可能提示比对到重复区域
MQ（Mapping Quality）：
- 反映reads在该区域的比对可靠性
- MQ < 40提示比对可能存在问题
FS（Fisher Strand Bias）：
- 检验变异是否被正负链均衡支持
- FS > 60提示严重的链偏倚，变异可能不可靠
SOR（Strand Odds Ratio）：
- FS的替代指标，对覆盖度不均更稳健
- SOR > 3提示链偏倚
MQRankSum：
- 检验支持变异的reads与支持参考的reads的比对质量是否存在显著差异
- 负值过大提示支持变异的reads比对质量较差
ReadPosRankSum：
- 检验变异是否主要出现在reads末端
- 末端碱基质量通常较低，可能导致假阳性

转换/颠换比率（Ti/Tv ratio）是评估SNP集合质量的重要指标：
- 全基因组范围的期望Ti/Tv：约2.0-2.2
- 外显子区域的期望Ti/Tv：约2.8-3.0（因为同义SNP只发生在特定密码子位置）
如果检测到的SNP集合Ti/Tv显著低于期望值，可能提示：
- 存在大量假阳性（特别是颠换类型的假阳性较多）
- 过滤条件不够严格
- 样本存在污染

（1）DNA聚合酶滑移（Replication Slippage/Strand Slippage）
- 在重复序列区域（尤其是微卫星），模板链和新生链可能发生相对滑动
- 导致重复单元的插入或缺失
- 这是Indel最主要的产生机制，解释了为什么Indel在重复序列区域富集
（2）不等交换（Unequal Crossing Over）
- 减数分裂过程中的同源重组发生错位
- 导致一个染色体上重复序列增加，另一个上减少
（3）DNA损伤修复错误
- 双链断裂修复过程中可能发生小片段的插入或缺失

Indel检测比SNP检测困难的主要原因：
（1）比对困难
- 含有Indel的reads与参考基因组比对时，需要在Indel位置引入"空位"（gap）
- 比对算法对空位有罚分，倾向于将Indel区域的碱基错配解释为SNP
- 例如，参考序列"ATCGATCG"，reads为"ATC---GATCG"（3bp缺失），比对软件可能错误地比对为"ATCGATCG"（多个错配）
（2）Indel周围的碱基质量下降
- Illumina测序在Indel附近容易产生碱基识别错误
- 需要专门的Indel重比对步骤来纠正
（3）Indel大小的不确定性
- 对于较大的Indel（>10 bp），单个read可能无法完全跨越，导致Indel边界不确定
- 需要多个reads的信息来确定Indel的精确边界

（1）基于局部重比对的Indel检测（GATK HaplotypeCaller）
- HaplotypeCaller在局部窗口内生成候选单倍型（包含可能的SNP和Indel组合）
- 将reads与所有候选单倍型比对，选择最优解释
- 优点：同时考虑SNP和Indel，避免"SNP vs Indel"的比对歧义
（2）基于split-read的Indel检测（Pindel）
- 识别"split reads"——一条read的两部分分别比对到Indel两侧
- 通过分析未比对上的read片段来推断Indel序列
- 优点：对较大Indel（>20 bp）的边界定位准确
- 缺点：只能检测被reads跨越的Indel
（3）基于de novo局部组装的Indel检测（Platypus, Dindel）
- 在候选Indel区域提取相关reads
- 进行局部de novo组装
- 将组装结果与参考基因组比对，识别Indel

Indel的功能影响通常比SNP更严重：
| Indel类型 | 功能影响 | 示例 |
|----------|---------|------|
| 编码区非3倍数Indel | 移码突变，通常导致功能丧失 | BRCA1移码突变导致乳腺癌风险 |
| 编码区3倍数Indel | 氨基酸插入/缺失，可能影响功能 | CFTR ΔF508导致囊性纤维化 |
| 剪接位点Indel | 影响pre-mRNA剪接 | 深度内含子Indel激活隐蔽剪接位点 |
| 调控区Indel | 影响转录因子结合 | 胰岛素基因启动子Indel影响表达 |
| UTR区Indel | 可能影响mRNA稳定性或翻译 | 5'UTR的uORF相关Indel |
经典病例：
- 囊性纤维化：CFTR基因最常见的突变是ΔF508（一个3bp缺失），导致苯丙氨酸缺失，约占所有囊性纤维化病例的70%
- 亨廷顿舞蹈症：HTT基因中CAG重复序列的扩增（可视为一种特殊Indel），当CAG重复数>36时导致疾病

（1）基于Read深度（Read Depth, RD）的方法
原理：测序深度与基因组拷贝数成正比。
- 缺失区域：reads深度显著低于平均值（约50%对于杂合缺失）
- 重复区域：reads深度显著高于平均值
代表软件：CNVnator、BIC-seq、CNV-seq、DELLY（结合RD和PE）
优点：
- 可以检测大尺度的CNV
- 计算相对简单
局限：
- 无法精确定位断点（只能确定区域）
- 无法检测平衡性SV（如倒位、易位不改变拷贝数）
- 受GC含量、比对偏好等因素影响
（2）基于Split Read（SR）的方法
原理：识别一条read的两部分分别比对到不同位置的"split reads"，从而精确定位断点。
代表软件：Pindel、DELLY、LUMPY、Manta
工作流程：
1. 提取未完全比对或部分比对的reads（soft-clipped reads）
2. 将未比对部分重新比对到基因组其他位置
3. 如果找到匹配，确定断点坐标
优点：
- 断点定位精确（可达bp级别）
- 可以检测小至50 bp的SV
局限：
- 只能检测被reads跨越的SV
- 需要足够长的reads（对Illumina短读长有限制）
（3）基于Paired-End（PE）的方法
原理：利用paired-end测序中insert size和方向的异常来检测SV。
正常配对特征：
- 两个read在同一染色体上
- 方向相对（read1正向，read2反向）
- 距离在预期insert size范围内
异常配对类型及对应的SV：
| 异常类型 | 提示的SV |
|---------|---------|
| 距离远大于预期 | 缺失（reads跨越了被删除的区域） |
| 距离远小于预期 | 插入（reads之间的区域被插入序列替代） |
| 方向相同 | 倒位 |
| 位于不同染色体 | 易位 |
代表软件：DELLY、LUMPY、BreakDancer、Manta
优点：
- 可以检测大尺度SV
- 对read长度要求不高
局限：
- 断点定位不精确（通常在数百bp范围内）
- 受insert size分布的影响
（4）基于de novo组装的方法
原理：对SV区域进行局部de novo组装，将组装结果与参考基因组比对来识别SV。
代表软件：TIGRA、MindTheGap、SMRT-SV（PacBio专用）、Sniffles（Nanopore专用）
优点：
- 可以检测最复杂的SV，包括串联重复扩增、倒位等
- 可以获得SV的完整序列
局限：
- 计算资源需求大
- 组装质量依赖于数据质量和覆盖度

现代SV检测工具通常整合多种信号：
DELLY：整合PE、SR和RD信号
LUMPY：整合PE和SR信号，使用概率框架
Manta：结合PE、SR和局部组装，速度快，准确性高
GRIDSS：基于基因组图（assembly graph）的SV检测

长读长测序（PacBio/Nanopore）对于SV检测具有革命性优势：
1. 跨越复杂SV：长读长可以直接跨越整个SV区域，获得完整的SV序列
2. 精确定位断点：无需依赖PE或SR推断，直接比对即可定位
3. 检测重复区域SV：短读长难以比对的重复区域，长读长可以跨越
专用工具：
- Sniffles：Nanopore SV检测
- PBSV：PacBio SV检测
- SVIM：同时支持PacBio和Nanopore

SV检测的假阳性率通常高于SNP/Indel，验证至关重要：
实验验证方法：
- PCR+Sanger测序：对断点区域进行PCR扩增和测序，是SV验证的金标准
- qPCR：验证拷贝数变化
- Southern blot：验证大尺度重复/缺失
- 光学图谱（Bionano）：验证大尺度SV
计算验证方法：
- 使用不同检测工具的结果取交集
- 检查SV区域的reads比对情况（IGV可视化）
- 与已知SV数据库（如gnomAD-SV、DGV）比对

6.4.1 微生物组概述

人体微生物组的规模：
- 总重量：约1-2 kg（主要在肠道）
- 细胞数量：约3×10¹³个（与人体自身细胞数量相当或略多）
- 基因数量：约300-500万个基因（是人类基因组基因数的150-250倍）
人体各部位的微生物组特征：
| 部位 | 主要门类 | 细菌密度 | 主要功能 |
|------|---------|---------|---------|
| 结肠 | Firmicutes, Bacteroidetes | 10¹¹-10¹²/g | 发酵膳食纤维、合成维生素 |
| 口腔 | Firmicutes, Proteobacteria | 10⁹/mL唾液 | 参与氮循环、免疫调节 |
| 皮肤 | Actinobacteria, Firmicutes | 10⁶-10⁷/cm² | 屏障保护、免疫训练 |
| 阴道 | Lactobacillus | 10⁸-10⁹/mL | 维持酸性环境、防御病原 |
| 鼻腔 | Firmicutes, Actinobacteria | 10⁶-10⁷/cm² | 免疫防御 |

传统微生物学估计，环境中可培养的微生物仅占1%左右（"Great Plate Count Anomaly"）。主要原因包括：
1. 营养需求未知：很多微生物需要特殊的生长因子或复杂的共生关系
2. 环境条件难以模拟：如极端pH、温度、压力、厌氧条件等
3. 生长缓慢：某些微生物代时可达数周甚至数月
4. 群体感应依赖：需要达到一定密度才能生长
5. 严格的共生关系：依赖其他微生物产生的代谢产物
宏基因组学通过直接提取环境DNA进行测序，彻底 bypass 了培养瓶颈。

16S rRNA（原核生物）和18S rRNA/ITS（真核微生物）是微生物分类鉴定的"条形码"。
16S rRNA的优势：
1. 普遍存在：所有细菌和古菌都含有16S rRNA基因
2. 功能保守：核糖体功能对生命至关重要，序列演化相对缓慢
3. 大小适中：约1,500 bp，既包含保守区（设计通用引物），又包含可变区（区分物种）
16S rRNA的可变区：
- V1-V9共9个可变区
- 常用测序区域：V3-V4（约465 bp，Illumina MiSeq 2×300可覆盖）
- 不同可变区的分辨能力不同：V1-V2适合属水平，V3-V4可区分到种/株水平
分类阈值：
- 16S rRNA序列相似性 >97%：通常认为是同一个种
- 16S rRNA序列相似性 >95%：通常认为是同一个属

原始测序数据 (FASTQ)
    ↓
质量控制 (FastQC, cutadapt)
    ↓
去接头、去低质量、去嵌合体 (DADA2/Deblur/QIIME2)
    ↓
生成ASV/OTU表
    ↓
物种注释 (SILVA/Greengenes/RDP数据库)
    ↓
多样性分析 (α, β)
    ↓
统计分析 (LEfSe, ANCOM, DESeq2)
    ↓
可视化 (phyloseq, ggplot2)

（1）DADA2去噪算法
DADA2（Divisive Amplicon Denoising Algorithm 2）通过以下步骤生成ASV：
1. 质控过滤：去除低质量和含N的reads
2. 去重复：将相同序列合并，记录丰度
3. 错误学习：从数据中估计测序错误模型
4. 去噪：将含错误的序列纠正回其真实序列
5. 去嵌合体：识别并去除PCR嵌合体
6. 合并双端reads
（2）OTU聚类 vs ASV
- OTU聚类以97%相似性为阈值，将序列分组。缺点：阈值是人为设定的；近缘种可能因高于阈值而无法区分；信息在聚类过程中丢失
- ASV保留了每个独特的真实序列，分辨率更高，且可重复性更好

原始测序数据 (FASTQ)
    ↓
质量控制 (FastQC, Trimmomatic)
    ↓
去宿主DNA (Bowtie2比对宿主参考基因组)
    ↓
宏基因组组装 (MEGAHIT, metaSPAdes)
    ↓
基因预测 (Prodigal, MetaGeneMark)
    ↓
基因 catalog构建 (CD-HIT聚类)
    ↓
功能注释 (eggNOG-mapper, KEGG)
    ↓
分箱 (MetaBAT2, MaxBin2, CONCOCT)
    ↓
MAG质量评估 (CheckM)
    ↓
物种分类和系统发育分析

分箱的目标是将来自同一基因组的contig聚类到一起。
（1）基于组成的分箱（Composition-based Binning）
- 原理：不同物种的基因组具有不同的碱基组成特征（如GC含量、四核苷酸频率）
- 方法：将contig表示为高维特征向量（如136维四核苷酸频率向量），使用聚类算法（如K-means、高斯混合模型）进行聚类
- 代表工具：MetaBAT2、MaxBin2
（2）基于覆盖度的分箱（Coverage-based Binning）
- 原理：来自同一物种的contig在不同样本中的覆盖度变化模式一致（因为它们来自同一个基因组）
- 方法：在多个样本中计算每个contig的覆盖度，构建覆盖度矩阵，进行聚类
- 代表工具：CONCOCT、MaxBin2
（3）混合分箱
- 同时使用组成和覆盖度信息进行分箱
- 代表工具：MetaBAT2（使用TNF和覆盖度）、VAMB（使用变分自编码器）
MAG质量评估：
- CheckM：基于保守单拷贝基因（SCG）评估MAG的完整性和污染度
- 高质量MAG：完整性>90%，污染度<5%
- 中等质量MAG：完整性≥50%，污染度<10%

α多样性指数：
| 指数 | 衡量 | 公式 | 特点 |
|------|------|------|------|
| 物种数（Richness） | 丰富度 | S = 物种数量 | 不考虑丰度 |
| Shannon | 丰富度+均匀度 | H' = -Σ(pi × ln pi) | 对丰富度敏感 |
| Simpson | 丰富度+均匀度 | D = 1 - Σ(pi²) | 对优势种敏感 |
| Chao1 | 估计总丰富度 | Chao1 = S_obs + n₁²/(2n₂) | 考虑稀有物种 |
| Pielou's J | 均匀度 | J' = H'/ln(S) | 纯度的衡量 |
β多样性距离度量：
| 距离 | 公式 | 特点 |
|------|------|------|
| Bray-Curtis | 1 - 2Σmin(xi, yi) / Σ(xi + yi) | 考虑丰度，最常用 |
| Jaccard | 1 - |A∩B| / |A∪B| | 仅考虑有无 |
| Unweighted UniFrac | 考虑系统发育距离 | 考虑进化关系 |
| Weighted UniFrac | 考虑丰度和系统发育 | 综合考量 |
β多样性可视化：
- PCoA（主坐标分析）
- NMDS（非度量多维尺度分析）
- UPGMA聚类树

识别在不同组间丰度显著差异的微生物或功能基因：
- LEfSe（Linear Discriminant Analysis Effect Size）：结合统计检验和LDA评分，筛选组间差异标志物
- ANCOM（Analysis of Composition of Microbiomes）：考虑微生物组数据的组成性特点
- DESeq2/edgeR：基于负二项分布的差异分析（适用于OTU/基因计数数据）
- ALDEx2：考虑组成性和稀疏性的差异分析

微生物组数据具有以下"大数据"特征：
1. 高维度：数千个OTU/ASV × 数百个样本
2. 稀疏性：很多物种只在少数样本中出现（零膨胀分布）
3. 组成性：相对丰度总和为100%，变量间存在约束
4. 高变异性：个体差异大，受饮食、环境等多种因素影响
5. 多层次：物种、功能、代谢物等多个层面的数据

（1）松散的整合（Loose Integration）
- 分别分析各组学数据
- 比较不同组学层面的结果，寻找一致性
- 例如：比较差异物种和差异代谢物，寻找相关性
（2）统计学整合（Statistical Integration）
- 使用多变量统计方法整合多个数据矩阵
- CCA（Canonical Correlation Analysis）：寻找两个数据集之间的最大相关
- sPLS（sparse Partial Least Squares）：寻找预测关系的同时进行变量选择
- MCIA（Multiple Co-Inertia Analysis）：整合多个组学数据集
- MOFA（Multi-Omics Factor Analysis）：分解共享和特异的变异源
（3）网络整合（Network Integration）
- 构建跨组学的关联网络
- 例如：物种-代谢物关联网络、基因-代谢物关联网络
- 识别网络中的关键节点（hub）
（4）机器学习整合
- 将多组学特征作为输入，训练预测模型
- 使用特征选择识别最重要的预测因子

（1）分类任务：疾病预测
- 输入：微生物组组成特征（物种丰度、功能通路丰度）
- 输出：疾病状态（如IBD vs 健康）
- 常用算法：
- 随机森林（Random Forest）：对高维稀疏数据表现好，可输出特征重要性
- 支持向量机（SVM）：适合小样本高维数据
- 深度学习（神经网络）：可学习复杂的非线性关系
案例：使用随机森林预测CRC（结直肠癌）
- 训练集：数百个粪便样本的16S/宏基因组数据
- 特征：物种丰度（或MetaCyc通路丰度）
- 模型：随机森林分类器
- 结果：AUC可达0.8-0.9，某些研究报道准确率>90%
- 关键标志物：Fusobacterium nucleatum、Peptostreptococcus anaerobius等
（2）回归任务：表型预测
- 输入：微生物组特征 + 宿主因素
- 输出：连续变量（如BMI、药物响应、血糖水平）
- 常用算法：LASSO回归、弹性网络、随机森林回归
（3）聚类任务：亚型鉴定
- 输入：微生物组组成
- 输出：患者亚型（enterotypes，肠型）
- 肠型概念：人类肠道微生物组可分为几种主要的组成模式
- 肠型1：以Bacteroides为主
- 肠型2：以Prevotella为主
- 肠型3：以Ruminococcus为主

微生物组研究面临的一个核心挑战是因果关系的确定：
- 相关性不等于因果性：微生物组变化可能是疾病的原因，也可能是疾病的结果
- 混杂因素：饮食、药物、生活方式等同时影响微生物组和疾病状态
- 策略：
- 纵向研究：追踪疾病发生前后的微生物组变化
- 动物模型：无菌小鼠或抗生素处理后的定植实验
- 孟德尔随机化：利用宿主基因型作为工具变量
- 干预研究：益生菌、FMT等干预后的效果评估

（1）疾病诊断
- 结直肠癌（CRC）筛查：F. nucleatum等细菌标志物+粪便免疫化学检测（FIT）联合筛查
- 艰难梭菌感染（CDI）诊断：结合临床表现和微生物组特征
- 非酒精性脂肪性肝病（NAFLD）评估：肠道菌群组成与肝纤维化程度相关
（2）疾病治疗
- 复发性CDI：FMT治愈率>90%，已被FDA批准
- 溃疡性结肠炎（UC）：多项RCT显示FMT对UC有一定疗效
- 自闭症谱系障碍（ASD）：小规模临床试验显示FMT可改善胃肠道症状和行为
- 黑色素瘤免疫治疗：肠道菌群组成影响PD-1抑制剂的疗效
（3）药物代谢
- 肠道微生物组参与多种药物的代谢，影响药效和毒性
- 例如：伊立替康（irinotecan）的毒性由肠道细菌的β-葡萄糖醛酸酶介导
- 通过抑制特定细菌酶，可以降低药物毒性

（1）植物微生物组与农业
- 根际微生物组：影响植物养分吸收、抗病性和抗逆性
- 生物肥料：固氮菌、解磷菌等促进植物生长
- 生物防治：利用拮抗微生物防治植物病原菌
- 案例：印度矮生小麦与根际促生菌（PGPR）联合使用，可减少化肥用量30%
（2）动物微生物组与畜牧业
- 反刍动物瘤胃微生物组：决定饲料转化效率
- 益生菌：改善家禽和家畜的生长性能和免疫力
- 替代抗生素：欧盟禁用抗生素促生长剂后，微生物组产品成为替代方案

（1）生物修复
- 利用环境微生物组降解污染物
- 石油污染修复：海洋微生物组中的烃降解菌可降解石油烃
- 重金属污染修复：某些微生物可将有毒重金属转化为低毒形态
（2）功能基因挖掘
- 从环境宏基因组中发现新的酶和生物活性物质
- 纤维素酶：从瘤胃、白蚁肠道等环境中发现高效纤维素降解酶
- 新型抗生素：从土壤宏基因组中发现新结构抗生素（如Teixobactin）
（3）碳循环和气候变化
- 土壤微生物组在全球碳循环中发挥关键作用
- 冻土融化释放古老微生物，可能加速温室气体排放
- 海洋微生物组（特别是浮游植物共生菌）影响全球碳固定

（1）益生菌设计
- 传统益生菌：Lactobacillus, Bifidobacterium等
- 下一代益生菌（Next-Generation Probiotics）：Akkermansia muciniphila, Faecalibacterium prausnitzii等
- 工程化益生菌：通过基因工程赋予益生菌特定功能（如产丁酸盐、分泌治疗性蛋白）
（2）合成微生物组
- 目标：构建具有稳定功能的最小微生物群落
- 方法：
- 自上而下：从复杂群落中简化出核心功能群
- 自下而上：从单一菌株开始，逐步构建稳定群落
- 挑战：微生物间互作复杂，群落稳定性难以预测

（1）样本采集和保存
- 微生物组在采样后迅速变化（如氧气暴露导致厌氧菌死亡）
- 保存方法影响结果：冷冻、RNAlater、FTA卡等各有优缺点
- 建议：采样后尽快冷冻（-80°C），或使用稳定化试剂
（2）DNA提取偏差
- 不同试剂盒对不同类型细菌的提取效率差异很大
- 革兰氏阳性菌（厚壁菌门）细胞壁含大量肽聚糖，比革兰氏阴性菌更难裂解
- 策略：在同一研究中统一使用同一种提取方法；使用裂解对照（spike-in controls）
（3）PCR扩增偏差
- 通用引物对不同物种的扩增效率不同
- 嵌合体（chimera）和异源双链（heteroduplex）的形成
- 策略：使用高保真酶、优化PCR循环数、使用去嵌合体算法
（4）测序深度不足
- 复杂样本（如土壤）可能需要>10 Gb数据才能充分覆盖稀有物种
- 深度不足会导致稀有物种被遗漏，影响多样性估计
（5）生物信息学分析的异质性
- 不同分析流程（QIIME vs mothur）、不同数据库（SILVA vs Greengenes）会产生不同结果
- 缺乏统一的金标准

（1）样本量计算
- 微生物组研究通常需要比传统临床研究更大的样本量
- 因为个体差异大，效应量通常较小
- 建议：初步探索性研究每组至少20-30个样本；验证性研究每组50-100个
（2）对照设置
- 阴性对照（空白对照）：检测试剂污染
- 阳性对照（mock community）：评估技术流程的准确性
- 裂解对照（spike-in）：评估DNA提取效率
（3）批次控制
- 随机化：将不同组的样本随机分配到不同批次
- 平衡设计：确保每批次包含所有组的样本
- 记录所有批次信息，以便后续统计校正
- 统计校正：使用ComBat等方法校正批次效应
（4）纵向设计
- 对于疾病研究，纵向采样（疾病前-中-后）比横断面设计更能揭示因果关系
- 建议：基线（T0）、急性期（T1）、恢复期（T2）

标准化倡议：
- MIxS标准（Minimum Information about any (x) Sequence）：包括MIMARKS、MIMS等，规定了宏基因组研究需要报告的元数据
- Earth Microbiome Project（EMP）：制定了标准的采样和测序方案
数据共享平台：
- NCBI SRA/EBI ENA：原始测序数据
- Qiita/MG-RAST：处理后的微生物组数据
- 数据共享是提高研究可重复性和促进元分析的关键

（1）方法学进展
- 长读长宏基因组学：Nanopore/PacBio长读长有助于解决重复序列和分箱问题
- 单细胞宏基因组学：在单细胞分辨率上解析微生物群落
- 空间宏基因组学：保留微生物的空间分布信息
（2）人工智能与微生物组
- 深度学习用于物种分类和功能预测
- 生成式AI用于设计合成微生物组
- 大语言模型整合文献知识和微生物组数据
（3）精准微生物组医学
- 基于个体微生物组特征的个性化益生菌
- 微生物组作为疾病早期诊断的生物标志物
- 微生物组靶向治疗（如噬菌体疗法、CRISPR编辑肠道菌）
（4）全球微生物组计划
- 类似于人类基因组计划，全球微生物组计划旨在系统性地绘制地球微生物多样性图谱
- 应用：生物勘探、环境监测、气候变化预测

基因组学作为生物信息学的核心领域，在过去二十年中经历了前所未有的技术革命和方法学进步。从人类基因组计划的完成到个人基因组测序进入千元时代，从短读长测序主导到长读长技术崛起，基因组学的发展深刻改变了生命科学研究的面貌。

本章核心要点回顾

测序技术层面：Illumina短读长测序以其高准确性和高通量继续占据主流地位，而PacBio HiFi和Nanopore长读长测序则在复杂基因组组装和结构变异检测中展现出独特优势。混合测序策略正在成为高质量基因组项目的标准配置。
组装与注释层面：de Bruijn图算法是二代测序组装的核心方法，而OLC算法在长读长时代重新焕发活力。Hi-C技术使染色体级别组装成为可能。基因预测已从单纯的从头预测发展为整合多种证据的共识预测策略。功能注释依赖于同源性比对、结构域识别和通路映射等多种方法。
变异检测层面：SNP检测已达到极高的准确性，Indel检测需要特殊的比对处理，而SV检测仍是挑战性最大的领域。长读长测序为SV检测带来了突破性进展。
宏基因组学层面：宏基因组学通过bypass培养瓶颈，揭示了微生物世界的巨大多样性。从16S扩增子测序到宏基因组shotgun测序，从多样性分析到功能挖掘，宏基因组学正在从描述性研究向机制性研究和应用转化迈进。

未来展望

完整基因组时代：端粒到端粒（T2T）组装技术的成熟将使每个物种都能获得真正完整的基因组参考序列，包括着丝粒、端粒和核糖体DNA等 previously intractable 区域。
泛基因组学：单一参考基因组已不足以代表物种内的遗传多样性。泛基因组（pan-genome）研究将揭示核心基因组和可变基因组的动态演化，为作物改良和疾病研究提供新的视角。
单细胞基因组学：单细胞测序技术的进步将使研究者能够在单细胞分辨率上解析基因组变异，揭示肿瘤异质性、发育谱系和微生物群落结构。
人工智能与基因组学：深度学习已经在变异检测（DeepVariant）、蛋白质结构预测（AlphaFold）等领域取得突破。未来，AI将进一步整合基因组、转录组、蛋白质组和表观遗传组数据，实现从序列到功能的端到端预测。
精准基因组医学：随着测序成本的进一步下降和数据分析方法的成熟，全基因组测序有望成为临床常规检测手段，为罕见病诊断、肿瘤精准治疗和药物基因组学提供个体化的遗传信息。
基因组学的发展远未止步。正如人类基因组计划的领导者之一Francis Collins所言："基因组测序只是开始，真正的挑战和机遇在于解读这本生命之书。"作为生物信息学工作者，我们正处于这一伟大事业的核心，肩负着从海量基因组数据中提取生物学洞见、推动精准医学发展的历史使命。

本章思考题

1. 回顾基因组学发展的历史，你认为哪三项技术突破对领域的影响最大？为什么？
2. 如果你要启动一个全新物种的基因组项目，请设计一个完整的技术路线，包括测序策略、组装方案和注释流程。
3. 比较SNP、Indel和SV在检测难度、功能影响和进化意义方面的异同。
4. 宏基因组学如何改变了我们对微生物世界和宿主-微生物互作的理解？
5. 展望基因组学的未来，你认为下一个重大突破可能出现在哪个方向？

推荐阅读
1. Goodwin S, McPherson J D, McCombie W R. Coming of age: ten years of next-generation sequencing technologies [J]. Nature reviews genetics, 2016, 17: 333-351.
2. Nurk S, Koren S, Rhie A, et al. The complete sequence of a human genome [J]. Science, 2022, 376: 44-53. （T2T-CHM13）
3. Quince C, Walker A W, Simpson J T, et al. Shotgun metagenomics, from sampling to analysis [J]. Nature biotechnology, 2017, 35: 833-844.
4. Turner T N, Hormozdiari F, Duyzend M H, et al. Genome sequencing of autism-affected families reveals disruption of putative noncoding regulatory DNA [J]. American journal of human genetics, 2016, 98: 58-74.

（1）测序技术的三次革命
第一代测序（Sanger, 1977）：
- 原理：双脱氧链终止法
- 特点：准确性极高（>99.99%），读长可达1000 bp，但通量极低
- 贡献：完成首批模式生物基因组（噬菌体φX174、流感嗜血杆菌、酿酒酵母）
- 局限：人类基因组计划耗时13年、耗资27亿美元
第二代测序（Illumina, 2005）：
- 原理：桥式PCR + 可逆终止子边合成边测序
- 特点：通量提升百万倍，读长100-300 bp，准确性>99.9%，成本降至$1000/人基因组
- 贡献：千人基因组计划、肿瘤基因组图谱（TCGA）、宏基因组学兴起
- 局限：短读长无法跨越重复序列和复杂区域
第三代测序（PacBio/Nanopore, 2011-2015）：
- 原理：单分子实时测序（ZMW荧光检测 / 纳米孔电流检测）
- 特点：读长可达数Mb，无需PCR扩增，可检测碱基修饰
- 贡献：T2T基因组、复杂结构变异解析、直接RNA测序
- 局限：原始错误率较高（PacBio HiFi已解决），成本仍高于二代
（2）组装算法的演进
- 2001-2008：基于OLC的Sanger组装（Celera Assembler, Arachne）
- 2008-2015：基于DBG的二代短读长组装（Velvet, SOAPdenovo, SPAdes）
- 2015-2020：长读长OLC组装（Canu, FALCON）
- 2020至今：HiFi精准长读长组装（Hifiasm, Verkko）+ T2T整合策略
（3）变异检测的精度跃升
- 早期：基于pileup的统计方法（Samtools mpileup）
- 中期：单倍型模型（GATK HaplotypeCaller）
- 近期：深度学习（DeepVariant, Google Brain）——在多个基准测试中达到或超越人类专家水平

（1）数据层面的挑战
- 数据量爆炸：全球测序数据已超100 EB（Exabyte），存储和传输成本急剧上升
- 数据异质性：不同平台（Illumina/PacBio/Nanopore）、不同版本（hg19/hg38/T2T-CHM13）的数据整合困难
- 数据质量不均：公共数据库中大量低质量组装和错误注释存在
（2）分析层面的挑战
- 计算资源需求：T2T组装、泛基因组构建、单细胞分析需要PB级内存和万核级计算
- 算法瓶颈：高度重复序列、复杂结构变异、多倍体基因组仍是组装的难点
- 标准化缺失：不同分析流程的结果可比性差，缺乏统一的金标准
（3）解读层面的挑战
- 功能未知基因：即使在人类基因组中，仍有数千个基因功能不明
- 非编码区解读：增强子、沉默子、绝缘子等调控元件的精确预测和功能验证困难
- 变异临床意义：大量罕见变异的致病性难以判定（ACMG指南仍依赖专家判断）
- 多基因复杂疾病：大多数常见疾病（糖尿病、高血压、精神疾病）受数千个微效变异影响，单一基因解释力有限

（1）完整基因组时代（T2T for All）
- 目标：为地球上所有物种（>1000万种真核生物）提供T2T参考基因组
- 驱动力：地球生物基因组计划（EBP）、Zoonomia项目（哺乳类）、达尔文树计划（中国）
- 技术：HiFi + ONT + Hi-C的标准化流程，自动化组装和审校
- 意义：理解生物多样性、保护濒危物种、挖掘新资源
（2）泛基因组与图基因组
- 单一参考基因组已不能满足研究和临床应用的需求
- 图基因组将线性参考升级为包含群体变异的图结构，减少比对偏差
- 人类泛基因组参考联盟（HPRC）正在构建包含>40个多样化单倍型的泛基因组
- 应用：改进变异检测、发现参考基因组中缺失的序列、群体遗传学研究
（3）单细胞基因组学
- 从"bulk平均信号"到"单细胞分辨率"
- 单细胞基因组测序揭示肿瘤异质性、发育谱系、微生物群落结构
- 技术挑战：全基因组扩增的偏倚、低覆盖度下的变异检测、海量数据的分析
- 前沿方向：单细胞多组学（同时测基因组+转录组+表观组+蛋白质组）
（4）人工智能与基因组学
- 变异检测：DeepVariant、Clair3等深度学习工具已达到金标准水平
- 蛋白质结构预测：AlphaFold2革命性地解决了蛋白质折叠问题
- 调控元件预测：Enformer、Basenji等模型从序列预测基因表达和调控
- 端到端预测：从DNA序列直接预测功能（如SpliceAI预测剪接位点）
- 大语言模型：基于Transformer的基因组语言模型（如DNABERT、Nucleotide Transformer）正在学习DNA序列的"语义"
（5）精准基因组医学
- 罕见病诊断：全基因组测序（WGS）已成为未确诊遗传病的首选诊断工具，诊断率约25-40%
- 肿瘤精准治疗：基于ctDNA的液体活检、MRD（微小残留病灶）监测、个性化疫苗
- 药物基因组学：基于CYP2D6、HLA-B57:01等基因型指导用药，避免严重不良反应
- 新生儿筛查：扩展性基因组筛查（NICUSeq）在重症新生儿中快速诊断遗传病
- 预防医学：基于多基因风险评分（PRS）的疾病风险预测和早期干预
（6）合成基因组学*
- 从"读取"到"编写"：基因组合成成本从2003年的$4/bp降至2020年的$0.001/bp以下
- 最小基因组：Mycoplasma laboratorium（2010年，首个合成细胞）
- 酵母基因组合成计划（Sc2.0）：重新设计并合成16条酵母染色体
- 人类基因组的合成规划（GP-write）正在讨论中
- 应用：设计新型生物制造底盘、创造抗病毒/抗逆作物、开发活体疗法

中国在基因组学领域已从"跟跑"进入"并跑"甚至部分"领跑"阶段：
测序技术自主化：
- 华大智造（MGI）的DNBSEQ技术成为全球第二大测序平台
- 真迈生物、齐碳科技等在纳米孔测序领域取得突破
大科学计划：
- 中国十万人基因组计划
- 百万微生态计划
- 万种鸟类基因组计划（B10K中国部分）
- 地球生物基因组计划中国节点
数据基础设施建设：
- 国家基因组科学数据中心（NGDC）
- 中国国家生物信息中心（CNCB）
- 基因组数据安全存储和共享的国家标准

6.1.1 基因组学的概念

基因组学研究的范畴可以从三个层次理解：
（1）结构基因组学（Structural Genomics）
结构基因组学关注基因组的物理结构，包括DNA序列的测定、基因组图谱的构建、基因的定位和结构分析等。其核心目标是获得基因组的完整序列，并理解基因在基因组中的排列方式。
（2）功能基因组学（Functional Genomics）
功能基因组学研究基因组中基因的功能、基因之间的相互作用、基因表达调控等。与结构基因组学回答"是什么"不同，功能基因组学回答"做什么"和"怎么做"。主要技术包括RNA-Seq、ChIP-Seq、ATAC-Seq等。
（3）比较基因组学（Comparative Genomics）
比较基因组学通过比较不同物种或不同个体基因组的异同，揭示基因组的进化规律、基因的起源和功能分化、物种之间的亲缘关系等。

基因组学中所研究的问题主要分为三类：
第一类：如何获得基因组序列
这是基因组学最基础的问题。随着DNA测序技术的发展，从一代Sanger测序到二代高通量测序（NGS），再到三代单分子长读长测序，获得基因组序列的成本急剧下降、速度大幅提升。一个物种的基因组项目通常包括：基因组大小估计、测序策略设计、序列组装、质量评估等步骤。
第二类：如何解读/解码基因组
获得序列后，需要解读这本"生命之书"。这包括：基因预测与结构注释（识别编码区、内含子、外显子、调控元件等）、功能注释（确定基因的功能、参与的通路）、比较分析（与近缘物种比较，理解进化关系）、变异分析（识别SNP、Indel、SV等）。
第三类：如何重写/编写新的基因组
这是合成生物学的前沿方向。通过基因编辑技术（如CRISPR-Cas9）对现有基因组进行精确修饰，甚至从头设计和合成全新的基因组。2010年，Craig Venter团队首次合成了支原体的人工基因组，标志着人类开始具备"编写生命"的能力。

人类基因组计划采用了层次测序策略：通过构建DNA大片段克隆文库、重叠克隆群和物理图谱，再进行逐个克隆鸟枪法测序。中国科学家承担了1%的测序任务（3号染色体端部约3000万个碱基对），成功跟进世界先进测序技术，推动了国内测序仪研发和生物信息学软件开发。
HGP的完成不仅提供了人类基因组的参考序列，更重要的是：
1. 技术推动：测序技术从Sanger法发展到高通量测序
2. 数据共享：建立了基因组数据开放共享的文化
3. 疾病研究：为疾病基因定位、药物基因组学奠定基础
4. 进化研究：为理解人类起源和进化提供基础数据

Sanger测序的核心是利用ddNTP终止DNA链延伸：
1. 反应体系：模板DNA、引物、DNA聚合酶、四种dNTP（其中一种带有放射性或荧光标记）和少量ddNTP
2. 链延伸与终止：DNA聚合酶沿模板合成新链，遇到dNTP时正常延伸，遇到ddNTP时因缺少3'-OH而终止
3. 四个独立反应：分别加入ddATP、ddCTP、ddGTP、ddTTP，产生四组不同长度的片段
4. 电泳分离：通过聚丙烯酰胺凝胶电泳（PAGE）或毛细管电泳按大小分离片段
5. 读取序列：从短到长读取末端碱基，获得互补链序列
自动化Sanger测序使用四种不同荧光标记的ddNTP，四个反应在同一管中进行，通过激光检测荧光颜色确定碱基类型。
局限性：通量低（每次反应只能读取一条序列），成本高，不适合大规模基因组项目。

Illumina测序流程包括四个关键步骤：
（1）文库制备
- 将基因组DNA随机打断成200-500 bp片段
- 末端修复、加A尾、连接测序接头（adapter）
- PCR扩增文库（可选）
（2）Flow Cell杂交与桥式PCR
- Flow Cell表面固定有两种与接头互补的引物
- 单链DNA文库与引物杂交
- 通过桥式PCR扩增：DNA片段弯曲成桥状，以另一端为模板进行扩增
- 经过多轮扩增，每个原始分子形成一个约1000-2000个克隆的"簇"（cluster）
- 变性后得到单链簇，作为测序模板
（3）边合成边测序（SBS）
- 加入DNA聚合酶和四种带有不同荧光标记的dNTP（每种dNTP的3'-OH被可切除的化学基团阻断）
- 每个循环：
a. 四种dNTP竞争掺入，每次只有一个碱基被加入
b. 激光扫描，根据荧光颜色确定掺入的碱基
c. 切除荧光基团和3'阻断基团，使链可以继续延伸
- 重复上述循环，逐碱基读取序列
（4）双端测序（Paired-End Sequencing）
- 从簇的两端分别测序，获得两段各约150 bp的序列
- 两段序列之间距离已知（约200-500 bp）
- 双端信息对于基因组组装和结构变异检测至关重要
数据格式与质量控制
测序数据以FASTQ格式存储，每条记录包含四行：序列标识符、碱基序列、"+"分隔符、碱基质量值（Phred分数）。Phred质量值Q与碱基错误率ε的关系为：Q = -10log₁₀(ε)。Q30表示碱基错误率约为0.1%，即准确性为99.9%。

（1）PacBio SMRT测序
PacBio测序的核心是零模波导孔（ZMW）技术：
- 在芯片上蚀刻数十万个直径约70 nm、深度约100 nm的纳米孔
- 每个孔底部固定一个DNA聚合酶分子
- 待测DNA与引物结合，被聚合酶捕获
- 加入四种带有不同荧光标记的dNTP（标记在磷酸基团上，而非碱基上）
- 当dNTP被掺入新链时，荧光标记被释放到溶液中，产生闪光信号
- 记录每个ZMW中的闪光序列和时间，即可获得DNA序列
PacBio Sequel II可产生约15-20 kb的长读长，最新的Revio系统进一步提高了通量。通过环化一致性测序（CCS，Circular Consensus Sequencing），同一个分子多次循环测序，可获得准确性>99.9%的HiFi reads。
（2）Nanopore测序
Nanopore测序的核心原理：
- 蛋白质纳米孔嵌入脂质双分子层膜中
- 在膜两侧施加电压，产生离子电流
- DNA/RNA分子在动力蛋白（helicase）驱动下穿过纳米孔
- 不同碱基（或k-mer）穿过纳米孔时引起特征性的电流变化
- 通过深度学习模型（如Bonito、Dorado）将电流信号解码为碱基序列
Nanopore的优势：
- 超长读长：目前最长记录超过4 Mb
- 便携性：MinION设备仅U盘大小，可在野外、太空等极端环境使用
- 直接测序RNA：无需逆转录，可检测RNA修饰
- 实时测序：数据边测边出，可即时分析
局限性：原始读取准确性较低（约92-97%），但R10.4.1 flow cell配合Dorado basecaller的准确性已大幅提升。

由于二代测序准确性高但读长短，三代测序读长长但原始准确性低，混合测序（Hybrid-Seq）成为当前基因组项目的标准配置：
- 用三代长读长进行基因组骨架组装
- 用二代短读长对组装结果进行纠错和填补缺口
- Hi-C数据辅助染色体级别组装

传统的基因组组装在高度重复序列区域（如着丝粒、端粒、rDNA阵列）存在大量缺口，主要原因是：
- 二代短读长无法跨越长重复单元
- 重复序列使组装算法难以确定正确的路径
T2T联盟实现人类完整基因组的关键技术组合：
（1）PacBio HiFi长读长
- 读长约15-20 kb，准确性>99.9%
- 足以跨越大多数重复单元（如着丝粒卫星重复单元约171 bp）
- 高准确性使得区分高度相似的重复序列成为可能
（2）Oxford Nanopore超长读长
- 读长可达100 kb至数Mb
- 直接跨越最大的重复阵列（如5S rDNA阵列可长达数Mb）
- 提供长距离的序列连续性信息
（3）Hi-C辅助scaffold连接
- 将HiFi组装产生的contig提升到染色体级别
- 确定contig的顺序、方向和间隙大小
（4）手动审校与验证
- 对复杂区域（如着丝粒）进行手动检查
- 使用荧光原位杂交（FISH）和光学图谱进行验证

Hi-C实验流程：
1. 甲醛交联：将空间上接近的DNA片段共价交联
2. 酶切：用限制性内切酶（如HindIII、DpnII）切割DNA
3. 末端修复和生物素标记：补平粘性末端并掺入生物素标记的dNTP
4. 连接：稀释后连接，使空间上接近的片段连接在一起
5. 去交联和测序：打断DNA，富集生物素标记的嵌合片段，进行双端测序
Hi-C数据分析用于基因组组装：
- 构建接触矩阵（contact matrix）：染色体各区域之间的相互作用频率
- 同一染色体上的区域相互作用频率高（呈现对角线富集）
- 相邻区域相互作用频率高于远端区域
- 使用算法（如3D-DNA、SALSA、ALLHiC）将contig聚类到染色体组、确定顺序和方向

泛基因组分为两部分：
- 核心基因组（Core Genome）：物种所有个体共有的基因，通常与基本生命功能相关
- 可变基因组/附属基因组（Accessory/Dispensable Genome）：仅在部分个体中存在的基因，通常与环境适应、致病性、代谢多样性等相关
泛基因组构建方法：
- 迭代式：将新个体的基因组与已有泛基因组比较，添加新序列
- 图基因组（Graph Genome）：将所有序列变异表示为图结构，而非线性序列
泛基因组的应用：
- 捕获参考基因组中缺失的遗传多样性
- 改进变异检测（减少比对偏差）
- 识别与重要性状相关的基因（如作物抗病基因）

一个典型的T2T基因组项目测序策略：
| 数据类型 | 测序深度 | 用途 |
|----------|---------|------|
| PacBio HiFi | 30-50× | 主要组装，产生高质量contig |
| Nanopore Ultra-long | 10-20× | 跨越超长重复区域，解决复杂区域 |
| Illumina PE150 | 50-100× | 纠错、质量验证 |
| Hi-C | 60-100× | 染色体级别scaffold |
| 光学图谱 | 5-10× | 验证组装结构 |

（1）线性序列展示
- 以染色体为坐标轴，从左到右展示DNA序列
- 支持不同尺度的缩放（从单碱基到全染色体）
- 显示正反链的序列和注释
（2）多轨道（Track）整合
- 每个track代表一类数据（基因注释、测序reads、变异、保守性等）
- 用户可自定义加载、排列和配置track
- 支持track的叠加和比较
（3）交互式查询
- 通过基因名、染色体坐标、rsID等定位基因组区域
- 支持序列搜索（如查找特定基序）
- 支持不同物种间的基因组比对（synteny view）
（4）数据导出
- 导出特定区域的序列（FASTA）
- 导出track数据（BED、GTF等）
- 生成高质量图片用于发表

（1）UCSC Genome Browser
- 最经典的基因组浏览器之一
- 整合了大量注释数据和工具
- 支持自定义track加载（通过URL或本地文件）
- 强大的Table Browser用于数据查询和下载
（2）Ensembl Genome Browser
- 欧洲生物信息研究所（EBI）维护
- 注释质量高，更新及时
- 强大的比较基因组学工具（如Ensembl Compara）
- 提供API和BioMart数据查询工具
（3）IGV（Integrative Genomics Viewer）
- 本地运行的桌面软件
- 支持大规模数据（如BAM、VCF、Hi-C矩阵）
- 特别适合查看测序reads比对情况和变异位点
- 支持基因组重排和结构变异的可视化
（4）JBrowse/GBrowse
- 可部署在本地服务器的网页浏览器
- 适合构建物种特异性的基因组数据库
- 支持插件扩展

BED格式
BED文件至少包含3列，最多12列：
- chrom：染色体名称
- chromStart：起始位置（0-based）
- chromEnd：终止位置（1-based，不包含）
- name：名称（可选）
- score：分数（0-1000，可选）
- strand：方向（+/-/.）
- thickStart：编码区起始（可选，用于显示）
- thickEnd：编码区终止（可选）
- itemRgb：颜色（可选）
- blockCount：外显子数量（BED12）
- blockSizes：外显子大小列表（逗号分隔）
- blockStarts：外显子相对起始位置列表
WIG格式
variableStep格式示例：

variableStep chrom=chr1 span=50
1001 12.5
1051 13.2
1101 11.8

表示chr1上1001-1050位置的值为12.5，1051-1100位置的值为13.2。
GTF格式
GTF是GFF2的变体，每行9列：
1. seqname：染色体或scaffold
2. source：注释来源
3. feature：特征类型（gene, transcript, exon, CDS等）
4. start：起始位置（1-based）
5. end：终止位置（1-based，包含）
6. score：分数（通常为".")
7. strand：方向（+/-/.)
8. frame：阅读框（0/1/2/.）
9. attributes：属性（gene_id, transcript_id等）

现代基因组研究往往涉及多组学数据的整合：
- 基因组层：基因结构、变异位点
- 转录组层：RNA-Seq信号、剪接异构体
- 表观组层：ChIP-Seq峰（组蛋白修饰、转录因子结合）、DNA甲基化、ATAC-Seq开放染色质
- 三维结构层：Hi-C相互作用矩阵、TAD边界
- 进化层：跨物种保守性分数（phyloP, phastCons）
整合可视化可以帮助发现跨层面的关联，例如：
- 某个增强子区域（H3K27ac峰+ATAC-Seq开放区）与远处基因的启动子有Hi-C相互作用
- 一个SNP位于保守非编码区，且与附近基因的表达量显著相关

6.2.1 基因组的组装

（1）重复序列
基因组中普遍存在重复序列，包括串联重复（如微卫星）和散在重复（如转座子、segmental duplication）。重复序列导致组装歧义——相同的reads可能来源于基因组的不同位置，组装算法无法确定其真实来源。
（2）测序错误
Illumina测序错误率约0.1-1%（Q30-Q20），Nanopore原始错误率可达5-15%。测序错误会干扰k-mer的正确连接，在DBG中产生"气泡"（bubble）和"死胡同"（tip）。
（3）覆盖度不均
基因组不同区域的测序深度差异很大。高GC区域、AT富集区域、重复序列区域的覆盖度往往偏低，导致这些区域组装困难。
（4）杂合度
二倍体生物的两条同源染色体存在差异。组装软件需要区分同源序列并可能分别组装（phased assembly），否则杂合位点会导致组装碎片化。

DBG算法是二代短读长组装的核心方法，其工作流程如下：
（1）k-mer分解
将所有reads拆分为长度为k的短片段（k-mer）。例如，序列"ATGCGGAT"在k=5时产生：ATGCG、TGCGG、GCGGA、CGGAT。通常k取奇数（如21、31、51、71），且不超过read长度的2/3。
（2）构建DBG
- 每个k-mer作为图中的一个顶点（node）
- 若两个k-mer之间有k-1个碱基重叠，则在它们之间建立一条有向边
- 例如，ATGCG和TGCGG共享TGCG（4个碱基），建立从ATGCG到TGCGG的边
- 将不同reads产生的相同k-mer合并为一个顶点，记录覆盖度
（3）图简化
测序错误会在DBG中产生两类错误结构：
- Tips（死胡同）：由于测序错误产生的低频k-mer连接到主路径后又很快终止。处理：去除长度和覆盖度低于阈值的tips。
- Bubbles（气泡）：由于杂合位点或测序错误，主路径在某处分叉为两条路径，在下游又汇合。处理：保留覆盖度较高的路径，或在杂合组装中保留两条路径。
（4）contig提取
在简化后的DBG中，一条连续的路径（从没有入度的顶点开始到没有出度的顶点结束）对应一个contig。
k-mer大小的选择：
- k值较小：对重复序列的分辨率低（不同重复单元可能共享k-mer），但覆盖度高（每个read产生更多k-mer）
- k值较大：对重复序列的分辨率高，但覆盖度低（要求更大的测序深度才能保证每个k-mer被覆盖）
- 多k-mer策略（如SPAdes）：从较小的k开始组装，逐步增大k值，利用小k的连续性和大k的准确性优势

OLC算法更适合长读长数据（PacBio/Nanopore），其工作流程：
（1）Overlap（重叠检测）
计算所有reads两两之间的重叠关系。对于N条reads，需要O(N²)次比对，计算量巨大。对于长读长数据，通常使用minimap2等快速比对工具进行all-vs-all比对。
（2）Layout（布局构建）
根据重叠关系构建布局图（string graph）：
- 每个read是图中的一个节点
- 如果两个reads有显著重叠，则在它们之间建立边
- 边的方向和长度由重叠的位置和长度确定
- 简化图结构（去除 transitive edges，压缩 unambiguous paths）
（3）Consensus（一致性序列生成）
对布局图中的每个unitig（unambiguous path），将对应的reads进行多序列比对，生成一致性序列作为contig。

利用二代和三代数据的互补性：
（1）二代纠错三代
- 用高准确性的二代reads对三代reads进行纠错
- 代表工具：PBcR、LoRDEC、Jabba
- 纠错后的三代reads准确性可提升至>99%
（2）三代搭建骨架，二代填充
- 用三代长读长构建初步的contig骨架
- 用二代reads对三代contig进行 polish（填补gap、纠正碱基错误）
- 代表工具：Pilon、Racon
（3）graph-based混合组装
- 将二代和三代数据同时输入组装算法
- 在组装图中同时利用短读长的准确性和长读长的连续性
- 代表工具：HybridSPAdes、MaSuRCA

连续性指标：
- N50/L50：衡量contig/scaffold长度的分布
- 最大contig长度：最长的连续序列
- gap数量：scaffold中N碱基的数量和总长度
完整性指标：
- BUSCO（Benchmarking Universal Single-Copy Orthologs）：检测基因组中保守单拷贝基因的存在情况。高BUSCO完整性（如>90%）表示基因组组装较完整。
- CEGMA（Core Eukaryotic Genes Mapping Approach）：类似BUSCO，基于核心真核基因集。
准确性指标：
- QUAST：计算N50、错误组装（misassembly）数量、基因组覆盖度等
- Pilon polish后的二代比对一致性：检测碱基错误
- k-mer spectrum分析：与原始测序数据的k-mer谱比较，评估完整性

原核生物基因特征：
- 基因密度高（大肠杆菌基因组约85%为编码区）
- 无内含子（少数古菌有内含子）
- 无5'端帽子和3'端polyA尾
n- 操纵子结构：多个基因共用一个启动子，转录为多顺反子mRNA
- 起始密码子通常为ATG（少数GTG、TTG）
- 启动子元件：-10区（Pribnow框：TATAAT）和-35区（TTGACA）
真核生物基因特征：
- 基因密度低（人类基因组仅约1-2%为编码区）
- 含内含子（人类基因平均含8个外显子、7个内含子）
- 5'端有帽子结构，3'端有polyA尾
- 单顺反子mRNA（一个转录本对应一个蛋白质）
- 复杂的剪接机制：内含子5'端为GT（供体位点），3'端为AG（受体位点）
- 核心启动子元件：TATA框（约-30区）、CAAT框（约-80区）
- 存在大量的可变剪接（人类约95%的多外显子基因发生可变剪接）

从头预测的核心是利用基因组的统计特征和信号特征：
（1）开放阅读框（ORF）检测
- 寻找起始于ATG、终止于TAA/TAG/TGA的连续密码子序列
- 计算ORF长度：真核生物外显子通常>300 bp（100个氨基酸），短ORF可能是假阳性
- 计算密码子使用频率（Codon Usage Frequency, CUF）：真实基因通常偏好使用特定密码子
（2）隐马尔可夫模型（HMM）
HMM是基因预测中最常用的统计框架：
- 隐藏状态：基因间区、外显子（起始/中间/终止）、内含子（起始/中间/终止）、UTR等
- 观测符号：A、C、G、T
- 状态转移概率：从一个状态转移到另一个状态的概率（如外显子→内含子的概率）
- 发射概率：在某个状态下观察到特定碱基的概率
- Viterbi算法：寻找最可能的隐藏状态序列（即基因结构）
经典HMM基因预测软件：
- Genscan：针对人类和脊椎动物，基于广义HMM
- Augustus：支持真核和原核生物，可整合外部证据
- Glimmer/GlimmerHMM：Glimmer用于原核，GlimmerHMM用于真核
- GeneMark/GeneMark-ES：自训练HMM，适合新物种
- SNAP：可训练参数，适合非模式生物
（3）支持向量机（SVM）和深度学习方法
近年来，机器学习方法也被用于基因预测：
- 将序列特征（k-mer频率、GC含量、密码子偏好等）编码为特征向量
- 用已知基因训练分类器（SVM、随机森林、神经网络）
- 对新序列进行预测

（1）蛋白质同源比对
- 将已知物种的蛋白质序列与目标基因组进行比对（如tblastn）
- 根据蛋白质比对结果推断基因的编码区
- 考虑剪接位点信号，将基因组DNA与外显子对应
- 代表工具：GeneWise、Exonerate（protein2genome模式）
（2）cDNA/EST比对
- 将已知的cDNA或EST序列与基因组比对（如BLAT、GMAP）
- 直接确定外显子-内含子边界
- EST数据质量参差不齐，需要过滤和聚类
（3）跨物种比较基因组学
- 利用进化保守性：编码区通常比非编码区更保守
- 利用密码子替换模式：同义替换率（Ks）和非同义替换率（Ka）
- 代表方法：Twinscan（整合人类-小鼠比较信息）、CONTRAST

（1）RNA-Seq数据比对与组装
- 将RNA-Seq reads比对到基因组（如HISAT2、STAR）
- 使用StringTie、Cufflinks等工具组装转录本
- 获得GTF格式的转录本结构注释
（2）全长转录本测序
- PacBio Iso-Seq或Nanopore direct RNA测序获得全长cDNA
- 无需组装，直接获得完整的外显子-内含子结构
- 可发现新的转录本异构体
（3）整合策略（Evidence Modeler, EVM）
- 将从不同来源（从头预测、同源比对、转录本组装）获得的基因模型进行整合
- 根据证据的可靠性和一致性进行加权评分
- 输出一致性最高的基因模型集合
- 代表工具：EVM（Evidence Modeler）、PASA（Program to Assemble Spliced Alignments）、BRAKER

非编码RNA（ncRNA）是不编码蛋白质的RNA分子，包括：
- miRNA（微小RNA）：~22 nt，调控基因表达
- lncRNA（长链非编码RNA）：>200 nt，多种调控功能
- circRNA（环状RNA）：共价闭合环，miRNA海绵等
- tRNA、rRNA、snRNA、snoRNA等
预测方法：
- miRNA：基于发夹结构、种子区域保守性、热力学稳定性（RNAfold、miRDeep2）
- lncRNA：基于开放阅读框长度、编码潜能（CPAT、CPC2）、序列保守性
- circRNA：基于反向剪接位点（backsplice junction）的reads支持（CIRI、find_circ）

GO数据库的核心是一个有向无环图（DAG）结构，其中：
- 每个节点（GO term）代表一个功能概念
- 边代表"is_a"（是一种）或"part_of"（是...的一部分）关系
- 子节点比父节点更具体（如"蛋白质激酶活性" is_a "转移酶活性"）
GO注释的常用方法：
（1）基于序列相似性的GO注释
- 用BLAST将查询蛋白与GO注释数据库（如UniProt-GOA）比对
- 将相似蛋白的GO term转移给查询蛋白
- 常用工具：Blast2GO、eggNOG-mapper
（2）基于结构域的GO注释
- 通过InterProScan识别蛋白质结构域
- 利用结构域与GO term的映射关系进行注释
- InterPro数据库提供了Pfam、SUPERFAMILY等结构与GO的映射
（3）基于机器学习的方法
- 使用蛋白序列特征（氨基酸组成、物理化学性质、亚细胞定位信号等）训练分类器
- 预测GO term
- 代表方法：DeepGO（基于深度学习）

KEGG数据库的核心是通路图（Pathway Map），包括：
- 代谢通路（Metabolism）：碳水化合物代谢、能量代谢、脂质代谢、核苷酸代谢、氨基酸代谢等
- 遗传信息处理（Genetic Information Processing）：转录、翻译、折叠/分选/降解、复制与修复
- 环境信息处理（Environmental Information Processing）：膜转运、信号转导、配体-受体相互作用
- 细胞过程（Cellular Processes）：细胞生长与死亡、细胞群落、运输与分解代谢
- 生物体系统（Organismal Systems）：免疫系统、神经系统、内分泌系统等
- 人类疾病（Human Diseases）：癌症、神经退行性疾病、心血管疾病等
KEGG注释方法：
- 序列比对法：用BLAST/KAAS（KEGG Automatic Annotation Server）将基因与KEGG GENES数据库比对
- KO（KEGG Orthology）注释：将基因映射到KO编号（功能单位），再根据KO编号映射到通路
- 工具：KAAS、KOFAM、eggNOG-mapper

InterPro整合了多个蛋白质特征数据库：
- Pfam：基于HMM的蛋白质家族和结构域数据库
- PRINTS：基于指纹（保守基序组）的数据库
- PANTHER：基于系统发育树的蛋白质家族分类
- ProSite：基于模式（pattern）和谱（profile）的功能位点数据库
- SUPERFAMILY：基于SCOP结构分类的远缘同源识别
- SMART：信号肽、跨膜区、结构域识别
InterProScan工作流程：
1. 用多种方法（HMM、模式匹配、profile等）扫描蛋白质序列
2. 整合不同数据库的预测结果
3. 为每个预测的domain/feature提供GO term映射
4. 输出GFF3或TSV格式的注释结果

（1）COG/KOG/EggNOG
- COG（Clusters of Orthologous Groups）：原核生物蛋白质直系同源簇
- KOG：真核生物版本
- EggNOG：扩展版本，覆盖更多物种，提供GO、KEGG、CAZy等注释
- 将蛋白质分为功能类别（如J: Translation, K: Transcription, S: Unknown）
（2）NR/UniProt/Swiss-Prot
- NR（Non-Redundant）：NCBI的非冗余蛋白数据库，用于序列相似性搜索
- UniProt：最全面的蛋白质数据库，包含Swiss-Prot（人工审编）和TrEMBL（自动注释）
- 注释流程：BLASTp比对NR/UniProt → 取top hit → 继承功能描述
（3）专用数据库
- CAZy：糖类活性酶分类（糖苷水解酶、糖基转移酶等）
- TCDB：膜转运蛋白分类
- VFDB：细菌毒力因子
- CARD：抗生素耐药基因
- antiSMASH：生物合成基因簇（次级代谢产物）

现代基因组项目通常使用自动化流程进行功能注释：

预测蛋白质序列
    ↓
1. InterProScan（结构域识别 + GO注释）
    ↓
2. KAAS/KOFAM（KEGG通路注释）
    ↓
3. eggNOG-mapper（COG/GO/KEGG综合注释）
    ↓
4. BLASTp against NR/Swiss-Prot（功能描述）
    ↓
5. 整合所有注释结果

6.3.1 基本原理

基因组变异按尺度可分为：
点突变尺度：
- SNP（Single Nucleotide Polymorphism）：单个碱基的替换，是最常见的变异类型
- Indel（Insertion/Deletion）：小片段（通常<50 bp）的插入或缺失
结构变异尺度：
- 缺失（Deletion）：>50 bp的序列缺失
- 插入（Insertion）：>50 bp的序列插入
- 倒位（Inversion）：序列方向的翻转
- 易位（Translocation）：序列从一个染色体位置移动到另一个位置
- 重复（Duplication）：序列拷贝数的增加
- CNV（Copy Number Variation）：>1 kb的拷贝数变化

原始测序数据 (FASTQ)
    ↓
质量控制 (FastQC, Trimmomatic)
    ↓
比对到参考基因组 (BWA, Bowtie2)
    ↓
比对后处理 (排序, 标记重复, 碱基质量校正)
    ↓
变异检测 (GATK, Samtools, FreeBayes)
    ↓
变异过滤和质量控制
    ↓
变异注释 (功能影响预测)
    ↓
最终变异集合 (VCF格式)

（1）排序（Sorting）
比对后的SAM/BAM文件需要按染色体坐标排序，以便后续处理。
（2）标记PCR重复（Mark Duplicates）
PCR扩增会引入系统性偏差——来自同一模板分子的多个reads并非独立采样。使用Picard MarkDuplicates等工具识别并标记PCR重复reads，避免重复reads对变异检测的干扰。
（3）Indel区域重比对（Indel Realignment）
比对软件在Indel附近可能产生比对偏差（如将Indel附近的SNP错误比对）。GATK的IndelRealigner通过识别Indel附近不一致的比对，进行局部重比对，提高Indel和周围SNP的检测准确性。
（4）碱基质量校正（Base Quality Score Recalibration, BQSR）
测序仪给出的原始碱基质量值可能存在系统性偏差（如某些测序循环质量系统性地偏低）。BQSR通过将碱基质量与已知变异数据库（如dbSNP）进行比较，建立校正模型，重新校准碱基质量值。

变异检测的核心是判断某个位置的碱基差异是真实的遗传变异还是测序/比对错误。这本质上是一个统计假设检验问题。
零假设（H₀）：该位置没有变异，所有观察到的差异都是测序错误
备择假设（H₁）：该位置存在变异
测序错误通常被建模为二项分布或贝叶斯模型：
- 设测序错误率为ε（通常ε≈0.001，对应Q30）
- 在覆盖度为D的位置，观察到k个非参考碱基
- 如果H₀成立，k ~ Binomial(D, ε)
- 计算P(k | H₀)，如果该P值小于阈值（如0.05），则拒绝H₀，判定为变异
基因型似然模型：
更精确的方法是对三种可能的基因型（参考纯合RR、杂合RA、变异纯合AA）分别计算似然值：
- L(RR) = ∏ P(read_i | RR)
- L(RA) = ∏ P(read_i | RA)
- L(AA) = ∏ P(read_i | AA)
选择似然值最大的基因型作为该位点的基因型。
GATK HaplotypeCaller使用更先进的单倍型模型，在局部区域（通常数百bp）同时考虑所有可能的变异，通过局部de novo组装生成候选单倍型，然后计算每个单倍型的似然值。

即使使用了严格的统计模型，变异检测结果中仍包含大量假阳性。需要通过多重过滤来提高准确性：
基于测序质量的过滤：
- QD（Quality by Depth）：质量值除以覆盖度，过低提示质量不可靠
- FS（Fisher Strand）：检验正负链reads是否均衡支持变异，不均衡提示比对偏差
- MQ（Mapping Quality）：比对质量值，过低提示比对不可靠
- ReadPosRankSum：检验变异是否主要出现在reads末端（末端碱基质量通常较低）
基于覆盖度的过滤：
- DP（Depth）：覆盖度过低（<5-10×）的变异不可靠
- 覆盖度过高（>平均深度的2-3倍）可能提示重复区域
基于群体频率的过滤：
- 使用已知变异数据库（如gnomAD）过滤常见变异
- 在罕见病研究中，过滤在人群中频率>1%的变异

SNP主要来源于DNA复制过程中的自发突变：
（1）脱氨基作用
- 胞嘧啶（C）自发脱氨基转变为尿嘧啶（U），DNA修复系统将U识别为T，导致C→T突变
- 这是转换突变的主要来源，解释了为什么转换/颠换比（Ti/Tv ratio）通常大于2
- 5-甲基胞嘧啶（5mC）的脱氨基率更高，导致CpG岛中的C→T突变频率显著升高
（2）复制错误
- DNA聚合酶在复制过程中偶尔掺入错误的碱基
- 尽管有校对功能，错误率仍约为10⁻⁹/碱基/复制
（3）氧化损伤
- 活性氧（ROS）可导致碱基氧化修饰，如8-oxo-guanine可导致G→T颠换

按功能位置分类：
| 位置 | 比例 | 功能影响 |
|------|------|---------|
| 基因间区 | ~90% | 通常无直接功能影响 |
| 内含子 | ~5% | 可能影响剪接调控 |
| 外显子同义 | ~1% | 不改变氨基酸（密码子简并性） |
| 外显子错义 | ~1.5% | 改变氨基酸序列 |
| 外显子无义 | ~0.1% | 引入提前终止密码子 |
| 调控区 | ~2% | 可能影响基因表达 |
| 剪接位点 | <0.1% | 可能导致剪接异常 |
按群体频率分类：
- 常见变异（Common Variant）：群体频率>5%，通常与复杂疾病的易感风险相关
- 低频变异（Low-frequency Variant）：群体频率1-5%
- 稀有变异（Rare Variant）：群体频率<1%，可能与孟德尔遗传病相关

（1）基于pileup的方法（Samtools/bcftools）
- 统计每个基因组位置上比对reads的碱基组成
- 使用贝叶斯模型计算每个可能基因型的后验概率
- 优点：速度快，计算资源需求低
- 缺点：不处理Indel附近的比对偏差，对Indel不敏感
（2）基于单倍型的方法（GATK HaplotypeCaller）
- 在局部窗口（数百bp）内进行de novo组装，生成候选单倍型
- 将reads与候选单倍型进行比对，计算每个单倍型的似然值
- 选择最可能的单倍型组合作为基因型
- 优点：对Indel和复杂区域的SNP检测更准确
- 缺点：计算量大，耗时较长
（3）基于深度学习的方法（DeepVariant）
- 将每个候选变异位点转化为"图像"（read pileup的视觉表示）
- 使用卷积神经网络（CNN）进行变异识别
- 优点：准确性极高，在多个基准测试中优于传统方法
- 缺点：计算资源需求高，训练数据依赖

QD（Quality by Depth）：
QD = QUAL / DP
- QD过低（<2）提示变异可能被低质量reads或重复区域支持
- QD过高（>40）可能提示比对到重复区域
MQ（Mapping Quality）：
- 反映reads在该区域的比对可靠性
- MQ < 40提示比对可能存在问题
FS（Fisher Strand Bias）：
- 检验变异是否被正负链均衡支持
- FS > 60提示严重的链偏倚，变异可能不可靠
SOR（Strand Odds Ratio）：
- FS的替代指标，对覆盖度不均更稳健
- SOR > 3提示链偏倚
MQRankSum：
- 检验支持变异的reads与支持参考的reads的比对质量是否存在显著差异
- 负值过大提示支持变异的reads比对质量较差
ReadPosRankSum：
- 检验变异是否主要出现在reads末端
- 末端碱基质量通常较低，可能导致假阳性

转换/颠换比率（Ti/Tv ratio）是评估SNP集合质量的重要指标：
- 全基因组范围的期望Ti/Tv：约2.0-2.2
- 外显子区域的期望Ti/Tv：约2.8-3.0（因为同义SNP只发生在特定密码子位置）
如果检测到的SNP集合Ti/Tv显著低于期望值，可能提示：
- 存在大量假阳性（特别是颠换类型的假阳性较多）
- 过滤条件不够严格
- 样本存在污染

（1）DNA聚合酶滑移（Replication Slippage/Strand Slippage）
- 在重复序列区域（尤其是微卫星），模板链和新生链可能发生相对滑动
- 导致重复单元的插入或缺失
- 这是Indel最主要的产生机制，解释了为什么Indel在重复序列区域富集
（2）不等交换（Unequal Crossing Over）
- 减数分裂过程中的同源重组发生错位
- 导致一个染色体上重复序列增加，另一个上减少
（3）DNA损伤修复错误
- 双链断裂修复过程中可能发生小片段的插入或缺失

Indel检测比SNP检测困难的主要原因：
（1）比对困难
- 含有Indel的reads与参考基因组比对时，需要在Indel位置引入"空位"（gap）
- 比对算法对空位有罚分，倾向于将Indel区域的碱基错配解释为SNP
- 例如，参考序列"ATCGATCG"，reads为"ATC---GATCG"（3bp缺失），比对软件可能错误地比对为"ATCGATCG"（多个错配）
（2）Indel周围的碱基质量下降
- Illumina测序在Indel附近容易产生碱基识别错误
- 需要专门的Indel重比对步骤来纠正
（3）Indel大小的不确定性
- 对于较大的Indel（>10 bp），单个read可能无法完全跨越，导致Indel边界不确定
- 需要多个reads的信息来确定Indel的精确边界

（1）基于局部重比对的Indel检测（GATK HaplotypeCaller）
- HaplotypeCaller在局部窗口内生成候选单倍型（包含可能的SNP和Indel组合）
- 将reads与所有候选单倍型比对，选择最优解释
- 优点：同时考虑SNP和Indel，避免"SNP vs Indel"的比对歧义
（2）基于split-read的Indel检测（Pindel）
- 识别"split reads"——一条read的两部分分别比对到Indel两侧
- 通过分析未比对上的read片段来推断Indel序列
- 优点：对较大Indel（>20 bp）的边界定位准确
- 缺点：只能检测被reads跨越的Indel
（3）基于de novo局部组装的Indel检测（Platypus, Dindel）
- 在候选Indel区域提取相关reads
- 进行局部de novo组装
- 将组装结果与参考基因组比对，识别Indel

Indel的功能影响通常比SNP更严重：
| Indel类型 | 功能影响 | 示例 |
|----------|---------|------|
| 编码区非3倍数Indel | 移码突变，通常导致功能丧失 | BRCA1移码突变导致乳腺癌风险 |
| 编码区3倍数Indel | 氨基酸插入/缺失，可能影响功能 | CFTR ΔF508导致囊性纤维化 |
| 剪接位点Indel | 影响pre-mRNA剪接 | 深度内含子Indel激活隐蔽剪接位点 |
| 调控区Indel | 影响转录因子结合 | 胰岛素基因启动子Indel影响表达 |
| UTR区Indel | 可能影响mRNA稳定性或翻译 | 5'UTR的uORF相关Indel |
经典病例：
- 囊性纤维化：CFTR基因最常见的突变是ΔF508（一个3bp缺失），导致苯丙氨酸缺失，约占所有囊性纤维化病例的70%
- 亨廷顿舞蹈症：HTT基因中CAG重复序列的扩增（可视为一种特殊Indel），当CAG重复数>36时导致疾病

（1）基于Read深度（Read Depth, RD）的方法
原理：测序深度与基因组拷贝数成正比。
- 缺失区域：reads深度显著低于平均值（约50%对于杂合缺失）
- 重复区域：reads深度显著高于平均值
代表软件：CNVnator、BIC-seq、CNV-seq、DELLY（结合RD和PE）
优点：
- 可以检测大尺度的CNV
- 计算相对简单
局限：
- 无法精确定位断点（只能确定区域）
- 无法检测平衡性SV（如倒位、易位不改变拷贝数）
- 受GC含量、比对偏好等因素影响
（2）基于Split Read（SR）的方法
原理：识别一条read的两部分分别比对到不同位置的"split reads"，从而精确定位断点。
代表软件：Pindel、DELLY、LUMPY、Manta
工作流程：
1. 提取未完全比对或部分比对的reads（soft-clipped reads）
2. 将未比对部分重新比对到基因组其他位置
3. 如果找到匹配，确定断点坐标
优点：
- 断点定位精确（可达bp级别）
- 可以检测小至50 bp的SV
局限：
- 只能检测被reads跨越的SV
- 需要足够长的reads（对Illumina短读长有限制）
（3）基于Paired-End（PE）的方法
原理：利用paired-end测序中insert size和方向的异常来检测SV。
正常配对特征：
- 两个read在同一染色体上
- 方向相对（read1正向，read2反向）
- 距离在预期insert size范围内
异常配对类型及对应的SV：
| 异常类型 | 提示的SV |
|---------|---------|
| 距离远大于预期 | 缺失（reads跨越了被删除的区域） |
| 距离远小于预期 | 插入（reads之间的区域被插入序列替代） |
| 方向相同 | 倒位 |
| 位于不同染色体 | 易位 |
代表软件：DELLY、LUMPY、BreakDancer、Manta
优点：
- 可以检测大尺度SV
- 对read长度要求不高
局限：
- 断点定位不精确（通常在数百bp范围内）
- 受insert size分布的影响
（4）基于de novo组装的方法
原理：对SV区域进行局部de novo组装，将组装结果与参考基因组比对来识别SV。
代表软件：TIGRA、MindTheGap、SMRT-SV（PacBio专用）、Sniffles（Nanopore专用）
优点：
- 可以检测最复杂的SV，包括串联重复扩增、倒位等
- 可以获得SV的完整序列
局限：
- 计算资源需求大
- 组装质量依赖于数据质量和覆盖度

现代SV检测工具通常整合多种信号：
DELLY：整合PE、SR和RD信号
LUMPY：整合PE和SR信号，使用概率框架
Manta：结合PE、SR和局部组装，速度快，准确性高
GRIDSS：基于基因组图（assembly graph）的SV检测

长读长测序（PacBio/Nanopore）对于SV检测具有革命性优势：
1. 跨越复杂SV：长读长可以直接跨越整个SV区域，获得完整的SV序列
2. 精确定位断点：无需依赖PE或SR推断，直接比对即可定位
3. 检测重复区域SV：短读长难以比对的重复区域，长读长可以跨越
专用工具：
- Sniffles：Nanopore SV检测
- PBSV：PacBio SV检测
- SVIM：同时支持PacBio和Nanopore

SV检测的假阳性率通常高于SNP/Indel，验证至关重要：
实验验证方法：
- PCR+Sanger测序：对断点区域进行PCR扩增和测序，是SV验证的金标准
- qPCR：验证拷贝数变化
- Southern blot：验证大尺度重复/缺失
- 光学图谱（Bionano）：验证大尺度SV
计算验证方法：
- 使用不同检测工具的结果取交集
- 检查SV区域的reads比对情况（IGV可视化）
- 与已知SV数据库（如gnomAD-SV、DGV）比对

6.4.1 微生物组概述

人体微生物组的规模：
- 总重量：约1-2 kg（主要在肠道）
- 细胞数量：约3×10¹³个（与人体自身细胞数量相当或略多）
- 基因数量：约300-500万个基因（是人类基因组基因数的150-250倍）
人体各部位的微生物组特征：
| 部位 | 主要门类 | 细菌密度 | 主要功能 |
|------|---------|---------|---------|
| 结肠 | Firmicutes, Bacteroidetes | 10¹¹-10¹²/g | 发酵膳食纤维、合成维生素 |
| 口腔 | Firmicutes, Proteobacteria | 10⁹/mL唾液 | 参与氮循环、免疫调节 |
| 皮肤 | Actinobacteria, Firmicutes | 10⁶-10⁷/cm² | 屏障保护、免疫训练 |
| 阴道 | Lactobacillus | 10⁸-10⁹/mL | 维持酸性环境、防御病原 |
| 鼻腔 | Firmicutes, Actinobacteria | 10⁶-10⁷/cm² | 免疫防御 |

传统微生物学估计，环境中可培养的微生物仅占1%左右（"Great Plate Count Anomaly"）。主要原因包括：
1. 营养需求未知：很多微生物需要特殊的生长因子或复杂的共生关系
2. 环境条件难以模拟：如极端pH、温度、压力、厌氧条件等
3. 生长缓慢：某些微生物代时可达数周甚至数月
4. 群体感应依赖：需要达到一定密度才能生长
5. 严格的共生关系：依赖其他微生物产生的代谢产物
宏基因组学通过直接提取环境DNA进行测序，彻底 bypass 了培养瓶颈。

16S rRNA（原核生物）和18S rRNA/ITS（真核微生物）是微生物分类鉴定的"条形码"。
16S rRNA的优势：
1. 普遍存在：所有细菌和古菌都含有16S rRNA基因
2. 功能保守：核糖体功能对生命至关重要，序列演化相对缓慢
3. 大小适中：约1,500 bp，既包含保守区（设计通用引物），又包含可变区（区分物种）
16S rRNA的可变区：
- V1-V9共9个可变区
- 常用测序区域：V3-V4（约465 bp，Illumina MiSeq 2×300可覆盖）
- 不同可变区的分辨能力不同：V1-V2适合属水平，V3-V4可区分到种/株水平
分类阈值：
- 16S rRNA序列相似性 >97%：通常认为是同一个种
- 16S rRNA序列相似性 >95%：通常认为是同一个属

原始测序数据 (FASTQ)
    ↓
质量控制 (FastQC, cutadapt)
    ↓
去接头、去低质量、去嵌合体 (DADA2/Deblur/QIIME2)
    ↓
生成ASV/OTU表
    ↓
物种注释 (SILVA/Greengenes/RDP数据库)
    ↓
多样性分析 (α, β)
    ↓
统计分析 (LEfSe, ANCOM, DESeq2)
    ↓
可视化 (phyloseq, ggplot2)

（1）DADA2去噪算法
DADA2（Divisive Amplicon Denoising Algorithm 2）通过以下步骤生成ASV：
1. 质控过滤：去除低质量和含N的reads
2. 去重复：将相同序列合并，记录丰度
3. 错误学习：从数据中估计测序错误模型
4. 去噪：将含错误的序列纠正回其真实序列
5. 去嵌合体：识别并去除PCR嵌合体
6. 合并双端reads
（2）OTU聚类 vs ASV
- OTU聚类以97%相似性为阈值，将序列分组。缺点：阈值是人为设定的；近缘种可能因高于阈值而无法区分；信息在聚类过程中丢失
- ASV保留了每个独特的真实序列，分辨率更高，且可重复性更好

原始测序数据 (FASTQ)
    ↓
质量控制 (FastQC, Trimmomatic)
    ↓
去宿主DNA (Bowtie2比对宿主参考基因组)
    ↓
宏基因组组装 (MEGAHIT, metaSPAdes)
    ↓
基因预测 (Prodigal, MetaGeneMark)
    ↓
基因 catalog构建 (CD-HIT聚类)
    ↓
功能注释 (eggNOG-mapper, KEGG)
    ↓
分箱 (MetaBAT2, MaxBin2, CONCOCT)
    ↓
MAG质量评估 (CheckM)
    ↓
物种分类和系统发育分析

分箱的目标是将来自同一基因组的contig聚类到一起。
（1）基于组成的分箱（Composition-based Binning）
- 原理：不同物种的基因组具有不同的碱基组成特征（如GC含量、四核苷酸频率）
- 方法：将contig表示为高维特征向量（如136维四核苷酸频率向量），使用聚类算法（如K-means、高斯混合模型）进行聚类
- 代表工具：MetaBAT2、MaxBin2
（2）基于覆盖度的分箱（Coverage-based Binning）
- 原理：来自同一物种的contig在不同样本中的覆盖度变化模式一致（因为它们来自同一个基因组）
- 方法：在多个样本中计算每个contig的覆盖度，构建覆盖度矩阵，进行聚类
- 代表工具：CONCOCT、MaxBin2
（3）混合分箱
- 同时使用组成和覆盖度信息进行分箱
- 代表工具：MetaBAT2（使用TNF和覆盖度）、VAMB（使用变分自编码器）
MAG质量评估：
- CheckM：基于保守单拷贝基因（SCG）评估MAG的完整性和污染度
- 高质量MAG：完整性>90%，污染度<5%
- 中等质量MAG：完整性≥50%，污染度<10%

α多样性指数：
| 指数 | 衡量 | 公式 | 特点 |
|------|------|------|------|
| 物种数（Richness） | 丰富度 | S = 物种数量 | 不考虑丰度 |
| Shannon | 丰富度+均匀度 | H' = -Σ(pi × ln pi) | 对丰富度敏感 |
| Simpson | 丰富度+均匀度 | D = 1 - Σ(pi²) | 对优势种敏感 |
| Chao1 | 估计总丰富度 | Chao1 = S_obs + n₁²/(2n₂) | 考虑稀有物种 |
| Pielou's J | 均匀度 | J' = H'/ln(S) | 纯度的衡量 |
β多样性距离度量：
| 距离 | 公式 | 特点 |
|------|------|------|
| Bray-Curtis | 1 - 2Σmin(xi, yi) / Σ(xi + yi) | 考虑丰度，最常用 |
| Jaccard | 1 - |A∩B| / |A∪B| | 仅考虑有无 |
| Unweighted UniFrac | 考虑系统发育距离 | 考虑进化关系 |
| Weighted UniFrac | 考虑丰度和系统发育 | 综合考量 |
β多样性可视化：
- PCoA（主坐标分析）
- NMDS（非度量多维尺度分析）
- UPGMA聚类树

识别在不同组间丰度显著差异的微生物或功能基因：
- LEfSe（Linear Discriminant Analysis Effect Size）：结合统计检验和LDA评分，筛选组间差异标志物
- ANCOM（Analysis of Composition of Microbiomes）：考虑微生物组数据的组成性特点
- DESeq2/edgeR：基于负二项分布的差异分析（适用于OTU/基因计数数据）
- ALDEx2：考虑组成性和稀疏性的差异分析

微生物组数据具有以下"大数据"特征：
1. 高维度：数千个OTU/ASV × 数百个样本
2. 稀疏性：很多物种只在少数样本中出现（零膨胀分布）
3. 组成性：相对丰度总和为100%，变量间存在约束
4. 高变异性：个体差异大，受饮食、环境等多种因素影响
5. 多层次：物种、功能、代谢物等多个层面的数据

（1）松散的整合（Loose Integration）
- 分别分析各组学数据
- 比较不同组学层面的结果，寻找一致性
- 例如：比较差异物种和差异代谢物，寻找相关性
（2）统计学整合（Statistical Integration）
- 使用多变量统计方法整合多个数据矩阵
- CCA（Canonical Correlation Analysis）：寻找两个数据集之间的最大相关
- sPLS（sparse Partial Least Squares）：寻找预测关系的同时进行变量选择
- MCIA（Multiple Co-Inertia Analysis）：整合多个组学数据集
- MOFA（Multi-Omics Factor Analysis）：分解共享和特异的变异源
（3）网络整合（Network Integration）
- 构建跨组学的关联网络
- 例如：物种-代谢物关联网络、基因-代谢物关联网络
- 识别网络中的关键节点（hub）
（4）机器学习整合
- 将多组学特征作为输入，训练预测模型
- 使用特征选择识别最重要的预测因子

（1）分类任务：疾病预测
- 输入：微生物组组成特征（物种丰度、功能通路丰度）
- 输出：疾病状态（如IBD vs 健康）
- 常用算法：
- 随机森林（Random Forest）：对高维稀疏数据表现好，可输出特征重要性
- 支持向量机（SVM）：适合小样本高维数据
- 深度学习（神经网络）：可学习复杂的非线性关系
案例：使用随机森林预测CRC（结直肠癌）
- 训练集：数百个粪便样本的16S/宏基因组数据
- 特征：物种丰度（或MetaCyc通路丰度）
- 模型：随机森林分类器
- 结果：AUC可达0.8-0.9，某些研究报道准确率>90%
- 关键标志物：Fusobacterium nucleatum、Peptostreptococcus anaerobius等
（2）回归任务：表型预测
- 输入：微生物组特征 + 宿主因素
- 输出：连续变量（如BMI、药物响应、血糖水平）
- 常用算法：LASSO回归、弹性网络、随机森林回归
（3）聚类任务：亚型鉴定
- 输入：微生物组组成
- 输出：患者亚型（enterotypes，肠型）
- 肠型概念：人类肠道微生物组可分为几种主要的组成模式
- 肠型1：以Bacteroides为主
- 肠型2：以Prevotella为主
- 肠型3：以Ruminococcus为主

微生物组研究面临的一个核心挑战是因果关系的确定：
- 相关性不等于因果性：微生物组变化可能是疾病的原因，也可能是疾病的结果
- 混杂因素：饮食、药物、生活方式等同时影响微生物组和疾病状态
- 策略：
- 纵向研究：追踪疾病发生前后的微生物组变化
- 动物模型：无菌小鼠或抗生素处理后的定植实验
- 孟德尔随机化：利用宿主基因型作为工具变量
- 干预研究：益生菌、FMT等干预后的效果评估

（1）疾病诊断
- 结直肠癌（CRC）筛查：F. nucleatum等细菌标志物+粪便免疫化学检测（FIT）联合筛查
- 艰难梭菌感染（CDI）诊断：结合临床表现和微生物组特征
- 非酒精性脂肪性肝病（NAFLD）评估：肠道菌群组成与肝纤维化程度相关
（2）疾病治疗
- 复发性CDI：FMT治愈率>90%，已被FDA批准
- 溃疡性结肠炎（UC）：多项RCT显示FMT对UC有一定疗效
- 自闭症谱系障碍（ASD）：小规模临床试验显示FMT可改善胃肠道症状和行为
- 黑色素瘤免疫治疗：肠道菌群组成影响PD-1抑制剂的疗效
（3）药物代谢
- 肠道微生物组参与多种药物的代谢，影响药效和毒性
- 例如：伊立替康（irinotecan）的毒性由肠道细菌的β-葡萄糖醛酸酶介导
- 通过抑制特定细菌酶，可以降低药物毒性

（1）植物微生物组与农业
- 根际微生物组：影响植物养分吸收、抗病性和抗逆性
- 生物肥料：固氮菌、解磷菌等促进植物生长
- 生物防治：利用拮抗微生物防治植物病原菌
- 案例：印度矮生小麦与根际促生菌（PGPR）联合使用，可减少化肥用量30%
（2）动物微生物组与畜牧业
- 反刍动物瘤胃微生物组：决定饲料转化效率
- 益生菌：改善家禽和家畜的生长性能和免疫力
- 替代抗生素：欧盟禁用抗生素促生长剂后，微生物组产品成为替代方案

（1）生物修复
- 利用环境微生物组降解污染物
- 石油污染修复：海洋微生物组中的烃降解菌可降解石油烃
- 重金属污染修复：某些微生物可将有毒重金属转化为低毒形态
（2）功能基因挖掘
- 从环境宏基因组中发现新的酶和生物活性物质
- 纤维素酶：从瘤胃、白蚁肠道等环境中发现高效纤维素降解酶
- 新型抗生素：从土壤宏基因组中发现新结构抗生素（如Teixobactin）
（3）碳循环和气候变化
- 土壤微生物组在全球碳循环中发挥关键作用
- 冻土融化释放古老微生物，可能加速温室气体排放
- 海洋微生物组（特别是浮游植物共生菌）影响全球碳固定

（1）益生菌设计
- 传统益生菌：Lactobacillus, Bifidobacterium等
- 下一代益生菌（Next-Generation Probiotics）：Akkermansia muciniphila, Faecalibacterium prausnitzii等
- 工程化益生菌：通过基因工程赋予益生菌特定功能（如产丁酸盐、分泌治疗性蛋白）
（2）合成微生物组
- 目标：构建具有稳定功能的最小微生物群落
- 方法：
- 自上而下：从复杂群落中简化出核心功能群
- 自下而上：从单一菌株开始，逐步构建稳定群落
- 挑战：微生物间互作复杂，群落稳定性难以预测

（1）样本采集和保存
- 微生物组在采样后迅速变化（如氧气暴露导致厌氧菌死亡）
- 保存方法影响结果：冷冻、RNAlater、FTA卡等各有优缺点
- 建议：采样后尽快冷冻（-80°C），或使用稳定化试剂
（2）DNA提取偏差
- 不同试剂盒对不同类型细菌的提取效率差异很大
- 革兰氏阳性菌（厚壁菌门）细胞壁含大量肽聚糖，比革兰氏阴性菌更难裂解
- 策略：在同一研究中统一使用同一种提取方法；使用裂解对照（spike-in controls）
（3）PCR扩增偏差
- 通用引物对不同物种的扩增效率不同
- 嵌合体（chimera）和异源双链（heteroduplex）的形成
- 策略：使用高保真酶、优化PCR循环数、使用去嵌合体算法
（4）测序深度不足
- 复杂样本（如土壤）可能需要>10 Gb数据才能充分覆盖稀有物种
- 深度不足会导致稀有物种被遗漏，影响多样性估计
（5）生物信息学分析的异质性
- 不同分析流程（QIIME vs mothur）、不同数据库（SILVA vs Greengenes）会产生不同结果
- 缺乏统一的金标准

（1）样本量计算
- 微生物组研究通常需要比传统临床研究更大的样本量
- 因为个体差异大，效应量通常较小
- 建议：初步探索性研究每组至少20-30个样本；验证性研究每组50-100个
（2）对照设置
- 阴性对照（空白对照）：检测试剂污染
- 阳性对照（mock community）：评估技术流程的准确性
- 裂解对照（spike-in）：评估DNA提取效率
（3）批次控制
- 随机化：将不同组的样本随机分配到不同批次
- 平衡设计：确保每批次包含所有组的样本
- 记录所有批次信息，以便后续统计校正
- 统计校正：使用ComBat等方法校正批次效应
（4）纵向设计
- 对于疾病研究，纵向采样（疾病前-中-后）比横断面设计更能揭示因果关系
- 建议：基线（T0）、急性期（T1）、恢复期（T2）

标准化倡议：
- MIxS标准（Minimum Information about any (x) Sequence）：包括MIMARKS、MIMS等，规定了宏基因组研究需要报告的元数据
- Earth Microbiome Project（EMP）：制定了标准的采样和测序方案
数据共享平台：
- NCBI SRA/EBI ENA：原始测序数据
- Qiita/MG-RAST：处理后的微生物组数据
- 数据共享是提高研究可重复性和促进元分析的关键

（1）方法学进展
- 长读长宏基因组学：Nanopore/PacBio长读长有助于解决重复序列和分箱问题
- 单细胞宏基因组学：在单细胞分辨率上解析微生物群落
- 空间宏基因组学：保留微生物的空间分布信息
（2）人工智能与微生物组
- 深度学习用于物种分类和功能预测
- 生成式AI用于设计合成微生物组
- 大语言模型整合文献知识和微生物组数据
（3）精准微生物组医学
- 基于个体微生物组特征的个性化益生菌
- 微生物组作为疾病早期诊断的生物标志物
- 微生物组靶向治疗（如噬菌体疗法、CRISPR编辑肠道菌）
（4）全球微生物组计划
- 类似于人类基因组计划，全球微生物组计划旨在系统性地绘制地球微生物多样性图谱
- 应用：生物勘探、环境监测、气候变化预测

基因组学作为生物信息学的核心领域，在过去二十年中经历了前所未有的技术革命和方法学进步。从人类基因组计划的完成到个人基因组测序进入千元时代，从短读长测序主导到长读长技术崛起，基因组学的发展深刻改变了生命科学研究的面貌。

本章核心要点回顾

测序技术层面：Illumina短读长测序以其高准确性和高通量继续占据主流地位，而PacBio HiFi和Nanopore长读长测序则在复杂基因组组装和结构变异检测中展现出独特优势。混合测序策略正在成为高质量基因组项目的标准配置。
组装与注释层面：de Bruijn图算法是二代测序组装的核心方法，而OLC算法在长读长时代重新焕发活力。Hi-C技术使染色体级别组装成为可能。基因预测已从单纯的从头预测发展为整合多种证据的共识预测策略。功能注释依赖于同源性比对、结构域识别和通路映射等多种方法。
变异检测层面：SNP检测已达到极高的准确性，Indel检测需要特殊的比对处理，而SV检测仍是挑战性最大的领域。长读长测序为SV检测带来了突破性进展。
宏基因组学层面：宏基因组学通过bypass培养瓶颈，揭示了微生物世界的巨大多样性。从16S扩增子测序到宏基因组shotgun测序，从多样性分析到功能挖掘，宏基因组学正在从描述性研究向机制性研究和应用转化迈进。

未来展望

完整基因组时代：端粒到端粒（T2T）组装技术的成熟将使每个物种都能获得真正完整的基因组参考序列，包括着丝粒、端粒和核糖体DNA等 previously intractable 区域。
泛基因组学：单一参考基因组已不足以代表物种内的遗传多样性。泛基因组（pan-genome）研究将揭示核心基因组和可变基因组的动态演化，为作物改良和疾病研究提供新的视角。
单细胞基因组学：单细胞测序技术的进步将使研究者能够在单细胞分辨率上解析基因组变异，揭示肿瘤异质性、发育谱系和微生物群落结构。
人工智能与基因组学：深度学习已经在变异检测（DeepVariant）、蛋白质结构预测（AlphaFold）等领域取得突破。未来，AI将进一步整合基因组、转录组、蛋白质组和表观遗传组数据，实现从序列到功能的端到端预测。
精准基因组医学：随着测序成本的进一步下降和数据分析方法的成熟，全基因组测序有望成为临床常规检测手段，为罕见病诊断、肿瘤精准治疗和药物基因组学提供个体化的遗传信息。
基因组学的发展远未止步。正如人类基因组计划的领导者之一Francis Collins所言："基因组测序只是开始，真正的挑战和机遇在于解读这本生命之书。"作为生物信息学工作者，我们正处于这一伟大事业的核心，肩负着从海量基因组数据中提取生物学洞见、推动精准医学发展的历史使命。

本章思考题

1. 回顾基因组学发展的历史，你认为哪三项技术突破对领域的影响最大？为什么？
2. 如果你要启动一个全新物种的基因组项目，请设计一个完整的技术路线，包括测序策略、组装方案和注释流程。
3. 比较SNP、Indel和SV在检测难度、功能影响和进化意义方面的异同。
4. 宏基因组学如何改变了我们对微生物世界和宿主-微生物互作的理解？
5. 展望基因组学的未来，你认为下一个重大突破可能出现在哪个方向？

推荐阅读
1. Goodwin S, McPherson J D, McCombie W R. Coming of age: ten years of next-generation sequencing technologies [J]. Nature reviews genetics, 2016, 17: 333-351.
2. Nurk S, Koren S, Rhie A, et al. The complete sequence of a human genome [J]. Science, 2022, 376: 44-53. （T2T-CHM13）
3. Quince C, Walker A W, Simpson J T, et al. Shotgun metagenomics, from sampling to analysis [J]. Nature biotechnology, 2017, 35: 833-844.
4. Turner T N, Hormozdiari F, Duyzend M H, et al. Genome sequencing of autism-affected families reveals disruption of putative noncoding regulatory DNA [J]. American journal of human genetics, 2016, 98: 58-74.

（1）测序技术的三次革命
第一代测序（Sanger, 1977）：
- 原理：双脱氧链终止法
- 特点：准确性极高（>99.99%），读长可达1000 bp，但通量极低
- 贡献：完成首批模式生物基因组（噬菌体φX174、流感嗜血杆菌、酿酒酵母）
- 局限：人类基因组计划耗时13年、耗资27亿美元
第二代测序（Illumina, 2005）：
- 原理：桥式PCR + 可逆终止子边合成边测序
- 特点：通量提升百万倍，读长100-300 bp，准确性>99.9%，成本降至$1000/人基因组
- 贡献：千人基因组计划、肿瘤基因组图谱（TCGA）、宏基因组学兴起
- 局限：短读长无法跨越重复序列和复杂区域
第三代测序（PacBio/Nanopore, 2011-2015）：
- 原理：单分子实时测序（ZMW荧光检测 / 纳米孔电流检测）
- 特点：读长可达数Mb，无需PCR扩增，可检测碱基修饰
- 贡献：T2T基因组、复杂结构变异解析、直接RNA测序
- 局限：原始错误率较高（PacBio HiFi已解决），成本仍高于二代
（2）组装算法的演进
- 2001-2008：基于OLC的Sanger组装（Celera Assembler, Arachne）
- 2008-2015：基于DBG的二代短读长组装（Velvet, SOAPdenovo, SPAdes）
- 2015-2020：长读长OLC组装（Canu, FALCON）
- 2020至今：HiFi精准长读长组装（Hifiasm, Verkko）+ T2T整合策略
（3）变异检测的精度跃升
- 早期：基于pileup的统计方法（Samtools mpileup）
- 中期：单倍型模型（GATK HaplotypeCaller）
- 近期：深度学习（DeepVariant, Google Brain）——在多个基准测试中达到或超越人类专家水平

（1）数据层面的挑战
- 数据量爆炸：全球测序数据已超100 EB（Exabyte），存储和传输成本急剧上升
- 数据异质性：不同平台（Illumina/PacBio/Nanopore）、不同版本（hg19/hg38/T2T-CHM13）的数据整合困难
- 数据质量不均：公共数据库中大量低质量组装和错误注释存在
（2）分析层面的挑战
- 计算资源需求：T2T组装、泛基因组构建、单细胞分析需要PB级内存和万核级计算
- 算法瓶颈：高度重复序列、复杂结构变异、多倍体基因组仍是组装的难点
- 标准化缺失：不同分析流程的结果可比性差，缺乏统一的金标准
（3）解读层面的挑战
- 功能未知基因：即使在人类基因组中，仍有数千个基因功能不明
- 非编码区解读：增强子、沉默子、绝缘子等调控元件的精确预测和功能验证困难
- 变异临床意义：大量罕见变异的致病性难以判定（ACMG指南仍依赖专家判断）
- 多基因复杂疾病：大多数常见疾病（糖尿病、高血压、精神疾病）受数千个微效变异影响，单一基因解释力有限

（1）完整基因组时代（T2T for All）
- 目标：为地球上所有物种（>1000万种真核生物）提供T2T参考基因组
- 驱动力：地球生物基因组计划（EBP）、Zoonomia项目（哺乳类）、达尔文树计划（中国）
- 技术：HiFi + ONT + Hi-C的标准化流程，自动化组装和审校
- 意义：理解生物多样性、保护濒危物种、挖掘新资源
（2）泛基因组与图基因组
- 单一参考基因组已不能满足研究和临床应用的需求
- 图基因组将线性参考升级为包含群体变异的图结构，减少比对偏差
- 人类泛基因组参考联盟（HPRC）正在构建包含>40个多样化单倍型的泛基因组
- 应用：改进变异检测、发现参考基因组中缺失的序列、群体遗传学研究
（3）单细胞基因组学
- 从"bulk平均信号"到"单细胞分辨率"
- 单细胞基因组测序揭示肿瘤异质性、发育谱系、微生物群落结构
- 技术挑战：全基因组扩增的偏倚、低覆盖度下的变异检测、海量数据的分析
- 前沿方向：单细胞多组学（同时测基因组+转录组+表观组+蛋白质组）
（4）人工智能与基因组学
- 变异检测：DeepVariant、Clair3等深度学习工具已达到金标准水平
- 蛋白质结构预测：AlphaFold2革命性地解决了蛋白质折叠问题
- 调控元件预测：Enformer、Basenji等模型从序列预测基因表达和调控
- 端到端预测：从DNA序列直接预测功能（如SpliceAI预测剪接位点）
- 大语言模型：基于Transformer的基因组语言模型（如DNABERT、Nucleotide Transformer）正在学习DNA序列的"语义"
（5）精准基因组医学
- 罕见病诊断：全基因组测序（WGS）已成为未确诊遗传病的首选诊断工具，诊断率约25-40%
- 肿瘤精准治疗：基于ctDNA的液体活检、MRD（微小残留病灶）监测、个性化疫苗
- 药物基因组学：基于CYP2D6、HLA-B57:01等基因型指导用药，避免严重不良反应
- 新生儿筛查：扩展性基因组筛查（NICUSeq）在重症新生儿中快速诊断遗传病
- 预防医学：基于多基因风险评分（PRS）的疾病风险预测和早期干预
（6）合成基因组学*
- 从"读取"到"编写"：基因组合成成本从2003年的$4/bp降至2020年的$0.001/bp以下
- 最小基因组：Mycoplasma laboratorium（2010年，首个合成细胞）
- 酵母基因组合成计划（Sc2.0）：重新设计并合成16条酵母染色体
- 人类基因组的合成规划（GP-write）正在讨论中
- 应用：设计新型生物制造底盘、创造抗病毒/抗逆作物、开发活体疗法

中国在基因组学领域已从"跟跑"进入"并跑"甚至部分"领跑"阶段：
测序技术自主化：
- 华大智造（MGI）的DNBSEQ技术成为全球第二大测序平台
- 真迈生物、齐碳科技等在纳米孔测序领域取得突破
大科学计划：
- 中国十万人基因组计划
- 百万微生态计划
- 万种鸟类基因组计划（B10K中国部分）
- 地球生物基因组计划中国节点
数据基础设施建设：
- 国家基因组科学数据中心（NGDC）
- 中国国家生物信息中心（CNCB）
- 基因组数据安全存储和共享的国家标准

6.1.1 基因组学的概念

基因组学研究的范畴可以从三个层次理解：
（1）结构基因组学（Structural Genomics）
结构基因组学关注基因组的物理结构，包括DNA序列的测定、基因组图谱的构建、基因的定位和结构分析等。其核心目标是获得基因组的完整序列，并理解基因在基因组中的排列方式。
（2）功能基因组学（Functional Genomics）
功能基因组学研究基因组中基因的功能、基因之间的相互作用、基因表达调控等。与结构基因组学回答"是什么"不同，功能基因组学回答"做什么"和"怎么做"。主要技术包括RNA-Seq、ChIP-Seq、ATAC-Seq等。
（3）比较基因组学（Comparative Genomics）
比较基因组学通过比较不同物种或不同个体基因组的异同，揭示基因组的进化规律、基因的起源和功能分化、物种之间的亲缘关系等。

基因组学中所研究的问题主要分为三类：
第一类：如何获得基因组序列
这是基因组学最基础的问题。随着DNA测序技术的发展，从一代Sanger测序到二代高通量测序（NGS），再到三代单分子长读长测序，获得基因组序列的成本急剧下降、速度大幅提升。一个物种的基因组项目通常包括：基因组大小估计、测序策略设计、序列组装、质量评估等步骤。
第二类：如何解读/解码基因组
获得序列后，需要解读这本"生命之书"。这包括：基因预测与结构注释（识别编码区、内含子、外显子、调控元件等）、功能注释（确定基因的功能、参与的通路）、比较分析（与近缘物种比较，理解进化关系）、变异分析（识别SNP、Indel、SV等）。
第三类：如何重写/编写新的基因组
这是合成生物学的前沿方向。通过基因编辑技术（如CRISPR-Cas9）对现有基因组进行精确修饰，甚至从头设计和合成全新的基因组。2010年，Craig Venter团队首次合成了支原体的人工基因组，标志着人类开始具备"编写生命"的能力。

人类基因组计划采用了层次测序策略：通过构建DNA大片段克隆文库、重叠克隆群和物理图谱，再进行逐个克隆鸟枪法测序。中国科学家承担了1%的测序任务（3号染色体端部约3000万个碱基对），成功跟进世界先进测序技术，推动了国内测序仪研发和生物信息学软件开发。
HGP的完成不仅提供了人类基因组的参考序列，更重要的是：
1. 技术推动：测序技术从Sanger法发展到高通量测序
2. 数据共享：建立了基因组数据开放共享的文化
3. 疾病研究：为疾病基因定位、药物基因组学奠定基础
4. 进化研究：为理解人类起源和进化提供基础数据

Sanger测序的核心是利用ddNTP终止DNA链延伸：
1. 反应体系：模板DNA、引物、DNA聚合酶、四种dNTP（其中一种带有放射性或荧光标记）和少量ddNTP
2. 链延伸与终止：DNA聚合酶沿模板合成新链，遇到dNTP时正常延伸，遇到ddNTP时因缺少3'-OH而终止
3. 四个独立反应：分别加入ddATP、ddCTP、ddGTP、ddTTP，产生四组不同长度的片段
4. 电泳分离：通过聚丙烯酰胺凝胶电泳（PAGE）或毛细管电泳按大小分离片段
5. 读取序列：从短到长读取末端碱基，获得互补链序列
自动化Sanger测序使用四种不同荧光标记的ddNTP，四个反应在同一管中进行，通过激光检测荧光颜色确定碱基类型。
局限性：通量低（每次反应只能读取一条序列），成本高，不适合大规模基因组项目。

Illumina测序流程包括四个关键步骤：
（1）文库制备
- 将基因组DNA随机打断成200-500 bp片段
- 末端修复、加A尾、连接测序接头（adapter）
- PCR扩增文库（可选）
（2）Flow Cell杂交与桥式PCR
- Flow Cell表面固定有两种与接头互补的引物
- 单链DNA文库与引物杂交
- 通过桥式PCR扩增：DNA片段弯曲成桥状，以另一端为模板进行扩增
- 经过多轮扩增，每个原始分子形成一个约1000-2000个克隆的"簇"（cluster）
- 变性后得到单链簇，作为测序模板
（3）边合成边测序（SBS）
- 加入DNA聚合酶和四种带有不同荧光标记的dNTP（每种dNTP的3'-OH被可切除的化学基团阻断）
- 每个循环：
a. 四种dNTP竞争掺入，每次只有一个碱基被加入
b. 激光扫描，根据荧光颜色确定掺入的碱基
c. 切除荧光基团和3'阻断基团，使链可以继续延伸
- 重复上述循环，逐碱基读取序列
（4）双端测序（Paired-End Sequencing）
- 从簇的两端分别测序，获得两段各约150 bp的序列
- 两段序列之间距离已知（约200-500 bp）
- 双端信息对于基因组组装和结构变异检测至关重要
数据格式与质量控制
测序数据以FASTQ格式存储，每条记录包含四行：序列标识符、碱基序列、"+"分隔符、碱基质量值（Phred分数）。Phred质量值Q与碱基错误率ε的关系为：Q = -10log₁₀(ε)。Q30表示碱基错误率约为0.1%，即准确性为99.9%。

（1）PacBio SMRT测序
PacBio测序的核心是零模波导孔（ZMW）技术：
- 在芯片上蚀刻数十万个直径约70 nm、深度约100 nm的纳米孔
- 每个孔底部固定一个DNA聚合酶分子
- 待测DNA与引物结合，被聚合酶捕获
- 加入四种带有不同荧光标记的dNTP（标记在磷酸基团上，而非碱基上）
- 当dNTP被掺入新链时，荧光标记被释放到溶液中，产生闪光信号
- 记录每个ZMW中的闪光序列和时间，即可获得DNA序列
PacBio Sequel II可产生约15-20 kb的长读长，最新的Revio系统进一步提高了通量。通过环化一致性测序（CCS，Circular Consensus Sequencing），同一个分子多次循环测序，可获得准确性>99.9%的HiFi reads。
（2）Nanopore测序
Nanopore测序的核心原理：
- 蛋白质纳米孔嵌入脂质双分子层膜中
- 在膜两侧施加电压，产生离子电流
- DNA/RNA分子在动力蛋白（helicase）驱动下穿过纳米孔
- 不同碱基（或k-mer）穿过纳米孔时引起特征性的电流变化
- 通过深度学习模型（如Bonito、Dorado）将电流信号解码为碱基序列
Nanopore的优势：
- 超长读长：目前最长记录超过4 Mb
- 便携性：MinION设备仅U盘大小，可在野外、太空等极端环境使用
- 直接测序RNA：无需逆转录，可检测RNA修饰
- 实时测序：数据边测边出，可即时分析
局限性：原始读取准确性较低（约92-97%），但R10.4.1 flow cell配合Dorado basecaller的准确性已大幅提升。

由于二代测序准确性高但读长短，三代测序读长长但原始准确性低，混合测序（Hybrid-Seq）成为当前基因组项目的标准配置：
- 用三代长读长进行基因组骨架组装
- 用二代短读长对组装结果进行纠错和填补缺口
- Hi-C数据辅助染色体级别组装

传统的基因组组装在高度重复序列区域（如着丝粒、端粒、rDNA阵列）存在大量缺口，主要原因是：
- 二代短读长无法跨越长重复单元
- 重复序列使组装算法难以确定正确的路径
T2T联盟实现人类完整基因组的关键技术组合：
（1）PacBio HiFi长读长
- 读长约15-20 kb，准确性>99.9%
- 足以跨越大多数重复单元（如着丝粒卫星重复单元约171 bp）
- 高准确性使得区分高度相似的重复序列成为可能
（2）Oxford Nanopore超长读长
- 读长可达100 kb至数Mb
- 直接跨越最大的重复阵列（如5S rDNA阵列可长达数Mb）
- 提供长距离的序列连续性信息
（3）Hi-C辅助scaffold连接
- 将HiFi组装产生的contig提升到染色体级别
- 确定contig的顺序、方向和间隙大小
（4）手动审校与验证
- 对复杂区域（如着丝粒）进行手动检查
- 使用荧光原位杂交（FISH）和光学图谱进行验证

Hi-C实验流程：
1. 甲醛交联：将空间上接近的DNA片段共价交联
2. 酶切：用限制性内切酶（如HindIII、DpnII）切割DNA
3. 末端修复和生物素标记：补平粘性末端并掺入生物素标记的dNTP
4. 连接：稀释后连接，使空间上接近的片段连接在一起
5. 去交联和测序：打断DNA，富集生物素标记的嵌合片段，进行双端测序
Hi-C数据分析用于基因组组装：
- 构建接触矩阵（contact matrix）：染色体各区域之间的相互作用频率
- 同一染色体上的区域相互作用频率高（呈现对角线富集）
- 相邻区域相互作用频率高于远端区域
- 使用算法（如3D-DNA、SALSA、ALLHiC）将contig聚类到染色体组、确定顺序和方向

泛基因组分为两部分：
- 核心基因组（Core Genome）：物种所有个体共有的基因，通常与基本生命功能相关
- 可变基因组/附属基因组（Accessory/Dispensable Genome）：仅在部分个体中存在的基因，通常与环境适应、致病性、代谢多样性等相关
泛基因组构建方法：
- 迭代式：将新个体的基因组与已有泛基因组比较，添加新序列
- 图基因组（Graph Genome）：将所有序列变异表示为图结构，而非线性序列
泛基因组的应用：
- 捕获参考基因组中缺失的遗传多样性
- 改进变异检测（减少比对偏差）
- 识别与重要性状相关的基因（如作物抗病基因）

一个典型的T2T基因组项目测序策略：
| 数据类型 | 测序深度 | 用途 |
|----------|---------|------|
| PacBio HiFi | 30-50× | 主要组装，产生高质量contig |
| Nanopore Ultra-long | 10-20× | 跨越超长重复区域，解决复杂区域 |
| Illumina PE150 | 50-100× | 纠错、质量验证 |
| Hi-C | 60-100× | 染色体级别scaffold |
| 光学图谱 | 5-10× | 验证组装结构 |

（1）线性序列展示
- 以染色体为坐标轴，从左到右展示DNA序列
- 支持不同尺度的缩放（从单碱基到全染色体）
- 显示正反链的序列和注释
（2）多轨道（Track）整合
- 每个track代表一类数据（基因注释、测序reads、变异、保守性等）
- 用户可自定义加载、排列和配置track
- 支持track的叠加和比较
（3）交互式查询
- 通过基因名、染色体坐标、rsID等定位基因组区域
- 支持序列搜索（如查找特定基序）
- 支持不同物种间的基因组比对（synteny view）
（4）数据导出
- 导出特定区域的序列（FASTA）
- 导出track数据（BED、GTF等）
- 生成高质量图片用于发表

（1）UCSC Genome Browser
- 最经典的基因组浏览器之一
- 整合了大量注释数据和工具
- 支持自定义track加载（通过URL或本地文件）
- 强大的Table Browser用于数据查询和下载
（2）Ensembl Genome Browser
- 欧洲生物信息研究所（EBI）维护
- 注释质量高，更新及时
- 强大的比较基因组学工具（如Ensembl Compara）
- 提供API和BioMart数据查询工具
（3）IGV（Integrative Genomics Viewer）
- 本地运行的桌面软件
- 支持大规模数据（如BAM、VCF、Hi-C矩阵）
- 特别适合查看测序reads比对情况和变异位点
- 支持基因组重排和结构变异的可视化
（4）JBrowse/GBrowse
- 可部署在本地服务器的网页浏览器
- 适合构建物种特异性的基因组数据库
- 支持插件扩展

BED格式
BED文件至少包含3列，最多12列：
- chrom：染色体名称
- chromStart：起始位置（0-based）
- chromEnd：终止位置（1-based，不包含）
- name：名称（可选）
- score：分数（0-1000，可选）
- strand：方向（+/-/.）
- thickStart：编码区起始（可选，用于显示）
- thickEnd：编码区终止（可选）
- itemRgb：颜色（可选）
- blockCount：外显子数量（BED12）
- blockSizes：外显子大小列表（逗号分隔）
- blockStarts：外显子相对起始位置列表
WIG格式
variableStep格式示例：

variableStep chrom=chr1 span=50
1001 12.5
1051 13.2
1101 11.8

表示chr1上1001-1050位置的值为12.5，1051-1100位置的值为13.2。
GTF格式
GTF是GFF2的变体，每行9列：
1. seqname：染色体或scaffold
2. source：注释来源
3. feature：特征类型（gene, transcript, exon, CDS等）
4. start：起始位置（1-based）
5. end：终止位置（1-based，包含）
6. score：分数（通常为".")
7. strand：方向（+/-/.)
8. frame：阅读框（0/1/2/.）
9. attributes：属性（gene_id, transcript_id等）

现代基因组研究往往涉及多组学数据的整合：
- 基因组层：基因结构、变异位点
- 转录组层：RNA-Seq信号、剪接异构体
- 表观组层：ChIP-Seq峰（组蛋白修饰、转录因子结合）、DNA甲基化、ATAC-Seq开放染色质
- 三维结构层：Hi-C相互作用矩阵、TAD边界
- 进化层：跨物种保守性分数（phyloP, phastCons）
整合可视化可以帮助发现跨层面的关联，例如：
- 某个增强子区域（H3K27ac峰+ATAC-Seq开放区）与远处基因的启动子有Hi-C相互作用
- 一个SNP位于保守非编码区，且与附近基因的表达量显著相关

6.2.1 基因组的组装

（1）重复序列
基因组中普遍存在重复序列，包括串联重复（如微卫星）和散在重复（如转座子、segmental duplication）。重复序列导致组装歧义——相同的reads可能来源于基因组的不同位置，组装算法无法确定其真实来源。
（2）测序错误
Illumina测序错误率约0.1-1%（Q30-Q20），Nanopore原始错误率可达5-15%。测序错误会干扰k-mer的正确连接，在DBG中产生"气泡"（bubble）和"死胡同"（tip）。
（3）覆盖度不均
基因组不同区域的测序深度差异很大。高GC区域、AT富集区域、重复序列区域的覆盖度往往偏低，导致这些区域组装困难。
（4）杂合度
二倍体生物的两条同源染色体存在差异。组装软件需要区分同源序列并可能分别组装（phased assembly），否则杂合位点会导致组装碎片化。

DBG算法是二代短读长组装的核心方法，其工作流程如下：
（1）k-mer分解
将所有reads拆分为长度为k的短片段（k-mer）。例如，序列"ATGCGGAT"在k=5时产生：ATGCG、TGCGG、GCGGA、CGGAT。通常k取奇数（如21、31、51、71），且不超过read长度的2/3。
（2）构建DBG
- 每个k-mer作为图中的一个顶点（node）
- 若两个k-mer之间有k-1个碱基重叠，则在它们之间建立一条有向边
- 例如，ATGCG和TGCGG共享TGCG（4个碱基），建立从ATGCG到TGCGG的边
- 将不同reads产生的相同k-mer合并为一个顶点，记录覆盖度
（3）图简化
测序错误会在DBG中产生两类错误结构：
- Tips（死胡同）：由于测序错误产生的低频k-mer连接到主路径后又很快终止。处理：去除长度和覆盖度低于阈值的tips。
- Bubbles（气泡）：由于杂合位点或测序错误，主路径在某处分叉为两条路径，在下游又汇合。处理：保留覆盖度较高的路径，或在杂合组装中保留两条路径。
（4）contig提取
在简化后的DBG中，一条连续的路径（从没有入度的顶点开始到没有出度的顶点结束）对应一个contig。
k-mer大小的选择：
- k值较小：对重复序列的分辨率低（不同重复单元可能共享k-mer），但覆盖度高（每个read产生更多k-mer）
- k值较大：对重复序列的分辨率高，但覆盖度低（要求更大的测序深度才能保证每个k-mer被覆盖）
- 多k-mer策略（如SPAdes）：从较小的k开始组装，逐步增大k值，利用小k的连续性和大k的准确性优势

OLC算法更适合长读长数据（PacBio/Nanopore），其工作流程：
（1）Overlap（重叠检测）
计算所有reads两两之间的重叠关系。对于N条reads，需要O(N²)次比对，计算量巨大。对于长读长数据，通常使用minimap2等快速比对工具进行all-vs-all比对。
（2）Layout（布局构建）
根据重叠关系构建布局图（string graph）：
- 每个read是图中的一个节点
- 如果两个reads有显著重叠，则在它们之间建立边
- 边的方向和长度由重叠的位置和长度确定
- 简化图结构（去除 transitive edges，压缩 unambiguous paths）
（3）Consensus（一致性序列生成）
对布局图中的每个unitig（unambiguous path），将对应的reads进行多序列比对，生成一致性序列作为contig。

利用二代和三代数据的互补性：
（1）二代纠错三代
- 用高准确性的二代reads对三代reads进行纠错
- 代表工具：PBcR、LoRDEC、Jabba
- 纠错后的三代reads准确性可提升至>99%
（2）三代搭建骨架，二代填充
- 用三代长读长构建初步的contig骨架
- 用二代reads对三代contig进行 polish（填补gap、纠正碱基错误）
- 代表工具：Pilon、Racon
（3）graph-based混合组装
- 将二代和三代数据同时输入组装算法
- 在组装图中同时利用短读长的准确性和长读长的连续性
- 代表工具：HybridSPAdes、MaSuRCA

连续性指标：
- N50/L50：衡量contig/scaffold长度的分布
- 最大contig长度：最长的连续序列
- gap数量：scaffold中N碱基的数量和总长度
完整性指标：
- BUSCO（Benchmarking Universal Single-Copy Orthologs）：检测基因组中保守单拷贝基因的存在情况。高BUSCO完整性（如>90%）表示基因组组装较完整。
- CEGMA（Core Eukaryotic Genes Mapping Approach）：类似BUSCO，基于核心真核基因集。
准确性指标：
- QUAST：计算N50、错误组装（misassembly）数量、基因组覆盖度等
- Pilon polish后的二代比对一致性：检测碱基错误
- k-mer spectrum分析：与原始测序数据的k-mer谱比较，评估完整性

原核生物基因特征：
- 基因密度高（大肠杆菌基因组约85%为编码区）
- 无内含子（少数古菌有内含子）
- 无5'端帽子和3'端polyA尾
n- 操纵子结构：多个基因共用一个启动子，转录为多顺反子mRNA
- 起始密码子通常为ATG（少数GTG、TTG）
- 启动子元件：-10区（Pribnow框：TATAAT）和-35区（TTGACA）
真核生物基因特征：
- 基因密度低（人类基因组仅约1-2%为编码区）
- 含内含子（人类基因平均含8个外显子、7个内含子）
- 5'端有帽子结构，3'端有polyA尾
- 单顺反子mRNA（一个转录本对应一个蛋白质）
- 复杂的剪接机制：内含子5'端为GT（供体位点），3'端为AG（受体位点）
- 核心启动子元件：TATA框（约-30区）、CAAT框（约-80区）
- 存在大量的可变剪接（人类约95%的多外显子基因发生可变剪接）

从头预测的核心是利用基因组的统计特征和信号特征：
（1）开放阅读框（ORF）检测
- 寻找起始于ATG、终止于TAA/TAG/TGA的连续密码子序列
- 计算ORF长度：真核生物外显子通常>300 bp（100个氨基酸），短ORF可能是假阳性
- 计算密码子使用频率（Codon Usage Frequency, CUF）：真实基因通常偏好使用特定密码子
（2）隐马尔可夫模型（HMM）
HMM是基因预测中最常用的统计框架：
- 隐藏状态：基因间区、外显子（起始/中间/终止）、内含子（起始/中间/终止）、UTR等
- 观测符号：A、C、G、T
- 状态转移概率：从一个状态转移到另一个状态的概率（如外显子→内含子的概率）
- 发射概率：在某个状态下观察到特定碱基的概率
- Viterbi算法：寻找最可能的隐藏状态序列（即基因结构）
经典HMM基因预测软件：
- Genscan：针对人类和脊椎动物，基于广义HMM
- Augustus：支持真核和原核生物，可整合外部证据
- Glimmer/GlimmerHMM：Glimmer用于原核，GlimmerHMM用于真核
- GeneMark/GeneMark-ES：自训练HMM，适合新物种
- SNAP：可训练参数，适合非模式生物
（3）支持向量机（SVM）和深度学习方法
近年来，机器学习方法也被用于基因预测：
- 将序列特征（k-mer频率、GC含量、密码子偏好等）编码为特征向量
- 用已知基因训练分类器（SVM、随机森林、神经网络）
- 对新序列进行预测

（1）蛋白质同源比对
- 将已知物种的蛋白质序列与目标基因组进行比对（如tblastn）
- 根据蛋白质比对结果推断基因的编码区
- 考虑剪接位点信号，将基因组DNA与外显子对应
- 代表工具：GeneWise、Exonerate（protein2genome模式）
（2）cDNA/EST比对
- 将已知的cDNA或EST序列与基因组比对（如BLAT、GMAP）
- 直接确定外显子-内含子边界
- EST数据质量参差不齐，需要过滤和聚类
（3）跨物种比较基因组学
- 利用进化保守性：编码区通常比非编码区更保守
- 利用密码子替换模式：同义替换率（Ks）和非同义替换率（Ka）
- 代表方法：Twinscan（整合人类-小鼠比较信息）、CONTRAST

（1）RNA-Seq数据比对与组装
- 将RNA-Seq reads比对到基因组（如HISAT2、STAR）
- 使用StringTie、Cufflinks等工具组装转录本
- 获得GTF格式的转录本结构注释
（2）全长转录本测序
- PacBio Iso-Seq或Nanopore direct RNA测序获得全长cDNA
- 无需组装，直接获得完整的外显子-内含子结构
- 可发现新的转录本异构体
（3）整合策略（Evidence Modeler, EVM）
- 将从不同来源（从头预测、同源比对、转录本组装）获得的基因模型进行整合
- 根据证据的可靠性和一致性进行加权评分
- 输出一致性最高的基因模型集合
- 代表工具：EVM（Evidence Modeler）、PASA（Program to Assemble Spliced Alignments）、BRAKER

非编码RNA（ncRNA）是不编码蛋白质的RNA分子，包括：
- miRNA（微小RNA）：~22 nt，调控基因表达
- lncRNA（长链非编码RNA）：>200 nt，多种调控功能
- circRNA（环状RNA）：共价闭合环，miRNA海绵等
- tRNA、rRNA、snRNA、snoRNA等
预测方法：
- miRNA：基于发夹结构、种子区域保守性、热力学稳定性（RNAfold、miRDeep2）
- lncRNA：基于开放阅读框长度、编码潜能（CPAT、CPC2）、序列保守性
- circRNA：基于反向剪接位点（backsplice junction）的reads支持（CIRI、find_circ）

GO数据库的核心是一个有向无环图（DAG）结构，其中：
- 每个节点（GO term）代表一个功能概念
- 边代表"is_a"（是一种）或"part_of"（是...的一部分）关系
- 子节点比父节点更具体（如"蛋白质激酶活性" is_a "转移酶活性"）
GO注释的常用方法：
（1）基于序列相似性的GO注释
- 用BLAST将查询蛋白与GO注释数据库（如UniProt-GOA）比对
- 将相似蛋白的GO term转移给查询蛋白
- 常用工具：Blast2GO、eggNOG-mapper
（2）基于结构域的GO注释
- 通过InterProScan识别蛋白质结构域
- 利用结构域与GO term的映射关系进行注释
- InterPro数据库提供了Pfam、SUPERFAMILY等结构与GO的映射
（3）基于机器学习的方法
- 使用蛋白序列特征（氨基酸组成、物理化学性质、亚细胞定位信号等）训练分类器
- 预测GO term
- 代表方法：DeepGO（基于深度学习）

KEGG数据库的核心是通路图（Pathway Map），包括：
- 代谢通路（Metabolism）：碳水化合物代谢、能量代谢、脂质代谢、核苷酸代谢、氨基酸代谢等
- 遗传信息处理（Genetic Information Processing）：转录、翻译、折叠/分选/降解、复制与修复
- 环境信息处理（Environmental Information Processing）：膜转运、信号转导、配体-受体相互作用
- 细胞过程（Cellular Processes）：细胞生长与死亡、细胞群落、运输与分解代谢
- 生物体系统（Organismal Systems）：免疫系统、神经系统、内分泌系统等
- 人类疾病（Human Diseases）：癌症、神经退行性疾病、心血管疾病等
KEGG注释方法：
- 序列比对法：用BLAST/KAAS（KEGG Automatic Annotation Server）将基因与KEGG GENES数据库比对
- KO（KEGG Orthology）注释：将基因映射到KO编号（功能单位），再根据KO编号映射到通路
- 工具：KAAS、KOFAM、eggNOG-mapper

InterPro整合了多个蛋白质特征数据库：
- Pfam：基于HMM的蛋白质家族和结构域数据库
- PRINTS：基于指纹（保守基序组）的数据库
- PANTHER：基于系统发育树的蛋白质家族分类
- ProSite：基于模式（pattern）和谱（profile）的功能位点数据库
- SUPERFAMILY：基于SCOP结构分类的远缘同源识别
- SMART：信号肽、跨膜区、结构域识别
InterProScan工作流程：
1. 用多种方法（HMM、模式匹配、profile等）扫描蛋白质序列
2. 整合不同数据库的预测结果
3. 为每个预测的domain/feature提供GO term映射
4. 输出GFF3或TSV格式的注释结果

（1）COG/KOG/EggNOG
- COG（Clusters of Orthologous Groups）：原核生物蛋白质直系同源簇
- KOG：真核生物版本
- EggNOG：扩展版本，覆盖更多物种，提供GO、KEGG、CAZy等注释
- 将蛋白质分为功能类别（如J: Translation, K: Transcription, S: Unknown）
（2）NR/UniProt/Swiss-Prot
- NR（Non-Redundant）：NCBI的非冗余蛋白数据库，用于序列相似性搜索
- UniProt：最全面的蛋白质数据库，包含Swiss-Prot（人工审编）和TrEMBL（自动注释）
- 注释流程：BLASTp比对NR/UniProt → 取top hit → 继承功能描述
（3）专用数据库
- CAZy：糖类活性酶分类（糖苷水解酶、糖基转移酶等）
- TCDB：膜转运蛋白分类
- VFDB：细菌毒力因子
- CARD：抗生素耐药基因
- antiSMASH：生物合成基因簇（次级代谢产物）

现代基因组项目通常使用自动化流程进行功能注释：

预测蛋白质序列
    ↓
1. InterProScan（结构域识别 + GO注释）
    ↓
2. KAAS/KOFAM（KEGG通路注释）
    ↓
3. eggNOG-mapper（COG/GO/KEGG综合注释）
    ↓
4. BLASTp against NR/Swiss-Prot（功能描述）
    ↓
5. 整合所有注释结果

6.3.1 基本原理

基因组变异按尺度可分为：
点突变尺度：
- SNP（Single Nucleotide Polymorphism）：单个碱基的替换，是最常见的变异类型
- Indel（Insertion/Deletion）：小片段（通常<50 bp）的插入或缺失
结构变异尺度：
- 缺失（Deletion）：>50 bp的序列缺失
- 插入（Insertion）：>50 bp的序列插入
- 倒位（Inversion）：序列方向的翻转
- 易位（Translocation）：序列从一个染色体位置移动到另一个位置
- 重复（Duplication）：序列拷贝数的增加
- CNV（Copy Number Variation）：>1 kb的拷贝数变化

原始测序数据 (FASTQ)
    ↓
质量控制 (FastQC, Trimmomatic)
    ↓
比对到参考基因组 (BWA, Bowtie2)
    ↓
比对后处理 (排序, 标记重复, 碱基质量校正)
    ↓
变异检测 (GATK, Samtools, FreeBayes)
    ↓
变异过滤和质量控制
    ↓
变异注释 (功能影响预测)
    ↓
最终变异集合 (VCF格式)

（1）排序（Sorting）
比对后的SAM/BAM文件需要按染色体坐标排序，以便后续处理。
（2）标记PCR重复（Mark Duplicates）
PCR扩增会引入系统性偏差——来自同一模板分子的多个reads并非独立采样。使用Picard MarkDuplicates等工具识别并标记PCR重复reads，避免重复reads对变异检测的干扰。
（3）Indel区域重比对（Indel Realignment）
比对软件在Indel附近可能产生比对偏差（如将Indel附近的SNP错误比对）。GATK的IndelRealigner通过识别Indel附近不一致的比对，进行局部重比对，提高Indel和周围SNP的检测准确性。
（4）碱基质量校正（Base Quality Score Recalibration, BQSR）
测序仪给出的原始碱基质量值可能存在系统性偏差（如某些测序循环质量系统性地偏低）。BQSR通过将碱基质量与已知变异数据库（如dbSNP）进行比较，建立校正模型，重新校准碱基质量值。

变异检测的核心是判断某个位置的碱基差异是真实的遗传变异还是测序/比对错误。这本质上是一个统计假设检验问题。
零假设（H₀）：该位置没有变异，所有观察到的差异都是测序错误
备择假设（H₁）：该位置存在变异
测序错误通常被建模为二项分布或贝叶斯模型：
- 设测序错误率为ε（通常ε≈0.001，对应Q30）
- 在覆盖度为D的位置，观察到k个非参考碱基
- 如果H₀成立，k ~ Binomial(D, ε)
- 计算P(k | H₀)，如果该P值小于阈值（如0.05），则拒绝H₀，判定为变异
基因型似然模型：
更精确的方法是对三种可能的基因型（参考纯合RR、杂合RA、变异纯合AA）分别计算似然值：
- L(RR) = ∏ P(read_i | RR)
- L(RA) = ∏ P(read_i | RA)
- L(AA) = ∏ P(read_i | AA)
选择似然值最大的基因型作为该位点的基因型。
GATK HaplotypeCaller使用更先进的单倍型模型，在局部区域（通常数百bp）同时考虑所有可能的变异，通过局部de novo组装生成候选单倍型，然后计算每个单倍型的似然值。

即使使用了严格的统计模型，变异检测结果中仍包含大量假阳性。需要通过多重过滤来提高准确性：
基于测序质量的过滤：
- QD（Quality by Depth）：质量值除以覆盖度，过低提示质量不可靠
- FS（Fisher Strand）：检验正负链reads是否均衡支持变异，不均衡提示比对偏差
- MQ（Mapping Quality）：比对质量值，过低提示比对不可靠
- ReadPosRankSum：检验变异是否主要出现在reads末端（末端碱基质量通常较低）
基于覆盖度的过滤：
- DP（Depth）：覆盖度过低（<5-10×）的变异不可靠
- 覆盖度过高（>平均深度的2-3倍）可能提示重复区域
基于群体频率的过滤：
- 使用已知变异数据库（如gnomAD）过滤常见变异
- 在罕见病研究中，过滤在人群中频率>1%的变异

SNP主要来源于DNA复制过程中的自发突变：
（1）脱氨基作用
- 胞嘧啶（C）自发脱氨基转变为尿嘧啶（U），DNA修复系统将U识别为T，导致C→T突变
- 这是转换突变的主要来源，解释了为什么转换/颠换比（Ti/Tv ratio）通常大于2
- 5-甲基胞嘧啶（5mC）的脱氨基率更高，导致CpG岛中的C→T突变频率显著升高
（2）复制错误
- DNA聚合酶在复制过程中偶尔掺入错误的碱基
- 尽管有校对功能，错误率仍约为10⁻⁹/碱基/复制
（3）氧化损伤
- 活性氧（ROS）可导致碱基氧化修饰，如8-oxo-guanine可导致G→T颠换

按功能位置分类：
| 位置 | 比例 | 功能影响 |
|------|------|---------|
| 基因间区 | ~90% | 通常无直接功能影响 |
| 内含子 | ~5% | 可能影响剪接调控 |
| 外显子同义 | ~1% | 不改变氨基酸（密码子简并性） |
| 外显子错义 | ~1.5% | 改变氨基酸序列 |
| 外显子无义 | ~0.1% | 引入提前终止密码子 |
| 调控区 | ~2% | 可能影响基因表达 |
| 剪接位点 | <0.1% | 可能导致剪接异常 |
按群体频率分类：
- 常见变异（Common Variant）：群体频率>5%，通常与复杂疾病的易感风险相关
- 低频变异（Low-frequency Variant）：群体频率1-5%
- 稀有变异（Rare Variant）：群体频率<1%，可能与孟德尔遗传病相关

（1）基于pileup的方法（Samtools/bcftools）
- 统计每个基因组位置上比对reads的碱基组成
- 使用贝叶斯模型计算每个可能基因型的后验概率
- 优点：速度快，计算资源需求低
- 缺点：不处理Indel附近的比对偏差，对Indel不敏感
（2）基于单倍型的方法（GATK HaplotypeCaller）
- 在局部窗口（数百bp）内进行de novo组装，生成候选单倍型
- 将reads与候选单倍型进行比对，计算每个单倍型的似然值
- 选择最可能的单倍型组合作为基因型
- 优点：对Indel和复杂区域的SNP检测更准确
- 缺点：计算量大，耗时较长
（3）基于深度学习的方法（DeepVariant）
- 将每个候选变异位点转化为"图像"（read pileup的视觉表示）
- 使用卷积神经网络（CNN）进行变异识别
- 优点：准确性极高，在多个基准测试中优于传统方法
- 缺点：计算资源需求高，训练数据依赖

QD（Quality by Depth）：
QD = QUAL / DP
- QD过低（<2）提示变异可能被低质量reads或重复区域支持
- QD过高（>40）可能提示比对到重复区域
MQ（Mapping Quality）：
- 反映reads在该区域的比对可靠性
- MQ < 40提示比对可能存在问题
FS（Fisher Strand Bias）：
- 检验变异是否被正负链均衡支持
- FS > 60提示严重的链偏倚，变异可能不可靠
SOR（Strand Odds Ratio）：
- FS的替代指标，对覆盖度不均更稳健
- SOR > 3提示链偏倚
MQRankSum：
- 检验支持变异的reads与支持参考的reads的比对质量是否存在显著差异
- 负值过大提示支持变异的reads比对质量较差
ReadPosRankSum：
- 检验变异是否主要出现在reads末端
- 末端碱基质量通常较低，可能导致假阳性

转换/颠换比率（Ti/Tv ratio）是评估SNP集合质量的重要指标：
- 全基因组范围的期望Ti/Tv：约2.0-2.2
- 外显子区域的期望Ti/Tv：约2.8-3.0（因为同义SNP只发生在特定密码子位置）
如果检测到的SNP集合Ti/Tv显著低于期望值，可能提示：
- 存在大量假阳性（特别是颠换类型的假阳性较多）
- 过滤条件不够严格
- 样本存在污染

（1）DNA聚合酶滑移（Replication Slippage/Strand Slippage）
- 在重复序列区域（尤其是微卫星），模板链和新生链可能发生相对滑动
- 导致重复单元的插入或缺失
- 这是Indel最主要的产生机制，解释了为什么Indel在重复序列区域富集
（2）不等交换（Unequal Crossing Over）
- 减数分裂过程中的同源重组发生错位
- 导致一个染色体上重复序列增加，另一个上减少
（3）DNA损伤修复错误
- 双链断裂修复过程中可能发生小片段的插入或缺失

Indel检测比SNP检测困难的主要原因：
（1）比对困难
- 含有Indel的reads与参考基因组比对时，需要在Indel位置引入"空位"（gap）
- 比对算法对空位有罚分，倾向于将Indel区域的碱基错配解释为SNP
- 例如，参考序列"ATCGATCG"，reads为"ATC---GATCG"（3bp缺失），比对软件可能错误地比对为"ATCGATCG"（多个错配）
（2）Indel周围的碱基质量下降
- Illumina测序在Indel附近容易产生碱基识别错误
- 需要专门的Indel重比对步骤来纠正
（3）Indel大小的不确定性
- 对于较大的Indel（>10 bp），单个read可能无法完全跨越，导致Indel边界不确定
- 需要多个reads的信息来确定Indel的精确边界

（1）基于局部重比对的Indel检测（GATK HaplotypeCaller）
- HaplotypeCaller在局部窗口内生成候选单倍型（包含可能的SNP和Indel组合）
- 将reads与所有候选单倍型比对，选择最优解释
- 优点：同时考虑SNP和Indel，避免"SNP vs Indel"的比对歧义
（2）基于split-read的Indel检测（Pindel）
- 识别"split reads"——一条read的两部分分别比对到Indel两侧
- 通过分析未比对上的read片段来推断Indel序列
- 优点：对较大Indel（>20 bp）的边界定位准确
- 缺点：只能检测被reads跨越的Indel
（3）基于de novo局部组装的Indel检测（Platypus, Dindel）
- 在候选Indel区域提取相关reads
- 进行局部de novo组装
- 将组装结果与参考基因组比对，识别Indel

Indel的功能影响通常比SNP更严重：
| Indel类型 | 功能影响 | 示例 |
|----------|---------|------|
| 编码区非3倍数Indel | 移码突变，通常导致功能丧失 | BRCA1移码突变导致乳腺癌风险 |
| 编码区3倍数Indel | 氨基酸插入/缺失，可能影响功能 | CFTR ΔF508导致囊性纤维化 |
| 剪接位点Indel | 影响pre-mRNA剪接 | 深度内含子Indel激活隐蔽剪接位点 |
| 调控区Indel | 影响转录因子结合 | 胰岛素基因启动子Indel影响表达 |
| UTR区Indel | 可能影响mRNA稳定性或翻译 | 5'UTR的uORF相关Indel |
经典病例：
- 囊性纤维化：CFTR基因最常见的突变是ΔF508（一个3bp缺失），导致苯丙氨酸缺失，约占所有囊性纤维化病例的70%
- 亨廷顿舞蹈症：HTT基因中CAG重复序列的扩增（可视为一种特殊Indel），当CAG重复数>36时导致疾病

（1）基于Read深度（Read Depth, RD）的方法
原理：测序深度与基因组拷贝数成正比。
- 缺失区域：reads深度显著低于平均值（约50%对于杂合缺失）
- 重复区域：reads深度显著高于平均值
代表软件：CNVnator、BIC-seq、CNV-seq、DELLY（结合RD和PE）
优点：
- 可以检测大尺度的CNV
- 计算相对简单
局限：
- 无法精确定位断点（只能确定区域）
- 无法检测平衡性SV（如倒位、易位不改变拷贝数）
- 受GC含量、比对偏好等因素影响
（2）基于Split Read（SR）的方法
原理：识别一条read的两部分分别比对到不同位置的"split reads"，从而精确定位断点。
代表软件：Pindel、DELLY、LUMPY、Manta
工作流程：
1. 提取未完全比对或部分比对的reads（soft-clipped reads）
2. 将未比对部分重新比对到基因组其他位置
3. 如果找到匹配，确定断点坐标
优点：
- 断点定位精确（可达bp级别）
- 可以检测小至50 bp的SV
局限：
- 只能检测被reads跨越的SV
- 需要足够长的reads（对Illumina短读长有限制）
（3）基于Paired-End（PE）的方法
原理：利用paired-end测序中insert size和方向的异常来检测SV。
正常配对特征：
- 两个read在同一染色体上
- 方向相对（read1正向，read2反向）
- 距离在预期insert size范围内
异常配对类型及对应的SV：
| 异常类型 | 提示的SV |
|---------|---------|
| 距离远大于预期 | 缺失（reads跨越了被删除的区域） |
| 距离远小于预期 | 插入（reads之间的区域被插入序列替代） |
| 方向相同 | 倒位 |
| 位于不同染色体 | 易位 |
代表软件：DELLY、LUMPY、BreakDancer、Manta
优点：
- 可以检测大尺度SV
- 对read长度要求不高
局限：
- 断点定位不精确（通常在数百bp范围内）
- 受insert size分布的影响
（4）基于de novo组装的方法
原理：对SV区域进行局部de novo组装，将组装结果与参考基因组比对来识别SV。
代表软件：TIGRA、MindTheGap、SMRT-SV（PacBio专用）、Sniffles（Nanopore专用）
优点：
- 可以检测最复杂的SV，包括串联重复扩增、倒位等
- 可以获得SV的完整序列
局限：
- 计算资源需求大
- 组装质量依赖于数据质量和覆盖度

现代SV检测工具通常整合多种信号：
DELLY：整合PE、SR和RD信号
LUMPY：整合PE和SR信号，使用概率框架
Manta：结合PE、SR和局部组装，速度快，准确性高
GRIDSS：基于基因组图（assembly graph）的SV检测

长读长测序（PacBio/Nanopore）对于SV检测具有革命性优势：
1. 跨越复杂SV：长读长可以直接跨越整个SV区域，获得完整的SV序列
2. 精确定位断点：无需依赖PE或SR推断，直接比对即可定位
3. 检测重复区域SV：短读长难以比对的重复区域，长读长可以跨越
专用工具：
- Sniffles：Nanopore SV检测
- PBSV：PacBio SV检测
- SVIM：同时支持PacBio和Nanopore

SV检测的假阳性率通常高于SNP/Indel，验证至关重要：
实验验证方法：
- PCR+Sanger测序：对断点区域进行PCR扩增和测序，是SV验证的金标准
- qPCR：验证拷贝数变化
- Southern blot：验证大尺度重复/缺失
- 光学图谱（Bionano）：验证大尺度SV
计算验证方法：
- 使用不同检测工具的结果取交集
- 检查SV区域的reads比对情况（IGV可视化）
- 与已知SV数据库（如gnomAD-SV、DGV）比对

6.4.1 微生物组概述

人体微生物组的规模：
- 总重量：约1-2 kg（主要在肠道）
- 细胞数量：约3×10¹³个（与人体自身细胞数量相当或略多）
- 基因数量：约300-500万个基因（是人类基因组基因数的150-250倍）
人体各部位的微生物组特征：
| 部位 | 主要门类 | 细菌密度 | 主要功能 |
|------|---------|---------|---------|
| 结肠 | Firmicutes, Bacteroidetes | 10¹¹-10¹²/g | 发酵膳食纤维、合成维生素 |
| 口腔 | Firmicutes, Proteobacteria | 10⁹/mL唾液 | 参与氮循环、免疫调节 |
| 皮肤 | Actinobacteria, Firmicutes | 10⁶-10⁷/cm² | 屏障保护、免疫训练 |
| 阴道 | Lactobacillus | 10⁸-10⁹/mL | 维持酸性环境、防御病原 |
| 鼻腔 | Firmicutes, Actinobacteria | 10⁶-10⁷/cm² | 免疫防御 |

传统微生物学估计，环境中可培养的微生物仅占1%左右（"Great Plate Count Anomaly"）。主要原因包括：
1. 营养需求未知：很多微生物需要特殊的生长因子或复杂的共生关系
2. 环境条件难以模拟：如极端pH、温度、压力、厌氧条件等
3. 生长缓慢：某些微生物代时可达数周甚至数月
4. 群体感应依赖：需要达到一定密度才能生长
5. 严格的共生关系：依赖其他微生物产生的代谢产物
宏基因组学通过直接提取环境DNA进行测序，彻底 bypass 了培养瓶颈。

16S rRNA（原核生物）和18S rRNA/ITS（真核微生物）是微生物分类鉴定的"条形码"。
16S rRNA的优势：
1. 普遍存在：所有细菌和古菌都含有16S rRNA基因
2. 功能保守：核糖体功能对生命至关重要，序列演化相对缓慢
3. 大小适中：约1,500 bp，既包含保守区（设计通用引物），又包含可变区（区分物种）
16S rRNA的可变区：
- V1-V9共9个可变区
- 常用测序区域：V3-V4（约465 bp，Illumina MiSeq 2×300可覆盖）
- 不同可变区的分辨能力不同：V1-V2适合属水平，V3-V4可区分到种/株水平
分类阈值：
- 16S rRNA序列相似性 >97%：通常认为是同一个种
- 16S rRNA序列相似性 >95%：通常认为是同一个属

原始测序数据 (FASTQ)
    ↓
质量控制 (FastQC, cutadapt)
    ↓
去接头、去低质量、去嵌合体 (DADA2/Deblur/QIIME2)
    ↓
生成ASV/OTU表
    ↓
物种注释 (SILVA/Greengenes/RDP数据库)
    ↓
多样性分析 (α, β)
    ↓
统计分析 (LEfSe, ANCOM, DESeq2)
    ↓
可视化 (phyloseq, ggplot2)

（1）DADA2去噪算法
DADA2（Divisive Amplicon Denoising Algorithm 2）通过以下步骤生成ASV：
1. 质控过滤：去除低质量和含N的reads
2. 去重复：将相同序列合并，记录丰度
3. 错误学习：从数据中估计测序错误模型
4. 去噪：将含错误的序列纠正回其真实序列
5. 去嵌合体：识别并去除PCR嵌合体
6. 合并双端reads
（2）OTU聚类 vs ASV
- OTU聚类以97%相似性为阈值，将序列分组。缺点：阈值是人为设定的；近缘种可能因高于阈值而无法区分；信息在聚类过程中丢失
- ASV保留了每个独特的真实序列，分辨率更高，且可重复性更好

原始测序数据 (FASTQ)
    ↓
质量控制 (FastQC, Trimmomatic)
    ↓
去宿主DNA (Bowtie2比对宿主参考基因组)
    ↓
宏基因组组装 (MEGAHIT, metaSPAdes)
    ↓
基因预测 (Prodigal, MetaGeneMark)
    ↓
基因 catalog构建 (CD-HIT聚类)
    ↓
功能注释 (eggNOG-mapper, KEGG)
    ↓
分箱 (MetaBAT2, MaxBin2, CONCOCT)
    ↓
MAG质量评估 (CheckM)
    ↓
物种分类和系统发育分析

分箱的目标是将来自同一基因组的contig聚类到一起。
（1）基于组成的分箱（Composition-based Binning）
- 原理：不同物种的基因组具有不同的碱基组成特征（如GC含量、四核苷酸频率）
- 方法：将contig表示为高维特征向量（如136维四核苷酸频率向量），使用聚类算法（如K-means、高斯混合模型）进行聚类
- 代表工具：MetaBAT2、MaxBin2
（2）基于覆盖度的分箱（Coverage-based Binning）
- 原理：来自同一物种的contig在不同样本中的覆盖度变化模式一致（因为它们来自同一个基因组）
- 方法：在多个样本中计算每个contig的覆盖度，构建覆盖度矩阵，进行聚类
- 代表工具：CONCOCT、MaxBin2
（3）混合分箱
- 同时使用组成和覆盖度信息进行分箱
- 代表工具：MetaBAT2（使用TNF和覆盖度）、VAMB（使用变分自编码器）
MAG质量评估：
- CheckM：基于保守单拷贝基因（SCG）评估MAG的完整性和污染度
- 高质量MAG：完整性>90%，污染度<5%
- 中等质量MAG：完整性≥50%，污染度<10%

α多样性指数：
| 指数 | 衡量 | 公式 | 特点 |
|------|------|------|------|
| 物种数（Richness） | 丰富度 | S = 物种数量 | 不考虑丰度 |
| Shannon | 丰富度+均匀度 | H' = -Σ(pi × ln pi) | 对丰富度敏感 |
| Simpson | 丰富度+均匀度 | D = 1 - Σ(pi²) | 对优势种敏感 |
| Chao1 | 估计总丰富度 | Chao1 = S_obs + n₁²/(2n₂) | 考虑稀有物种 |
| Pielou's J | 均匀度 | J' = H'/ln(S) | 纯度的衡量 |
β多样性距离度量：
| 距离 | 公式 | 特点 |
|------|------|------|
| Bray-Curtis | 1 - 2Σmin(xi, yi) / Σ(xi + yi) | 考虑丰度，最常用 |
| Jaccard | 1 - |A∩B| / |A∪B| | 仅考虑有无 |
| Unweighted UniFrac | 考虑系统发育距离 | 考虑进化关系 |
| Weighted UniFrac | 考虑丰度和系统发育 | 综合考量 |
β多样性可视化：
- PCoA（主坐标分析）
- NMDS（非度量多维尺度分析）
- UPGMA聚类树

识别在不同组间丰度显著差异的微生物或功能基因：
- LEfSe（Linear Discriminant Analysis Effect Size）：结合统计检验和LDA评分，筛选组间差异标志物
- ANCOM（Analysis of Composition of Microbiomes）：考虑微生物组数据的组成性特点
- DESeq2/edgeR：基于负二项分布的差异分析（适用于OTU/基因计数数据）
- ALDEx2：考虑组成性和稀疏性的差异分析

微生物组数据具有以下"大数据"特征：
1. 高维度：数千个OTU/ASV × 数百个样本
2. 稀疏性：很多物种只在少数样本中出现（零膨胀分布）
3. 组成性：相对丰度总和为100%，变量间存在约束
4. 高变异性：个体差异大，受饮食、环境等多种因素影响
5. 多层次：物种、功能、代谢物等多个层面的数据

（1）松散的整合（Loose Integration）
- 分别分析各组学数据
- 比较不同组学层面的结果，寻找一致性
- 例如：比较差异物种和差异代谢物，寻找相关性
（2）统计学整合（Statistical Integration）
- 使用多变量统计方法整合多个数据矩阵
- CCA（Canonical Correlation Analysis）：寻找两个数据集之间的最大相关
- sPLS（sparse Partial Least Squares）：寻找预测关系的同时进行变量选择
- MCIA（Multiple Co-Inertia Analysis）：整合多个组学数据集
- MOFA（Multi-Omics Factor Analysis）：分解共享和特异的变异源
（3）网络整合（Network Integration）
- 构建跨组学的关联网络
- 例如：物种-代谢物关联网络、基因-代谢物关联网络
- 识别网络中的关键节点（hub）
（4）机器学习整合
- 将多组学特征作为输入，训练预测模型
- 使用特征选择识别最重要的预测因子

（1）分类任务：疾病预测
- 输入：微生物组组成特征（物种丰度、功能通路丰度）
- 输出：疾病状态（如IBD vs 健康）
- 常用算法：
- 随机森林（Random Forest）：对高维稀疏数据表现好，可输出特征重要性
- 支持向量机（SVM）：适合小样本高维数据
- 深度学习（神经网络）：可学习复杂的非线性关系
案例：使用随机森林预测CRC（结直肠癌）
- 训练集：数百个粪便样本的16S/宏基因组数据
- 特征：物种丰度（或MetaCyc通路丰度）
- 模型：随机森林分类器
- 结果：AUC可达0.8-0.9，某些研究报道准确率>90%
- 关键标志物：Fusobacterium nucleatum、Peptostreptococcus anaerobius等
（2）回归任务：表型预测
- 输入：微生物组特征 + 宿主因素
- 输出：连续变量（如BMI、药物响应、血糖水平）
- 常用算法：LASSO回归、弹性网络、随机森林回归
（3）聚类任务：亚型鉴定
- 输入：微生物组组成
- 输出：患者亚型（enterotypes，肠型）
- 肠型概念：人类肠道微生物组可分为几种主要的组成模式
- 肠型1：以Bacteroides为主
- 肠型2：以Prevotella为主
- 肠型3：以Ruminococcus为主