1. 背景介绍
病原体基因检测，为各种严重感染的诊断提供了基础。病原体检测流程分成五个步骤：（1）采集病人的样本，比方说静脉血，痰液，肺泡灌洗液，或者脑脊髓液等。（2）对样本进行培养，提取样本中的核酸组织。（3）通过高通量基因测序仪，对核酸序列进行测序。高通量测序为了保证精度，会将很长的核苷酸序列，切成小的分片，分别进行测序，在读取基因序列上，一般是 50pb 到 200bp 不等。（4）高通量测序之后，需要查询相关的病原体标准序列，找到匹配的基因序列。（5）通过将这些小的片段进行分析，得到全部基因片段的组成成分（即检测结果），从而为病人的疾病确诊以及精准治疗提供有力的支持。
图 1. 病原体基因检测流程
对于生物分析检测来说，通常情况下，一次病原体检测，大约会生成 5 亿左右的 75bp 的基因片段。过滤掉一些人的基因组织序列之后，还需要查询 1 亿左右的基因片段。通常情况下，会使用 nBlast [1] 工具来进行基因匹配，这部分在整个病原体检测过程中，也是非常耗时的部分，大约需要 2 - 3 个小时。阿里云 AnalyticDB 向量版提供了一个高效的基因检索工具，大大提升了基因分析的性能，能够在几十分钟内完成整个病原体的查询检测过程。
2. 基因检索应用
2.1 基因检索功能
图 1 展示了病原体基因检索的界面。当前演示包含了 12182 个病毒的碱基序列，我们将病毒切分成 150bp 的小片段（总共 1590804 个片段），转化成向量之后，存储到 AnalyticDB 中。在检索框中，用户可以输入一段基因序列，到我们当前的系统中进行检索。为了方便大家使用，我们挑了新冠病毒，艾滋病毒，埃博拉病毒和中东呼吸综合症的基因序列来进行演示，用户可以拷贝相关的序列，来检测查询的性能。
图 2. 核酸查询演示
图 3 用户输入了一段新型冠状病毒的基因序列，可以看到排到前面的序列片段就是我们要找的新冠的病毒的序列。因为当前 AnalyticDB 提供高效的向量索引，系统会在毫秒级的时间内，返回相关的基因片段。
图 3. 基因检索结果
2.2 端到端的基因分析
我们模拟了人体的基因采样，将新型冠病毒基因（塞尔维亚 MT450872 [2]，美国 MT450873 [3]）和中东呼吸症 MERS 基因（NC_019843.3 [10]），三株病毒混合在一起，打散成 75bp 的序列，当做测试集合。我们希望，通过当前病毒库的检索分析，能够识别出当前测试集合中包含新型冠状病毒和 MERS 病毒。当前系统通过匹配检测，我们生成图 4。
图 4. 基因匹配结果
可以看到经过检测，病毒库检索系统返回了三个基因组（NC_045512.2，NC_019843.3 和 NC_038294.1）。NC_045512.2（65%）是武汉海鲜市场的新型冠状病毒的基因；NC_019843.3（20%）是 MERS 病毒的基因；而 NC_038294.1（13%）是 beta 型英国冠状病毒，经过查询这个是 MERS 病毒的另外一个名字 [8]，也属于 MERS 病毒。因此，通过分析，当前混合测试集合中包含了新型冠状病毒和 MERS 病毒。
2.3 应用架构总体设计
阿里云基因检索系统的总体架构如图 5 所示，AnalyticDB 负责整个应用的全部的结构化数据（比方说，基因序列的长度，基因的名称，基因的种类，以及基因的详细介绍，DNA 或者 RNA 等）和基因序列产生的特征向量的存储和查询。在查询的时候，我们使用基因向量抽取模型，将基因转化成向量，在 AnalyticDB 库中进行粗排检索。在向量匹配的结果集中，我们使用经典的 Needleman-Wunsch [4] 算法进行精排，返回最相似的基因序列。
图 5. 基因检索系统框架
3. 基因训练和查询模型
3.1 基因查询过程
基因模型的训练已经在上一篇文章中 [5]，进行了详细的讲解。通过训练好的 DNA K-Mer 模型，我们可以得到每个 k -mer 的向量。给定一段 12bp 的基因序列（如图 6），我们在这段基因序列中抽取出 5 个 8 -mers。我们将这 5 个 8 -mers 转成对应的向量，求和归一化之后，就是这段 12bp 的基因序列的最终的向量。当然，为了提升精度，我们也可以使用 doc2vec [6] 等学习模型来对整段基因片段进行转化。
图 6.DNA 序列转向量
3.2 基因精度分析
我们训练了两个模型，全部病毒模型和 21 个病原体细菌模型（痤疮丙酸杆菌，金黄色葡萄球菌，表皮葡萄球菌，溶血葡萄球菌，大肠埃希氏菌，鲍曼不动杆菌，结核分枝杆菌，肺炎链球菌，肺炎克雷伯氏菌，流感嗜血杆菌，副流感嗜血杆菌，嗜麦芽窄食单胞菌，铜绿假单胞菌，屎肠球菌，纹带棒状杆菌，人疱疹病毒 4 型 (EB 病毒)，细环病毒，人腺病毒 B 组，黄曲霉，白色假丝酵母，耶氏肺孢子菌）。我们将一个基因，每隔 150 个 bp，做一下切分。然后将 150bp 的小的分段，转化成向量存在库里面，进行检索。因此病毒数据集包括 12182 个病毒, 1590804 个分段；21 个细菌共 275 个基因，1521807 个分段。
实验 1（见表 1）系统随机的在当前的基因库里面，取出 75bp 的小的片段。我们知道这 75bp 的基因片段是在哪个基因的哪个片段中提取的。我们将这 75bp 的基因段到库里面进行检索，查找返回前 N 个结果集中，查看是否包含这 75bp 的基因段对应的基因片段。Top- n 的精度（Precision(n)），用公式（1）进行计算。
其中，n 表示查询返回的列表的长度。u 表示查询的次数，在实验中 u 取的是 1000 次。 表示在第 i 次查询中，序列 si 是否出现在 Top- n 的列表中，出现为 1，不出现为 0。如果 n 越小，精度越高，说明我们的方法在实际中非常有效。我们可以看到，针对两个模型来说，top20 的精度，都在 99% 以上。精度在 0.99 以上，对于实际检验基因片段包含物种基因是足够了的。
表 1. 基因查询精度检测
数据集
top 1
top 2
top 3
top 4
top 5
top 10
top 20
病毒
0.866
0.965
0.983
0.986
0.99
0.992
0.994
21 个细菌
0.901
0.975
0.987
0.987
0.993
0.994
1.0
实验 2（见表 2）系统随机的在当前的基因库里面，取出 75bp 的小的片段。我们将这 75bp 的基因段进行了 2% 的随机的突变（自然界中的基因突变的概率会更低，比方说人的 30 亿个碱基，新生儿会有 30 个基因发生突变。病毒 RNA 的突变概率会高一些，一般也都小于 1%），然后到库里面进行检索，查找返回前 N 个结果集中，查看是否包含这 75bp 的基因段对应的基因片段。基因突变之后，虽然查询的精度有所下降，但是 top20 的精度也都达到了 0.99。
表 2. 基因突变查询精度检测
数据集
top1
top 2
top3
top4
top5
top6
top7
病毒
0.846
0.954
0.960
0.976
0.98
0.982
0.99
21 个细菌
0.884
0.961
0.968
0.973
0.973
0.989
1.0
实验 3（见表 3）对基因检索的速度进行了比较。我们下载了病毒序列，菌类的基因序列，以及部分植物基因序列 [7]，总共 9.7G。我们分别将相关数据，导入到 AnalyticDB 数据库中和 Blast 库中。我们跑了 100 次不同的查询，对实验结果取了平均。Blast 需要 3.22 秒才能返回结果，我们算法精度在 top30 的情况下，保证精度在 0.95 的准确性下，测试端到端的查询（包括查询基因转向量，向量粗排和 Needleman-Wunsch 算法的精排）只需要 0.257s（提升了 12.5 倍）。
表 3. 检索时间
4. 结尾
详细的基因模型以及相关系统信息请加入我们的钉钉群，欢迎大家讨论和使用。
[1] blast+ https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/
[2] https://www.ncbi.nlm.nih.gov/nuccore/MT450872
[3] https://www.ncbi.nlm.nih.gov/nuccore/MT450873
[4] Needleman, Saul B. & Wunsch, Christian D. (1970). “A general method applicable to the search for similarities in the amino acid sequence of two proteins”. Journal of Molecular Biology. 48 (3): 443–53. doi:10.1016/0022-2836(70)90057-4. PMID 5420325.
[5] 汉朝. “ 阿里云提供高效基因序列检索功能，助力冠状病毒序列快速分析 ”, https://developer.aliyun.com/article/753097?utm_content=g_1000111278
[6] Mikolov Tomas; et al. (2013). “Efficient Estimation of Word Representations in Vector Space”. arXiv:1301.3781
[7] 基因数据集 https://www.ncbi.nlm.nih.gov/genome/viruses/variation/help/flu-help-center/ftp/
[8] de Groot RJ Baker SC Baric RS et al. Middle East respiratory syndrome coronavirus (MERS-CoV): announcement of the Coronavirus Study Group. J Virol. 2013; 87: 7790-7792
[9] https://www.ncbi.nlm.nih.gov/nuccore/NC_045512.2
[10] https://www.ncbi.nlm.nih.gov/nuccore/NC_019843.3
[11] https://www.ncbi.nlm.nih.gov/nuccore/NC_038294.1
往期文献：
[1] 戴口罩也能刷门禁？疫情下 AnalyticDB 亮出社区管理的宝藏神器！https://developer.aliyun.com/article/745160
[2] 阿里云提供高效基因序列检索功能，助力冠状病毒序列快速分析
https://developer.aliyun.com/article/753097