【分享】Nature一作深度解读:人类基因组单倍体组装重要技术细节(三—七)


【编者的话】续上篇。5月25日,《Nature Biotechnology》 发表了华大基因炎黄(YH)基因组单倍体组装的结果(De novo assembly of a haplotype-resolved human genome),这是迄今最好的人类单倍体组装结果,本文同样由共同一作所写,纯技术活还在继续。
 
这两天越来越多的同事、同学在询问文章中的诸多细节,一一回复起来着实挺累的,这迫使我再次感到把剩下的关键细节点写出来的重要性和紧迫性。
在上一篇博文中我们说到了fosmid clone的定位和每个fosmid pool中含有overlap fosmid clone的概率的计算问题。接下来我将重点说一下在phasing之前我们是如何预先对变异进行纠错以及如何评估最后的phasing结果的,这也是被问得最多的地方,另外我也想说一下我们是如何通过模拟出来的fosmid区域计算出理想状态下,使用fosmid 测序策略所能达到的phasing期望长度等其它几个重要的点。
那么以下就进入正文吧。

3. Phasing之前我们是如何对变异进行纠错的?

这里我想先指明一个重要的点,在获得了杂合SNP集之后(关于SNP集合的获取原文有详细说明,即通过合并WGS reads和fosmid clone reads各自检测到的结果所得到),还必须要明确地知道每一个fosmid在对应的这个杂合SNP位点中具体是哪一个碱基才能进行下一步的phasing。举个具体的例子,比如某一个位置上是A/T杂合,那你还必须知道FAH在这个位置上究竟是A还是T!

那关于这一点我们是如何办到的呢?原文章在这里没有作太多详细的描述。其实,我们采用了一个简单而直接的办法——fosmid组装出来的序列(FAH,见原文注释)在这个位置上是啥碱基我们就认为是啥!如果刚好和该位置上杂合SNP的两个碱基都不同的话则直接忽略fosmid组装序列在这一个位置上的变异。这样做的合理性(或者说是前提)是建立在所组装的序列有足够高的单碱基准确性的基础之上,但YH基因组的这个数据确实有这个自信,因为为了能够顺利进行组装,我们为每个fosmid clone测了将近50x,这是很夸张的数据量!!因此这个准确度的基础确实是有的。

不过,我们只为整个YH基因组构建了~8x的fosmid clone(注意不要和上面的50x混淆了,50x指的是对每个fosmid clone的测序深度),也就是说fosmid组装序列对整个YH基因组的覆盖是7x~8x,我们最终的phasing对象也是这些fosmid组装序列。对于这个深度的fosmid而言,即便每个fosmid内部的序列准确度能够得到保证,但还是太脆弱了,特别是对于那些fosmid与相邻fosmid之间相互overlap的杂合位点不够多的情况就更是如此了。而这些overlap的地方正是真正能影响phasing准确性的核心区域,这些地方就犹如大陆板块之间的结合地带一样,对出错碱基的响应是非常敏感的。对于这么低的fosmid clone深度,这个时候哪怕只有一个overlap杂合位点上的碱基call得不对,都会对其phasing的结果带来很大的影响!因此我们需要进行二次纠错。

刚刚也说了,fosmid组装序列只有7x~8x,那其实意味着每个单倍体只有3x~4x的fosmid覆盖,因此直接通过简单地计算overlap位点上占多数的碱基来判断的策略就行不通了!只能另寻其他更合适的办法,后来我们注意到了一个地方,同个fosmid上多个不同的杂合位点同时出错的概率是很低的!因此既然不能通过单个位点上的碱基覆盖情况来纠错,那么就通过与附近的其他杂合碱基来共同判别,这就是原文章的Methods中提到的linkage information!具体过程是这样的:(1)在每一个fosmid中,我们逐一计算每个杂合位点与相邻的3个其他杂合位点之间两两相连的碱基配对信息(这里简称“配对基”);(2)通过相同两个位点之间配对基的多寡来判断对应位点上的杂合碱基是否存在错误。因为,一个很明显的事实是,如果fosmid上某个位置的杂合碱基是错的,那么它的配对基一定都是少数派,也就是说conflict linkage information必是低频!这个时候,我们就只需要将其纠正为另一个杂合碱基就行了,当然必须保证纠正后它的配对基变成高频了!举个例子来说明吧,比如fosmid上有两个杂合位点 A/T和C/G, 在我们遍历了所有的fosmid之后发现这两个位点A-C/T-G的配对基连接支持数是最高的,那么当我们发现有一个fosmid在这两个位点的配对基是A-C/C-G之后,就能马上发现C-G连接中的C将是错的,而应是T,从而将其纠正为T-G连接!这看起来可能有点绕,但其实就是参考了基因组组装中使用k-mer频率对组装进行纠错的策略。
 
4. 在Phasing前的硬链接.
 
硬连接,就是直接连接,不使用任何推断过程。这是一个文章中没有明确提到的步骤,我想也将会是最易被忘掉的步骤。之所以有这一步其实是因为fosmid clone的深度不够高,如果把所有区域都一股脑儿扔去做phasing推断是不合适的。一开始我们就是这样操作,结果却发现有一部分fosmid clone之间的连接本来是非常明显而直接的,即相互之间overlap的位点数量很多,一致性也非常高,它们连接在一起应该是绑在板上的事情。但是在使用RefHap之后却发现被分到了另一边,后来分析这一情况的时候觉得低深度应该是导致这个情形的主要因素。因此后来在用RefHap phasing之前,我们都对那些能明确合并在一起的fosmid clone进行了硬连接!这个过程也是比较简单的,但为了尽可能确保正确性,这个硬连接其实进行得很保守,相邻两个fosmid clone区域要被直接合并在一起必须同时满足:overlap的杂合位点个数超过5个,并且一致性率超过80%才能合并。

