用 GCModeller 做泛基因组分析（一）：算法原理与香农信息熵

面向读者：具有一定 R 语言编程基础的生物信息学初学者、实验室里的微生物研究者。
本篇讲清楚两件事：GCModeller 的泛基因组分析是怎么算的，以及贯穿整个分析体系的"香农信息熵"到底在度量什么。
全文不依赖任何编程知识，算法原理一律用公式和自然语言伪代码讲解。实战操作请看系列第二篇。

1. 什么是泛基因组？

把大肠杆菌的两个菌株的基因组摆在一起看，你会发现一个 surprising 的事实：这两个"同一个物种"的基因组，基因组成可能相差 20% 以上。K-12 实验室菌株约 4600 个基因，而某些致病菌株可以超过 5900 个基因，其中相当一部分在 K-12 里根本不存在。

单一"参考基因组"无法代表一个物种的全部遗传潜力。于是有了泛基因组（pan-genome）的概念：

泛基因组 = 该物种全部个体基因组的并集，由三部分构成：

• 核心基因组（core genome）：所有个体都有的基因，通常是管家基因、中心代谢、遗传信息处理相关功能；
• 附属基因组（accessory genome）：只在部分个体中存在的基因，包括"壳基因"（中等频率出现）和"云基因"（只在极少数个体中出现），往往与血清型、致病性、耐药性、环境适应直接相关。

泛基因组分析要回答的核心问题就是：哪些基因是 everyone 都有的？哪些是部分人有的？哪些是个别人独有的？某个基因组离"平均"有多远？

2. GCModeller 泛基因组分析流水线总览

GCModeller（https://github.com/xieguigang/GCModeller）是一个 .NET 平台上的基因组学分析引擎。它的泛基因组分析模块由 R# 脚本语言（GCModeller 自带的 R 语言方言）驱动，整体数据流如下：

GenBank 文件库 (*.gbff)
      │  ① 导出蛋白序列 + 基因位置表     [pack_genes.R]
      ▼
proteins.faa  +  genes.csv
      │  ② CD-HIT 式序列聚类 → 基因家族   [kmers 包]
      │  ③ 并查集合并 → 家族分组
      ▼
基因家族 (gene families)
      │  ④ 构建 PAV 矩阵 + 家族四分类     [cdhit_family.R]
      │  ⑤ 泛基因组曲线 / 遗传距离 / 共线性
      │  ⑥ 结构变异检测
      │  ⑦ 香农信息熵 + PCA + KMeans
      ▼
result.html 报告 + 一组 CSV 统计表格

分析引擎由 GCModeller 源码库中的几个核心模块组成：泛基因组分析模块（含主分析器 GenomeAnalyzer、SV 信息熵计算模块 SVDomainEntropy、统计模块 PanGenomeStats）、提供序列聚类算法底层的序列比对工具包（SequenceAlignment 模块），以及对外暴露 R# 函数接口的比较基因组学工具包（comparative_toolkit 包的 pangenome 模块）。下面按流水线顺序拆解每一步的算法。

3. 第一步：从百万条序列到基因家族

3.1 为什么不能两两比对？

一个中等规模的泛基因组项目（比如 800 个大肠杆菌基因组）大约包含 380 万条蛋白序列。如果对每两条序列都做一次序列比对，需要的比对次数是千万亿级别的——完全不可行。

CD-HIT 算法给出了一条务实的捷径，GCModeller 内置的 CD-HIT 式聚类（kmers 包的 cdhit_clusters 函数）遵循同样的思路：

长序列优先、贪心占坑：先把序列按长度从长到短排序，然后从最长的序列开始，每条序列要么加入一个已有的家族（当它与该家族代表序列的相似度足够高时），要么自立门户成为一个新家族的代表。

用自然语言伪代码描述：

输入: 全部蛋白序列集合, 相似度阈值 t (本系列实验取 0.6)

