未来组又一篇Nature Communications || 高质量苹果基因组撩动你心
2019年4月2日,中国农业科学院果树研究所与武汉未来组等单位共同合作的研究成果“A high-quality apple genome assembly reveals the association of a retrotransposon and red fruit colour”发表于国际顶级期刊Nature Communications。中国农业科学院果树研究所张利义副研究员与未来组生物信息研发总监胡江为本文并列第一作者,中国农业科学院果树研究所丛佩华研究员与未来组创始人汪德鹏为本文共同通讯作者,武汉未来组的李净净,田玉等四人为本文共同作者。本研究以来源于花药的苹果纯合系HFTH1为材料,采用三代PacBio测序技术,结合光学图谱Bionano和Hi-C技术辅助组装,获得高质量苹果基因组图谱(contigN50为6.99Mb,2017年发表的苹果基因组contigN50为620kb)。有趣的是,进一步研究发现花青素生物合成的核心转录激活因子MdMYB1上游UTR区域的一个LTR反转录转座子的插入与果实红色相关联。本研究为探讨红色果实着色的分子机制提供了重要见解,同时突出了高质量的基因组组装在破解苹果重要农业性状中的实用性。
苹果基因组是苹果遗传研究和分子育种的基础,可推动苹果的可持续生产。尽管早前发布的高质量金冠(Golden Delicious)基因组使苹果育种研究取得了一些进展,但仅仅利用这些数据在发现新基因和描述基因组结构变异方面仍存在一些局限性。
在中国北方地区广泛栽培的寒富(Hanfu)品种是高度杂合二倍体,经过花药体外培养的纯合子植株HFTH1(三倍体)大大降低了杂合度(Fig.1);通过对相应Illumina reads进行SNP calling,发现与先前发表的双单倍体金冠(GDDH13)相比,HFTH1的纯合度更高(Fig.2),高度纯合的基因组有利于提升基因组组装质量。
Fig.1花药培养再生植株HFTH1表型与杂合子供体(HFP)基因型的显著变化
Fig.2 GDDH13及HFTH1基因组中的SNP数量
对~117× PacBio数据(平均读长13.1 kb)进行初步组装,获得了一个初步组装结果:基因组大小656.52 Mb,Contig N50 4.63 Mb;随后使用PacBio长reads及Illumina 数据对组装结果进行polish;应用~224× Bionano光学图谱数据将polish后的contigs组装成scaffold;随后,应用145× Hi-C数据进行精确聚类并排序,最终组装出一个高质量的苹果基因组,包含17条模拟染色体(Fig.3a,3b),基因组大小658.90 Mb,contig N50 of 6.99 Mb;此外,还将160,068 bp叶绿体基因组和396,939 bp线粒体基因组组装成完成图。
Fig.3 (a)17条染色体间的Hi-C互作图(100-kb resolution);(b)HFTH1基因组特征
本研究组装的HFTH1基因组在7条染色体的两端和8条染色体的一端共捕获22条长端粒序列(拷贝数从294到1073,Fig.3b)。BUSCO评估基因完整度达到~97.0%,覆盖公共数据库Malus domestica基因组中98.02%的ESTs。将组装结果进一步与水稻基因组(RGAP7)和拟南芥基因组(TAIR10)比较,发现HFTH1基因组BUSCO评估与两者都很接近,表明我们组装的基因组完整度达到模式植物的水平。
Fig.4 结构变异特征
转座子在物种适应性进化中扮演重要角色,而在苹果基因组中,超过75%的转座子都是LTR-RTs,为了追溯苹果基因组中LTR-RTs的动态演化历程,研究者在HFTH1基因组中鉴定了7,313个完整的LTR-RTs,平均读长7,868 bp(Fig.5a)。其中约62%的LTR-RTs在金冠基因组中同样存在(Fig.5b),说明这些LTR-RTs可能早已插入到寒富和金冠的共同祖先中。在蔷薇科基因组中,梨和苹果有共同相对保守的基因组,但由于梨基因组较片段化,在梨基因组中只发现了40个完整的LTR-RTs,进一步的分析表明,苹果基因组中的大多数LTR-RTs可能是在梨与苹果分化之后插入苹果基因组的,这一结果与这些LTR-RTs平均插入时间(0.8 MYA)显著低于苹果和梨估算的分化时间(8.1 MYA)相一致(Fig.5c)。比较HFTH1和GDDH13基因组,其中一半以上共有的LTR-RTs高度相似(identity≥0.99,Fig.5b),且共有的LTR-RTs及其侧翼区域的平均累积核苷酸替换(CNS)低于基因组的平均水平(Fig.5d)。
Fig.5 HFTH1基因组完整LTR反转录转座子演化
进一步对LTR-RTs进行分析,发现两侧有2个TSDs的特异LTR-RTs(占总LTR-RTs的31.9%)只在HFTH1基因组中发现(Fig.5b),在GDDH13基因组的对应位置却只有1个TSD,表明特异的LTR-RTs是在寒富和金冠分化之后插入的。这些特殊的插入事件,特别是在基因内或靠近基因区的,可能是在较强选择下固定下来的,因为蛋白编码基因占总基因组的22.22%,但只有12.05%的插入事件是在基因内发生。在基因区附近观察到的插入事件也比预期的少(Fig.5e)。而且,靠近这些插入位点的基因其平均表达水平也比总基因集低。但是,当插入在大于1kb的基因上时,则选择压力很弱或者丧失(Fig.5e)。进一步分析发现,与64.39%共有的LTR-RTs相比,82.54%的特异LTR-RTs至少在一个受测样本中表达,暗示大部分的特异LTR-RTs都年轻且富有活性。插入事件不仅影响附近基因的表达,而且增加了插入位点附近的突变率。本研究数据还发现随着与LTR-RTs距离的增加, CNSs逐渐减少(Fig.5d)。同时,研究者发现随着LTR-RTs分化时间的增长平均的CNSs也增加(Fig.5f),但是当转座子逐渐失活变成非功能性序列时,平均CNSs速率呈下降趋势并逐渐达到瓶颈(Fig.5g)。此外,研究者发现HFTH1基因组中2.9%的LTR-RTs在GDDH13基因组中可能已经被选择性清除。清除的LTR-RTs周围序列的平均CNSs在基因组中最高(Fig.5d),暗示在TEs演化过程中,清除事件可能比插入产生更大的影响。在HFTH1基因组中TEs的动态演化可能创造了独特的遗传和表型特征。
比较 HFTH1 和 GDDH13基因组MdMYB1序列,结果表明MdMYB1的编码序列是一致的,但是,在内含子区域发现了1个SNP,在上游区域发现15个SNPs和5个Indels。其中,GDDH13有1个501bp的插入,HFTH1有1个4097bp的插入,为gypsy-like LTR 反转录转座子,被称为redTE(Fig.6a)。尽管在HFTH1基因组中发现有3913个gypsy-like LTR 反转录转座子,但只有1个与redTE具有96.26%的相似性(被称为redTE-like),其它的相似度均低于75.10%。有意思的是,研究者发现redTE-like只存在于在HFTH1基因组中,在GDDH13基因组中没有发现,其两侧的LTR序列突变更多,这一结果暗示redTE-like比redTE更古老。
为研究上述的不同是否与苹果果实的颜色相关,研究者从数据库中调取这些有差别的序列进行比较。结果发现HFTH1基因组中的16个SNPs和2个Indels在非红皮品种也存在,而GDDH13的501bp插入在红皮品种中也被发现。随后,研究者研究了redTE与苹果红色的关系。对已知的112个红皮品种和非红皮品种进行redTE的PCR研究,结果发现在所有的红皮品种中均发现有redTE,而非红皮品种中均没有redTE(Fig.6b),暗示redTE插入可能与苹果的红果皮相关。接着,研究者筛选了75个非红皮品种HuaYue与红皮品种Honeycrisp的杂交后代(包括41个红皮品种和34个非红皮品种),同样证实了红皮表型与redTE的共分离(Fig.6c)。由此可见,redTE与苹果的红皮表型相关。后续研究者还通过Transient方法证明redTE可作为增强子,且redTE可调节苹果光响应基因MdMYB1的高表达。
Fig.6 苹果红色表型与一个LTR反转录转座子的关联性
由于苹果育种困难且周期长,研究者希望组装一个可用于指导苹果MAS(marker-assisted selection)育种的参考基因组。本研究中,研究者在基因组中发现了大量的结构变异,这将在苹果育种的分子标签开发中发挥重要作用。如,基于redTE的特定标记,是预选目标颜色苹果杂交苗特别有价值的工具,因为它比以前的标记更有效,更精确。这种标记可以通过消除大量的非靶标杂交苗而大大降低苹果育种成本。此外,来自于抗性亲本的HFTH1基因组将为苹果抗寒性和抗病性的标记开发奠定良好的基础,如苹果育种中的分枝环腐病和黑斑病。总之,这个近完成基因组图谱和其他基因组资源将有助于基因和功能标记的挖掘,并支持将研究成果转化为遗传改良,实现可持续的苹果生产。
参考文献
Zhang, L. et al. A high-quality apple genome assembly reveals the association of a retrotransposon and red fruit colour. Nature Communications 10, 1494 (2019).
论文链接
https://www.nature.com/articles/s41467-019-09518-x
未来组部分高分合作文章
海岛棉、陆地棉的基因组进阶之路
近期,海岛棉(Gossypium barbadense, AD2)和陆地棉(Gossypium hirsutum,AD1)基因组的文章继2018年之后再一次发表在Nature Genetics杂志上。这两种棉de novo基因组文章一波又一波,还能是Nature系列这样的顶级杂志,这是为什么呢?今天来介绍海岛棉、陆地棉基因组的进阶之路及其背后的科学故事~
Fig.1 海岛棉(左)与陆地棉(右)

