Haplotype-resolved assembly of a tetraploid potato genome

编辑 / 发布于2024-01-26 / 更新于2024-03-02 / 阅读 27

Haplotype-resolved assembly of a tetraploid potato genome using long reads and low-depth offspring data

Background

基于参考序列的单倍型组装方法 HapTreeand H-PoP

本文使用 hifiasm 组装hifi read 获得一个最初的组装,然后使用193个两个马铃薯品种的后代的低coverage测序数据和一个新的基于k-mer的方法识别四种单倍型,从获得的组装图上组装单倍型。

Results

Overall assembly strategy

  • a. Altus 测序的 hifi read,和193 个 Altus × Colomba杂交的Illumina测序数据

  • b. 用hifiasm组装pacbio hifi read 组装成一个组装图,这样就产生一个解决了部分单倍型的图,其中气泡代表不同的单倍型。对于组装图中的每个的unitig ,检测独特的k-mers。

  • c. 把hifi read 比对到contig上,给每个contig估计dosage(1-4,有几个单倍型是这个contig)

  • d. 用illumina reads统计193个后代上独特的k-mer数量,每个unitig由包括193个值的k-mer计数表示,具有相似计数模式的节点,意味着节点由相同的后代样本子集继承,因此很可能是同一单倍型的一部分

  • e. 对于组装图组件中的所有节点,计算 k 聚体计数模式的成对相关性,并对组件进行聚类以表示染色体

  • f. 在每个染色体簇中,估计剂量为 1 的节点首先被聚类为四个单倍型,同样基于成对相关性。

  • g. 剂量 > 1 的重叠群被添加到包含 k 聚体计数模式相关性方面最匹配节点的簇中

  • h. 这个过程产生了包含每个单倍型子簇的染色体簇

Initial assembly


每个单倍型进行24×PacBio HiFi 测序,平均coverage = 24x4=96。对103个Altus × Colomba杂交后代 进行 Illumina short-read测序,获得2 × 150 bp paired-end reads ,测序数据中包括两个altus的单倍型,两个colomba的单倍型。平均每个单倍型coverage 为1.5×,平均coverage 为 6 。

用hifiasm组装hifi reads输出了一个包含所有组装的,和未处理的unitigs,部分解决了4个单倍型,变异以泡泡结构在图上表示,既一个unitig分支两个而或多个unitigs

Dosage analysis

dosage (单倍型数量),可以是{1,2,3,4}中的任意值(4倍体)。通过分析比对到unitigs上的reads覆盖深度获得(analyzing the coverage of reads aligned to the unitigs)。

由于hifiasm的图包含overlap,只计算每个节点不与其相邻节点重叠的区域,仅计算这些独特区域的深度作为平均深度。

b图就是node长度超过100k的覆盖深度图,峰值23,46,69,分别表示dosage 1,2,3。92几乎不可见,可能表示只有很少一部分是纯和区域,并且完全不存在超过100kb的长纯和延伸。

Analysis of k-mers

统计unitigs中所有可能的 k-mers (k=71) 的数量,确定那些在整个组装图中只出现一次的k-mer,另外丢掉在 Colomba genome 中找到的那些,保证这些 k-mers对 altus 是独一无二的,把这些 k-mer 记录为 unique k-mer,基于unique k-mer 的方法促进了 二倍体后代三倍的单倍型组装和人类88号染色体的挑战性区域的组装。第二个例子中,作者制造了一个独一无二的核苷酸k-mer库,对 long read 编码并组装成 scaffolds。对于每个unitig,使用相应的的 unique k-mer 集合作为node标志,k-mer计数基于Jellyfish。

在后代193个样本中记录 unique k-mer的数量,每个四倍体后代从altus和colomba中分别继承了2个单倍型,通过评估unique k-mer的数量推理一个unitig中继承的copy的数量。对于亲本中dosage 1 ,代表unitig的k-mer在后代基因组可能不存在或存在。对于dosage 2,k-mer可能在后代中不存在,存在一次或存在两次,对于亲本 dosage 3 ,后代必须存在1或者2,对于dosage 4,两个遗传的单倍型都必须源自该unitig。

基于上述,可以将带有unique k-mer的标记为片段有效(phase informative),没有的为片段无效,图c表示node长度分析,一般无效片段都是更短的那些。

Correlation analysis


利用来自同一单倍型的unitigs有相似度k-mer计数,因为相应的单倍型上下文被遗传到后代样本的相同子集。使用 a. 具有有效信息,b. hifi数据和ont数据中有相同dosage估计 (将read比对到unitig上如果hifi数据中dosage估计是2,ont数据要求是0,1或者2)的unitig 进行分析。在候选unitig 的k-mer计数模式之间计算 Spearman correlation coefficients, 距离相近的基因组位点很可能一起传递给后代,因此由于重组,后代的单倍型信号随着距离的增加而变弱。下图的结果与预期一致。距离小于10Mb的contig之间有强烈的相关性,并随着距离的增加相关性降低。

