没经费也能干大项目-50块钱测序一个人的基因组

专题介绍：数据分析那些事儿是瑞典农业科学大学的一名科研工作者。在生物大数据分析领域已经有七年的科研经验。希望通过分享过去7年在国内外硕士、博士、博士后申请、求学、科研、求职、统计学、R&Python编程、生物大数据分析、国内外高校企业求职生命科学发展对我们当下社会的影响中收到的个人启发来让1000个人受益让求学求职生活变得容易一点点。如果你也想持续学习，克服懒惰的话。关注他参与“100天生信/数据科学自我挑战”的话题吧。

建库（10 RMB） + 测序到0.5X + 填充到全基因组密度（大于95%准确率）= 重测序人基因组成本低于50RMB。这也就是说，如果研究的物种已经被测序，不管是家系还是野生资源，一个1000个体（基因组大小3Gb以内）的重测序课题预算，只有5万RMB。本来以为这个事儿，打开了新世界的大门，结果进了门，发现院儿里已经有好多人。可能对于太美好的事儿，大家都警惕性比较高，站在门口使劲儿喊进来的人也有限。所以我写下这篇推文，再喊一嗓子。时间有限，可能漏掉了很多前辈的工作，请多多包涵。

大概17年底的时候，导师开会认识了一位普林斯顿研究果蝇的AP（助理教授），这位AP说可以用1美元建立全基因组重测序文库，并且可以无偿帮我们构建4000只鸡的文库。这些文库测序到0.5X，就可以填充到全基因组重测序的SNP密度，每个样本建库加测序大概5美元左右。

当时我们建库比较省钱的办法是订购illumina的Nextera试剂盒，用1/10的体系去做，这样可以把50个容量的试剂盒变成500个，建库成本是大约50RMB一个。当时瑞典公司商业建库报价是大概一个600RMB。这里除了建库，通过低深度测序（0.5X）后做基因型填充省下的钱，也非常可观。说到基因填充大家都不陌生，不管是从低丰度（2-4X）填充到重测序水平，还是从低密度芯片填充到重测序水平都有很多文章和例子。低到1X以下测序，2010年最早水稻的重组自交系有过类似的报道，但是自交系毕竟是inbred，只有一个单倍型，而且重组不多。对于野生2倍体生物，尤其是野生群体，每一个基因组上具体的区间上有很多个单倍型，比自交系高很多，测序到0.5X，意味着一个个体两个染色体（haplotype）都没有被同时捕捉到的情况下，填充到0.95的精度，乍一看还是非常吃惊的。

太美好的事情，来的又这么突然，本着绝对有诈的想法，课题组开始了自己摸索的过程。

一、10 RMB建库

2014年的时候，有一位卡洛琳斯卡医学院的博士生Picelli，优化了一个转座酶建库的体系，这个转座酶叫Tn5。这个转座酶是可以携带外源的一段DNA接头(图一A，蓝色和绿色部分)，在切碎长链基因组的同时把接头插入到切碎的片段中去(图一B，灰色为基因组，蓝色和绿色部分为加上的接头)。通过设计合适的barcode+接头互补序列(图一C，蓝色和绿色部分为加上的接头，橘色和红色为文库末端与flowcell cluster的序列，紫色为barcode)，我们可以通过PCR扩增切割产物的方式，得到全基因组文库。这个就是illumina nextera 建库的原理。而Picelli 的纯化了Tn5，（载体可以在Addgene订购，https://www.addgene.org/60240/），并且摸索了buffer的配方，并且把这一系列过程公开发表了。这个protocol Illumina有专利，所以不可以商业化，但是可以自己实验室用。

我们课题组仅仅是重复了Picelli的过程，把它搬到了96孔板上，做了一些微不足道的优化，使得通量上有所提高。当时峰值的时候与G师姐一天之内，两个人建了1000个库，而且QC几乎都合格。当时觉得如果一个样品500RMB，我们不久就要登上人生巅峰。完整的protocol在这里（https://www.protocols.io/view/a-high-throughput-cost-efficient-library-preparati-rt8d6rw）。这个protocol价格的关键是自己纯化Tn5，我们当时找了卡洛琳斯卡医学院的蛋白平台，一次纯化5000RMB得到了20ml酶,足够做几万个样本，但是最近听说国内1500RMB才16ul。纯化酶，是个很基础的工作，我没有亲自做过这个酶，就不误导大家了，但是相信大家都或多或少有同事熟悉。今天这个流程已经被很多同行继续优化，我知道的看中国农业大学的一篇文章，https://www.biorxiv.org/content/10.1101/754671v2

图一 Tn5 建库流程

二、低丰度测序（0.5X）基因型填充

2-5X测序深度的基因型填充，有很多很多报道。尽管如此，刚听说0.5X可以填充时候，还是觉得不可思议。