5. RefHap算法.

RefHap的算法我不打算在这里写下来,原因有3个:(I)篇幅和主题所限;(II)RefHap并非我们的原创;(III)RefHap自己的文章已经说的很明白了。所以我这里仅仅把它的大意梳列在这里:
(1)获得所有存在overlap关系的fragments集合;
(2)构建二分图,计算集合中所有fragments两两之间的打分值,并取出分值最大的两条fragment,将其分为两类,一类认为是haplotype I(记为I集),另一类为haplotype II(记为II集),这是RefHap中最具重要性的单体型初始划分办法;
(3)取出任意一条不属于I和II集的fragment,分别计算其在这两个集中的分值,按打分值分类,分到分值大的那一类。
注意,这里的fregment就是指单个的fosmid clone区域以及区域上的杂合变异。

6. Phasing之后如何对haplotype的结果进行评价的?

一般情况下,对phasing结果进行评估都是需要有一个由实验验证过的单体型序列来作为金标准的,而且这个序列还必须要有足够的长度,如GIAB(Genome in a Bottle Consortium)项目中提到的那样。Sanger测序序列在我们这个情形下是不适合的。后来我们也发现fosmid测序的结果是可以直接作为金标准去评估单倍体的phasing结果的,而且也已经有人这么做了。而我们本来就是大量的fosmid测序,并且质量也是够高的,这也就意味着我们已经没必要再去做其他的测序了,直接使用自身的fosmid来评估就行了!

我们的具体做法是这样的:从这~8x的fosmid clone中随机挑选出1x,并直接将其视为“金标准”。然后将最后phasing的结果与这个金标准比较,以此来评估phasing的准确性。文章也是这么说的 —— “We randomly selected fosmids that represent ~1× physical coverage of YHref and calculated the consistency of the marker’s linkage information between the selected fosmids and the HDG phasing”,不过,我觉得一定会有许多人仍然觉得糊涂,而且很难明白这里所说的“linkage information”到底指的是什么。下面就让我来解释一下吧。

这里的linkage information和上文“在Phasing前是如何对变异进行纠错的”中所提到的是一样的概念,也都是指杂合位点之间两两连接的碱基对(配对基)。以下是详细的说明和解释:

我们在取出了1x fosmid作为金标准之后,考虑到真正能够影响phasing连接准确性的地方,其实只有fosmid之间相互overlap的那些区域!每个fosmid内部的结果是根本不需要进行评价的,因为它们不会影响phasing连接的准确度,最多就是影响phasing之后碱基序列的准确性。因此,在评估的时候,我们选择了fosmid之间相互overlap并且能够被“金标准”fosmid完全覆盖到的区域作为评估对象,其他的区域则一概忽略。接着我们选择overlap区域中间的一个杂合位点作为分割点,将phasing结果分为左右两边,然后把左边和右边的杂合位点上的碱基分别两两配对(配对基),这些配对基的连接就是lingkage information。这样来进行左右划分的意义是为了让左边和右边的都有来自不同fosmid的杂合位点,这就有利于避免评估的时候只使用了落在同一个fosmid上中的位点,从而也就尽可能地避免了bias。最后,就是通过比较这些两两连接的碱基对和“金标准”fosmid上相同位置的碱基对的一致性率来评估phasing结果的准确性。

7. 我们是如何模拟fosmid区域并计算出phasing的期望长度的?
 
这也是我在这篇文章中最后一个想说的问题了。很多时候我们可能会以为fosmid clone测得越多,那么基因组phasing的长度就一定会越大,但实际上,这是错的!它总有一个顶部,决定的因素是测序序列的读长和物种基因组本身杂合marker(比如杂合SNP)的密度(见下图)。YH(中国人)基因组与其它的几个人种相比有着更多大范围的纯合区域,这导致了很多区域的杂合SNP密度要更低,因此即便是在相同的fosmid clone深度下,YH所能达到的phasing长度也将是最小的,而且对于YH基因组利用fosmid的策略最高能phasing的N50长度也只有1Mbp左右,现代人的祖先——非洲人(NA18507)则在各个深度上都是最大的,但关于这个问题我不打算展开来说了,以免跑偏了主题。我还是说一下我们是如何通过模拟fosmid 区域从而计算出了在不同的fosmid clone深度下,YH基因组所能达到的phasing长度的。
QQ20150603-4@2x.png

我们的fosmid区域模拟方法比较简单,具体是这样的:(1)规定每个fosmid区域的长度都统一为36.8kbp(这个值是我们的fosmid clone平均长度);(2)计算特定深度下应该模拟的fosmid区域总数N;(3)在整个基因组上随机选取N个起点,获得起点之后一律往后延伸36.8kbp,ok,这就获得了总共的这N个fosmid区域了。最后再与杂合的SNP集进行Overlap,确定出每个区域中都有哪些杂合位点,然后再把含有相同杂合位点的fosmid区域按顺序逐个合并直到不能找到下一个相互overlap的区域为止(注意这里没有任何纠错的必要),这样便可以得到特定fosmid clone深度下期望的phasing长度。

好了,这篇单倍体组装文章中的主要细节点就写到这了,如无必要,接下来我也暂不打算在这方面继续做更多的拓展了,虽然写了很多,但恐怕还是没能真正描述清楚,欢迎一起留言讨论。最后容我再show一下原文章中的主流程图,这个图也是修改了几十遍,总得show多几下 :)
 
QQ20150603-3@2x.png

De novo assembly of a haplotype-resolved human genome. Nat. Biotechnol. (2015). doi:10.1038/nbt.3200


1 个评论

后面几个微信上的图挂了呢。

要回复文章请先登录注册