1. 把所有序列按长度从长到短排序
2. 为每条序列制作一份"k-mer 指纹"(序列切成长度为 k 的小片段，
   用这些片段的哈希集合作为这条序列的快速"草图")
3. 依次处理每条序列 S:
     a. 用 k-mer 指纹粗筛：哪些已有家族"有可能"与 S 相似？
        (指纹差异大的家族直接排除，不做昂贵的真实比对)
     b. 对粗筛出来的少数候选家族的代表序列，计算 S 与它的
        真实序列同一性 (identity = 比对上一致位点的比例)
     c. 若 identity >= t  →  S 加入该家族
        否则              →  S 成为新家族的代表序列
4. 输出: 所有家族

k-mer 指纹在这里的作用类似于图书馆的索书号：不需要逐页比对两本书，只需要比较索引卡片就能快速排除绝大多数无关的书目。这一步把"每条序列对全部家族做比对"变成了"每条序列只对极少数候选做比对"，聚类才可能在普通工作站上完成。

3.2 从簇到家族：并查集

聚类得到的每个簇就是候选基因家族。在 GCModeller 的实现里，家族的维护用的是数据结构中经典的并查集（Union-Find）。可以把它想象成社交网络里的"朋友圈"：开始时每个人自成一座孤岛，每当确认两个基因同源，就把他们所在的两个朋友圈合并成一个。一个包含 k 个基因的家族，只需要 k-1 次"合并"操作就能建好。

4. 第二步：PAV 矩阵——泛基因组的数据基石

有了基因家族，我们把每个家族在每个基因组中的基因个数统计出来，就得到了整个泛基因组分析中最核心的数据结构——PAV 矩阵（Presence/Absence Variation，存在/缺失变异矩阵）：

矩阵元素是拷贝数：0 表示该基因组没有这个家族，1 表示单拷贝，2 以上表示多拷贝。行是基因家族，列是基因组。

对每一个家族，计算它的存在比例：

$r = \frac{\text{存在该家族的基因组个数}}{\text{基因组总数 } N}$

按存在比例把所有家族分成四类（GCModeller 默认阈值，可通过 R# 接口调整）：

此外分析器还会顺带打上两个标签：特有基因（只在 1 个基因组中出现的家族）和单拷贝直系同源（在所有基因组中都恰好只有一份的家族，做进化分析的黄金材料）。

为什么软核心阈值很重要？ 100% 核心的判定极其严格：只要有一个基因组漏注释了一个基因，该家族就会跌出核心集合。软核心（比如"80% 的基因组都有"）对注释噪声和个别异常基因组更宽容，在基因组数量很大时更能反映"稳定遗传的核心功能"的真实规模。

5. 泛基因组曲线：物种的基因池还在扩张吗？

泛基因组曲线回答一个问题：随着不断加入新的基因组，还会不断发现新的基因吗？

做法是蒙特卡洛模拟：把 N 个基因组的加入顺序随机打乱，逐个加入，记录每加入一个新基因组后"已发现的家族总数"（泛基因组大小）和"在已加入基因组中全部存在的家族数"（核心基因数）；这样的随机顺序重复 100 次，取平均值抹平随机波动。

重复 100 次:
    随机打乱基因组的加入顺序
    家族计数器清零
    依次加入每个基因组:
        对该基因组携带的每个家族:
            若该家族此前从未见过   → 泛基因组大小 +1
            若该家族此前每个基因组都有 → 核心基因数 +1（若它在这个新
                                      基因组中也存在，核心基因数保持；
                                      否则下一次不再计入）
取 100 次模拟的平均值 → 泛基因组曲线

曲线形态本身就是结论：

• 封闭泛基因组（closed）：曲线早早走平——这个物种的基因库已被采样殆尽；
• 开放泛基因组（open）：曲线持续攀升——每加一个基因组都在冒出新基因，说明该物种的基因水平转移极其活跃。

微生物学里的经验规律（Heaps 律）告诉我们：大多数自由生活的细菌物种都是开放泛基因组。第二篇博客的大肠杆菌实测数据会印证这一点。