Fig.2 棉花属系统发育和进化.不同棉花种基因组大小和核型(zhang ,Peng et al 2008.)
2012年8月,棉花重要种质雷蒙德氏棉(Gossypium raimondii)基因组问世[1],为棉花基因组学研究领域提供了有效的参考信息。雷蒙德氏棉的祖先被认为是海岛棉和陆地棉的D基因组供体。同年12月,又一篇棉花基因组文章发表在Nature杂志,揭示棉花基因组多倍化事件及棉花纤维的演化历程[2]。2014年,海岛棉和陆地棉的另一个供体A基因组亚洲棉 (Gossypium arboreum) 全基因组草图也被绘制出来[3]。棉花A组和D组基因组测序的完成, 为异源四倍体海岛棉和陆地棉的全基因组测序奠定了基础。
Fig.3 棉花纤维的演化
海岛棉和陆地棉是棉花的主要栽培品种,陆地棉产量高,适应性强,海岛棉产量低,栽培区域性强,仅能在新疆、埃及、美国亚利桑那州等少数干旱地区种植,但是其纤维品质比陆地棉优。基于高质量、高精度的基因组认识棉花的复杂性,进化的多样性,挖掘性状背后的分子机制培育出综合海、陆棉花高产、纤维优良、适应性强的新品种一直是棉花工作者努力的目标。因此,海岛棉和陆地棉基因组备受人们的关注。
Table 1 海岛棉与陆地棉基因组文章列表
2015年4月,两篇陆地棉基因组文章同时发表在Nature biotechnology期刊:
以中国农科院棉花研究所为主的科研人员对陆地棉遗传标准系—TM-1进行了全基因组测序和组装,最终获得了26条染色体[4]。利用比较基因组学方法,首次从全基因组水平揭示了四倍体棉花是由A基因组的祖先和D基因组的祖先通过染色体融合而形成,并初步揭示了四倍体棉花基因组的进化规律。
Fig.4 陆地棉的演化及基因组共线性分析[5]
另一篇文章同样以陆地棉为研究主体,利用Illumina测序平台,结合BAC末端序列和高密度的SNP遗传图谱,成功构建了高质量的陆地棉全基因组图谱[5];解析了四倍体棉花两个亚基因组的非对称进化机制;发现了MYB基因家族中的一个分支在纤维发育中起重要的作用,陆地棉中多个纤维素合酶等糖代谢基因在驯化过程中受到了显著的正选择作用,与棉纤维品质改良有直接关系。

