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有更高的一致性。

史上最全的长读长数据校错方法大比拼

相较于第二代测序技术(the Second Generation Sequencing,SGS)的短读长(short reads,SRs),第三代测序技术(the Third Generation Sequencing,TGS)因其长读长、无GC偏好性等优势能够很好的检测到结构变异、碱基修饰从而解析复杂的基因组和组合的基因组事件,随着数据质量的不断改善,相应的文章产出也逐年攀升(上图以PacBio为例),但其信噪比低而彰显的错误率问题也不可回避。自2012年研发出来的一系列校错方式可以分为两大阵营:自校正和混合校正。自校正通过PacBio或ONT内置的自校正模块进行换装一致性序列(PacBio)或2D reads(ONT)的构建,建立在一定测序深度的覆盖基础上混合校正

使用高精度、低成本的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]

———-   评 估 结 果   ———-

灵 敏 度
上图显示的是10种校正方法对4个PacBio数据集在5种SR覆盖度时的灵敏度,从结果可以看出大多数校正方法的灵敏度会随着SR覆盖度的提升而升高。FMLRC、Jabba、LoRDEC、HALC、CoLoRMap、ECTools和proovread在达到一定灵敏度时趋于饱和,而当SR覆盖度从5×到20×时,LSC、Nanocorr和pacBioToCA仅呈现轻度的灵敏度提升,可能是因为在SR和LR比对后,LSC等三种方法能在低SR覆盖度校正LR区域,但FMLRC、Jabba、LoRDEC、HALC、CoLoRMap这些图像类策略的方法则需要足够的SR覆盖度以构建图像。LoRDEC、HALC和proovread于20×SR覆盖度时达到灵敏度饱和,而FMLRC和ECTools则是在50×时饱和。在LR数据集的作用下,CoLoRMap与Jabba分别在20-50×和75-100×时饱和。在饱和灵敏度方面,FMLRC在所有LR数据集中的表现相较于其他为最佳,不过有意思的是,当使用仅有5×的SR覆盖度的小数据集时,LoRDEC、 HALC、CoLoRMap、LSC、Nanocorr、 pacBioToCA和 proovread的灵敏度均超过了FMLRC。pacBioToCA的表现则在这次比较中基本上处于最下风。整体来看,除了ECTools在S. cerevisae 数据上的饱和灵敏度表现略胜于E. coli的PacBio数据,其他方法都是E. coli数据表现最佳。再者,对于相同大小的LRs,除了FMLRC, LoRDEC, HALC和 proovread例外,校正ONT数据的灵敏度峰值普遍较PacBio数据要低。与PacBio数据相似,E. coli的ONT数据灵敏度基本上都高于S. cerevisae,但也有ECTools例外。 

准 确 度

和灵敏度的整体趋势相似,准确度基本上也是随着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

科研速递|Nature Communications最新文章连连看

单分子测序组装首条非裔人类Y染色体

2019年1月2日,Nature Communications上收录了一篇短小精悍的文章,基于目前仅有欧裔人类参考基因组级别的Y染色体现状,来自西班牙进化生物学研究所的学者们使用ONT平台装配出了首条非裔人类Y染色体,且提升了800%的连续性!原始RNA测序还便于研究者从Nanopore电信号数据中识别甲基化修饰,该研究提供的方法还可以推广到其他任意大型重复基因组的简化组装中。

简单来说,研究从一个淋巴母细胞系(HG02982)中分离出约9000000条染色体,依照自主设计方案纯化后在4个Nanopore芯片上生成了305528条reads,然后进行Illumina双端测序,在使用Canu进行组装后使用Nanopolish进行校正,再通过pilon对组装结果进行校错,最终组装出的Y染色体含35条contigs,N50为1.46Mb,序列总长为21.5Mb。

研究者将其与GRCh38以及大猩猩的Y染色体进行了比较,相对于大猩猩组装的contig N50为17.95kb,本次组装的连续性提升了两个数量级。

下表是基于本研究组装的Y染色体(HG02982)和GRCh38以及NA23385进行序列覆盖的比较。

组装统计

此外,研究者还将基因注释结果与GRCh38进行了比较,并鉴定了其中的结构变异和CpG甲基化的分析等。

人类Y染色体因其高度重复的特性而使其组装充满挑战,但是通过Oxford Nanopore组装出的结果明显优于之前的长读长WGS组装的结果,且更为经济,这预示着将来完整的Y染色体组装将成为可能。

文章来源:

https://www.nature.com/articles/s41467-018-07885-5#Sec9

3000份基因组数据描绘亚洲水稻逆转录转座全景图

2019年1月3日,NC刊登了一份从物种水平研究水稻遗传多样性的成果,来自法国佩皮尼昂大学的Olivier Panaud研究团队通过开发一款名为TRACKPOSON 的工具,从3000份水稻重测序项目数据集中识别出了32种不同家族的逆转录转座子的超过50000个TIPs(转座元件插入多态性),描绘了亚洲稻逆转录转座的全景图。

研究人员使用TRACKPOSON对传统的基于成对末端mapping和reads分离的方法进行优化,先把成对的reads 比对到特定转座子的一致序列上再把提取后的reads比对到参考基因组上,从而省去了传统的将全部reads比对到参考基因组上的繁琐计算以快速、准确地鉴定已知转座子家族的TIPs。

TRACKPOSON工具分析流程

通过这个方法检测了32个家族的逆转录转座子在3000个水稻基因组中的位点,识别出53262个TIPs并且发现识别到的TIPs在群体中大部分是低频率出现的,推测这些TIPs可能是在进入到农业时代后于近期形成。

研究人员还进行了GWAS全基因组关联性分析,探查到相比于转座元件沉默通路上相关遗传因子的改变,这些TIPs可能由外部刺激引起个别转座子的跳跃而诱导激发。另外,TIPs数据还被用来进行水稻起源的分析,研究使用古生物基因组学分析手段发现籼稻(Indica)、粳稻(Japonica)和Aus/Boro水稻的基因组在水稻驯化之前已经具有明显分化,可以推断水稻起源于三个不同的驯化事件。

文章来源:

https://www.nature.com/articles/s41467-018-07974-5#Sec8

参考微生物中m6A的单分子测序检测

DNA碱基修饰N6-甲基腺嘌呤(m6A)参与许多与细菌生存及其与宿主交互作用相关的路径。本研究通过与之前的方法(衡量在每个腺嘌呤通过一个纳米孔时测量值和预期电流值之间的偏差)进行比较,展示了一种新的便携式神经网络算法来检测碱基修饰。该模型通过McAller软件包实现,可以扩展为基于未经模拟甲基化预测的候选甲基转移酶靶模体检测工具。

研究使用PacBio、Oxford Nanopore、甲基化DNA免疫沉淀测序(medip-seq)和全基因组亚硫酸氢盐测序数据来生成、验证8种微生物参考物种甲基化组。这些特征描述详尽的微生物参考基因组可用来研发、评估未来从单分子测序数据中识别碱基修饰的方法。

当m6A周围的DNA穿过nanopore时picoampere电流偏离模型值

研究成果于2019年2月4日发表在NC上,总之,和传统方法相比,该方法能够提升碱基修饰的检测。近期发生在国际空间站、西非埃博拉病毒以及美国寨卡病毒的爆发期间都使用到了Nanopore测序方法,证实了远程基因组研究和表观基因组学研究的无限价值。

文章来源:

https://www.nature.com/articles/s41467-019-08289-9

Direct RNA-seq重新定义病毒病原体转录复杂性

传统的RNA测序方法鉴定复杂的病毒转录组时会因高基因密度、重叠阅读框以及复合的剪切产物而变得困难重重,而direct RNA-seq不失为一种可靠的替代方法。本研究对HSV-1病毒转录组进行直接RNA测序,展示了其在定义与所有多聚腺苷酸化病毒RNA相关的转录起始和RNA裂解位点上的优势,并证明低水平的read-through转录产物产生了一类新的嵌合HSV-1转录本,包括一种编码病毒E3泛素连接酶ICP0和病毒膜糖蛋白L融合的功能性mRNA。因此,direct RNA-seq为表征复杂基因组病毒转录格局变化提供了强有力的方法。

有意思的是,研究将direct RNA-seq和Illumina RNA-seq进行了比较,通过分析发现无论是Illumina还是Nanopore测得的基因区的平均read深度较基因间区要高,这种差异反映了Illumina标准流程中在poly(A)选择阶段富含腺苷区域的mis-priming,相反,在direct RNA-seq中,接头耦合启动蛋白的连接要求意味着只有多聚腺苷酸化的RNA才可以通过nanopore。

图a是poly(A) RNA对HSV-1基因组范围100nt滑动窗口的覆盖情况,黑色线条代表Nanopore测序结果,红色虚线为Illumina数据,红色实线为转化为与Nanopore数据相同覆盖度后的Illumina结果;图b是Nanopore和Illumina测序产生的HSV-1基因组覆盖度关联分析;图c表示direct RNA-seq和标准化后的Illumina数据集的基因和基因间区域的read深度值(100 nt间隔)的点图。基因区和基因间区的read深度平均相差12.08倍(Nanopore)和6.82倍(Illumina)。图d和e是通过和两个版本的HSV-1转录组比对计算Nanopore和Illumina数据集的转录丰度。