6. 基因组间的遗传距离：Jaccard 距离

衡量两个基因组"差多远"，GCModeller 采用基于基因家族集合的 Jaccard 距离：

$d(A, B) = 1 - \frac{|F_A \cap F_B|}{|F_A \cup F_B|}$

其中 $F_A$、$F_B$ 分别是两个基因组所拥有的基因家族集合。交集越大（共有的家族越多），距离越接近 0；两个基因组各玩各的，距离趋近 1。

这个定义完全基于基因的有无，不需要序列比对，因此对泛基因组这种"基因组成差异巨大"的场景特别合适。算出全部基因组对之间的距离后，就得到了一张 N × N 的遗传距离矩阵，后续可以做聚类树、多维缩放（MDS）等下游分析，直观呈现菌株群体按血清型/致病型/生态来源的分化格局。

工程细节（供好奇者）：GCModeller 把每个基因组的家族集合编码成一个二进制位图，两个基因组的交集大小用位图按位与（AND）配合 popcount 指令计算，单对基因组的距离计算只需"家族数 ÷ 64"次机器指令，几百个基因组两两比较也只是一眨眼的功夫。

7. 共线性与结构变异

共线性（collinearity）指的是：两个基因组上，一串同源基因不仅同源，而且保持相同的排列顺序。找出共线性区块，就能看清基因组的大尺度结构是保守的——大肠杆菌的不同菌株之间，染色体的基因顺序通常高度共线。

在共线性的基础上，GCModeller 逐家族扫描 PAV 矩阵和拷贝数，检测五类结构变异（SV）事件：

这些 SV 事件正是下一节 SV 信息熵分析的原料。

8. 核心章节：泛基因组中的香农信息熵

信息熵是本系列的重头戏。GCModeller 在两个层面使用香农信息熵：基因组层面（一个基因组压缩成一个熵值）和基因家族层面（一个家族压缩成两个熵值）。

8.1 一分钟看懂香农信息熵

香农（Claude Shannon）1948 年给出的定义：一个离散随机变量的熵是它所有可能取值带来的"平均意外程度"：

$H = -\sum_{i} p_i \log p_i$

其中 $p_i$ 是第 i 种取值出现的概率。约定 $0 \log 0 = 0$（不可能发生的事件不贡献信息）。

直觉版本：

• 抛一枚均匀硬币（正/反各 50%），结果最难猜，熵最大；
• 抛一枚两面都是正面的硬币，结果毫无悬念，熵为 0；
• 取值越"摊得平"、越均匀，熵越大；越"扎堆"在少数取值上，熵越小。

（以 2 为底时熵的单位是比特；GCModeller 的实现使用自然对数，单位是 nat，数值上乘以换算系数 1/ln2 即得比特数。）

把这套度量搬进泛基因组：把"某个基因组有哪些家族"或"某个家族出现在哪些基因组"看作一个概率分布，熵就成了一把度量基因多样性的尺子。

8.2 基因组层面的熵（一）：存在/缺失均衡度 —— GCModeller 的实现

设泛基因组共 N 个基因家族，某个基因组存在其中 K 个、缺失 N-K 个。以"存在/缺失"两种状态定义概率分布：

$p = \frac{K}{N}, \qquad H = -\Big(p \log p + (1-p) \log (1-p)\Big)$

这就是把 PAV 矩阵的一整列（一个基因组对全部 N 个家族的存在状态）压缩成一个数字。由于只有两种状态，这个熵在 p = 0.5 时达到最大值 $\log 2 \approx 0.693$ nat，p 越偏离 0.5（无论是趋近 0 还是趋近 1）熵越低。

生物学解读：