Fig.5 异源四倍体棉花进化的示意图
2018年12月在Nature Genetics上发表的棉花基因组文章[7]又再次提升了海岛棉和陆地棉的基因组组装质量,该研究利用第三代测序技术(PacBioRS II),BioNano光学图谱技术和染色质高级结构捕获技术(Hi-C)进行联合组装。与以前发表的基因组草图相比,该项研究组装的基因组序列在连续性和完整性上有极大提高。通过比较分析,研究者还发现了棉花基因组多倍化后的复杂结构变异,同时定位了与优质纤维相关的基因座。
Fig.6 海岛棉和陆地棉基因组染色体特征
2019年3月,又一篇棉花基因组文章发表于Nature Genetics,其中应用llumina、10Xgenomics及Hi-C技术获得的海岛棉Hai7124及陆地棉TM-1基因组再次提升参考基因组的精度[8]。研究者指出,本次组装在在重复DNA富集的着丝粒区有明显的改善。全基因组比较分析显示,基因表达、结构变异和基因家族扩张是这些物种形成和进化的原因。这些发现有助于阐明棉花基因组的进化及其驯化历史。
Fig.7 海岛棉Hai7124及陆地棉TM-1基因组特征
参考文献:
[1] Wang, K. et al. The draft genome of a diploid cotton Gossypium raimondii. Nature Genetics 44, 1098 (2012).
[2]Paterson, A.H. et al. Repeated polyploidization of Gossypium genomes and the evolution of spinnable cotton fibres. Nature 492, 423 (2012).
[3] Li, F. etal. Genome sequence of the cultivated cotton Gossypium arboreum Nature Genetics 46,567(2014).
[4] Li, F. et al. Genome sequence of cultivated Upland cotton (Gossypium hirsutum TM-1) provides insights into genome evolution. Nature Biotechnology 33, 524 (2015).
[5] Zhang, T. et al. Sequencing of allotetraploid cotton (Gossypium hirsutum L. acc. TM-1) provides a resource for fiber improvement. Nature Biotechnology 33, 531 (2015).
[6] Liu X, Zhao B, Zheng HJ, et al. Gossypium barbadense genome sequence provides insight into the evolution of extra-long staple fiber and specialized metabolites. Sci Rep. 2015;5:14139. Published 2015 Sep 30. doi:10.1038/srep14139
[7] Wang, M. et al. Reference genome sequences of two cultivated allotetraploid cottons, Gossypium hirsutum and Gossypium barbadense. Nature Genetics 51, 224-229 (2019).
[8] Hu, Y. et al. Gossypium barbadense and Gossypium hirsutum genomes provide insights into the origin and evolution of allotetraploid cotton. Nature Genetics (2019).
Nature Genetics | 八倍体草莓基因组的起源和进化
种植园结出的草莓香甜美味,一直是广受世界各地消费者喜爱的水果,但是小伙伴们知道吗?和野生的四倍体东方草莓或者人工诱导的八倍体小黑麦不一样,栽培草莓中含有天然的八套异源染色体组!有意思的是,融合成为这八倍体基因组的四个二倍体物种,其中两个已经鉴定出来了,而另外两个至今仍然是谜,这就导致八倍体草莓形成的历史进化过程也扑朔迷离。
大家都知道蔷薇科是双子叶植物中的大类,其中的草莓属已有22个种为人们所熟知,因其倍性水平内部和之间的高度可杂交性,自然呈现出从二倍体乃至十倍体的多样性。相较于同源多倍体,异源多倍体更容易呈现多样化的农艺性状,且其中每个亚基因组都包含一个在单核内进化独立的遗传和表观遗传,但是当前研究对亚基因组优势的潜在机制和终极导向还不甚明朗。
由于大尺度的染色体改变和同源基因的交换导致了亲本染色体之间同源基因的乱序和替换,使得异源多倍体系统的分析无法很好的将亲本基因拷贝分配到每个亚基因组上。而选择八倍体草莓则有如下两个优势。一,其四套亚基因组仍包含完整的同源染色体子集,大大简化了同源染色体的分配;二,二倍体祖先物种现生近缘种的基因序列(可能仍存在与八倍体草莓中)能用于将同源染色体精确划分到各自亲本亚基因组中。同时,要充分利用草莓作为研究异源多倍体的模型系统,并为确定生物学和农业上重要的基因和应用基因组育种方法提供一个平台,需要一个高质量的八倍体参考基因组。
因此,本研究选取八倍体草莓通过PacBio、10×Genomics和Illumina等平台覆盖了615×的基因组数据,最终组装出0.81G的基因组、包含28条染色体水平的pseudomolecules,约占到预估基因组大小的99%,同时使用其遗传图谱进行校错并通过野草莓(Fragaria vesca)进行同源染色体鉴定。
研究者还注释了108087个编码蛋白的基因以及30703个编码长链非编码RNAs(LncRNAs)的基因,并在注释中鉴定到了高等植物数据库中99.17%的核心基因(1440个),证实了组装的优质。另外,研究也进行了重复元件如DNA转座子、LTR-RTs以及非LTR反转录转座子等的注释,分析发现转座元件(TE)占到了基因组约36%,而其中LTR-RTs丰度最高,占约28%。最后,对质体和线粒体基因组的组装注释也体现了组装的完整度。
图1 二倍体和八倍体草莓基因组的共线性(a. 本研究中的八倍体草莓和二倍体的F. vesca基因组宏观共线性比较,红色为F. vesca,紫色为F. nipponica,蓝色为F. iinumae,F. viridis 是绿色;b. 1号染色体四个同源拷贝的基因保留模式,颜色编码同a。c. 二倍体F. vesca和八倍体草莓的四个同源区在1号染色体某区域上的微观共线性比较)
起源与演化推断
系统发育分析最重要的切入点其一是物种的选取具有代表性,其二是包含足够的遗传信号。本研究从头组装了31个已描述的二倍体草莓转录组,预测出了19302个直系同源核基因用于鉴定祖先种,是草莓属目前遗传信号最丰富的分子系统发育分析。
①研究揭示了前人未知的两种二倍体祖先,结合地理分布、历史事件和遗传印记推断出了八倍体草莓形成的过程,详见图2。
图2 八倍体草莓的进化历史(图中标识了八倍体草莓的二倍体祖先现生亲缘种、推断的中间四倍体、六倍体祖先和北美现生八倍体野生种,每种二倍体祖先颜色和图1一致)
②系统发育分析表明Fragaria iinumae和Fragaria nipponica是四种现生二倍体祖先的其中两个种,为日本特有,在地理分布上毗邻中国的所有五个四倍体种。
③第三个为分布于欧亚的二倍体Fragaria viridis,与独有的六倍体种Fragaria moschata在分布上部分重叠,研究据此假设四倍体和六倍体是从二倍体到八倍体进化的中间产物,该假设也得到了之前研究的支持。
④研究鉴定到的第四个种为F. vesca subsp. Bracheata,仅分布于从墨西哥到不列颠哥伦比亚的北美西部。
⑤系统发育结合现存种的地理分布推断八倍体草莓起源于北美,而F. vesca subsp.Bracheata极有可能贡献了八倍体草莓形成的最后一个二倍体。该发现与之前研究一致,也得到了本研究卡麦罗莎(美国二十世纪九十年代育成草莓品种,长势健旺)质体基因组分析的证实。因此,可以推断六倍体祖先可能从亚洲传入北美并在约1.1个百万年前与F. vesca subsp. Bracheata的当地种群杂交。
异源多倍体的亚基因组优势
大部分古老的异源多倍化事件后,其中一个亚基因组通常会呈现主导优势,如基因含量更高、高表达以及更强的选择压力等,且和亚基因组含量的差异及TEs(转座元件)调控相关(基因表达水平与其附近的TEs密度呈负相关),因此可以基于TEs富集和分布来预测基因表达优势和单个同源染色体水平上最终的基因缺失。
基于上述鉴定的二倍体近亲,研究对四个亚基因组进行了动态分析,鉴定出F. vesca祖先种为优势亚基因组的供体(见图1),保留了多出20.2%的蛋白编码基因和多出14.2%的lncRNA基因,并比其他同源染色体少了19.5%的TEs。相较于其他亚基因组,F. vesca同源染色体基因附近的TEs密度也最低,同时其串联基因重复也多出约40.6%等。这些都表明了F. vesca亚基因组承受了更多的选择压力以保留基因,包括串联重复基因。
研究还分析了影响草莓产量的疾病抗性基因(R genes)。进来研究证实很多R蛋白通过整合的诱饵结构域(decoy domains)识别病原体效应物,而F. vesca基因组编码了20个这种蛋白模型。八倍体草莓扩张了105个融合到R蛋白结构上分化的结构域且具有潜在的整合诱捕功能。
另外,研究还发现不同于F. vesca祖先种,其他祖先种染色体部分区域保留了更多的祖先基因,这些区域是同源染色体交换(HEs)或基因转化事件的产物。值得注意的是,大部分八倍体草莓的HEs涉及显性F. vesca亚基因组对应区域的同源染色体替代。如系统发育和比较基因组分析显示相对F. iinumae,HEs呈现7.3×的偏向于F. vesca亚基因组,但它们并不是像以前报道的那样是单向的。这些结果都表明,F. iinumae亚基因组部分已被F. vesca亚基因组所取代。当然,结果也表明F. vesca亚基因组在草莓抗病性上起主要作用,同时其他三个二倍体祖先种也贡献了抗性机制的多样性。
最后研究还进行不同器官的基因表达分析,而结果证实显性F. vesca亚基因组确实有更高的表达,也支持了亚基因组表达优势受亚基因组之间TEs密度差异影响的观点(图3)。
图3 亚基因组表达优势(灰色柱状图是所有可测试的同源对的表达偏倚,即HEB,可测试的同源对标准基于共线性、>80%的系统发育自展支持率以及转录组数据中至少包含一条read来判定;红色表示同源对显著偏向F. vesca同源基因,偏向三个二倍体祖先种中的一个则用黑色表示)
此外,研究还揭示八倍体草莓中的大多数HEs最终都会导致显性F. vesca亚基因组替代其他亚基因组其中一个对应的同源区域。因此,图3中观测到的F. vesca亚基因组同源基因表达倾向低估了转录组范围的表达优势(占所有转录本的68.7%)。这种偏向导致某些生物学通路如本研究中包括草莓风味、颜色和香味等的代谢通路很大程度上由一个显性亚基因组控制。
参考文献
Edger P P, Poorten T J, VanBuren R, et al. Origin and evolution of the octoploid strawberry genome[J]. Nature genetics, 2019: 1.
参考级基因组如何组装?野生大豆告诉你

Fig.1野生大豆与栽培大豆种子
方法与结果

Fig.2野生大豆基因组de novo组装流程

Fig.3控制大豆种皮色素沉着的结构变异
转座元件(TEs)是一种重复的DNA序列,是植物基因组的重要组成部分,是自然选择和人工选择的重要变异源。基于同源分析和从头预测,研究者发现W05基因组中53.9%的序列为重复元件,其中最主要的重复类型为LTR。通过将基因组进行比较分析,研究者在W05、Wm82_v2和ZH13中分别鉴定出了~2300–3000个TE插入事件。其中,W05基因组中分别发现了419和400个有别于Wm82_v2 和 ZH13的转座插入事件。
高质量的基因组使精确的结构变异比较分析成为可能。研究者通过比较分析发现,与W05相比,Wm 82_v2和ZH13分别有32个和12个大的结构变异,其中最大的结构变异为发生在11号和13号染色体上的染色体间相互易位。研究者还发现控制蔓生长、单株种子数和单株荚数的QTL位点位于易位区。