1. Bi-parental cross-F2 + RIL

2010年，水稻里发表了一篇文章，在RIL 测序到0.03X，把基因组分割成没有重组的区段之后（Bin），每个Bin来自于哪个亲本可以通过找亲本特异性的SNP求均值而得到。原理就是，在一个俩个亲本建立的杂交群体中，重组非常有限，具体到某一个基因组上的区段上只有两个haplotype，而且序列还已知（有亲本的重测序数据），由于重组有限，基因组上一个个没有重组的独立区段都很大很大，这样就可以通过挑选能够区分开两个haplotype的informative SNP然后求均值，来得到这个区段的基因型。有了这个Bin的基因型，就可以去做QTL定位。文章里面有Perl语言里面写的代码。类似的在玉米里也有很多很多文章，欲知细节请关注各位玉米界的中外大佬们。

2. Multi-parental cross

由于我们课题组研究的群体是多个outbred的parent所构建的一个20代的AIL，具体到基因组上一个区段上，有很多个haplotype，20个世代，重组也相对多了一点。但是在一个outbred的家系里，每个nuclear family都是一个Bi-parental cross，从上一代到下一代重组有限，因此尽管这个家系里包括20个世代，很多个F0亲本，我们仍然可以通过在每一个世代中iterate nuclear family的方法，把它变成一个Bi-parental cross，从而用同样的方法得到每个bin的基因型。这里我做了一个R版本的包https://github.com/yanjunzan/GBSA。

但是这个case里有一个棘手问题，由于群体世代很高，所以我们必须要能够精确的定位到recombination break point。就是说希望第一个bin和第二个bin之间物理距离很短，否则这里会出现一个gap，F2的gap会变成F3的gap，越往高世代，gap随着世代的积累会越来越大，可以impute的区段随之会越来越少。所以我们采用了马尔科夫链，来尽可能的缩小Gap，文章在这里https://link.springer.com/article/10.1186/s12711-019-0487-1, 软件在这里https://github.com/yanjunzan/Stripes，文章看不看没关系，结论就是在任何一个有系谱的复杂家系里，我们可以通过低丰度测序来完美填充到重测序水平。

插播一个前人在MAGIC群体里低丰度基因型填充的例子，http://mtweb.cs.ucl.ac.uk/mus/www/19genomes/magic.html，这里用了类似的逻辑，可以在没有系谱的条件下，impute MAGIC Design群体。

3. 野生群体里0.5X的基因型填充

读了上面的例子，我相信大家看出来规律了，这里几个点特别重要，一、被填充的群体，具体到某一个位置haplotype diversity多高？二、没有重组的区段有多大，也就是LD衰减有多快? 三、关于1和2 我们知道多少。

在野生群体里，我们测序到0.5X时候，对1和2 几乎不知道。但是我们可以去根据我们测到的结果去推测，只要我们测得个体足够多（目前看到的是依据群体LD结构和群体历史，500-1000精度已经很高了），我们可以用每个个体的reads合起来去做某一个小区段的de novo assembly,一旦assembly成功，我们就大致知道了1。我们还可以用个体的0.5X测序结果和assembly结果去无限逼近1的真实情况. 知道了1我们就可以用关键的informative SNP和马尔科夫链去推测个体的基因型。这个NG的文章在这里https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4966640/pdf/emss-68776.pdf，软件在这里https://github.com/rwdavies/STITCH。亲测在鸡里有效，据我所知道的在人、猪(近期中国农业大学一篇文章https://www.biorxiv.org/content/10.1101/754671v2)、狐狸里完美重复，目前应该可以说只要基因组不是特别复杂（比如针叶树20GB，几百万个contig），assembly不是稀烂，应该没有物种的限制。

至此可以说DNA重测序的成本已经非常低，技术难度也不是很大。这个建库方法也可以改良到RNA测序以及一些研究转录活性的DNA测序中。由于转座酶比较温和，不经过机械破摔，很少产生nick。未来可能会被拓展用于三代测序长链DNA文库制备。

虽然数据分析和存储也会有一定的花费，但是有了这个低成本测序方案，理论上一个中等项目就可以做很多东西，尤其是不依赖于表型的群体遗传学研究。但是没有一个高通的表型监测方案，需要耗费很大的人力去做数量遗传学课题。所以下一篇我打算写写如何用1000RMB解决植物表型组的DIY，来相对高通量拿到株型，花型，花数目，果夹，叶片侵染情况等通过RGB图像可以获取的表型，解放繁琐的人工统计过程。

未来故事，请扫码关注。拜托各位老师、同学和家长，把它传播给需要的师弟师妹。也欢迎各位师弟师妹，如果想知道或者了解什么方面，请留言互动。我会尽力和身边大佬同仁学习后，努力给大家一个全面的总结。