该研究阐述了direct RNA-seq在复杂病毒上的价值,确定RNA多样化在病毒感染中的机制将大大提高人们对宿主-病毒之间的交互作用的理解。

文章来源:

https://www.nature.com/articles/s41467-019-08734-9#Abs1

PacBio技术破译首个木兰类植物——牛樟基因组

木兰亚纲(Magnoliidae)隶属于种子植物中的被子植物,含超过9000个种和四个目(白桂皮目、胡椒目、樟目和木兰目)。来自我国宝岛台湾生物多样性研究中心和医学科技学院的多名研究人员首次聚焦木兰亚纲樟目樟科樟属的牛樟,结合PacBio、Hi-C、Illumina等技术平台组装出了其reference级别的基因组并进行注释。13种代表性的种子植物系统发育分析表明,相较于单子叶植物,木兰亚纲和双子叶植物之间有着更近的祖先;在木兰亚纲谱系间发生了两次全基因组复制事件(WGD):一次早于樟目和木兰目之间的分化,另一次在樟科的分化之间。樟属演化史中也有小片段复制和串联重复事件的作用,如月桂中萜合成酶基因亚家族的扩张造就了肉桂单萜和倍半萜的多样性。文章见刊Nature Plants,该期刊影响因子为11.471。芳香药用植物在人类历史上有着重要的工业和药用价值,而樟属更是分布于亚洲和南美在经济和生态上重要的常绿芳香树种。但是木兰类植物之间的谱系关系却一直有待明确。牛樟作为台湾特有的植物因其巨大的树干、芳香、耐腐的特性而有着不可估量的价值,同时也因为砍伐过度、种子发芽率低而濒临灭绝。正因如此,研究者组装出了染色体水平的牛樟基因组并结合其他10种被子植物和2种裸子植物指明了木兰亚纲的系统发育地位还解析了对开花植物基因组进化问题的疑惑。那么研究采取什么样的策略又是如何进行分析的呢?请随组学君一探究竟。
方法策略
采样(12龄牛樟)
基因组DNA提取(叶片)并测序①Illumina双端测序②PacBio测序(20kb SMRT文库)③构建Chicago和Hi-C文库并进行Illumina测序

RNA提取(花、二期花芽、未成熟叶、幼叶、成熟叶、幼茎和果实)并进行Illumina测序
染色体数目评估
基因组大小评估
基因组从头组装
基因预测和功能注释
基因组杂合度分析
重复元件鉴定
基因家族和同源类群推断以及蛋白域分析
系统发育分析
分歧时间估算
基因组共线性和WGD分析
R基因和TPS基因鉴定
结果概述
牛樟基因组的组装和注释

牛樟预估基因组大小为823.7± 58.2 Mb/1 C, 2n=24。本研究先用PacBio测序覆盖85×深度,reads N50为11.1kb,contig N50为0.9Mb,组装基因组728.3Mb。然后使用141×Illumina reads进行校正,接着用207בChicago’重组染色质和204×Hi-C双端reads进行染色体挂载。最终,整合基因组大小730.7Mb,含2153条scaffolds,占流式细胞仪估算基因组大小的91.3%,scaffold N50达到50.4Mb,占12条对应于牛樟染色体超过90%的pseudomolecules,且注释蛋白质BUSCO评估不低于89%。

基因组特征

杂合位点的空间分布高度可变,有23.9%的基因组每1kb仅有不到1个SNP位点,而有10%的基因组1kb有至少12.6个。杂合区域呈现随机分布,在第11号scaffold上达到最多20.2Mb(参见图1a),但在序列覆盖上较其他杂合区域更均衡,推断与选择性清除、同系繁殖或最近的种群瓶颈效应相关。同时,位于这些杂合区域上的基因在木质素生物合成过程与半乳糖代谢中尤为富集,可能在木质素-碳水化合物复合物的合成中起到作用(见图1b)。

图1 牛樟基因组杂合度(a.每100-kb非重叠窗口的杂合双等位基因的SNPs数目对最大的12条scaffolds作图,Indels被排除;b.使用PSMC方法推断的有效种群大小的历史;c.对每个100-kb非重叠窗口,分布从顶端到末端为基因密度、转录组和三类重复序列。红色T字母表示scaffold末端存在端粒重复簇,LINE表示长散在重复序列 )

转座元件和散在重复序列占到基因组的48%,最大的类别为长末端重复转座子(LTR RT)(25.53%),其次是DNA转座子占12.67%。在LTR RT中,又分别有40.75%为Ty3/Gypsy、23.88%为Ty1/Copia(图1c)。另外,LTR富集的区域较基因组其他区域平均多覆盖了35%,正是其在组装中被折叠导致了流式细胞仪和k-mer分析对基因组大小评估的差异。

再者,染色体水平的scaffolds在染色体中心区域呈现蛋白编码基因密度低、转座元件分布密度高的特点(图1c)。同时研究还发现了687kb核质体类DNA序列(NUPTs),其中96%是小于500bp的短片段。

系统发育分析

研究从上述13个物种中搜寻到了211个单拷贝直系同源数据集,整合为超级矩阵后使用最大似然法进行系统发育树的构建,结果显示木兰亚纲和双子叶植物支系形成姐妹类群关系(见图2)。该系统树在加入额外22个木兰类物种的转录组数据后仍然保持原来的拓扑结构,只是自展支持率较之前稍低。

最后,通过化石年代的标定,研究估算出木兰亚纲和真双子叶谱系于136.0-209.4Ma分化出来,这一结果也得到了最近其他研究结果的支持(图2)。

图2 基于13个物种中211个单拷贝直系同源基因构建的物种系统发育树(+和-旁的数字分别表示基因家族的扩张和收缩;括弧里的绿色数字表示分歧时间的估算;如无额外说明,所有自展支持率为100%)

共线性分析和全基因组复制

研究在72.7%的基因组中鉴定到16498个基因对分布于992个共线性区上,这些共线性区又有72.3%与基因组上超过一个位点呈共线性,也就是说:牛樟祖先发生了不止一次WGD事件(参见图3a)。染色体区域成对的广泛共线性以及每个区域与另外两个基因组片段显著、但非共线性的配对和两轮古代WGD事件息息相关(图3a)。

无油樟(Amborella trichopoda )是其它现生所有被子植物姐妹类群的唯一物种,自现生开花植物谱系最近的共同祖先分化以来没有发生过WGD事件的证据。因此,研究做了一个假设:两轮WGD事件是牛樟祖先在牛樟和无油樟分化之后发生的,两者基因组的共线性分析发现牛樟1-4个片段匹配上了无油樟基因组的单个区域,证实了上述假设(图3b)。

