转载自yizhuo
基因、测序、分析
基因,生命的基本因素,是人类和其他生物的基础遗传物质。人有 23 对染色体,总共记录了大约 3Gb 个碱基(这里的 b 是 base,即碱基,可不是 bit,参考这里),每个位置上的碱基可能是 ATCG 中的一个。简单理解起来,就是有了这 3Gb 长的字符串,就能克隆一个你。基因测序,就是用化学和物理的方法,把你身体里这 3Gb 字符串检测出来。
当然,由于受当前测序技术的限制,我们并不能一次性测得一个完整的 3Gb 字符串,而是无数个 150bp 左右长度的小碎片。把这无数个小碎片重新组合还原成 3Gb 的长字符串的过程,叫全基因组组装。人类基因组计划干的就是这个组装拼图的事情,到了 2003 年,基本上算是拼完了。于是就有了一个标准的 3Gb 人类字符串,业界称其为『人类基因组参考序列』,也常被简称为『参考序列』。
基因测序分析,就是检测每个个体身上的这 3Gb 字符串,然后再跟标准的字符串(参考序列)做比对,来看是不是哪里发生了变异,是否有已知的遗传疾病风险。
现在随着基因测序技术的成熟,成本在飞速下降,所谓『旧时王谢堂前燕,飞入寻常百姓家』简直就是指日可待的事情。
于是越来越多的基因数据需要计算分析。
原来的分析计算流程
测序出来的数据是一大堆长度 150bp 的小碎片,但由于已经有了完整的人类参考序列的拼图,那么在这个拼图上寻找位置要比还原拼图的过程容易很多。多年来,从输入碎片数据比对到标准的 3Gb 参考序列,再到变异的检测,形成了如 BWA、Samtools/Picard、GATK 等为单机运行优化良好的业界公认软件。
将这些软件组合使用起来可以形成一个基因组数据分析的流程。华大基因基于阿里云 ECS 的 BGI Online 平台、七桥在 Amazon 和 Google 云上架构的 IGOR 分析平台,都有这样的业界金标准流程(下图为七桥,1.7G 左右的输入数据跑了 1 小时 48 分)。
但是上面的分析只是外显子测序分析(WES),所谓外显子就是那些会直接影响我们外在性状和健康与否的序列,它的长度只占了全基因组总长的 1%。因此需要测序的数据量比全基因组测序分析(WGS)少得多。使用 WES 的原因和意义,一方面是因为我们对这 3Gb 的人类基因组所蕴含的全面信息所知甚少,这 1% 的序列是我们目前能够很有效进行解读的位点,同时,也是由于它直接影响着我们的性状和健康,所以在很多疾病的研究中,通常只检测和分析这些位点和区域;另一方面也是因为全基因组测序目前成本相对较高,数据量巨大,单机计算时间基本要几天。
基于 MaxCompute 的方案
数据表达
来看基因数据的几个重要表达。
测序出来的原始输入数据 FastQ 格式(截取了两个 record),原本就是四行一个 record 的文本:
@ERR194147.1 HSQ1004:134:C0D8DACXX:1:1104:3874:86238/1
GGTTCCTACTTCAGGGTCATAAAGCCTAAATAGCCCACACGTTCCCCTTAAATAAGACATCACGATGGATCACAGGTCTATCACCCTATTAACCACTCACG
+
CC@FFFFFHHHHHJJJFHIIJJJJJJIHJIIJJJJJJJJIIGIJJIJJJIJJJIJIJJJJJJJJJJIJHHHHFFFDEEEEEEEEDDDCDDEEDDDDDDDDD
@ERR194147.2 HSQ1004:134:C0D8DACXX:2:2104:2852:75174/1
ACTTCAGGGTCATAAAGCCTAAATAGCCCACACGTTCCCCTTAAATAAGACATCACGATGGATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTC
+
@BBFDFFFHFFHHIHIJJJJFIHHFHFHJCIHFHIJJJJJJJIJIJIJJIIHIJJJJJJJBEGIGHIHGHHHEFCDFFEDEEDEEDDD?CCCDDDDDDDDC
测序产生的中间及最终结果的 VCF 格式(截取两个 record),也是一个多行表头,然后每行一个数据 record 的形式。
##fileformat=VCFv4.1
##...
##other headers
#CHROM POS ID REF ALT QUAL FILTER INFO
1 10019 rs376643643 TA T . . OTHERKG;R5;RS=376643643;RSPOS=10020;SAO=0;SSR=0;VC=DIV;VP=0x050000020001000002000200;WGT=1;dbSNPBuildID=138
1 10054 rs373328635 CAA C,CA . . NOC;OTHERKG;R5;RS=373328635;RSPOS=10055;SAO=0;SSR=0;VC=DIV;VP=0x050000020001000002000210;WGT=1;dbSNPBuildID=138
用表来表达描述完全没问题。当然,PairEnd 的 FastQ 总是成对出现,那么就编排成 8 列的表就好了。VCF 格式的多行表头也要稍微处理一下,不过这都不是事儿。
问题可分割
这块涉及到生物学的专业知识,要描述清楚比较困难。基本的想法,是将 3Gb 标准数据切成若干份,然后分布式执行每一份的计算。理论上,切割肯定会对计算的精度带来影响。但是影响多大要试过才知道。实践中,也有些办法来弥补这些影响:比如适当扩大 reference 的 interval,或者通过数据重叠来缓解。
简言之,跟单机计算结果对比,一致率 98.8%,这个数字可以接受。
MaxCompute Streaming 作业
如果数据如何在分布式系统里流动都想得很清楚的话,那么最大的障碍其实已经不在,剩下的就都是工程问题。
MaxCompute 支持 Streaming 作业,允许用户使用任何语言开发 MR 作业,开启隔离环境之后甚至可以用在公共云环境。用 Streaming 作业,可以非常方便的将 BWA、samtools、GATK 集成在一个 MR 作业当中:
这个作业,运行时间不到 3 个小时。
结语
或许很快,癌症患者可以当天拿到靶向药物是否有疗效的根据;新生病儿的早期诊断当天就可以排除是否遗传疾病;基因测序分析就像今天验血一样成为普通人司空见惯的体检手段。
或许将来,MaxCompute 可以更加发挥计算的能力,辅助挖掘 DNA 上目前 99% 还未知的领域。
我们期待这一天早些到来。
欢迎加入“数加·MaxCompute购买咨询”钉钉群(群号: 11782920)进行咨询,群二维码如下: