Blog Post
Velvet 基因组组装器详解：从论文到 C 代码实现
系统讲解 Velvet 基因组组装算法的原理、设计和实现，包括 de Bruijn 图构建、错误校正和图简化等核心模块
引言
基因组组装（Genome Assembly）是生物信息学中最基础也最具挑战性的问题之一。2008 年，Daniel Zerbino 和 Ewan Birney 在《Genome Research》上发表了 Velvet 算法，开创了短读长基因组组装的新纪元。
本文是我学习基因组组装算法的笔记，旨在从论文原理、算法设计和C 代码实现三个层面，系统梳理 Velvet 的来龙去脉。
作者注：本文所有内容均为本人学习和理解的总结，如有错误，欢迎指正。
一、背景与挑战
1.1 基因组组装问题
基因组组装的核心问题是：给定大量短的 DNA 序列片段（reads），重建原始的基因组序列。
Reads: ATGCGA CGATGC GCTAGC TACGCT
|||||| |||||| |||||| ||||||
ATGCGA CGATGC GCTAGC TACGCT
Assembly: ATGCGATGCGCTAGCTACGCT对于短读长测序（如 Illumina），我们面临三大挑战：
- 数据量大：数百万到数十亿条 reads
- 测序错误：错误率约 0.1%-1%
- 重复序列：基因组中存在大量重复区域
1.2 为什么选择 de Bruijn 图？
在 Velvet 之前，主流的组装算法基于重叠 - 布局 - 共识（Overlap-Layout-Consensus, OLC）方法。但 OLC 需要计算所有 reads 之间的两两比对，时间复杂度为 O(n²)，无法处理海量短读长数据。
Velvet 采用 de Bruijn 图（de Bruijn Graph, DBG）方法，将序列分解为固定长度的 k-mer，时间复杂度降为 O(n)。
二、de Bruijn 图的理论基础
2.1 k-mer 与 de Bruijn 图
定义：对于一个 DNA 序列，将其拆分为所有长度为 k 的子串，称为 k-mer。
序列：ATGCGAT
k = 3 时的 k-mers:
ATG → TGC → GCG → CGA → GAT
每个 k-mer 是图中的一个节点，相邻 k-mer 之间有边连接de Bruijn 图的构建：
- 节点（Node）：每个唯一的 k-mer
- 边（Edge）：如果两个 k-mer 有 k-1 个碱基重叠，则存在一条有向边
2.2 双向链表示
由于 DNA 是双链结构，Velvet 同时存储 正向链 和 反向互补链：
// DNA 碱基互补配对
A ↔ T
C ↔ G
// 反向互补示例
// 正向：5'-ATGCGA-3'
// 反向：3'-TACGCT-5'
// 反向互补：5'-TCGCAT-3'在 de Bruijn 图中，每个 k-mer 节点都有对应的反向互补节点，两者通过特殊边连接。
三、Velvet 算法设计
3.1 算法流程概览
┌─────────────────────────────────────────────────────────────┐
│ Velvet Assembly Pipeline │
├─────────────────────────────────────────────────────────────┤
│ │
│ 1. 预处理 │
│ └─→ 读取 FASTQ 文件，提取 k-mers │
│ │
│ 2. 构建 de Bruijn 图 │
│ └─→ 创建节点和边，建立哈希表索引 │
│ │
│ 3. 错误校正 │
│ └─→ 去除测序错误产生的"气泡"和"尖端" │
│ │
│ 4. 图简化 │
│ └─→ 合并线性路径，解析简单分支 │
│ │
│ 5. 配对末端信息整合 │
│ └─→ 利用 paired-end reads 解决重复区域 │
│ │
│ 6. 路径搜索与 contig 输出 │
│ └─→ 寻找最优路径，输出组装结果 │
│ │
└─────────────────────────────────────────────────────────────┘3.2 核心数据结构
Velvet 的核心数据结构设计非常精巧：
// 简化版的节点结构（基于 Velvet 源码）
typedef struct _Node {
char *sequence; // k-mer 序列
long coverage; // 覆盖度（支持该 k-mer 的 reads 数）
struct _Arc *arcs; // 出边列表
struct _Node *next; // 哈希表链地址法的下一个节点
unsigned int flags; // 节点标记（用于图简化）
} Node;
// 边的结构
typedef struct _Arc {
struct _Node *destination; // 目标节点
int overlap; // 重叠长度（通常为 k-1）
struct _Arc *next; // 下一条边
} Arc;
// 图的全局结构
typedef struct _Graph {
Node **hashTable; // 哈希表，用于快速查找 k-mer
unsigned int kmerSize; // k 值
long totalNodes; // 节点总数
long totalArcs; // 边总数
} Graph;四、C 代码实现详解
4.1 k-mer 编码
Velvet 使用 2-bit 编码 将 DNA 序列压缩存储：
// 2-bit 编码表
// A = 00, C = 01, G = 10, T = 11
// 将 k-mer 转换为整数（用于哈希表索引）
unsigned long encodeKmer(const char *sequence, int k) {
unsigned long hash = 0;
for (int i = 0; i < k; i++) {
hash <<= 2; // 左移 2 位
switch (sequence[i]) {
case 'A': case 'a': hash |= 0; break; // 00
case 'C': case 'c': hash |= 1; break; // 01
case 'G': case 'g': hash |= 2; break; // 10
case 'T': case 't': hash |= 3; break; // 11
default:
// 遇到 N 或其他模糊碱基，返回特殊值
return INVALID_KMER;
}
}
return hash;
}
// 反向互补 k-mer
unsigned long reverseComplementKmer(unsigned long kmer, int k) {
unsigned long rc = 0;
for (int i = 0; i < k; i++) {
// 取出最低 2 位，取反，放到对应位置
unsigned int base = kmer & 3; // 取出最低 2 位
base = 3 - base; // A↔T, C↔G
rc = (rc << 2) | base; // 左移并添加新碱基
kmer >>= 2; // 右移处理下一个碱基
}
return rc;
}4.2 哈希表构建
// 初始化哈希表
Graph* initGraph(unsigned int kmerSize, unsigned int tableSize) {
Graph *graph = (Graph*)malloc(sizeof(Graph));
graph->kmerSize = kmerSize;
graph->totalNodes = 0;
graph->totalArcs = 0;
// 分配哈希表
graph->hashTable = (Node**)calloc(tableSize, sizeof(Node*));
return graph;
}
// 向图中添加 k-mer
Node* addKmerToGraph(Graph *graph, const char *sequence, long readId) {
unsigned long hash = encodeKmer(sequence, graph->kmerSize);
unsigned int index = hash % HASH_TABLE_SIZE;
// 检查是否已存在
Node *node = graph->hashTable[index];
while (node != NULL) {
if (strcmp(node->sequence, sequence) == 0) {
node->coverage++;
return node;
}
node = node->next;
}
// 创建新节点
Node *newNode = (Node*)malloc(sizeof(Node));
newNode->sequence = strdup(sequence);
newNode->coverage = 1;
newNode->arcs = NULL;
newNode->next = graph->hashTable[index];
newNode->flags = 0;
graph->hashTable[index] = newNode;
graph->totalNodes++;
return newNode;
}4.3 图的构建
// 从 reads 构建 de Bruijn 图
Graph* buildDeBruijnGraph(char **reads, int numReads, int k) {
Graph *graph = initGraph(k, HASH_TABLE_SIZE);
for (int r = 0; r < numReads; r++) {
int readLen = strlen(reads[r]);
// 滑动窗口提取 k-mers
for (int i = 0; i <= readLen - k; i++) {
char *kmer = strndup(reads[r] + i, k);
// 添加节点
Node *currentNode = addKmerToGraph(graph, kmer, r);
// 添加边（与前一个 k-mer 连接）
if (i > 0) {
char *prevKmer = strndup(reads[r] + i - 1, k);
unsigned long prevHash = encodeKmer(prevKmer, k);
unsigned int prevIndex = prevHash % HASH_TABLE_SIZE;
Node *prevNode = graph->hashTable[prevIndex];
while (prevNode != NULL) {
if (strcmp(prevNode->sequence, prevKmer) == 0) {
addArc(prevNode, currentNode, k - 1);
break;
}
prevNode = prevNode->next;
}
free(prevKmer);
}
free(kmer);
}
}
return graph;
}
// 添加边
void addArc(Node *from, Node *to, int overlap) {
Arc *arc = (Arc*)malloc(sizeof(Arc));
arc->destination = to;
arc->overlap = overlap;
arc->next = from->arcs;
from->arcs = arc;
from->outDegree++;
}4.4 错误校正：去除尖端（Tips）
测序错误会在图中产生短的”尖端”结构：
A → T → G → C
/
ATG → C
\
G → A → T// 去除尖端：短的线性路径
void removeTips(Graph *graph, int maxTipLength) {
for (unsigned int i = 0; i < HASH_TABLE_SIZE; i++) {
Node *node = graph->hashTable[i];
while (node != NULL) {
// 检查是否为尖端起点（入度为 0，出度为 1）
if (node->inDegree == 0 && node->outDegree == 1) {
int tipLength = 0;
Node *current = node;
// 遍历尖端路径
while (current != NULL &&
current->inDegree <= 1 &&
current->outDegree == 1 &&
tipLength < maxTipLength) {
tipLength++;
current = current->arcs->destination;
}
// 如果尖端足够短，删除整个路径
if (tipLength < maxTipLength) {
removePath(graph, node);
}
}
node = node->next;
}
}
}4.5 错误校正：去除气泡（Bubbles）
测序错误还会产生”气泡”结构：
A → T → G
/ \
ATG → C C → GAT
\ /
G → C → A// 去除气泡：并行路径
void removeBubbles(Graph *graph, int maxBubbleLength) {
for (unsigned int i = 0; i < HASH_TABLE_SIZE; i++) {
Node *node = graph->hashTable[i];
while (node != NULL) {
// 检查是否有多个出边（气泡起点）
if (node->outDegree >= 2) {
Arc *arc1 = node->arcs;
Arc *arc2 = arc1->next;
while (arc1 != NULL && arc2 != NULL) {
// 寻找两条路径的汇合点
Node *mergePoint = findMergePoint(
arc1->destination,
arc2->destination,
maxBubbleLength
);
if (mergePoint != NULL) {
// 保留覆盖度更高的路径
if (arc1->destination->coverage >
arc2->destination->coverage) {
removePath(arc2->destination, mergePoint);
} else {
removePath(arc1->destination, mergePoint);
}
}
arc2 = arc2->next;
}
}
node = node->next;
}
}
}4.6 图简化：合并线性路径
// 合并非分支路径
void simplifyGraph(Graph *graph) {
int changed = 1;
while (changed) {
changed = 0;
for (unsigned int i = 0; i < HASH_TABLE_SIZE; i++) {
Node *node = graph->hashTable[i];
while (node != NULL) {
// 入度=1 且 出度=1 的节点可以合并
if (node->inDegree == 1 && node->outDegree == 1) {
Node *nextNode = node->arcs->destination;
// 如果下一个节点也是入度=1 且 出度=1
if (nextNode->inDegree == 1 && nextNode->outDegree == 1) {
mergeNodes(node, nextNode);
changed = 1;
}
}
node = node->next;
}
}
}
}
// 合并两个节点
void mergeNodes(Node *node1, Node *node2) {
// 拼接序列（重叠部分只保留一次）
int overlap = node1->arcs->overlap;
int newLen = strlen(node1->sequence) + strlen(node2->sequence) - overlap;
char *newSeq = (char*)malloc(newLen + 1);
strcpy(newSeq, node1->sequence);
strcat(newSeq, node2->sequence + overlap);
free(node1->sequence);
node1->sequence = newSeq;
// 更新边
node1->arcs = node2->arcs;
node1->outDegree = node2->outDegree;
// 标记 node2 为已删除
node2->flags |= NODE_DELETED;
graph->totalNodes--;
}五、配对末端信息整合
Velvet 的核心创新之一是有效利用 paired-end reads 信息解决重复区域。
5.1 配对末端约束
Read1: 5'────────────3'
←── insert size ──→
5'────────────3'
Read2
已知：Read1 和 Read2 之间的距离约为 insert size（±标准差）5.2 约束传播算法
// 使用配对末端信息解析重复区域
void resolveRepeatsWithPairedEnds(Graph *graph, PairedEndReads *pairs) {
for (int p = 0; p < pairs->numPairs; p++) {
Node *node1 = findNodeByKmer(graph, pairs[p].read1);
Node *node2 = findNodeByKmer(graph, pairs[p].read2);
if (node1 != NULL && node2 != NULL) {
// 添加距离约束
addDistanceConstraint(node1, node2,
pairs[p].insertSize,
pairs[p].stdDev);
}
}
// 传播约束，解决冲突
propagateConstraints(graph);
}六、路径搜索与 Contig 输出
6.1 寻找最优路径
// 深度优先搜索寻找最长路径
void findContigs(Graph *graph, FILE *output) {
for (unsigned int i = 0; i < HASH_TABLE_SIZE; i++) {
Node *node = graph->hashTable[i];
while (node != NULL) {
// 从起始节点（入度=0）开始搜索
if (node->inDegree == 0 && !(node->flags & NODE_VISITED)) {
List *path = createList();
dfsTraverse(node, path, output);
freeList(path);
}
node = node->next;
}
}
}
// DFS 遍历
void dfsTraverse(Node *node, List *path, FILE *output) {
node->flags |= NODE_VISITED;
appendToList(path, node->sequence);
if (node->outDegree == 0) {
// 到达终点，输出 contig
writeContig(output, path);
} else {
// 选择覆盖度最高的边继续搜索
Arc *bestArc = selectBestArc(node);
if (bestArc != NULL) {
dfsTraverse(bestArc->destination, path, output);
}
}
}七、总结与思考
7.1 Velvet 的核心贡献
| 特性 | 描述 |
|---|---|
| de Bruijn 图 | 将 O(n²) 复杂度降为 O(n) |
| 错误校正 | 自动去除测序错误产生的气泡和尖端 |
| 配对末端整合 | 利用距离约束解决重复区域 |
| 内存优化 | 2-bit 编码，紧凑存储 |
7.2 局限性
- k 值选择：需要手动指定，不同 k 值影响组装质量
- 高重复基因组：对于重复序列比例高的基因组效果有限
- 杂合度：对高杂合度基因组的处理不够理想
7.3 后续发展
Velvet 之后，基因组组装算法继续演进：
- SPAdes：多 k-mer 策略，自动选择最优 k 值
- MEGAHIT：succinct de Bruijn 图，更低内存占用
- Flye/Canu：针对长读长（PacBio/Nanopore）的组装算法
参考文献
-
Zerbino DR, Birney E. Velvet: algorithms for de novo short read assembly using de Bruijn graphs. Genome Research. 2008;18(5):821-829. doi:10.1101/gr.074492.107
-
Compeau PEC, Pevzner PA, Tesler G. How to apply de Bruijn graphs to genome assembly. Nature Biotechnology. 2011;29(11):987-991. doi:10.1038/nbt.2023
-
Velvet GitHub Repository. https://github.com/dzerbino/velvet
-
Pevzner PA, Tang H, Waterman MS. An Eulerian path approach to DNA fragment assembly. PNAS. 2001;98(17):9748-9753. doi:10.1073/pnas.171285098
附录：完整代码示例
完整的 Velvet 风格 de Bruijn 图实现可以在我的 GitHub 上找到：