倒排索引在基因测序中的应用

[郭行健]

最近在听MIT的一门定量生物学慕课（courses.edx.org/courses/MITx/7.QBWx/2T2014/），在基因测序一节中居然也提到了倒排索引，还是挺有趣的，在这里给大家简要介绍一下：

人类基因组总共约有30亿个碱基对，顺序读取显然太慢而且稳定性会有问题，因此通常的做法是：把DNA样本切成众多长约300碱基对的微小片段，同时准备一个植于玻璃片上的测序阵列（每一个测序节点的尺度只有微米级，因此集成度是很大的），如此一来这些片段的测序便可以很快完成。

假如我们手头有100份全同的DNA样本，每一个都经由上述步骤处理，那么从不同样本切下的片段间或多或少会有重叠，这一信息给我们提供了拼接的依据。如何从这些片段中复原整个基因组呢？这本质上是一个算法问题：

输入：30亿÷300×100=10亿个DNA片段，譬如说[(1, ATGTATTCCATCA...ATGGG), (2, GTATTCCATCA...GGGCC), ...]

输出：总DNA序列

要求：利用不同片段间的重叠信息完成拼接。

很朴素的一个贪心算法是：先用第一个片段与余下所有片段作比较，找到重叠最好的那一个拼在一起，重复这一步骤直至所有片段都拼在一起——这显然是个O(N^2)的算法，对N=10亿的数据集显然是不可行的！倒排索引便在这时出场了。

仔细想想，上述算法问题与全文搜索有着许多类似之处：我们都需要在海量的数据集中找到与搜索关键词（DNA片段）相关度高（重叠片段长）的那些文档（其他片段）——与文档（片段）总数相比，它们通常是很少的，倒排索引正是利用了这一语义上“稀疏性”，省下了遍历查找上浪费的无用功。分布式计算则进一步提高了构建倒排索引的效率。

为了在测序问题上利用倒排索引，我们需要定义基因序列的term是什么。为了抓住“重叠”这一核心，我们可以将片段以10个碱基为单位有重复地做划分——例如片段ATGTATTCCATCA...ATGGG的terms就是(1, ATGTATTCCA), (2, TGTATTCCAT), (3, GTATTCCATC), ...实际上一个碱基用2位表示就足够了，因此存储整个term只需3个字节。另外所有可能出现的term只有4^10≈100万个，整个索引完全可以放在内存里，带来了许多方便。

倒排索引建好之后该怎么做课程里并没有说，我想贪心算法应该可以，或者可以把问题转化成最小生成树？实际的做法要来得更复杂一些，一方面需要考虑测序中出现差错的可能，另一方面“转座子（transposable elements）”的大量存在会使某些posting list格外长而需要额外处理——但利用倒排索引这一核心idea是不变的。

上述是在未知总基因序列的情形下测序的算法——一旦我们测定了第一个个体的基因组，此后其他个体的测序就可以抄捷径了：我们可以利用先前测得的序列作为“样板”确定某一片段最可能出现的位置，所借助的当然还是倒排索引。类似的方法还可以用来确定某一mRNA片段所属的基因。

倒排索引在生物学上或多或少出人意料的应用，充分说明了大数据手段与方法的广泛应用空间——这也是cs402开设的意义所在吧。

附：Eric Lander教授关于基因测序的授课视频（挂在YouTube上；6，7，8，9与本文内容直接有关）

1 2 3 4 5 6 7 8 9 10 11

[Krasus C]

赞

在 2014年7月18日星期五UTC+8下午11时52分08秒，郭行健写道：

[杨博文]

膜拜腿腿……

这居然也可以啊

在 2014年7月18日星期五UTC+8下午11时52分08秒，郭行健写道：

[Yq Peng]

edx的慕课是翻墙听的么@@

在 2014年7月18日星期五UTC+8下午11时52分08秒，郭行健写道：

[郭行健]

不用……视频可以下载来看，虽然说翻墙看一般来得快一些

在 2014年7月22日星期二UTC+8上午9时20分44秒，Yq Peng写道：