Fig.4 I locus区种皮色素沉着的倒位现象
此外,研究者还利用光学图谱(OM)数据对发现的结构变异与其他的大豆品种进行比较分析及验证。对I locus 区的比较分析表明,有色种皮的大豆如W05在该区域未发生倒位。而浅色种皮的Wm 82和ZH13具有相同的倒位现象(Fig.4)。
野生种质的高质量参考基因组提高了复杂基因组群体遗传分析的精度。在本研究中,研究者展示了W05参考基因组在遗传变异分析、比较基因组和进化研究中的重要作用,例如大结构变异的鉴定、QTL、基因和等位基因的识别。同时还强调了高质量参考基因组与光学作图技术相结合在研究多个种质结构变异中的优势。
参考文献:Xie, M. et al. A reference-grade wild soybean genome. Nature Communications 10, 1216 (2019).
全面升级!未来组PacBio 宏基因组V3.0震撼来袭!
比起其它的地球生命体,人类对微生物的了解可谓“冰山一角”。宏基因组学作为新型研究工具,可以完整获得环境样本中全部微生物的物种信息和功能基因信息。但是,二代短reads 高通量测序很难有效获得微生物群落各方面的特性。基于单分子实时测序(SMRT)技术的PacBio 平台全面升级以后,序列长度、单cell 产出以及CCS reads 准确度明显提升,使在种甚至亚种水平深入解析物种组成、基因和功能注释、比较基因组分析、进化分析等成为可能,进而从根本上阐明微生物群落在生态系统中发挥作用的机制。未来组紧跟PacBio官方升级脚步,PacBio宏基因组测序项目也迎来了全面升级。
全面升级
系统升级
Sequel 系统软件升级至V6.0,配套试剂升级至V3.0
更长读长,平均读长达到30-40Kb
更高产出,单Cell 产出25-40Gb
更高准确度,CCS reads 准确度达99.999%
指标升级
直观反映群落组成和群落功能
高丰度菌株基因组完成图
DNA 甲基化修饰分析
服务升级
极速服务周期,从样品到标准分析报告仅需4周
研究策略
技术路线
未来组宏基因组数据展示
Chemistry V2.1 vs V3.0
与V2.1 试剂相比,V3.0 试剂下机有效CCS reads、单个cell 产出及高质量CCS 序列数量都明显增加。
NGS vs PacBio
Mock metagenomic sample ——20 株菌,总共基因组大小61,499,437bp
部分分析数据展示
某肠道微生物样本,基于Q20 以上的CCS 序列进行宏基因组组装:
基于有参分Bin 及无参分Bin 分别获得7 个和9 个单菌基因组草图,最高完整度分别高达98.29% 和95.82%。以Bin13 为例,进行DNA 甲基化修饰分析,统计结果如Table 2 所示。
分析图例展示
微生物是世界上存在范围最广泛的生物,从与人类生存息息相关的陆地、空气、水体,到特殊的极寒极热地区,再到工业生产、昆虫肠道、植物根际,以及人类皮肤、肠道等等,尚有许多关于微生物与环境互作的秘密等待人类去探索。未来组应用第三代长读长PacBio测序技术,将为您展示最全面最直观的微观世界全景图。还等什么,抓紧时间联系我们吧!
延伸阅读:
讲真,DNA甲基化多样性还需长读长技术来搞定
先睹为快!未来组Bionano DLS实测数据首发
武汉未来组基因组测序中心拥有新型Bionano Saphyr平台,已完成Bionano测序分析项目上百个,具有丰富的实际操作经验和超强的数据分析能力。目前基于DLS标记技术,实现了高效率、高质量的数据产出(Figure 1)。
Figure 1未来组某植物Bionano实测数据统计
未来组Bionano DLS辅助基因组组装案例展示
未来组成功在多个基因组组装项目中应用Bionano数据辅助提升基因组组装质量,辅助基因组组装提升能力可达百倍级别。如某植物基因组,辅助组装前contig N50为0.27Mb,以53X的有效覆盖数据辅助组装以后,scaffold N50达到38.9Mb,提高了144倍(Table 1)。
Table 1未来组Bionano DLS辅助组装项目案例
新技术DLS的推出,使得Bionano可以在基因组研究领域中得到更广泛的应用。在基因组结构变异研究方面,DLS标记技术超高的灵敏度和超强分辨率能够更加特异性的揭示基因组真实结构;在基因组图谱研究方面,DLS标记技术帮助研究者更快速、更高效地提高基因组组装的准确性和完整性。Bionano DLS将成为世界各国基因组研究者不可缺少的研究工具,未来组积累的丰富项目经验为您的基因组研究保驾护航。
附:近两年发表的三代+Bionano辅助组装基因组高分文章(“*”为未来组参与合作项目文章)