遗传过程中交叉互换,一个单倍型中相近位置的contig在遗传到后代中更有可能遗传到同一个单倍型中。所以相近位置的contig具有更高的相关性。

基于高相关性值,可以重新构建出在初始组装中没有连接的区域,像气泡结构和未连接的片段。ch3的重构表示如上图b所示,最初的组装,染色体由3个相连组件和2个更长未连接的contig组成。,通过连接高正相关的contig对,可以重构出染色体的片段结构,并相应的对组件和重叠群排序。

Graph traversal and final haplotype assembly


Methods

Data production

Dosage estimation in unitigs

计算每个contig i 去掉和其他contig有overlap部分序列的平均coverage c_i 整体平均coverage(每个单倍型的coverage or 四倍体的coverage/4)是 m,对每个contig估计dosage

d_i = d\ for\ (d-0.5)m <c_i\le(d+0.5)m,d\in\{1,2,3,4\}

c_i > 4.5m 设置 d_i=5 标识这是一个重复的contig

Connection of graph components

将unitigs聚合到单倍型-resolved染色体簇包括两个步骤。首先,尝试在染色体层次解决,染色体可能有多个相连的组件以及额外的单体,确定图中的那些组件属于应该在一起。然后将每个染色体cluster 分成4个不同的clusters,代表4个单倍型。

利用先前计算的k-mer计数(某个node(片段)上的unique k-mer 在193个后代测序数据中的k-mer计数),根据相似k-mer计数将unitigs分组。分组依据同一个单倍型的unitig具有高度相似的k-mer计数模式,而相反模式更有可能来自同一染色体的不同单倍型,模式不相关的unitig可能来自不同的染色体。
两个nodes之间的k-mer计数模式相似性根据seperman correlation coefficient p计算。高正相关代表两个nodes来自同一个单倍型,负相关代表node在不同的单倍型中,只有来自相同染色体会有高的相关性(且距离相近的contig),在两个不同染色体上node,k-mer计数应该是不相关的。

根据最高成对相关性(p>0.5)对所有节点分组,将聚合组件和单个unitigs聚集到染色体cluster里,当发现其中的contigs源自相同的图组件时,合并这些clusters。第一个聚类步骤产生了12个大簇(土豆是4倍体,单倍有12条染色体),定义为相关的染色体cluster或伪染色体。

Clustering of unitigs based on similar k-mer patterns

使用先前计算的doage估计,从dosage1 开始,确定每个染色体的单倍型,因为这些node只能被分到一个组里。首先构建最高相关性种子,然后合并成更大的组添加更多节点。

对每个dosage 为1 的node,创建一个cluster,只包含和它相关性>0.5的unitigs,生成种子cluster,然后根据这些种子cluster之间包含的公共节点合并这些cluster。为了分开不同的组,把高负相关的两个系欸但表示不同的单倍型。

在cluster c_i,c_j 之间如果存在 node pair n_i,n_j\ where\ n_i\in c_i\ and\ n_j \in c_j具有高负相关性( p_{i,j} < -0.3 ) 则构建负边,相反如果存在node pair 其中 p_{i,j} >0.5 则构建正边。如果不存在矛盾边,例如一个从 c_i 到另一个 c_k有正边,但 c_kc_j 之间连接有负边,则通过正边连接的两个 c_ic_j cluster可以合并。合并后,和其他节点有高相关性的节点都被分配到cluster里。对于剩余的node,如果3个最好的匹配节点都在 c_i 中,就把它分配到 c_i ,否则就不把它分配到任何一个簇中。

最终,分配具有更高dosage的unitigs到先前计算的单倍型clusters。分配一个dosage x(x\in\{2,3,4\}) 的cluster n ,计算 n 和其他所有cluster中所有nodes的相关性,把node添加到与 n 正相关节点的比率最高的 x个(2,3,4)个cluster

这是分配是保守的,可能有contig未被添加,采取后处理,重新把剩余的节点分配给具有最高相关性的cluster。

Assembly of clustered unitigs

我们从每个染色体的cluster开始,重新构建染色体的contig顺序,去找出它们之间所有的联系,创建更大更连续的单倍型。

如果 phased node 在任一方向只有一个邻居,该邻居也是 phased 的,对于简单气泡结构(源节点,两个分支,汇合节点)。如果源节点,汇合节点是phased的,则两个分支节点中的一个在phased路径上,如果两个分支都不是phased,说明不存在可用的信息选正确的,任选一个,相应的序列被填充为占位符而不是节点序列,表示缺乏正确的单倍型序列信息。

然后考虑剩余部分没有被连接的所有phased节点。这些形成了一个线性的块状结构,对于每个块,我们能够识别两个终端节点,并通过块重新创建节点路径(以及序列,块内部的序列)。为了重建这些单倍型块的顺序,我们随后搜索了仅包含unphased节点的不同块的末端节点之间的路径。对于可以唯一的连接到一个额外块的块,我们将两个块序列连接起来,并再次使用占位符字符来表示中间unphased片段的长度。最后,我们解决了扩展节点路径之间剩余的任何重叠,从而产生了最终的输出序列。