那么两轮WGD事件什么演化节点上发生的呢?研究者又估算了基因组内和种间的同源染色体Ks(同义替换)分布。牛樟基因组内重复片段峰值在0.46和0.76附近(图4a),正好和两轮WGD事件对应,研究据此来推断核型进化事件。洛矶山耧斗菜(Aquilegia coerulea隶属于其它所有现生真双子叶植物姐妹类群的毛茛目,其与牛樟直系同源的分析显示主峰在Ks=1.41左右(图4a),而耧斗菜基因内重复片段Ks=1,由此推断在牛樟和耧斗菜谱系分化后发生了独立的WGD事件。通过挖掘17个樟目和木兰目转录组数据的信息,研究发现樟科Ks分布的两个峰与牛樟Ks分布一致,对应上了牛樟祖先WGD事件的两个基于共线性的推断(图4b和图3)。而其他樟目和木兰目物种仅发现一个Ks峰,推断WGD事件在其演化过程中仅发生一次。在耧斗菜数据中观测到的Ks峰可能对应的是真双子叶植物和木兰亚纲分化后毛茛目内部的WGD事件(图4a)。

图3 牛樟基因组的进化分析(a.牛樟基因组中的637个共线性区基因组内关系示意图;共线性区由五色线条表示,代表祖先染色体组型;紫色区块表示分配到第一种颜色区域上的共线性区;b.牛樟基因组内第一种颜色区域和无油樟对应的共线性关系)

图4 牛樟和其他物种之间同义替换的密度图(a.牛樟与洛杉矶耧斗菜及两者间共线性区鉴定到的成对直系同源重复序列;b.樟科和木兰目基因组内成对重复序列Ks,虚线为在牛樟上观察到的Ks峰,棕色和灰色的线条分别为牛樟和其他樟科Ks分布)

特化木兰亚纲蛋白质组学

研究试图通过对蛋白家族(pfam)域进行注释、评估其在13个用于系统发育分析的种子植物基因组中的分布来鉴定牛樟特有的基因和蛋白域。分析发现,牛樟、真双子叶植物和单子叶植物之间有相当大的重合,说明三个谱系在分化后功能发生了显著的多样化。包括萜烯合成酶(tps)羧基末端结构域在内的蛋白质结构域的获得涉及植物蒸腾效率中的防御反应和富含亮氨酸的重复序列。有意思的是,研究者发现牛樟拥有21个EIL转录因子的拷贝,比之前报道的拥有最多拷贝数的香蕉基因组还多4个;同时EIL是通过激活乙烯应答因子(ERF)来启动乙烯信号应答的,牛樟中的ERF同样高度扩张。通过EILs的扩张来刺激ERF从而实现对下游效应的调控形成了牛樟特有的属性。

接下来,研究还评估了种子植物系统发育中直系同源类群的扩张和收缩。基因家族大小的演化在系统发育过程中是动态的,进化为牛樟的支系没有呈现出显著的扩张或收缩变化。GO富集分析揭示出牛樟不同的基因家族享有共同的功能或者单个基因家族经历了大规模的扩张。

R基因和TPS基因家族

牛樟基因组注释包含了387个R基因模型,其中82%都属于核苷酸结合位点上亮氨酸富集的重复序列(NBS-LRR)或者卷曲螺旋结构的NBS-LRR类型。在13个研究基因组中,牛樟在非栽培植物中含有最多数目的R基因。2465个NBS结构域的系统发育分析也显示基因家族内的分支在真双子叶植物、单子叶植物和木兰亚纲中是独立分化的。引人注意的一点是,牛樟NBS基因最分化的一支与双子叶植物NBS基因最保守的一支形成了姐妹类群。

图5 101个CkTPS基因的系统发育位置

对于牛樟基因组,最显著的特征是其庞大数量的TPS基因(CkTPS)。本研究在牛樟基因组中预测、注释了101个CkTPS,是迄今基因组中发现数目最多的。在加了两个木兰亚纲物种转录组数据并进行系统发育分析后在之前已描述过的7个种子植物TPS基因亚族中明确了6个亚族CkTPS的系统位置(图5和表1)。101个CkTPS基因有7个与催化形成20-碳异戊二烯类化合物的关键酶有关,另外94个很可能编码10-碳单萜合成酶、15-碳倍半萜烯合成酶以及其他20-碳双萜合成酶(表1)。

值得引起注意的是CkTPS基因系统树还解决了TPS-aTPS-b TPS-fTPS-g基因亚族内部樟科特异性的TPS基因分支。总之,分析表明不论是樟科起源之前还是之后,木兰亚纲TPS基因都处于不断分化之中。

最后研究还分析了不同染色体上TPS基因的分布密度。研究提到,TPS基因在染色体上是不均匀分布的,并且独立亚族中的成员聚类以串联重复的形式存在。在牛樟基因组最大的12条scaffolds上观察到了76个TPS基因,而7号scaffold包含了隶属于多个亚族的29个TPS基因,相比之下,1号scaffold仅含CKTPS-c的两个成员,再有24个CKTPS基因定位在其他小一些的scaffolds上等等。

结 语
本研究涉及的牛樟自19世纪以来遭到过度砍伐,追溯到900万年前,它的有效种群规模其实就在持续下降,这些不仅由于牛樟复杂的种群历史也与台湾的诞生、发展休戚相关,而且还和中新世晚期以及5-6Ma的造山运动脱不开关系。那么本研究最大的贡献在哪里呢?①首次进行了木兰亚纲代表物种的测序,并纳入13个种子植物代表物种的系统发育研究中。虽然金粟兰科和Ceratophyllacae (疑似金鱼藻科)尚无基因组数据可用,但是将牛樟和真双子叶植物归为姐妹类群对于真双子叶植物演化的比较基因组学分析有着重大的意义。②与之前同工酶分析结果相一致,本研究推断出了木兰亚纲发生的两次独立WGD事件的节点,其对病原体、食草动物及其共生交互作用的基因家族的扩张和多样性演化起到了举足轻重的推进作用。③对于核心真双子叶植物六倍体祖先起源的二倍体事件重建,牛樟基因组扮演重要的参考外类群的角色。

④六个被子植物TPS亚族分别进行的基因树拓扑结构揭示了牛樟祖先TPS基因的分化和基因功能,例如证实TPS-f基因亚族在所有被子植物最近的共同祖先中存在但在禾本科中是缺失的等。

⑤研究鉴定了牛樟基因组中101个CkTPS基因并在12条染色体对应的scaffolds上是不均匀分布的且包含来自于多个亚族的基因簇。

简言之,牛樟基因组的问世为揭示其他木兰亚纲物种的遗传多样性和演化打下了坚实的基因组学基础,并对开花植物基因组进化和分化提供了更完善的见解,同时对这种具有重大文化、经济价值的阔叶林物种的多样性保护也起到了作用。

史上最全的ONT碱基识别软件大比拼

大家都知道Oxford Nanopore Technologies(ONT)是基于单条DNA链通过纳米孔引起电流改变来进行测序的,但是这个过程可不简单!因为噪音信号和随机数据都会干扰碱基读取;而纳米孔的电阻受位于孔最狭窄处的5个核苷酸中的碱基(R9.4)决定,从而有多达45即1024种电波模式,一旦出现碱基修饰如5mc,那么就会有55即3125种电波,因此碱基识别软件(basecaller)在这个过程中至关重要。

现在的basecaller都使用的神经网络算法,且需要用原始的、即包含碱基修饰的DNA进行模拟学习;识别的准确性则由read精确度(单条read对参考数据的相似度)或一致性序列准确度(重叠reads构建的一致性序列对基因组同样位点的序列相似度)来评估,一致性准确度通过增加读取深度来提高。

事实上,read精确度和一致性准确度没有直接的相关性,在错误产生是随机的且读取深度足够深的情况下,低精度的reads也可以产生完美的一致性序列,但如果包含系统错误而不考虑读取深度的话,即使是高精度的reads也可能产生一致性准确度差的序列。

一致性准确度对于基因组组装等需要高深度数据支持的应用来说是关键的影响因素,对读取深度需求较低的应用,则是read准确度显得更重要。

 

鉴于R9.4型纳米孔自2016年十月份延用至今,来自澳大利亚莫纳什大学和伦敦卫生与热带医学院的研究人员想到通过量化R9.4适配的各种碱基识别软件的表现来帮助学者们充分利用ONT电信号,特别是给不知使用自定义数据模拟软件还是新碱基识别软件的用户以参考。研究成果最近发表在了bioRxiv上,和组学君先睹为快吧!

//方法流程//

研究者先是测试了ONT旗下四款碱基识别程序——Albacore、Guppy、Scrappie和Flappie并运行和R9.4reads配套的所有可用版本,另外还测试了一款仍处于研发状态的、基于深度神经网络算法的碱基识别软件Chiron等。

接下来是进行自定义模式模拟,研究者使用Sloika工具包——一种可生成用于Guppy的模式神经网络模拟工具,模拟的数据来自30个不同谱系的肺炎克雷伯菌、10种肠杆菌以及10种蛋白菌。对收集到的5629714条原始reads进行过滤、处理,如去掉短的和低质量的reads等,然后将这些数据修正为不同批次的定量模拟数据,再只用上述软件进行不同模式下的模拟运算。

研究者使用了MinION R9.4芯片来测序原始的肺炎克雷伯菌DNA以测试basecallers的性能,同时也准备了同样样本的高质量Illumina数据以作参考。研究还提取了肺炎克雷伯菌的ONT数据,在使用Guppy进行识别后比对到参考基因组上并排除极低质量和外源reads。除此之外,为了拓宽数据集的范围,研究用到了R9.4和R9.4.1芯片读取肺炎克雷伯菌之外的六种菌种作为测试数据集。

最后就是进行read精确度和一致性准确度的分析了,为评估read精确度,研究使用minimap2将获取到的reads识别数据集比对到参考基因组上计算“BLAST identity”,然后使用Rebaler从这些数据集中生成一致性序列再进行分析。还使用了NUCmer将组装结果比对到参考基因组上对不同的一致性序列错误类型进行分类。

考虑到在组装下游一般都会使用Nanopolish进行数据打磨,那么basecaller的选择是否还那么重要呢?为了得到这个问题的答案,研究者对每一个最终的组装结果都进行了Nanopolish打磨和一致性评估。

//结果概述//

默认模式性能

①Albacore

2017年4月释放的Albacore v1.0.1在均聚物识别上表现良好,而2017年8月释放的Albacore v2.0.1则省去了事件切割(event-segmentation )步骤直接识别原始信号(图1),这之后的版本表现相对稳定,read精确度为Q9.2,一致性准确度为Q21.9。运行速度约为120000bp/s。

②Guppy

表现与Albacore相当,不过最新版本Guppy v2.2.3在默认模式下的read精确度(Q8.9)上不如一致性准确度(Q22.8)。在使用flip-flop模式时,一致性准确度相当(Q23.0),但read精确度更佳(Q9.7)。相较其他软件,其在GPU加持下运行速度最快(约1500000bp/s)

③Scrappie和Flappie

Scrappie表现最不理想,Flappie的一致性准确度不如read精确度且速度慢(约14000bp/s)。

④Chiron

Chiron在read精确度上性能不佳(见图1),但Chiron v0.3在默认模式下得到了最高的一致性准确度。最新版v0.4.2亦不如预期。另外,其运行速度为最慢,约2500bp/s。

图1 不同日期释放的每种basecaller版本评测的一致性准确度、read精确度和识别速度。准确度以qscore(Phred quality scores)的对数刻度来表示,Q10=90%,Q20=99%,Q30=99.9%等。

自定义模式性能

和默认的模式相比,在custom-Kp模式(使用了50种菌种的模拟数据)的基准设置下运行Guppy v2.2.3生成的read精确度有了一定的提升(Q9.5),而一致性准确度提升更加显著(Q28.5)(参见图1)。该结果证实了使用分类阶元特异性模拟数据的优势。另外,两种模式下的运行速度相似。

而custom-Kp-big-net在两个方面均有更好的表现,尤其是一致性准确度达到了Q31.6,表明更复杂的神经网络也能够升级结果,只不过速度要略逊一筹。为观察结果对其他基因组是否也适用,研究将可用的Guppy模式也运用于其他的数据集上(参见图2)。结果显示,flip-flop模式相较于默认模式在所有的基因组中都表现更好。而在所有的案例中,custom-Kp-big-net模式得到了准确度最高的read和序列一致性。

custom-Kp模式和custom-Kp-big-net模式都没有使用Guppy的flip-flop模式中存在的新神经网络结构,可以推测:基于flip-flop模式且在自定义模拟数据的模拟下既可以从flip-flop模式也可以从custom-Kp模式中体现优势。

图2 使用Guppy v2.2.3在默认RGRGR模式、flip-flop模式和两种自定义模式下运行各基因组得到的read准确度和一致性准确度。

一致性错误描述

为了探索不同basecallers对各种一致性识别错误的影响,研究者量化了肺炎克雷伯菌基准基因组在Dcm甲基化位点、均聚物以及其他位点上的错误数(详见图3)。结果发现,由于模拟数据中缺少同类信息,ONT basecallers在使用默认模式下对Dcm甲基化位点的识别都不太理想,一致性错误率达到了约0.4%;相反,因模拟数据集中包含了Dcm信息,Guppy v2.2.3在自定义模拟模式下几乎没有Dcm错误(约0.002%)。

除了Dcm模体,错误的均聚物长度占据了大部分的错误(图3)。而ONT的升级则在Albacore上充分体现,不仅一致性准确度有所提升,均聚物的错误也从v.8.4的0.53%降到了v2.3.4的0.13%。同时,最新的Guppy v2.2.3也展现了较好的性能,将均聚物错误降到0.07%。最后,结果还显示对于均聚物,其在custom-Kp模式下表现略逊于默认模式,而custom-Kp-big-net模式则表现最好。

图3 肺炎克雷伯菌基准基因组在不同basecallers下的一致性错误

Nanopolish性能

Nanopolish包含对Dcm甲基化和均聚物特异的逻辑算法体系,能够修正原始信号数据中的一致性序列错误。研究者发现,Nanopolish除了custom-Kp-big-net模式之外,几乎在所有案例中都可以提升一致性准确度(见图4)。分析表明,前Nanopolish一致性准确度和后Nanopolish准确度相关,R2=0.580,这表示即使使用了Nanopolish,basecaller的一致性准确度仍然非常关键。

图4 对肺炎克雷伯菌基准基因组进行Nanopolish之前(红)和之后(蓝)的一致性准确度

因ONT测序平台的便携、经济,其在食源性或其他爆发性病原体DNA碱基替换的检测中有着潜在用途。在本研究中,研究者使用Guppy v2.2.3在custom-Kp-big-net模式下对组装基因组检出的最小替代错误是337个,这差不多是细菌和病原体基因组之间预期的真实SNPs数目的十倍,因此在这类应用中会造成不可估量的高假阳性率。该问题可以通过调用SNP识别策略得到控制,不过基于组装的序列比较则需要识别准确性更高才行。

//研究要点//

1. 从肺炎克雷伯菌基准数据集的read精确度和一致性准确度来看,碱基识别软件Guppy v2.2.3在自定义模式下进行数据模拟得到的结果最佳。这种优越的性能主要取决于对Dcm甲基化的正确处理。因不同物种之间的差异,本研究结论更多的代表一种趋势,那就是当使用相同或亲缘关系足够近(含相似的DNA修饰)的物种DNA进行模拟时,原始的DNA碱基识别准确度是最高的。

2. 对大多数basecallers而言,在默认模式下用于模拟的物种及DNA类型的信息并不是公开的,研究建议软件开发者提供多种模拟模式以供用户选择与研究对象最匹配的一款,同时在条件允许的情况下推荐自定义模式以将碱基识别准确度最大化。以本研究为案例,如神经网络更全面,custom-kp-big-net模式能够获取更精确的结果,但是可能会要耗时一些。

3. ONT自诞生以来在产量和精确度方面均在不断升级,但是仍有上升空间。本研究的模拟结果中最佳的一致性准确度为Q32.2(99.94%的一致性),可换算为5Mbp的基因组中约有3000个错误,其中很多属于会导致SNP识别假阳性的替代错误。对于完美的细菌基因组而言,标准可预期为Q70即10Mbp中一个错误。目前可采纳的升级方式可以从技术、试剂、basecaller升级或组装后polish工具等入手。因此,在达到预期目标之前,使用Illumina数据混合组装或polish可能对于高精度序列仍是不可或缺的。

未来组2019成果首秀!五谷“元老”庆丰年

路人甲:组学君组学君,今天和大家聊点啥?组学君:硕鼠硕鼠,无食我黍……路人甲:讲诗经?

组学君:错,今天我们聊“黍”!稻、黍、稷、麦、菽五谷中的黍!

路人甲:那不就是糜子,又叫黄米嘛!

组学君:对!民以食为天,从古诗歌中我们了解到在西周末年到春秋初期糜子就是关乎生计的重要作物,而其最早的栽培史可以追溯到距今七千到一万年前,五谷“元老”名副其实呀!

路人甲:

组学君:最近未来组参与合作的糜子基因组近完成图发表于Nature Communications,当然要细说一番啦!

路人甲:快进入正题吧!

多技术结合打造糜子基因组近完成图[1]

Chromosome conformation capture resolved near complete genome assembly of broomcorn millet (Nature Communications, IF: 12.353)

糜子(Panicum miliaceum L., 2n = 4 × = 36)作为一种古老的作物,具有生长周期短(~60–90d),耐盐碱的特性,尤其是极端抗旱,其耐旱能力甚至比蒸腾系数低于高粱、玉米、小麦的谷子还要强。这种对非生物胁迫的超强耐受性机制一直为研究者所好奇。近日,糜子高质量参考基因组终于问世,该研究由中国农业大学赖锦盛教授团队主导,相关研究文章发表在Nature Communications杂志,从分子角度挖掘了潜藏这一机能的秘密。武汉未来组为该研究提供PacBio测序及组装

技术方案 
 
实验材料:Longmi4种子播种14天后(25 °C,暗处理),采集叶片,液氮速冻,提取高质量基因组DNA。

测序方法: 

PacBio、Illumina、BioNano、Hi-C

组装策略:使用Falcon对~170×PacBio长读长数据进行基因组组装,用~116×Illumina短读长数据进行polish,结合~235×BioNano数据辅助构建scaffold,再利用~140.2×Hi-C数据进行模拟染色体构建。

研究结果 
01 糜子基因组组装分析

研究首先用K-mer分析预估了Longmi4基因组大小约为谷子的两倍,证实了其四倍体化事件,然后基于PacBio和Bionano数据组装出了约848.4Mb的基因组,组装指标Contig N50 达2.55Mb,scaffold N50 达到8.24 Mb,包含了18条super scaffolds,覆盖了~95.6%的基因组(见Table 1)。接着研究通过全基因组测序和RNA-seq数据分别对组装好的基因组进行比对以评估组装的质量,比对率达到了99.6%和91.5%,说明组装质量非常好;经BUSCO评估发现基因完整度高达98.0%,标明组装完整度和准确度都非常高!

最后,精益求精的研究者还将Hi-C数据纳入其中以将scaffolds锚定、排序,通过深度约140.2×的覆盖将444条scaffolds装出了19条超长scaffolds以进一步分析。研究发现Longmi4染色体内Hi-C互作矩阵中出现明显的反对角图案(Fig. 1a),反映出所谓的染色体Rabl构象(间期染色体长短臂平行折叠)。基因组比较分析揭示了糜子和谷子之间的2对1的共线关系(Fig. 1b)。与玉米基因组不同,糜子的染色体在基因组四倍化后并未发生明显的融合现象,但染色体内的重排尤其是倒置现象在糜子基因组中广泛存在。例如在pm5和pm8两条染色体上,就有两个大小约为11.8 Mb和8.9 Mb的区域发生了倒置(Fig. 1b)。同时,研究者还发现pm3、pm10、pm14和pm15号染色体末端的一致性,它们之间是两对同源染色体,由此推测这对染色体末端之间的交换应是发生在糜子基因组四倍化之前。

Fig. 1 Hi-C技术辅助构建Longmi4模拟染色体

(a)200kb分辨率下的Hi-C互作热图

(b)糜子(pm1–pm18)和小米(si1–si9)基因组比较

02 基因注释和基因家族分析

通过结合~68.6 Gb的RNA-seq数据及一些近缘物种的蛋白序列,研究者在糜子基因组中预测出63,671个蛋白质编码基因,其基因数量几乎是谷子基因数量(34,584)的2倍。基因组比较分析发现,糜子和谷子共有19,609个基因,在这些基因中,有16,884(~86.2%)个基因在糜子基因组中都拥有两个同源拷贝,表明在发生全基因组复制事件(WGD)之后糜子基因组中丢失的基因较少(Fig.2 )。

Fig.2 以小米基因组为参考分析糜子基因组中的基因丢失与保留

通过对糜子基因组的深度挖掘,研究者揭示了糜子基因组中与抗生物胁迫和非生物胁迫的可能相关基因。共鉴定出493个含有NB-ARC结构域(可能与抗病性有关)的基因,其中20个基因(7个基因家族)在糜子中具有特有。与珍珠稷类似,糜子中的含NB-ARC结构域的基因具有一定的偏向性,主要集中在pm13和pm18染色体末端(Fig.3)。同时,糜子作为典型的耐旱植物,在其基因组中鉴定出了15 个ABA 或WDS应答基因。

更重要的是,与一般ABA基因的诱导表达不同,这15个ABA基因中有四个基因存在持续表达的表达谱,充分暗示了这些基因可能参与了糜子的极端非生物胁迫抗性。

Fig.3 NB-ARC结构域基因在糜子的18条

模拟染色体上的分布情况

03 Gypsy元件的近期大爆发

研究者发现在糜子基因组中~54.1%的序列为重复序列,其重复度在谷子和珍珠稷之间(二者基因组重复序列含量分别为~46.8%和~68.0%)。与其他谷物类似,LTR反转座子尤其是Gypsy超家族是糜子基因组中最主要的重复元件(Gypsy超家族占比~31.4%)。研究者对包括糜子在内的5种作物中的GypsyCopia家族的插入时间进行了分析,发现玉米和高粱中Gypsy元件在近期有一次爆发,其他三个物种的爆发时间较早,但都发生在距今~100万年以内(Fig.4a)。相比之下,糜子基因组中的Copia家族的爆发时间要早得多,大约发生在距今200万年前(Fig.4b)。

Fig.4 糜子 (P. miliaceum)、谷子 (S. italica)、珍珠稷 (P. glaucum)、高粱 (S. bicolor)、 玉米(Z. mays)基因组中的Gypsy (a) Copia (b)的插入时间分析

04 亚科大类进化分析

研究利用组装的高质量的糜子基因组、新发表的珍珠稷基因组和Dichanthelium oligosanthe基因组进行联合分析,构建了黍族系统发育树(Fig.5)。系统发育分析揭示了糜子的异源四倍体化发生在591万年之内,而糜子和谷子的分化大约发生在1310万年前。

Fig.5 黍族的系统发育分析

结 语

本研究所报道的高质量糜子基因组序列不仅对理解糜子基因组四倍体化后的动态进化具有重要意义,而且对今后的糜子分子育种也有一定的参考价值,并将促进黍属植物与其他作物的比较基因组研究。

同时,研究结合了三代长读长全基因组测序技术、二代短读长测序技术、Bionano光学图谱、Hi-C染色体构象捕获等技术优势,获得了高质量的糜子基因组,说明了多组学平台的结合对于复杂程度高的作物具有显著的优势。近年来,不断有高等植物通过该手段获得了良好的组装质量,如武汉未来组2018年参与的玉米基因组文章[2]借助PacBio辅以Illumina、Bionano技术登上了Nature Genetics,Bionano光学图谱已成为破解作物遗传密码的首选秘密武器。还有借助Hi-C技术登上Science的小麦[3]和其他作物如日本晴[4]、苦荞[5]等,可以说杂合、重复已不再是基因组难题,对于它们基因组的解析在分子育种、农业生产以及特性优化方面都有着不可估量的价值。

 

延伸阅读:

Nature Genetics | 三代测序揭示玉米种内存在广泛的基因结构变异

ONT + 光学图谱 = 染色体水平的植物基因组

Nanopore 和 BioNano DLS 双剑合璧实现染色体级别高粱基因组

未来组三代基因组项目再出新篇!“英雄树”木棉基因组草图首发

小麦基因组草图到精细图的利器:长读长测序+光学图谱

参考文献

[1] Junpeng Shi, Xuxu Ma, Jihong Zhang, et al. Chromosome conformation capture resolved near complete genome assembly of broomcorn millet[J]. Nature Communications, 2019.

[2] Sun, S. et al. Extensive intraspecific gene order and gene structural variations between Mo17 and other maize genomes. Nature Genetics (2018).

[3] Avni, R., Nave, M., Barad,et al.  (2017). Wild emmer genomearchitecture and diversity elucidate wheat evolution and domestication. Science, 357(6346), 93-97.

[4] Dong Q , Li N , Li X , et al. Genome-wide Hi-C analysis reveals extensive hierarchical chromatin interactions in rice[J]. The Plant Journal, 2018.

[5]Zhang L , Li X , Ma B , et al. The Tartary Buckwheat Genome Provides Insights into Rutin Biosynthesis and Abiotic Stress Tolerance[J]. 分子植物:英文版, 2017, 10(9):1224-1237

Cell| PacBio升级解析人类基因组结构变异

继登顶Nature后,近日,PacBio文章再度收录于顶级期刊Cell上!第三代测序技术的实力日渐凸显,这次又在哪个领域有所突破呢?原来,华盛顿大学医学院的Peter A. Audano等人瞄准了当前人类基因组注释中的缺陷,试图用长读长技术修复人们对SVs认知的误差,下面和组学君先睹为快吧!
 一句话搞定?
长读长测序技术助力人类SVs分类解析并促进短读长数据对其进行基因分型的算法研究,明确了SVs在人类基因组研究中的重要作用!
 一张图说明?

 一分钟看完? 

为了优化人类结构变异(SVs)信息,研究者对15个人类进行了长读长测序并且分析了SVs,最后找到了99604个插入、缺失和倒位,其中2238个(约合1.6Mb)为已揭示的人类基因组中所共有,另外还有13053个(约合6.9Mb)在大多数人类基因组中找得到,证实参考基因组中含次要等位基因或者错误。附加的440个基因组分型结果证实了特异染色质中最常见的SVs被解析出来。研究发现:人类染色体最末端的5Mb中所包含的SVs是其他位置的9倍之多,其中55%的可变数目串联重复序列都映射到该区域。研究者鉴定出影响编码和非编码调控位点的SVs,优化了注释和对功能变异的解析,为精细人类参考基因组图谱构建了框架并为捕获等位基因多样性提供了重要信息。

如果有时间慢慢读,请继续……

材料与数据

该研究使用PacBio长读长测序技术对11个人类基因组进行测序,然后添加了两个之前已由本研究测序过的葡萄胎(胎盘绒毛发生良性病变的胚胎)CHM1和CHM13。另外,还加入了由未来组参与的华夏一号HX1以及同样是16年发表的亚洲人AK1基因组数据(见表1)。

深度挖掘SVs
对每一个人类基因组样本使用SMRT-SV鉴定、组装出50bp或者和GRCh38关系更密切的插入、缺失以及倒位SVs,且排除不可靠的SVs calling结果,如具有密集串联重复或间隙结构的着丝粒周围区域。最后平均每一个样本中鉴定到22755个SVs,且将它们融合为一个99604个非冗余SVs数据集(见表1和图1A),并分为四大类别:共有的、主要的(≥50%但并非所有样本中存在)、多态的(多于一个但<50%样本中存在)以及特异的SVs。

图1A 使用非冗余策略融合每个样本的变异成一个数据集

和预期一致,非洲样本多态性最丰富,平均每一个非洲人样本贡献了11.1%的特异性SVs,而非非洲人样本平均是5.6%,可以推断加入非洲样本可以将SVs识别翻倍(见图1B)。

图1B 每个样本中每一个分类类别的变异数量

非冗余数据集起初增长急剧,但是随着样本增加,增速放缓,这也说明这15个人中共有的SVs比例较高。同样,共有SVs数据集一开始降低急剧,随着样本增加逐渐平缓,所有样本中共有的SVs是2238个,且携带每一个SV的样本比例也呈现类似模式,共有的SVs增加了100%。同时,研究鉴定出15291个主要的SVs,表明当前的人类参考基因组在这些位点上中也含有次要等位基因或者错误。相较于多态SVs,主要变异数量更多且倾向于在重复DNA(约占80%)中富集。当与GRCh38基因组进行缺失比较的时候,如分析样本中共有SVs比例更高则定义其为插入SVs(见图1C)。

图1C 插入(INS)、缺失(DEL)和倒位(INV)三种变异在各类别中发现的频率

有意思的是,研究还和Illumina测序的人类基因组数据SVs识别结果进行了比较,发现本研究使用长读长测序技术挖掘出了87.3%的SVs在之前的二代测序数据中没有找到,尤其是插入SVs,有93.5%是之前不曾鉴别到的,其次是缺失SVs,新发现的比率也很高。这再一次证明:第三代长读长测序技术较之第二代短读长测序技术能够鉴定到更多的SVs。当然研究还将结果和人类基因组结构变异联盟做了比较,因篇幅限制就不详述了。

附加人类基因组的基因分型
为更好的理解SVs的群体分布,研究者优化了基因分型工具SV genotyper并应用于使用Illumina数据构建的440个人类基因组上(图1D)。结果显示,在至少95%的样本中,55.1%的SVs成功地进行了基因分型,92.6%的SVs成功地进行了一半或更多的基因分型。在那些能够成功进行基因分型的基因中,我们观察到至少一个附加人类基因组中有97.2%的SVs。这表明绝大多数SVs代表真正的人类多态性,而不是个体变异或体细胞伪影。

图1D 440个人类基因组样本中可分型的SVs其不同类别的出现频率

与预期的一样,在共有的和主要的SVs中分别找到了95.4%和66.7%的次要等位基因。在本研究中发现的507例(0.74%)共有和主要SVs中,研究者只观察到替代等位基因,却未观察到任何人类参考基因组序列。对于这些基因位点,人类参考基因组要么代表一个极端次要的等位基因(<0.2%),要么就是错误。

SV密度和染色体分布

SVs在基因组中是非随机分布的,研究者观察到在重复片段富集的染色体臂末端5Mb区域内,SVs呈现明显的偏倚(见图1E),并且推断亚端粒区域,SVs的密度是其他区域的9倍!共有SVs偏倚要小一些,但仍有3倍的密度差。

图1E  将SV划分为500Kb的bins并在不同染色体臂距离上进行聚类

当研究者在人类染色体中发现这一现象时,这些SVs却不是均匀分布,尤其是染色体长臂端更倾向于出现SVs亚端粒聚集,不过5号、19号以及X染色体是例外。研究者为了深入分析SVs偏倚这一现象,对其重复类型分类检验其亚端粒上的富集情况。他们观测到:SVs在VNTR上的富集密度是其他区域的4.8倍,其次是STRs(短串联重复区域),为2.9倍(见图2A)。

图2A SVs在STR和VNTR位点上的分布密度

虽然不同染色体富集情况不一,但是相较于短臂,人类染色体长臂上SVs普遍呈现出更宽区域的VNTR富集(图2B)。

图2B VNTR在不同的染色体臂距离上的聚集比率

另外,研究者还观察到双链断裂和VNTR密度之间的显著相关性,强烈的暗示了容易出现双链断裂的区域和VNTR形成之间的关系(图2C)。

图2C STR和VNTR数和双链断裂数相关性

基因的和潜在的调控SVs

接着,研究者又将共有的和主要的SVs和RefSeq注释结果交叉分析,解析了86个影响编码序列的事件、47个UTRs(非翻译区)事件、7417个内含子或任何基因2Kb空白区域中的事件。另外,还特意鉴定了1033个影响推断的非编码调控序列事件,本研究中定义为注释了的DNase I hypersensitive、H3K27Ac、H3K4Me1以及H3K4Me3位点的联合(见表2)。

这些事件中的许多嵌入了GCRICH或低复杂度DNA的区域,并可能影响基因结构。以图3A为例,在UBEQ2L1的5’端,研究者鉴定了一个1.6 kb的插入,主要由94 bp 富含GC的序列附近的二核苷酸和三核苷酸CACA重复单元组成。插入的断点精确地map到5’ UTR的第一个碱基,很可能扩展了UBEQ2L1启动子的长度。富含AT的序列照样可以被解析,例如458 bp重复元件在载脂蛋白APOOL 3’UTR内map上(见图3B)。

图3 缺失的基因或调控序列(部分)

优化mapping和SVs挖掘
基因组注释的优化以及对人类单倍型结构差异的理解深度对SVs的发现、解析具有重要的影响。在30个Illumina WGS样本中,研究发现如果将SV contigs添加到人类参考基因组及其替代contigs上,可以找到之前2.62% unmapped的reads,且有1.24% map到这些contigs上的reads提高了mapping质量。甚至,这些新map到的reads促使了插入SVs间SNVs和插入缺失的发现。例如,研究者使用GATK HaplotypeCaller鉴定到21969个特异变异,含68656个替代等位基因。通过短读长测序技术或者简单的线性参考基因组无法确定这些SVs,当缺失的序列映射到了编码序列时,这之间的差异直接影响对SVs的解析。举例说明,研究者鉴定到了FOXO6 exon 2上200bp的插入,而这200bp的片段与海马记忆加强和树突状脊柱密度紧密相关(见图4C)。恰恰这一片段在RefSeq以及Ensembl基因注释上都是缺失的,第二个(及最后一个)外显子在该插入的位点被分离了:RefSeq将外显子结合到一个0bp的内含子以及1391bp片段的第三位外显子上,而Ensembl则将它们接合到一个1bp的内含子以及477bp片段的第三位外显子上。本研究分析发现包含这200bp片段的序列形成了一个连续的编码外显子,加了67个氨基酸到ORF(开放性阅读框)上,相较于RefSeq注释,改变了基因终止密码子的位点。基因组突变频率数据库(gnomAD)报道了FOXO6上发现的7个功能缺失(LoF)的SVs,本研究通过纠正FOXO6阅读框,将两个推断的LoF变异正名为同义SNVs,还有一个更正为3’UTR中的SVs(见图4D)。

图4 纠正FOXO6阅读框

除了上述分析内容,本研究还分析了共有的和主要的等位基因SVs的特性、偏向性的GC组成并对人类参考基因组进行了补洞,对SVs进行了表达分析等,感兴趣的话可以阅读原文一探究竟。

总之,文章有如下几大亮点:

1. 测序注释了99604个常见的人类结构变异;

2. 发现了55%的可变数目串联重复序列(variable number of tandem repeats, VNTRS)映射到染色体末端,经分析其与双链断裂有着密切关联

3. 发现长读长测序技术能够鉴定到更多的SVs,尤其对于编码序列,SVs识别更加准确

4. 完善了参考基因组并为人类泛基因组研究丰富了多样性

原文内容博大精深,详情请点击原文链接

https://www.cell.com/cell/fulltext/S0092-8674(18)31633-7

相关阅读:

https://www.grandomics.com/research/h_x_w_r_sv/

SVs识别哪家强?PromethION为您揭晓

一滴蚊子血克隆出恐龙太扯?嗯,克隆出这只蚊子没毛病!

经典好莱坞大片《侏罗纪公园》相信小伙伴们还记忆犹新,斯皮尔伯格拿一只唯一不会吸血的还是只公蚊子当道具实在是功课没做好——不过也有技术宅质疑,如此微量的血液真的足够克隆出一个庞然大物?那么问题来了,基因组测序到底需要多少样本量?科学怪蜀黍要克隆出来自侏罗纪的蚊子,一只不够吧?
PacBio告诉你——完成全基因组测序,一只蚊子的DNA够够的了!

近期 Wellcome Sanger 研究所和 PacBio 公司的科学家公布了一种新的利用低起始量样本DNA进行全基因组测序的protocol,仅输入100ng DNA即获得了高质量的科氏按蚊基因组de novo组装效果。下面让我们通过预印版文章抢先来看看这一protocol到底特别在哪里吧!

材料:一只雌性按蚊(Anopheles coluzzii

HMW DNA提取:

在200μl PBS中捣碎虫体,使用Qiagen MagAttract HMW kit (PN-67653)试剂盒提取基因组DNA,并作如下调整:将Buffer ATL替换为200ul 1X PBS;在组织匀浆和孵育之前将PBS与RNAse A、蛋白酶K、Buffer AL混合;孵育时间缩短至2h;轻击管壁混匀;清洗1min;提取过程中转移DNA均使用宽口tips;

这些调整与10×genomics基因组提取方法一致,目的在于获取>50Kb的基因组DNA分子。结果显示,通过此法获得的gDNA总量约为250ng,且DNA长度主要分布在~150Kb左右。但经过低温运输过程之后,部分gDNA发生降解(Fig.1);使用Qubit fluorometer检测DNA浓度,随后使用其中的100ng DNA进行文库构建。

文库构建:

1 使用SMRTbell Express Prep kit v2.0构建SMRTbell文库(由于大部分基因组DNA已经片段 化至20Kb以上,因此省略gDNA剪切步骤)
2 对100ng DNA进行第一次酶促反应,除去单链突出端
3 随后用对DNA链进行损伤修复;
4 双链DNA末端加A;
5 于20°C条件下连接含T末端的SMRTbell接头(此过程耗时60min);
6 使用AMPure PB bead 进行SMRTbell文库纯化:
①首先使用0.45X AMPure②随后使用0.80X AMPure;
③最后使用FEMTO Pulse及Qubit Fluorometer对最终的文库进行浓度及大小检测(Fig.1)。

Fig.1 Anopheles coluzzii gDNA及最终文库

SMRT测序:

测序引物版本v4.0;聚合酶v3.0;使用Sequel系统测序,测序试剂版本v3.0,运行时间1200min;软件版本v6.0;共测序3 个SMRT Cells。

数据产出平均25 Gb /SMRT Cell (20 h movies)

组装指标:使用FALCON-Unzip软件进行de novo组装,contig N50 3.5 Mb,基因完整度大于98%

 操作要点: 
  • 基因组提取过程轻柔操作,保留大片段高质量DNA;
  • 建库过程不经历DNA剪切及片段筛选过程,避免剪切、纯化等步骤对DNA带来的损失,尽可能的保留DNA。

怎么样?这么省样品的建库测序方法你get了吗?据悉,这一新的官方protocol将于2019年2月正式发布,届时,PacBio必将在一些典型的涉及低DNA起始量的领域大展拳脚,如小生物膜的宏基因组分析、穿刺样本中分离的DNA样本、及需要有限扩增的靶向测序和单细胞测序等应用。

武汉未来组早在2016年就已经成功搭建基于PacBio的三代测序平台,并于2017年9月成功搭建Nanopore测序平台,一直致力于第三代测序技术的应用和推广,应用三代测序的合作研究成果多次登上Nature Genetics、Molecular Plant及Nature Communications等国际知名期刊。三代测序首选的合作者是谁?当然是未来组啦!

参考文献:

Kingan, Sarah, et al. “A High-Quality De Novo Genome Assembly from a Single Mosquito using PacBio Sequencing.” bioRxiv (2018): 499954.

Science首秀!看Nanopore如何杠上拉沙热病毒

过去三年,

有一种病毒疯狂肆虐

几内亚(科纳克里)、利比里亚、

塞拉利昂和尼日利亚部分地区!

2018年,

更是一发不可收拾!

截止2018年3月,仅尼日利亚

疑似病例就达1495例!

不排除其它西非国家存在的可能…

拉沙热听过?

一种主要通过动物传染的病毒性出血热疾病,

一旦感染,人畜共患!

Fig.1 2018年1月1日-3月18日拉沙热确诊病例分布图。(A) 受影响国家;(B)测序样本来源(橙色标记)

拉沙热病毒(LASV)是一种RNA病毒,且基因组高度可变,L区段(基因组大片段,编码RNA聚合酶和锌结合蛋白)和S区段(基因组小片段,编码糖蛋白和核蛋白)的种间核苷酸变异高达32%和25%,很难用短读长扩增子测序的方法准确检出。

怎么办?不忍看到水深火热中的西非人民继续遭罪,英格兰公共卫生署的Liana Eleni Kafetzopoulou及其同事们想到了纳米孔测序技术。对了,就是一种比U盘大不了多少的MinION纳米孔测序仪,通过对36个基因组及120份临床样本进行实时宏基因组测序分析,帮助他们揭示了LASV的多样化及其与早期发现的毒株的系统发育相关性,研究成果于2019年1月4日登上Science杂志——学者仁心,大大的赞!

Fig.2 样本测序时间轴

研究者调取了120份LASV-阳性临床样本,7周内搞定测序(Fig.2)。为了弄明白产生病毒间的亲缘关系,研究者将basecall reads和原始信号数据映射到参考序列上,使用Nanopolish进行突变体检测。对于非人源序列,研究者用canu进行了de novo组装,平均每个样本中包含LASV序列4.26%—42.9%,组装出了可以对91个样本中至少一个直系同源片段进行系统发育构建的矩阵。

研究者还留了个心眼,用Illumina测序对14个SISPA文库测序进行验证,发现Nanopore与Illumina序列高度一致,所以三代测序就是这么牛,靠谱,还走哪带到哪(Table 1)。

研究采用Centrifuge软件对110号样本进行宏基因组分类,鉴定出其中0.10%的reads源自甲型肝炎病毒,20×的测序深度达到了74%的基因组覆盖率。在同一个样本中,LASV reads占到了0.83%,达到了96%的基因组覆盖率。说明啥?嗯,哪怕多个长得很像的RNA排排站,Nanopore都能火眼金睛,一眼看穿——那些作为混合感染存在的病毒也无一例外!

那么,2018年尼日利亚拉沙热爆发的分子流行机制又是怎么回事?

别急,研究者使用生成的LASV序列及一些已有序列构建了拉沙热病毒系统发育树,真正的寻根溯源!基于S区段的最大似然法构树显示:2018年的LASV毒株都归属于尼日利亚LASV变种,尤其是Ⅱ型和Ⅲ型(Fig.2)。这个结果与基于L区段构建的系统发育树仅有7个毒株不一致。最后,真相大白:啮齿动物宿主污染是2018年拉萨热爆发的主要原因。所以奉劝小伙伴,没事零食别到处乱扔……

所以说,纳米孔测序技术在宏基因组学研究及疾病传播研究中的价值不容小觑,再加上实时测序、快速分析等其他技术无可比拟的优势,真的可以造福人类。这个研究缓解了人们对拉沙热在人际间广泛传播的恐惧,使公共卫生资源得到了合理分配,还指出LASV防治重点是加强社区鼠类控制、环境卫生和食品储存安全。

参考文献

Kafetzopoulou, L. E., Pullan, S. T., Lemey, P., Suchard, M. A., Ehichioya, D. U., Pahlmann, M., … Wozniak, D. M. (2019). Metagenomic sequencing at the epicenter of the Nigeria 2018 Lassa fever outbreak. Science, 363(6422), 74–77. doi:10.1126/science.aau9343 

辞旧迎新贺岁篇|走出非洲,人类基因组丢了10%?

人类基因组约3Gb,也就是30亿个碱基,如果将所有碱基排印成书,那么厚度将超过100米[1]。自从1977年测序技术诞生以来,人们已经推开了生命的奥秘之门,如今,又将好奇的目光投向人类自身。地球46亿年的漫长岁月,已有35亿年的生命史,而人的出现才几百万年。怎样定义人?人类是从地球孕育而来吗?人可以永生吗?人的终极演化形态是怎样的?也许有一天,这一切都可以从人类基因组中找到答案。

许多年以后,面对曾经科幻片[2]里才有的万能医疗舱,人类一定会想起他们第一次宣布HGP启动的1986年。

生命之书卷帙浩繁

——测序技术敲开人类基因组大门

1990年,在经过长达四年的争论和筹划后,人类基因组计划(HGP)终于获批启动,计划15年内完成绘制分析,投入资金30亿美元。

2000年,国际人类基因组测序联盟与Celera公司联合发布了基于全基因组鸟枪发测序的人类基因组草图,在2001年成果分别见刊Nature和Science杂志[3-4],发现人类基因数目仅3-3.5万个左右。值得一提的是,中国作为六个参与国家中唯一的发展中国家,测定3116Mb的序列,即完成了人类基因组的1%,精度达到了99.99%[5]。问题是,常染色质序列覆盖度只有90%,且序列之间存在近15万个空缺,导致了早期建立的很多基因模型是错误的。

2003年,中、美、英、德、日、法六国宣布比预期提前了两年完成了人类基因组序列图并于2004年发表在Nature上,进一步压缩人类编码蛋白的基因到2-2.5万个,精度达99.999%[6]。相较于2001年,常染色质的空缺只有341个,在这前后,研究者们也陆续将性染色体注解出来。

2005年,我国参与度达10%的人类单体型图谱问世[7]。

2006年,基因含量最多、解码难度最大的1号染色体登上Nature[8],标志着HGP的传奇乐章画上了休止符。

探索奥秘从未止步

——第二代测序掀起测序行业革命

2007年全球第一个白种人基因组图谱的公布标志着个体基因组时代的来临[9]。很快,深圳华大基因研究院就骄傲的宣布:第一个亚洲人基因组图谱“炎黄一号”发表于Nature[11],覆盖了36×的深度,拿到了一千一百七十七亿碱基对,比对了NCBI人类相关基因组,短reads序列达到99.97%覆盖率,而且根据参考的基因组,研究人员利用唯一的mapped reads获得了一个92%亚洲个体基因组的高质量序列集合。同时研究人员从中识别出了大约300万个SNPs,其中13.6%在dbSNP数据库中没有出现过,基因型分析证明这些SNP具有高精确性和一致性。研究人员还将这些序列与另外两个个体基因组(J. D. Watson and J. C. Venter)进行了比较,证明了第二代测序技术在个人基因组方面的应用潜力[12]。得益于第二代测序技术的高通量,整个项目不过一年时间,耗资1000万人民币!这项里程碑式的成果对中国以至整个亚洲人的治病基因、疾病预测等研究都有着非同寻常的意义。接着,第一张女性个人基因组图谱、第一张非洲人基因组图谱也相继出炉。但是第二代测序技术读长比重复元件要短,而人类基因组中已知的重复序列和片段化的重复元件占了近一半,这就导致在拼接的时候难免遗漏很多重要的信息。

 2008年炎黄一号首张中国人基因组 图谱登 《Nature》封面

技术革新再创辉煌

——第三代测序助力人类精细图谱问世

第三代测序技术以单分子测序且读长超长著称,因无需PCR,所以几无测序偏好性,由于荧光基团并不是附着于碱基而是磷酸键之上,大大降低了测序过程中的三维阻力,再加上ZMW孔锁定荧光检测区域,使读长远超二代测序。

通过PacBio SMRT,未来组助力暨南大学粤港澳中枢神经再生研究院主导的亚洲人参考基因组“华夏一号”收录于Nature Communications[10]。

“华夏一号”基因组组装策略结合了PacBio SMRT单分子实时测序技术和BioNano光学图谱分析技术,从头组装得到2.93G基因组,Contig N50为8.3Mb,Scaffold N50为22Mb得到一个中国人个体的基因组接近完成图。

图2 相较已发布人类参考基因组,“华夏一号”的Contig N50有将近10倍的提高

研究者还发现PacBio数据可以轻松跨越从5’末端到3’-Poly A tail的完整转录本,从而准确鉴定异构体,并对可变剪接、融合基因、等位基因表达等进行精确分析。在对Illumina和PacBio的测序数据的覆盖率比较后发现PacBio数据不受GC含量高低的影响,所以可以覆盖到多Illumina数据覆盖不到的区域,所以在基因组组装上优势就很明显。

图3 a.PacBio数据对GC含量异常区域覆盖更均匀,b.PacBio覆盖到Illumina覆盖不到的区域

“华夏一号”的发布填补了中国人群的疾病研究缺少精细参考基因组的不足,并将推进临床和科研大数据应用的重要基础性工作,大力推动中国的遗传疾病研究与诊断的发展。

第三代测序技术不仅有PacBio,更有人类历史上首次实现的纳米级别、也是唯一一种通过电信号的波动进行测序的技术Oxford Nanopore Technology(ONT),如2018年Nature Biotechnology上就发布了使用ONT首次高精度解析人类Y染色体着丝粒的研究[13]。

开拓视野解码升级

——非裔泛基因组补充10%人类DNA

最新的GRCh38基因组只有875个gaps[14],虽则如此,研究者们的目光多是聚焦在单个人身上,这样无疑滞碍了混血群体的研究,例如非裔人群。

近期,有约翰斯·霍普金斯大学的研究者使用全基因组鸟枪法测序法深度测序了910个非洲人种,构建了人类参考基因组中缺失但却在这910个非洲人中共有的DNA序列集,并鉴定出非洲人泛基因组在参考基因组中缺失的区域,最后发现了125715个特异contigs相当于比人类参考基因组多出超过约10%的DNA。研究者揭示出其中387个contigs来自315个特异的蛋白编码基因,余者来自基因间区域。这一研究成果发布在Nature genetics[15]上。

在这份研究中,研究人员收集了来自910个非洲人后裔群体的基因组,横跨全球20个地区包含美国、中非和加勒比等地的CAAPA(美洲非裔群体哮症协会)成员,使用Fig.1图示步骤去除了污染及冗余的contigs,最后鉴定到了GRCh38基因组中缺失的296.5Mb共125715条新序列,并且研究者将其中1548条序列(4.4Mb)锚定到了GRCh38基因组上的特异位点上,平均每个个体包含了859条插入序列,其中有一条序列同时可以在六个个体中找到。1548条序列中的302条完整定位到了基因组中的位置并解决了剩下1246条序列插入末端的断点。最长的定位到的序列为79938bp,存在于在197个样本中,而最长的未定位到的序列为152806bp,存在于11个样本。所有定位到参考基因组上的序列中的387条与已知基因相交,48个特异基因在外显子上,另外267个基因属于内含子区域;其中315个基因含插入序列,其中292个被命名(非“假设”或无意义的鉴定)。研究组装出的contigs中的31354079个碱基可以比对到GRCh38基因组上(一致性≥80%),组装成单个基因组后可以匹配上60202871个碱基(一致性≥80%)。

另外,研究还将该研究中的125715个泛基因组contigs比对到未来组参与完成的华夏一号(HX1)基因组和韩国人基因组(KOREF1.0)上,发现有42207个contigs共120.7Mb可以比对到韩国人或华夏一号基因组上且一致性≥90%,覆盖度≥80%,优于对GRCh38的匹配度,其中一个区段的示例见表1和图4,这个发现表明生成GRCh38基因组的个体缺失了一部分序列。

表1 非裔泛基因组contigs和华夏一号及韩国人基因组的比较

图4 将非裔泛基因组和华夏一号、GRCh38基因比对

Shi et al等人于2016年组装的华夏一号基因组报道了12.8Mb的新DNA,研究者则发现华夏一号和该研究中生成的特异序列共享68.1Mb的DNA。总之,该研究发现,亚洲人泛基因组产生的序列中有296.5Mb相当于10%的基因组大小在标准人类参考基因组中是缺失的,这其中有120.7Mb可以在韩国人或者华夏一号基因组中找得到,间接表明这些DNA代表的基因区域在GRCh38基因组所代表的群体中于更近的时期丢失了或者十分罕见,也可以说明单个参考基因组不适宜基于群体的人类遗传学研究,将来或许会有更好的方法获取综合性的人类泛基因组,捕获所有人类中的DNA。

随着“中国十万人基因组计划”、“地球生物基因组计划”的相继问世,人类仍在组学研究的大潮中寻找属于自己的那朵浪花,现有的成果距离揭示人类的奥秘还有很长的路要走。武汉未来组有幸也成为了一名弄潮儿,由未来组发起的“个人参考基因组服务计划”、“华夏万人结构变异计划”正在如火如荼的进行,并且进展顺利。人的智慧无穷尽,探索的脚步永不停,总有一天,人们只需要带着自己的基因图谱去看医生,扫描一下数据就可以直达病灶,然后躺进万能的医疗舱,出来的时候百病全消……

辞旧岁,迎新年,在这里,组学君默默祝祷,期望每一个善良的人都平安喜乐,百病不生。在新的一年,未来组将创造更多的成绩回馈社会,回馈组学领域,也祝愿每一位科研工作者硕果累累,万事如意!

已发表精细人类基因组图谱[10,16]

2018年NCBI上收录的人类基因组组装版本[17]

参考文献

[1]http://jiyongqing.blogchina.com/2427017.html

[2]《第五元素》、《极乐空间》、《普罗米修斯》等好莱坞科幻影片中均出现过能够复原生命或者治疗人类疾病的医疗舱。

[3]International Human Genome Sequencing Consortium. Initial sequencing and analysis of the human genome. Nature 409, 860–921 (2001).

[4]Venter, J. C. et al. Te sequence of the human genome. Science 291,1304–1351 (2001).

[5]骆建新, 郑崛村, 马用信, et al. 人类基因组计划与后基因组时代[J]. 中国生物工程杂志, 2003, 23(11):87-94.

[6]Finishing the euchromatic sequence of the human genome[J]. Nature:931-945.

[7]A haplotype map of the human genome : Article : Nature[J]. Nature, 2005.

[8]Gregory S G , Barlow K F , Mclay K E , et al. Corrigendum: The DNA sequence and biological annotation of human chromosome1[J]. Nature, 2006, 441(7091):315-321.

[9]高媛. 后基因组时代的生物信息学发展[J]. 中国科技信息, 2009(10):225-226.

[10]L. Shi, et al., Long-read sequencing and de novo assembly of a Chinese genome. Nature Communications (2016)

[11]The diploid genome sequence of an Asian individual[J]. Nature.

[12]https://m.antpedia.com/news/49844.html

[13]Jain M, Olsen H E, Turner D J, et al. Linear assembly of a human centromere on the Y chromosome[J]. Nature biotechnology, 2018.

[14]Schneider, V. A. et al. Evaluation of GRCh38 and de novo haploid genome assemblies demonstrates the enduring quality of the reference assembly. Genome Res. 27, 849–864 (2017)

[15]Sherman, R. M. at al. Assembly of a pan-genome from deep sequencing of 910 humans of African descent. Nature Genetics.51, pages30–35 (2019)

[16]De novo assembly and phasing of a Korean human genome. Nature 538,243–247 (13 October 2016) doi:10.1038/nature20098

[17]https://www.ncbi.nlm.nih.gov/assembly/organism/9606/latest/