讲真,DNA甲基化多样性还需长读长技术来搞定
DNA甲基化是一类重要的表观遗传修饰,除真核生物外在多种原核生物中也有发现,如保护宿主细胞免受噬菌体或细胞外DNA入侵的序列特异性限制性修饰(RM)系统,已经成为生物技术研究的重要工具。
此外,原核DNA甲基化还有调节基因表达、错配DNA修复以及细胞周期等功能。但是,目前大部分研究者仅局限于实验室可培养稀有原核生物的研究上,因此对原核生物甲基化系统多样性的了解知之甚少。基于单分子实时测序(SMRT)技术的PacBio平台可以说是DNA甲基化研究者的新宠,像原核菌株中的N6-methyladenine (m6A),5-methylcytosine (m5C)和 N4-methylcytosine (m4C)等修饰都已有使用该技术进行解析的案例。尤其是其环状一致性(CCS)测序模式能够对同一模板生成整合多条序列的单条超长read,从而无需进行无性系种系培养,目前CCS仅用于部分鸟枪法宏基因组研究,还没有应用到宏表观基因组或环境微生物群体的直接甲基化分析上。
近日,来自东京大学的Satoshi Hiraoka等人使用PacBio CCS技术研究日本最大的淡水湖——琵琶湖微生物群落的鸟枪法宏基因组和宏表观基因组,绘制了19个细菌和古细菌草图基因组,并揭示了22种甲基化motifs,其中9种是首次发现,还通过计算预测和实验验证了与之相关的4种甲基转移酶(MTases),其中2种被发现能够识别新的motif序列。该研究表明宏表观基因组作为一种强大的方法可用于鉴定自然界中未探索的多种原核DNA甲基化系统,研究成果发表在Nature Communications上。
采样测序和分类分析
采集琵琶湖不同深度水样,PacBio测序得到的数据和环状一致性分析统计见表1,其中90%的CCS reads都是高质量的。
通过Kaiju和NCBI数据库对CCS reads进行分类,结果见图1。其中,>88%鉴定到门,>56%鉴定到属,高于同样计算方法的基于Illumina测序的鸟枪法宏基因组。在门水平上,变形菌门在两种样本中均占主导,其次为放线菌门、疣微菌门和拟杆菌门(Fig.1)。绿弯菌门和奇古菌门在深水中尤其丰富,古生菌在浅水中尤其稀少。后鞭毛生物为最优势的真核生物,其次为囊泡虫总门和不等鞭毛类。而病毒中占优势的依次为有尾噬菌体目和藻类去氧核糖核酸病毒科。研究结果与之前淡水湖微生物群体研究一致,表明PacBio平台和短读长测序技术平台一样适合进行宏基因组分析。
图1 CCS reads的系统发育分布,a、b、c分别代表域、门、纲的富集推断,忽略真核生物和病毒reads
宏基因组组装和基因组binning
研究从浅、深水中分别组装出554、345条contigs,对应的N50为83kb和76kb,最长contig分别为481kb和740kb,远远超过了之前将CCS应用于鸟枪法宏基因组分析活性泥微生物群体的组装结果。然后使用MetaBAT进行binning,浅水和深水分别有52.3%和29%分配到15个和4个bins上,且有46.9%和44.8%的CCS reads可以map到草图基因组上(见图2和表2)。
图2 对组装的contigs进行基因组binning。每一种颜色和大小的圆圈代表所属的bin和序列大小。
对每个bin进行草图基因组组装,基因组完整性在17%-99%之间,污染均低于3%,基因组大小在1.0-5.6Mb之间,GC含量从29%-68%,平均N50为24kb,最大1.67Mb。19个草图基因组划为7个门,其中7个推断为放线菌门、4个为疣微菌门,另外,最丰富的变形菌门只在浅水中装出2个草图基因组,深水甚至没有装出来。总之,系统发育重建很可能反映出了琵琶湖中优势、但尚未培养的微生物谱系。
宏表观基因组分析
从10个草图基因组中检测出29种候选甲基化motifs,甲基化比例从19%-99%(可能低于真实的甲基化水平),见表3,其map上的subreads覆盖度从28.7×到297.3×。
可能由于单个甲基化motif检测不完整或包含在该基因组中近缘谱系间的异质motif序列,变形菌BS12基因组中的3个motifs包含了相似的序列,而拟杆菌BS15基因组中则观测到了回文motif和5个互补motif对,值得注意的是,绿弯菌门的3个草图基因组——BS1、BS3和BD1共享相同的motif序列集,这可能是由于进化共享的甲基化机制引起。大体上,每个草图基因组的contigs都呈现相似的甲基化模式,也为基因组binning提供了表观基因组上的支持。
当然,也有至少9个motifs没有匹配上任何现有的REBASE序列上,这也暗示环境原核生物DNA甲基化体系中还存在着诸多没有被挖掘的多样性,包括一些没有培养过的菌株。
对应检出甲基化motif上的已知MTase
研究还尝试鉴定能够催化检出甲基化motif甲基化反应的MTase,对MTase基因进行系统性的注释。首先通过相似性搜索从9个草图基因组中鉴定到20个MTase基因,找到了丰度最高的Type II MTase等,以及编码REases等蛋白的基因,20个MTase基因中有7个和通过宏表观基因组分析的结果一致,见表3和表4。如奇古菌BD3基因组包含了两个MTase,和识别AGCT和GATCmotif序列的宏表观基因组分析结果高度统一,说明对于环境原核生物甲基化系统的鉴定,宏表观基因组分析是行之有效的手段。
未挖掘的原核生物甲基化系统多样性
20个MTases中也有13个和宏表观分析鉴定的酶序列不相似(表3和表4),说明揭示环境原核生物多样性甲基化体系也需要借助直接观测。研究以拟杆菌BS15基因组、疣微菌BS8、BS10和BS6基因组以及消化螺旋菌BD2基因组、变形菌BS12基因组等为例,从中均鉴定出数目不等的MTases和甲基化motif,但也都出现和报道的MTase或近缘MTasemotif序列不一致或完全匹配不上的情况,甚至也有检出甲基化motif但未检出MTase基因的个例。研究者推断可能的原因是基因组完整度不够,也有可能是这些MTase基因和可培养菌株中的MTase存在分化,还有可能是它们属于新的类群。
含新甲基化motif的MTases实验验证
对于宏表观基因组分析与motif序列高度相似的MTases,研究实验验证了其中四种的甲基化特异性,构建每一个携带人工合成MTase基因的质粒,转导到大肠肝菌细胞中强制表达,然后通过REase消化观察独立质粒DNA的甲基化状态,甚至个别还进行了PacBio测序来检测两个MTase基因各自转化的大肠杆菌染色体DNA甲基化情况。最后,研究鉴定出一系列新的、具有甲基化特异性的MTase基因,并对其表达的蛋白进行命名。
探索自然界中原核生物甲基化系统的
宏表观基因组
一系列的分析证实,PacBio SMRT测序结合CCS模式进行宏表观基因组的分析是非常有效的,较基于序列一致性和基于培养的甲基化系统分析以及短读长宏基因组分析都有显著的优势。
CCS reads让宏基因组组装、binning和基于序列的蛋白分类都变得更容易,最重要的是,该方法揭示了多个甲基化motif,包括环境原核生物中的一些新的motif和本研究中鉴定出的4种MTases。当然,由于需要更深度的覆盖,SMRT测序目前应用于更多种、更复杂样本的宏表观基因组还略显不足,检测到的DNA甲基化类型尚有限,但是要捕获足够的reads以构建稀有样本的长contigs并检测甲基化motif本身难度还是很大。同时,Oxford Nanopore也可以提供基于另一种原理的检测技术。
本研究充分说明即使亲缘关系极近的菌株之间甲基化motif和MTase也会有很大的差异,而且MTase除对抗噬菌体外,还发挥着不为人知的适应性作用。值得注意的是,宏表观基因组分析可以用于多种生物信息分析上面,也能够增进人们对环境原核生物甲基化机制及其应用的了解,而这里新鉴定的MTase也有其他生物技术上的应用。
怎么样?阅览了这篇文章,是不是认可了DNA甲基化多样性还需长读长技术来搞定呢?武汉未来组是中国最早一批引入第三代长读长测序技术的生物科技公司,基于多年的摸索研发经验,已经基于Pacbio CCS reads搭建了成熟的分析流程(见下图),完成了多个环境微生物宏基因组的分析项目,积累了丰富的经验。
未来组宏基因组分析流程
13.5Kb CCS reads升级人类基因组变异识别和组装
短读长测序技术因准确度较高,常应用于单核苷酸变异(SNVs)和小型插入缺失(indels)的识别上,但是较少进行从头组装、单体型定相和结构变异(SVs)检测等远程应用上。基于单分子检测的长读长测序技术通过read-to read的校错也可以达到高准确度,只是计算量较大,并且在校正过程中错误映射的reads和混合单体型中仍存在错误,故较少应用于SNVs和indels识别。同时,人类基因组仍需要结合各种测序技术全面覆盖各类基因变异,而长读长测序技术尤其是基于SMRT测序的PacBio平台正好可以通过CCS来进行高精度的变异捕获。
来自PacBio的科学家Aaron M. Wenger联手多家机构的研究人员选取了理想的基准基因组样本HG002进行长CCS reads的表现分析以验证其在大、小片段变异检测、基因组组装和单体型定相等方面上的优势,研究结果已于1月发表于bioRxiv预印本期刊上,组学君再为大家细读好文。
CCS文库制备,测序及reads质量验证
长度紧密分布在15Kb的SMRTbell文库用于CCS测序,仅保留准确度高于Q20(99%)的reads,生成的89G数据平均读长13.5±1.2Kb,预估准确度中位数为Q30(99.9%),平均值Q27(99.8%),且与GRCh37一致性准确度高。
图1 HG002的高精度长读长测序(a. 多轮DNA模板测序生成CCS序列,b. 不同测序轮数得到的CCS序列质量预测,c. CCS reads的读长和预测质量)
研究统计了CCS reads和HG002参考基因组的不一致,其中3.4%为错配,4.6%为非均聚物背景下的indels,92%为均聚物背景下的indels。通过read-to-read比对的方式计算错误率发现与预估的CCS reads质量一致,平均reads精度高达99.8%。同时,研究还将其映射到GRCh37上,匹配度达97.5%,还发现对应的NGS短读长数据上报道的193个医学相关基因可以在CCS数据上映射152个,包含 CYP2D6、GBA、PMS2和STRC等基因。
图2 使用CCS reads的人类基因组可映射性(a.不同映射质量阈值下NGS和CCS数据对GRCh37人类基因组的覆盖度,b. NGS和CCS数据对先天性耳聋基因STRC分别在以10为映射质量阈值下的覆盖度,c. 13.5Kb CCS reads对193个人类基因的映射提升,这193个基因之前使用NGS数据进行映射发现与医学相关且存在问题)
小片段变异检测和定相
GAKT用于SNVs和小型indels识别,相较于基准基因组,SNVs精度达99.468%,召回了99.559%;而indels的精度和召回率分别为78.977%和81.248%。针对SNVs和NGS数据进行比较,识别的indels相对少一些。研究还通过一种叫Google DeepVariant的模型进行数据模拟,发现相较于Illumina数据,CCS数据可以大大提高SNVs和indels的准确度和召回率。
接下来,研究又使用WhatsHap来定相DeepVariant识别的结果进行定相以确定CCS reads能够生成单体型所需的高精度变异识别和远程信息。几乎所有(99.64%)的常染色体杂合变异都被定相到19215个区域上,N50为206Kb。定相区域大小分布几乎和理论上的极限完全匹配,而这个极限是通过生成变异之间超过平均CCS读长13.5Kb的间隔来评估的,这说明定相区域大小受HG002基因组读长和变异量所限,而非变异识别的覆盖度或质量。
图3 使用CCS reads进行变异识别和定相(a. 使用DeepVariant进行的SNV和indels识别和基准基因组的一致性,b. 杂合位点的DeepVariant变异识别定相和使用13.5Kb reads对HG002理论定相的比较,c. 整合的CCS SVs识别和基准基因组一致性,d. 通过变异大小进行上述一致性的比较)
用单体型定相升级小片段变异检测
GATK和DeepVariant在识别变异时都不直接并入远程单体型相位信息。研究对CCS reads基于GIAB的三相变异对进行单体型标记,然后使用DeepVariant模型对reads按单体型分选通过的顺序进行模拟,发现对于SNVs的识别表现和原始的DeepVariant模型相近,但是对于indels的识别却有着显著的升级:精度达97.835%且召回率达97.141%(见表1)
表1 CCS reads小片段变异识别的表现
CSS reads的SVs检测
≥50bp的插入和缺失SVs使用两款基于read映射的工具pbsv和Sniffles,使用paftools分析Falcon和Canu的从头组装来升级更大片段变异的识别,研究发现使用SURVIVOR生成整合的callset对于<1kb和≥1kb的变异以及插入和缺失的表现都相仿,凸显了基于映射的和基于组装的SVs识别之间的互补性(见图3c和d)。
另外,研究还对Illumina短reads和10×Genomics的连接reads(linked reads)进行SVs识别比较,发现所有的短reads和链接reads无论在准确度还是召回率上面都不及CCS callsets。
从头组装
使用了Falcon、Canu和wtdbg2对CCS reads分别进行从头组装,跳过原始的read-to-read校错步骤的组装结果连续性都很好,contig N50从15.43到28.95Mb,其中Canu因为将一些杂合位点的等位基因算作分离的contigs上而组装的基因组大于预期(表2)。同时,研究选择了HG002亲本的短reads来鉴定对一方亲本唯一的k-mers然后基于单倍型对CCS reads进行分区(trio bin),最后选择51-mer的binning进行组装。
表2 CCS reads从头组装统计
三款组装软件独立对父系和母系reads进行组装也都得到高连续性的近完整组装结果,N50从12.10到19.99Mb,基因组从2.67到3.04Gb。父母系的组装中都鉴定出从95.3%到98.2%的单拷贝人类基因(表2),并且无论是混合组装还是一方亲本的组装都与HG002基准基因组保持高度一致,甚至大大超过了之前发表的使用PacBio长reads和用Illumina数据打磨过的ONT数据组装结果(图4a)。
此外,研究还组装到了跨度超60Mb的大片段重复,较之前提升了20%的连续性,以及解析了一致性99%-99.5%的15kb重复片段(图4b)。
图4 read准确度对从头组装的影响
变异识别和从头组装需要的覆盖深度
研究者调取部分数据进行分析,发现对于SNVs,使用DeepVariant进行准确度和召回率超过99.5%的识别仅需15×的覆盖,而对于indels,超过90%的准确度和召回率需要17×。再者,使用混合单倍型的wtdbg2的组装,只要覆盖度高于15×,就可以保持高于Q42的一致性准确度等。
基准基因组校错和扩展
最后,研究对SVs仍处于草图形式的基准基因组进行校错,如鉴定到小片段变异中,31个基准基因组中的均聚物不一致处有29个被误为正确的;而在一项Illumina全基因组病例研究中,仅有很少的假阴性SNVs和indels以及假阳性indels和CCS callset相一致,通过人工处理后评估基准基因组中有2434个(95%的置信区间为1313-2611)错误可以通过CCS reads来校正。
对于SVs,基准基因组精确度比较高,但不一定完整。将CCS DeepVariant callset加入基准基因组小片段变异整合流程中可以将基准区域扩展1.3%相当于418875个变异。而18832个常染色体变异识别中仅有9232个和基准区域相重叠,也表明并入CCS变异识别可能可以将基准基因组中变异的数量翻2倍以上。
结语
读长和准确度宛如测序天平的两端,这份研究开发出了一种新的流程力图让天平达到完美的平衡:研究基于单分子环状一致性序列(CCS)生成准确度达99.8%、平均读长达13.5Kb的长reads,还将该流程运用到已精确解析的人类基因组(HG002/NA24385)上。
通过优化现有工具对HG002/NA24385进行变异检测,召回了超过99.91%的SNVs、95.98%的indels以及95.99%的结构变异。另外,研究还评估了和参考基因组中2434处不一致的可纠正性错误,且几乎所有(99.64%)检测到的的变异都可以定相到单体型上,而从头组装得到contig N50超过15Mb、一致性达99.998%的高连续性、高精度基因组。将CCS reads匹配上短reads以进行小片段变异检测,对于能检测到结构变异和相似组装连续性的区域,CCS reads明显比含噪音信号的长reads有更高的一致性。
史上最全的长读长数据校错方法大比拼