1. 基因组可塑性的代理指标。大多数细菌基因组中，泛基因组家族的绝大多数是附属基因，因此任一基因组的 p 通常远小于 0.5。如果一个基因组通过水平基因转移（HGT）大量获取附属基因，K 增大、p 向 0.5 靠拢，熵值升高；反之，经历基因大量丢失的基因组熵值走低。
2. 识别极端基因组。熵值极低的基因组往往对应生活在封闭稳定环境中的内生菌、专性寄生菌或长期实验室驯化菌株（基因组缩减）；熵值偏高的基因组则常见于土壤、肠道等复杂环境中的泛化型菌株，它们需要庞大的附属基因库来应对多变的环境。
3. 作为表型补充参与关联分析。可以把"基因组熵"本身当作一种量化表型：找出哪些家族的存在/缺失显著拉高或压低某个基因组的熵，这些家族往往就是驱动基因组分化的关键变异。

一个亲手算的例子（数字来自本系列第二篇的真实结果）：801 个大肠杆菌基因组共聚出 N = 252,301 个基因家族。实验室标准株 K-12 MG1655 存在 K = 2,421 个家族：

$p = \frac{2421}{252301} \approx 0.0096, \qquad H = -(0.0096 \log 0.0096 + 0.9904 \log 0.9904) \approx 0.054 \text{ nat}$

而致病血清型 O157:H7 的一个菌株 K = 5,645，算出 H ≈ 0.107 nat——几乎翻倍。实验室驯化株与携带大量毒力因子的致病株，在熵这一把尺子上被清晰地区分了开来。

8.3 基因组层面的熵（二）：进化独特性（延伸思路）

除了均衡度熵，泛基因组分析中还有第二种常见的基因组级熵定义，思想类似文本挖掘里的 TF-IDF——越罕见的东西信息量越大：

先对每个家族 i 计算它在全部 M 个基因组中的出现频率 $p_i$，其信息量定义为：

$I_i = -\log p_i$

（全物种都有的核心基因 $p_i \approx 1$，信息量趋近 0；只在 1 个基因组里出现的特有基因 $p_i = 1/M$，信息量最大。）再把基因组 j 所携带的所有家族的信息量加总：

$Hj = \sum{i \in G_j} -\log p_i$

这种"熵"度量的是基因组携带了多少罕见遗传信息：高值基因组拥有大量别人没有的基因，往往对应独特的生态功能（耐药、特殊代谢通路、宿主互作因子）或致病潜力。

GCModeller 当前默认实现的是 8.2 节的第一种定义（存在/缺失均衡度）；第二种定义可以直接基于导出的 PAV 矩阵在 R 里自行计算，适合作为深入分析的手段。论文写作时务必注明采用的是哪种概率分布定义。

8.4 家族层面的熵：SV 结构变异信息熵 —— GCModeller 的实现

熵的视角还可以从"基因组"翻转到"基因家族"：一个家族在群体中的变异模式，本身也携带进化信息。

先把结构变异事件组织成两个矩阵（行 = 有 SV 事件的家族，列 = 基因组）：

• CopyNumber 矩阵：元素 $c(f,g)$ = 家族 f 在基因组 g 中该 SV 的拷贝数；
• Median 矩阵：元素 $m(f,g)$ = 该家族 SV 的中位拷贝数。

① CopyNumber 熵——以"各基因组拷贝数占比"作为离散概率分布：

$H{cn}(f) = -\sum{g} \frac{c(f,g)}{\sum{g'} c(f,g')} \log \frac{c(f,g)}{\sum{g'} c(f,g')}$

• 高 $H_{cn}$：拷贝数在各基因组间千差万别（有的 1 份、有的 5 份、有的 10 份）→ 该家族经历了频繁的拷贝数扩张/收缩；
• 低 $H_{cn}$：所有基因组清一色相同拷贝数 → 剂量受到严格的纯化选择。
• 取值为 0 的基因组（无 SV 事件）自然不参与求和，不干扰结果。

② Median 熵——把每一行的取值当作离散类别统计频率：

$H{med}(f) = -\sum{v} q_v \log q_v$