使用高精度、低成本的SRs数据即使低深度也可以实现,基于算法设计分为比对类、图像类以及结合两种策略的双重(比对/图像)类三种
那么,针对长读长数据如何选择最合适的校正方法呢?
为了回答这个问题,美国爱荷华大学内科部区健辉副教授及其团队的Shuhua Fu、Anqi Wang等人选取了10款针对长读长数据性能最卓越的校正工具在常规基准下进行灵敏度、准确度、产率、比对率、读长、运行事件、内存占用以及从头组装和单体型序列解析等方面的比较评估,从而基于可用数据大小、计算资源和个体研究目的等给出了全方位的方法选择指南[1],研究成果已于2019年2月4日被Genome Biology在线收录,影响因子为13.214,下面请随组学君先睹为快吧!
———- 数 据 来 源 ———-
基于PacBio或ONT平台的Escherichia coli(大肠杆菌)和Saccharomyces cerevisiae(酿酒酵母)作为“小”数据集;
基于PacBio平台的Drosophila melanogaster(黑腹果蝇)和 Arabidopsis thaliana(拟南芥)作为“大”数据集。
———- 评 估 设 置 ———-
为考察SR覆盖度对错误校正性能的影响,研究对覆盖度为5×、20×、50×、75×以及100×的每一个SR数据集生成随机子数据集;所有的计算均在16核、256G内存的20台设备上进行。
———- 评 估 策 略 ———-
为评估不同校正方法的性能,原始的和校正的reads使用BLASR比对到对应的参考基因组上。真阳性(True Positive,TP)位点为由单个校正工具进行校错的位点,假阴性(False Negative,FN)位点则是没有经过校正的错误位点。TP和FN位点能够通过比较对参考基因组进行原始和错误数据的比较来使用Error Correction Evaluation Toolkit(校错评估工具包)进行计算;错误率通过比对域(alignment)上的插入、缺失和替代总数除以每条read比对上的总长度来计算。统计参数的计算方式如下:
- 灵敏度:TP/(TP+FN),TP为校错工具校正的错误位点数目,FN为未校正的错误位点
- 准确度:1 – error rate
- 输出率:原始数据中reads的比例
- 比对率:输出数据的reads比对到对应参考基因组上的比例
- 输出读长:输出reads的长度
- 运行时间:校错工具消耗的运行时间
- 内存使用:校错工具占用的峰值存储空间
- 对输出拆分reads的Jabba、ECTools、pacBioToCA和Nanocorr等,研究使用如TH×TP/(TP+FN)这样的灵敏度计算策略,这里TH是输出碱基的数目相对于原始reads全部碱基数目的比值
- 运算时间超过20天的工具不计入评估(如比对类的ECTools和LSC)
- 选用方法:selective method,已知Jabba、ECTools和pacBioToCA丢弃未校正的碱基和剪切片段并输出部分LR数据,因此产量较低,将其称作选用方法,在图中会标注下划线。
十款本研究进行评估的校错工具信息如下表,参数设置还可以参考附件中的Note 1[2]
———- 评 估 结 果 ———-