在 GCModeller 当前数据模型中，同一家族的 SV 中位拷贝数是家族级常量，因此每行的取值实际上只有 {0, 中位数} 两类，$H_{med}$ 刻画的是"该家族的 SV 状态在基因组之间分布的均衡程度"：

• $H_{med}$ 越高 → 有 SV 的基因组与没有的大约各占一半（高度多态）；
• $H_{med}$ 越低 → 要么几乎人人都有、要么几乎人人没有（状态一致）。

（若上游检出流程把每个基因组各自的 SV 尺寸/断裂点都记录下来，这个熵便自然升级为"SV 结构特征的离散程度"——直接度量各基因组间断裂点与尺寸形态的多样性——而算法模块本身无需任何改动。）

③ 两个熵 → 一张散点图 → KMeans 聚类

把每个家族表示成平面上的点 $(H{cn}, H{med})$，再用 KMeans 聚类（默认 k = 4）自动划出具有不同演化特征的家族类群。为了让聚类结果稳定可靠，GCModeller 在工程实现上做了三步处理：

1. 过滤稀疏数据：SV 出现频率低于 5% 基因组的家族不参与聚类——这种家族的熵容易被噪声主导；
2. Z-score 标准化：两个熵的量纲可能不同，聚类前都变换成"距均值若干个标准差"，防止方差大的维度霸占距离计算；
3. 确定性簇编号：KMeans 的初始中心是随机的，聚类完成后按簇中心坐标重排簇编号，保证同一数据重跑得到完全一致的结果。

四类典型进化模式（象限解读）：

把聚类得到的各簇家族拿去做 GO/KEGG 富集，常能发现"某类功能基因偏爱某种 SV 进化策略"的规律——这是从泛基因组数据里挖掘进化故事的富矿。

8.5 小结：熵把 0/1 海洋压缩成一条连续的尺子

• 基因组层面：25 万 × 800 的 0/1 矩阵 → 每个基因组一个熵值 → 直接比较基因组可塑性；
• 家族层面：两个熵值（$H{cn}$，$H{med}$）→ 二维特征空间 → 自动划出保守型/剂量调谐型/混沌型/结构微调型四类进化模式。

相比"基因总数""基因组大小"这类粗放计数，香农熵同时考虑了数量的分布与稀有度，是更稳健的数学度量。

9. GCModeller 泛基因组分析 API 速查表

9.1 数据准备（seqtoolkit 包）

9.2 分析与导出（comparative_toolkit 包的 pangenome 模块）

所有结果表格都实现了 as.data.frame 泛型，可以无缝接入熟悉的 R 风格数据处理流程。

10. 小结

• 泛基因组分析的本质是：序列聚类建家族 → PAV 矩阵 → 从矩阵中榨取统计量；
• GCModeller 用 CD-HIT 式贪心聚类 + 并查集解决规模问题，用蒙特卡洛模拟画泛基因组曲线，用 Jaccard 距离刻画基因组分化；
• 香农信息熵是贯穿这套体系的统一语言：基因组层面的熵度量可塑性，家族层面的双熵度量结构变异的进化模式；
• 全部分析通过两条 R# 脚本即可完成——具体怎么跑、结果怎么读，请看系列第二篇：《用 GCModeller 做泛基因组分析（二）：801 个大肠杆菌基因组实战》。

• Author
• Recent Posts

谢桂纲

高级数据科学家 at 苏州帕诺米克

Working on Engineered bacteria CAD design on its genome from scratch. Writing scientific computing software for Tianhe & Sunway TaihuLight supercomputer. Do scientific computing programming in R/R# language, he is also the programming language designer of the R# language on the .NET runtime.

Latest posts by 谢桂纲 (see all)

• 用 GCModeller 做泛基因组分析（二）：801 个大肠杆菌基因组实战 - 2026年9月26日
• 用 GCModeller 做泛基因组分析（一）：算法原理与香农信息熵 - 2026年9月26日
• 给果蝇通电，让它替你玩贪吃蛇——兼谈”缸中脑”与数字永生 - 2026年9月25日