和灵敏度的整体趋势相似,准确度基本上也是随着SR覆盖度增高而提升,但是也有Jabba、Nanocorr和pacBioToCA这三种方法例外。当Jabba输出选定比例的LRs,准确度在所有数据集中是最高的,几乎达到了Illumina数据的质量。因此,对于Jabba而言,虽然更多的SRs可以有更好的灵敏度,但是增加SR覆盖度几乎没有其他额外的作用。
同样的原因,pacBioToCA也表现出了高准确度。有趣的是,虽在≥20×的SR覆盖度时,另一个选用方法ECTools的准确度与Jabba近似,但其在5×时却相对要低。总而言之,即便这些选用方法都没有输出全长LRs的全部数据集,它们在足够的SR覆盖度下还是输出了高准确度的LR片段。
研究还对小数据集的准确度饱和作了分析,但相对最引入注意的还是大部分方法都在相同SR覆盖度下达到了灵敏度和准确度的饱和值这一现象,当然也有Jabba、ECTools和Nanocorr例外。
非比对类校正方法能应用于所有四个PacBio数据集中,其在小数据集上的准确度比大数据集要更高或相当。校正前,原始的ONT LRs相较于PacBio的准确度要低,这一劣势在校正后仍复如是(由选用方法校正的除外)。不过在处理小数据集时,不同选用方法输出的PacBio和ONT LRs在准确度上并无本质的差异,因其只剪切并输出高质量的LRs片段。应用了其他方法后发现校正后的ONT和PacBio LRs之间的准确度差异降低,LSC除外。
LRs的错误分布随机,还有一些易出错的LRs因难以校正而使校正工具无法输出数据,因此,评估校错后的数据量很重要。对于PacBio数据,除了Jabba、ECTools和pacBioToCA,其他方法产出都超过输入数据的90%。尤其要指出,FMLRC、HALC和CoLoRMap产出了100%的数据集。LoRDEC和proovread仅剔除了少许零碎的reads,LSC和Nanocorr因比对类方法可能无法校正而只输出含比对SRs的LRs仅丢失小部分的数据,尤其是这类方法的输出率并不取决于SR覆盖度和LR数据大小。高输出率便于用户继续依照自己的研究目的执行分析,而通过改变方法设计和执行方式以保留校正和未校正的LRs可以实现100%的输出率。
相反,Jabba、ECTools和pacBioToCA这三种选用方法只输出校正区域的剪切reads, pacBioToCA在所有测试方法和数据中输出率是最低的。当然也不绝对,如在除E.coli数据外、5×SR覆盖度时的Jabba以及小数据集在5×SR覆盖度时的ECTools,只不过这两个方法多少还是损失了些关键数据。
FMLRC和CoLoRMap在校正ONT数据集时输出率也达到100%,LoRDEC、HALC和proovread同样接近98%。需要引起注意的是,LoRDEC和HALC都没有成功校正长度>100kb的LRs。除开这些高输出率的方法,针对同样的物种,其他方法相较基于pacBio的数据集输出均要低一筹——基本规律是:原始错误率越高的数据,得到校正的LRs就会越少。
因准确度会影响比对率而比对工具能够接受一定数量的错误,所以这两个指标独立评估。更重要的是,SVs及可变剪切位点的识别对准确度要求严格,而融合基因检测和富集评估等则对高比对率有更高的需求。对基于PacBio的小数据集而言,比对率和输出率形成了鲜明的对比:ECTools、pacBioToCA、LSC和Nanocorr比对率均相对较高,而输出率却十分低。FMLRC、LoRDEC、HALC、CoLoRMap和proovread 则恰好相反。导致这种情况的原因可能是Jabba、ECTools和pacBioToCA等选择性输出容易比对的高质量LR片段,而那些高输出率的方法则包含了未校正的LRs,从而导致比对率较低。总的来说,比对类方法(除proovread外)倾向于获得更高的比对率,因其输出的LRs已经由SRs校正过,反过来,它可以作为种子文件来比对输出的LRs。
另外,在处理小数据集时,LoRDEC、HALC、CoLoRMap和proovread以及LSC、Nanocorr、pacBioToCA等都随着SR覆盖度的增加表现出不同程度的比对率降低,而FMLRC、Jabba和ECTools则没有明显的受SR覆盖度的影响。同时,图像类方法在处理大数据集时的比对率相对于小数据集要低得多。 Jabba在大数据集中得到了最高的比对率。
再者,LoRDEC和CoLoRMap的比对率在SR覆盖度增加时表现出轻度的提升,不过仍只在30%左右。与之相对,FMLRC、Jabba和HALC在SR覆盖度增加时其比对率显著提升,Jabba比对率仍是最高的。
因为准确度相对PacBio要低,ONT原始LRs的比对率也较其低一些。 在校正之后,对于同样的物种,大多数方法输出的ONT数据比对率要明显低于PacBio数据,不过选用方法和Nanocorr除外,它们只输出部分数据集。
除了选用方法和Nanocorr,其他所有方法的输出读长和输入读长都相近。整体看来,由于PacBio数据的主要错误类型是插入错误,所以随着SR覆盖度的加深,输出读长普遍会有轻微的减小。与此相对,选用方法输出读长则有着显著的特征:要么远长于、要么远短于输入LRs。这可以用两个原因去解释:①选用方法只输出一小部分LRs;②其输出了含校正区域的剪切后的reads。
对于小数据集,Jabba在5×SR覆盖度处输出读长非常短,随着SR覆盖度升高,读长显著提升,这与输出率的表现一致。虽ECTools在5×SR覆盖度处输出读长中值远小于输入读长中值,但自5×到50×时,读长有明显提升随后达到饱和,同时,50×处的读长中值远超过输入读长。当SR覆盖度足够时,ECTools就会选择性的去除一些短LRs。
PabBioToCA则不一样,输出读长短,即使增加SR覆盖度也没有本质的提升,其在运行时会将LRs剪切成很多小的校正片段。而Nanocorr输出读长就更短了,不过和Jabba相似,随着SR覆盖加深会有提升。只是在低SR覆盖时,Nanocorr输出读长普遍较Jabba长些。尤其是Nanocorr输出率始终在95%左右,而选用方法则从0.40%到88.35%不等,这可能的原因是Nanocorr将每一个LR剪切为小片段如5’和3’末端而丢失了SR覆盖度,选用方法是将单条read剪切为多个校正的片段,而其他的方法输出包含校正和未校正区域的完整reads。
和E. coli的PacBio数据集一样,输出的ONT读长随着SR覆盖加深会有轻度的降低,选用方法和LSC以及Nanocorr例外。LSC输出读长较原始的ONT LRs要长的多,因更长的ONT LRs更有可能进行用于LSC校正的SR比对。
对于PacBio LRs而言,图像类方法处理小数据集普遍耗时较比对类方法要少,随着SR覆盖加深,差异更加鲜明,因图像类方法使用SR构建图像而非直接使用校正的SRs,且在SR数据足够时,图像大小也没有明显增加,同时构图阶段比较省时——所以SR覆盖不像对比对类方法影响那么大。
综合比较,处理小数据集时不同SR覆盖度的Jabba表现完胜其他方法,仅32-1041秒,但在处理大数据集时,FMLRC和LoRDEC用时最短。不出预料,结合两种策略的双重类方法HALC和CoLoRMap耗时居于两类方法之间。
对于比对类方法,pacBioToCA耗时最短,但在高深度SR覆盖和处理大数据集时崩溃了;而同为比对类的ECTools、LSC、Nanocorr和proovread方法在处理小数据集时则耗时要久的多——这主要是因为将LRs分离进行SR比对较为耗时。
ONT LRs相对PacBio数据就省时一些,因原理上的差异,比对类和图像类方法处理错误区域的校正方式不一样,因此除了选用方法外,其他方法中图像类均比比对类要省时,双重类方法居中,该趋势和PacBio数据一致。
在校正Pacbio LRs时,LoRDEC不受SR覆盖影响成为内存使用最小的方法,无论大、小数据集都只要1-2G,可以理解为其使用GATB内核构建SRs de Bruijn图和图中的遍历路径,这个过程相当节省空间。但Jabba和部分FMLRC的内存使用在处理大数据时就显得受SR覆盖影响较大,如Jabba处理黑腹果蝇(D. melanogaster)数据时由于增强型稀疏后缀数组的存储影响,在100×SR覆盖处占用超过180G内存。
双重类方法JALC和CoLoRMap较图像类和部分比对类方法的存储效率要低。而比对类的ECTools、LSC和pacBioToCA比较节省内存,另两种比对类方法Nanocorr和proovread则比较耗内存,主要由LR数据大小或SR覆盖度的影响。
考虑到实际应用,常规的32G内存可以满足FMLRC、Jabba、 LoRDEC、 ECTools和LSC处理小数据集的需求,CoLoRMap、Nanocorr和proovread则需要更高的配置以应对SR覆盖度的提升。在处理大数据时,除LoRDEC仅需2G内存,FMLRC、Jabba、HALC和CoLoRMap都需要高配置。
总体上,ONT LRs需要的内存较PacBio小。比对类较耗内存的proovread在处理ONT数据时没有像处理PacBio数据时一样崩溃,而图像类方法也比较节约内存。
综合评估上述指标,以E. coli的PacBio数据为例,十款方法可以分割为三个大类:图像/双重类,比对/双重类和选用方法。普遍来讲,图像/双重类的综合表现最佳,而FMLRC较其他方法更为突出。
对于第一大类,图像类的FMLRC和LoRDEC与双重类的HALC表现相近,无论是灵敏度、准确度、输出率和比对率还是内存占用都比较好,尤其是效率高,仅次于Jabba,虽然Jabba只输出部分数据。不过唯一需要警剔的一点是,它们的灵敏度相当依赖于SR覆盖度。
第二大类由比对类的LSC、Nanocorr和proovread以及双重类的CoLoRMap构成,它们的灵敏度、准确度和输出率较第一大类略低,但比对率要稍高一些,尤其是和LoRDEC及HALC比较而言,不过问题也很明显,那就是计算效率低(较耗内存的LSC除外),从而这一大类都比较耗时,其中的Nanocorr还有输出读长较短的劣势。
第三大类由于仅输出原始数据的部分子集,因而输出读长和输出率都是值得关注的问题,不过胜在内存占用低、准确度和比对率高,且运行时间长(除Jabba外)。
在处理测评用到的两个小数据集时,所有方法表现类似,但在运行大数据集时,图像类和双重类方法都成功,比对类方法却失败了,尤其是,第一大类仍表现出高准确度和高输出率。
很可能是由于计算设计的差异,整体上看,图像类方法相较比对类方法在计算效率上有明显的优势。但另一方面,比对类方法则在低深度的SR覆盖时表现更正常。另外,大部分方法在处理ONT数据时的表现和处理PacBio数据的表现相近,只不过比对率和准确度稍低。
校错的一个重要目的是为从头组装提供更高质量的LRs,研究人员通过提升从头组装的比较来评估十种方法的校错性能。所有组装统一由Miniasm完成。
1.组装contigs的数目评估
理想状态下,congtig数应该和染色体数目相近或相当。当组装基于PacBio数据的E. coli基因组时,原始的或者校正的LRs都组装出和染色体数目相近的contigs数,但对于复杂的基因组组装,校正后的LRs比原始数据组装出的contigs要多很多,只有选用方法Jabba例外。同时,ONT LRs组装出的contigs数和用对应的PacBio数据组装出来的相近。
2.Contig N50评估
因只输出LRs的小片段,选用方法校正过的LRs装出来的N50比其他方法要小很多,甚至比一些原始数据装出来的还短。对于复杂一些的基因组,选用方法还是不敌其他方法对于提升N50的表现。比较有意思的一点是,原始LRs在校正后组装出的N50反而比不校正还短,有可能是因为混合校正引入了错误导致组装错误。
3.基因组分数评估
基因组分数反映组装完整度。使用包括FMLRC、HALC、CoLoRMap、ECTools以及pacBioToCA等方法进行处理PacBio小数据集后的组装拥有差不多或完整的基因组片段(即基因组分数接近或等于1),而Jabba在SR覆盖达100×时基因组分数仍偏低,即使相对另两种选用方法其输出率要高一些。但对于大数据集而言,不论是原始的还是校正过的LRs,基因组分数都比小数据集要高。另外,ONT LRs校正后组装的基因组分数较PacBio LRs要低。
4.Contig序列准确度
所有的方法都是在校正后会对contig序列准确度有所提升,且都是随SR覆盖加深而提升并在饱和后略有改变。虽Jabba校正后拿到了最高的contig序列准确度,不过其基因组分数较低的事实仍无法忽略。
总的来说,除选用方法外,其他混合校正的LRs组装大部分情况下都改进了contig数、contig N50和基因组分数,且随SR覆盖加深也提升了contig序列准确度。选用方法尤其是Jabba虽拥有最高的contig序列准确度,但是牺牲了组装的完整度和连续性。在SRs足够的情况下,ECTools和pacBioToCA输出数据的组装能达到近1的基因组分数。迄今仅有为数不多的算法开发出内置错误校正模块以组装基因组。
杂合度分析对二倍体或多倍体物种非常关键,研究选取父系和母系基因组LRs通过来自这两个单体型、随机混合的SRs进行校正以进行模拟人类基因组数据的杂合位点(314307个)碱基校正并据此评估TGS校错方法的性能。
研究主要展示了资源利用率高并省时的FMLRC和LoRDEC结果——模拟PacBio LRs杂合位点很多错误未被校正,还有一些正确碱基被修改为错误的。研究得到的结论是:在校正过程中,浅层SR深度相对有助于FMLRC维持单体型信息,而LoRDEC则更适合高深度的SR覆盖。
FPR和FNR两个参数在模拟ONT LRs数据中的值要高于PacBio LRs,也就表明校正ONT数据在保有单体型信息上面临更大的挑战。分析还发现,虽然FMLRC和LoRDEC在其他评估中表现卓越,但对于PacBio或ONT LRs,两者都不能复原或维持真实的杂合碱基。因此和SGS基因组数据一样,LRs的杂合位点校错也面临着巨大的困难。
分析表明,自校正和混合校正都能提升LRs的准确度。通道数(pass number)表示生成校准一致性序列的原始reads条数,如CCS reads至少有两条reads生成,因而准确度通常较subreads要高,且随着通过通道数的增加而提升,并在通过通道数为5的时候达到饱和。相似的,ONT 2D reads也较原始reads有更高的准确度,但由于原始ONT数据准确度相对低些且仅含两轮通道数,故而2D LRs准确度没有那么高。
混合校正方法下,PacBio reads在≥50×SRs达到准确度的饱和,相当于PacBio CCS reads在≥5轮通道时饱和的准确度。而ONT 2D reads的准确度通常在混合校正的5×和20×SR覆盖度之间远低于≥50×SRs时饱和的准确度。考虑到PacBio CCS reads在读长和通道数之间的平衡,混合校错非常有助于获得高精度的LRs。
1. 长读长数据目前还非常依赖校错工具,自校正和混合校正策略都可以提升原始LRs准确度,但自校正的局限更大,如就PacBio CCS reads而言,加大通道数代价就是损失读长,对于ONT数据,通道数最多也就2轮。相较下,混合校错方法能获得高准确的LRs,尤其是使用SRs进行校错能便于LRs比对、组装,稀释TGS数据成本高的劣势。
2. 本研究在灵敏度、准确度、输出率、比对率、输出读长、运行时间以及内存占用几个方面全方位的比较评估了十种校错方法的性能。综合测评显示,图像类方法优于比对类方法,其中FMLRC综合性能较其他图像类方法又更胜一筹。但需要指出的是,仅有低SR覆盖的数据,图像类并不如比对类有更强的鲁棒性;除了选用方法(指本研究中基于比对的ECTools和pacBioToCA)外,比对类方法通常具有相当甚至略低的灵敏度、准确率、输出率和比对率,主要缺点是运行时间长和内存占用高。ECTools和pacBioToCA与图像类方法内存占用相当,其与Jabba较适于对准确度要求高但对读长和数据损失不太关注的研究。
3. 与SR覆盖度不同,LR覆盖度不是影响混合校错方法性能的关键因素,一些依赖LR覆盖度的进程如布置图生成或一致性推断,混合校错方法几乎不会进行,但自校正方法就需要考虑到。普遍来说,十种方法在小数据集的表现比大数据集要好,推断原因是基因组的复杂性所致,尤其是高重复的区域对校错造成负担。再者,参数设置也是影响校错工具运行的重要因素,并且线程数是影响运行时间的关键。最后就是安装和启用操作,虽无法量化,但在实际应用中也是至关重要的环节。
4. 基于研究目的的差异,用户可以权衡各种方法的不同性能指标,如在单碱基分辨率或者序列分析的研究中,高准确度是首要关注指标,而在检测基因亚型、融合基因和丰度估计上面,即使准确度相对低一些,但是高的比对率也非常有用。因此,在选择校错方法时,用户应始终综合考虑数据大小、计算资源和研究兴趣。另一方面,本研究的性能评估还确定了在今后优化现有纠错方法或开发新的纠错方法具有一定意义的关键因素。
文章内容博大精深,需要更进一步的了解,请参考原文吧!
参考文献:
[1] Fu S , Wang A , Au K F . A comparative evaluation of hybrid error correction
methods for error-prone long reads[J]. 2019.
[2] https://static-content.springer.com/esm/art%3A10.1186%2Fs13059-018-
1605-z/MediaObjects/13059_2018_1605_MOESM2_ESM.pdf