A. 三代测序入门
移步github
共有的特点:
10X Genomics,是常规Illumina二代测序的升级版,由于开发出了一套巧妙的Barcoding建库方案,使得Illumina这种短读长二代测序能够得到跨度在30-100Kb的linked reads信息,与二代测序数据相结合,在Scaffold的组装上能够得到媲美三代测序的组装结果
其GC偏好性如何?
10X Genomics技术相对于Illumina来说,有改进,但依旧是个拱形,而PacBio则是无偏倚的均一分布。10X的技术,其Coverage一样是受GC含量影响较大的,那么如果真要应用10X技术,那么必须注意目标DNA的GC含量分布最好能控制在30~70%。
真正的单分子测序(Helicos True Single Molecule Sequencing)
待测DNA 被随机打断成小片段,在每个小片段( 200bp)的末端加上poly-dA,并于玻璃芯片上随机固定多个 poly-dT 引物,其末端皆带有荧光标记,以利于精确定位。
首先,将小片段 DNA 模板与检测芯片上的poly-dT 引物进行杂交并精确定位,然后逐一加入荧光标记的末端终止子。这个终止子与 Illumina 的终止子可不一样,不是四色的,是单色的,也就是说所有终止子都标有同一种染料。
在掺入了单个荧光标记的核苷酸后,洗涤,单色成像,之后切开荧光染料和抑制基团,洗涤,加帽,允许下一个核苷酸的掺入。通过掺入、检测和切除的反复循环,即可实时读取大量序列。最后以软件系统辅助,可分析出完整的核酸序列。
缺点 :Heliscope 在面对同聚物时也会遇到一些困难,但可以通过二次测序提高准确度;由于在合成中可能掺有未标记的碱基,因此其最主要的错误来源是缺失。
PacBio SMRT(single molecule real time sequencing)技术也应用了边合成边测序的思想,并以SMRT 芯片为测序载体。
基本原理是:DNA 聚合酶和模板结合,4 色荧光标记4 种碱基(即是dNTP),在碱基配对阶段,不同碱基的加入,会发出不同光,根据光的波长与峰值可判断进入的碱基类型。
DNA 聚合酶是实现超长读长的关键之一,读长主要跟酶的活性保持有关,它主要受激光对其造成的损伤所影响。
PacBio SMRT 技术的一个关键是怎样 将反应信号与周围游离碱基的强大荧光背景区别出来 :
优缺点:
该技术的关键之一是,它们设计了一种特殊的纳米孔,孔内共价结合有分子接头。当DNA 碱基通过纳米孔时,它们使电荷发生变化,从而短暂地影响流过纳米孔的电流强度(每种碱基所影响的电流变化幅度是不同的),灵敏的电子设备检测到这些变化从而鉴定所通过的碱基。
测序原理:
特点:
Nanopore 测序仪 MinION 的一些特征:
ONT公司目前推出的几款测序仪:
在analysis文件夹中,下机的数据被分割为三个文件进行存储
数据的命名:
Pacbio 数据的文库模型是两端加接头的哑铃型结构,测序时会环绕着文库进行持续的进行,由此得到的测序片段称为 polymerase reads ,即一条含接头的测序序列,其长度由反应酶的活性和上机时间决定。目前,采用最新的 P6-C4 酶,最长的读长可达到 60kb 以上。
polymerase reads 是需要进行一定的处理才能获得用于后续分析的。这个过程首先是去除低质量序列和接头序列:
处理后得到的序列称为 subreads ,根据不同文库的插入片段长度,subreads 的类型也有所不同。
对长插入片段文库的测序基本是少于2 passes的(pass即环绕测序的次数),得到的reads也称为 Continuous Long Reads (CLR) ,这样的reads测序错误率等同于原始的测序错误率。
而对于全长转录组或全长16s测序,构建的文库插入片段较短,测序会产生多个passes,这时会对多个reads进行一致性校正,得到一个唯一的read,也称为 Circular Consensus Sequencing(CCS)Reads ,这样的reads测序准确率会有显着的提升。
不同于二代测序的碱基质量标准Q20/Q30,三代测序由于其随机分布的碱基错误率,其单碱基的准确性不能直接用于衡量数据质量。那么,怎么判断三代测序的数据好不好呢?
需要关注的是两个比例:
目前采用的组装策略:
这四种组装策略并不是完全孤立的,在一个组装任务的不同阶段会用到不同的方法
不同的组装策略可以选用的工具:
基因组的组装问题,实际上就是从序列得到的图中搜寻遍历路径的问题,有两种构建图的方法:
可以看到,随着reads长度的增加,基于OLC算法的组装工具组装出的contigs的长度几乎在线性增长,而基于de Bruijn图算法的组装效果并没有随着reads长度的增加而提高
三代单分子测序会产生较高的随机错误,平均正确率在82.1%-84.6%。这么高的错误率显然不能直接用于后续的分析,需要进行错误校正:
校正过程中会将short reads未覆盖到的Gap进行裁剪,short reads在PacBio long reads上的覆盖情况:
这样做的其中一个考虑是去除adapter
那么是什么原因导致了低覆盖度区域的产生的呢?
Base-calling做的就是从测序仪输出的电流信号波形图中将碱基解码 (decoding) 出来
第一步就是就是对波形图进行分段 (segmentation),即检测每个current shift的边界,这一步由ONT公司提供的 MinKNOW 完成,但是分段基于的假设是ssDNA分子匀速穿过nanopores,但是由于ssDNA穿过nanopore的速度很快,很容易产生一两个碱基的速度差异,这样就容易在decoding时造成insert和delete
接着就基于current shift进行base calling,ONT公司提供的base caller为Metrichor,其底层算法基于HMM,将可能的k-tuple(由k个碱基组成的序列)作为隐藏状态,将current signals作为观测状态。ONT公司最新开发出的Metrichor用RNN取代了HMM,并将其整合到其开发出的新的生物信息数据分析平台EPI2ME中
随后,科研圈又开发出了开源的base calling工具,Nanocall 和 DeepNano。
ONT后来又在github上开源了一个RNN base-caller —— Nanonet
测序时,测序仪 MinION 连接上主机,安装在主机上的软件 MinKNOW 控制测序仪,对于每条reads,其 signal segmentation 结果(包括segment mean, variance and ration)以及测序过程中的 metadata 会被保存成FAST5格式的二进制文件(基于 HDF5标准 的变种)。
保存在FAST5文件中的原始数据会经过云端的Metrichor的处理,产生的解码的序列会被保存在另外的以 .FAST5 为后缀的HDF5文件中,包含一条template read和一条complement read或只有一条 2D read 。
MAP (MinION Access Programme) community 开发出的用于处理FAST5文件的工具,它们均能从FAST5文件中解析出FASTA/FASTQ文件,除此之外还有各自特色的质量统计功能:
参考资料:
(1) 生物技能树论坛:PacBio sequence error correction amd assemble via pacBioToCA
(2) 天津医科大学,伊现富《系统生物学-chapter2》
(3) Nanopore 第四代测序技术简介
(4) Magi A, Semeraro R, Mingrino A, et al. Nanopore sequencing data analysis: state of the art, applications and challenges.[J]. Briefings in Bioinformatics, 2017.
(5) 细节曝光!Oxford Nanopore真机还原,听听圈内人怎么说
(6) 三代测序--QC篇
(7) PacBio Training: Large Genome Assembly with PacBio Long Reads
(8) Koren S, Schatz M C, Walenz B P, et al. Hybrid error correction and de novo assembly of single-molecule sequencing reads[J]. Nature Biotechnology, 2012, 30(7):693-700.
(9) 冷泉港ppt:Hybrid De Novo Assembly of Eukaryo6c Genomes
(10) Leggett R M, Darren H, Mario C, et al. NanoOK: multi-reference alignment analysis of nanopore sequencing data, quality and error profiles[J]. Bioinformatics, 2016, 32(1):142-144.
B. 什么是人类基因组计划科学家们使用怎么样的技术策略来完成它
人类基因组计划(human genome project, HGP)是由美国科学家于1985年率先提出,于1990年正式启动的。美国、英国、法兰西共和国、德意志联邦共和国、日本和我国科学家共同参与了这一价值达30亿美元的人类基因组计划。这一计划旨在为30多亿个碱基对构成的人类基因组精确测序,发现所有人类基因并搞清其在染色体上的位置,破译人类全部遗传信息。
HGP的主要任务是人类的DNA测序,包括的四张谱图[遗传图谱(genetic map),物理图谱(physical map),序列图谱 ,基因图谱],此外还有测序技术、人类基因组序列变异、功能基因组技术、比较基因组学、社会、法律、伦理研究、生物信息学和计算生物学、教育培训等目的。
大规模测序基本策略
逐个克隆法:对连续克隆系中排定的BAC克隆逐个进行亚克隆测序并进行组装(公共领域测序计划)。
全基因组鸟枪法:在一定作图信息基础上,绕过大片段连续克隆系的构建而直接将基因组分解成小片段随机测序,利用超级计算机进行组装(美国Celera公司)。
C. 睡莲基因组与开花植物的早期进化
睡莲属于被子植物睡莲目(Nymphaeales)。无油樟目(Amborellales)、睡莲目(Nymphaeales)和木兰藤目(Austroleyales)共同组成了早期被子植物类群(ANA被子植物类群),它们是现存被子植物的代表,从谱系分化的最早期发展到现存的被子植物。在此,我们公布了蓝星睡莲( Nymphaea colorata )的基因组(409Mb)。系统发育组显示,睡莲和无油樟属于早期被子植物类群。通过蓝星睡莲基因组及其他19个睡莲的转录组分析,显示睡莲科祖先发生了一次全基因组复制事件,这次复制事件可能由睡莲科(Nymphaeaceae)和莼菜科(Cabombaceae)所共享。在全基因组复制事件保留的基因中,有调节花期转变和花期发育的同源基因。蓝星睡莲中花ABCE同源基因的广泛表达,可能揭示了在早期被子植物花器官中可能具有类似广泛活跃的ABCE祖先模型。睡莲进化出了迷人的花香和颜色,这是被子植物共有的特征,在蓝星睡莲中我们推测出了它们的生物合成基因。花香味背后的化合物和生物合成基因表明,它们的进化与被子植物是平行的。由于其独特的系统发育位置,蓝星睡莲基因组揭示了被子植物的早期进化。
许多睡莲属( Nymphaea )物种,特别是睡莲科(Nymphaeaceae)的睡莲,花朵大而艳丽,属于被子植物(也称为开花植物)。它们的美深深吸引着许多着名的艺术家,例如法国印象派画家莫奈(Claude Monet)。睡莲的花被(外部花器官)分化有限,但它们同时具有雄性和雌性器官,并且具有多种不同的气味和颜色,与许多被子植物(核心被子植物,包括双子叶植物,单子叶植物和木兰科植物)相似(Supplementary Note 1)。此外,一些睡莲的生命周期短,种子数量多,这增加了它们作为早期被子植物类群(ANA被子植物类群)模式植物代表并研究被子植物内部早期进化事件的潜力。特别是,蓝星睡莲( N. colorata )的基因组相对较小(2n = 28;约400 Mb),且蓝色的花瓣使它在育种中很受欢迎(Supplementary Note 1)。
在此,我们利用PacBio RSII单分子实时(SMRT)测序技术获得的蓝星睡莲( N. colorata )基因组序列。基因组组装成1429条contigs (contig N50为2.1Mb),总长度409 Mb, 804个scaffolds,其中770个scaffolds锚定在14条染色体上(Extended Data Fig. 1 and Extended Data Table 1)。基因组完整性评估为94.4%(Supplementary Note 2)。我们对31580个编码蛋白基因进行了注释,并预测了总长度为160.4 Mb的重复元件,占总基因组的39.2%(Supplementary Note 3)。
蓝星睡莲( N. colorata )基因组为解决无油樟目(Amborellales)、睡莲目(Nymphaeales)和所有现存被子植物之间的关系提供了一个机会( Fig. 1a )。使用六个真双子叶植物,六个单子叶植物,蓝星睡莲及无油樟属,每三个裸子植物(银杏( Ginkgo biloba )、云杉( Picea abies )和火炬松( Pinus taeda ))作为又一个类群,我们分别鉴定了2169、1535和1515个 直系同源低拷贝核基因(****LCN****) ( Fig. 1b )。当使用银杏( Ginkgo biloba )作为外群时,从核苷酸序列推断出的LCN基因树中,62%(475中的294)将无油樟(Amborella)作为所有现存被子植物的姐妹系,且自展支持度(bootstrap support)大于80%(type II, Fig. 1c )。而当使用云杉( Picea abies )和火炬松( Pinus taeda )作为外群时,在LCN基因树中,分别有57%和54%将无油樟(Amborella)作为所有现存被子植物的姐妹系,即支持无油樟(Amborella)是最早的被子植物类群。且利用氨基酸序列推断出的LCN基因树具有相似的系统发育模式(Supplementary Note 4.1)。
为了使稀疏分类单元采样的潜在缺陷最小化,我们还使用来自44个基因组和71个转录组的序列,包括ANA被子植物类群,双子叶植物,木兰类植物,单子叶植物和裸子植物外群(买麻藤( Gnetum montanum )、银杏( Ginkgo biloba )、云杉( Picea abies )和火炬松( Pinus taeda ))的代表来推论被子植物的物种进化树。为了对这115个物种进行进一步的系统发育推断,我们根据不同的标准选择了5种不同的LCN基因集,包括1167、834、683、602和445个基因。对这五个数据集的分析均得出与无油樟(Amborella)相似的树形拓扑,睡莲目(Nymphaeales)作为所有其他现存被子植物的连续姐妹系。
使用101个严格的LCN基因以及基于21个化石的年龄对被子植物谱系的分子年代测定进行校准。推断出被子植物的冠龄为2.34-2.63亿年前(Ma)( Fig. 1d )。单子叶植物和双子叶植物之间的分界估计在1.71-2.02亿年之间,而睡莲科(Nymphaeaceae)和莼菜科(Cabombaceae)之间的分化在1.47-1.85亿年之间。
基因组共线性揭示了蓝星睡莲( N. colorata )发生全基因组复制(WGD)事件的证据(Extended Data Figs. 1f, 2a and Supplementary Note 5.1)。蓝星睡莲(N. colorata)旁系同源基因的每个同义位点上的同义替换( Ks )分布的数量进一步表明,有一个 Ks 约为0.9的显着峰值( Fig. 2a ),而在其他睡莲科(Nymphaeaceae)物种中也鉴定到了类似的 Ks 峰值(Supplementary Note 5.2)。这表明,一个古老的单一的全基因组复制事件(WGD)可能是睡莲科成员所共有的。通过比较蓝星睡莲( N. colorata )旁系同源与蓝星睡莲( N. colorata )和其他睡莲目世系(Nymphaeales lineages)、红茴香( Illicium henryi )、无油樟( Amborella )之间的直系同源(代表物种形成事件) Ks 分布,发现全基因组复制事件(WGD)发生在睡莲科(Nymphaeaceae)与莼菜科(Cabombaceae)分化之后( Fig. 2a )。相比之下,对至少包含一个来自蓝星睡莲( N. colorata )共线区域的旁系同源基因家族的系统基因组学分析表明,全基因组复制事件(WGD)在睡莲科(Nymphaeaceae)和莼菜科(Cabombaceae)之间共享( Fig. 2b , Supplementary Note 5.4)。如果属实,那么莼菜科水盾草( Cabomba caroliniana )似乎保留了很少的重复( Fig. 2b, c ),这也可以解释水盾草( Cabomba caroliniana )旁系同源Ks分布中没有明显的峰(Supplementary Note 5.2)。考虑到Nymphaealean谱系中可变替换率( Fig. 2a****, b , Extended Data Fig. 2c),对蓝星睡莲( N. colorata )的绝对年代测定确实表明,全基因组复制事件(WGD)可能发生在睡莲科(Nymphaeaceae)与莼菜科(Cabombaceae)分化之前或接近于它们的分化(Extended Data Fig. 2d, Supplementary Note 5.3)。对上述结果的另一种解释可能是,全基因组复制事件来自于发生在睡莲科祖先和莼菜科系谱之间的异源多倍事件,在它们分化后不久,睡莲科(但不是莼菜科)的主干分支得以兴起( Fig. 2d , Supplementary Note 5.4)。
睡莲起源于被子植物早期分化的一个分支,早于被子植物大范围的辐射扩张。因此,睡莲家族为了解被子植物,特别是开花植物的早期进化,提供了一个独特的窗口。我们鉴定了70个MADS-box基因,包括参与花器官发育ABCE模型的同源基因: AP1 (还有 FUL) 及 AGL6 (A参与萼片和花瓣发育), AP3 和 PI (B参与花瓣和雄蕊发育), AG (C参与雄蕊和心皮发育), 以及 SEP1 (E与ABC功能蛋白相互作用)。对MADS-box基因及其基因组邻域的系统发育和共线性分析表明,在种子植物分化之前就存在古老的串联重复,产生了A功能基因( FUL )和E功能基因( SEP )的祖先(Extended Data Fig. 3, Supplementary Note 6.1)。此外,由于睡莲(Nymphaealean)全基因组复制事件(WGD),蓝星睡莲( N. colorata )具有两个旁系同源基因,即C功能基因AG的AGa和AGb(Extended Data Fig. 4)。类似地,由睡莲(Nymphaealean)WGD衍生的重复序列同与心皮和雄蕊发育相关的其他基因、以及调控开花时间及生长素调控花的昼夜开合的基因是同源的(Extended Data Figs. 4–6, Supplementary Note 6.2–6.4)。
蓝星睡莲( N. colorata )ABCE同源基因的表达谱与它们在花器官中推测的作用基本一致( Fig. 3a )。值得注意的是,蓝星睡莲( N. colorata ) AGL6 同源基因主要在萼片和花瓣中表达,而 FUL 同源基因主要在心皮中表达,说明 AGL6 在蓝星睡莲( N. colorata )中起A功能基因的作用。两种C功能同源基因 AGa 和 AGb 分别在雄蕊和心皮中高表达,而 AGb 也在萼片和花瓣中表达,表明它们可能在睡莲(Nymphaealean)WGD后经历了花发育的亚功能化和可能的新功能化。此外,与双子叶模型系统相比,蓝星睡莲( N. colorata )的ABCE同源基因在花器官中的表达范围更广( Fig. 3b )。这种更广泛的表达模式,与至少一些ABCE基因在一些双子叶植物中更广泛的表达相结合,代表了一个早期分化谱系,一些单子叶植物和木兰类植物,提出了一种古老的ABCE花发育模型,在被子植物,特别是核心双子叶植物的进化过程中,随后渠限化基因的表达和功能受到更特异的ABCE基因的调控。这也可以解释为什么在睡莲属植物中萼片和花瓣的分化是有限的,这与被子植物祖先花中花被器官的单一类型是一致的。
花香为昆虫传粉者提供嗅觉线索。然而无油樟属的花是无香味的,蓝星睡莲的花释放11种不同的挥发性化合物,包括萜类化合物(倍半萜烯)、脂肪酸衍生物(甲基癸酸酯)及苯环型化合物( Fig. 4a )。蓝星睡莲基因组包含92个假定的萜烯合酶( TPS )基因,这些基因归属于被子植物中4个已知的TPS亚家族:TPS-b, TPS-c, TPS-e/f 及TPS-g( Fig. 4b ),但是在被子植物中没有发现负责倍半萜生物合成的TPS-a。值得注意的是,在蓝星睡莲中,TPS-b亚家族含有80多个基因;其中NC11G0123420在花中高表达(Extended Data Fig. 7);这一结果表明,该基因可能是蓝星睡莲倍半萜烯生物合成酶的候选基因。此外,并未在单子叶和双子叶挥发性化合物中检测到癸酸甲酯,其被认为是由蓝星睡莲( N. colorata )SABATH甲基转移酶家族合成的。蓝星睡莲( N. colorata )基因组包含13个 SABATH 同源基因,其中12个形成睡莲目特异性家族(Supplementary Fig. 41)。在这12个成员中,NC11G0120830在花瓣中表达最高( Fig. 4c ),并且其相应的重组蛋白被证明是脂肪酸甲基转移酶,其以癸酸为底物具有最高的活性( Fig. 4d , Supplementary Note 7.1)。这些结果表明,蓝星睡莲( N. colorata )的花香生物合成是通过酶的功能完成的,而酶的功能是独立于被子植物的功能而进化的( Fig. 4e )。
睡莲( Nymphaea colorata )美丽迷人的蓝色花瓣被认为是很有价值的,这在观赏植物中是较为罕见的特征。为了理解蓝色的分子基础,我们鉴定到翠雀素(3′- O -(2″- O -galloyl-6″- O -acetyl-β-galactopyranoside))为主要蓝色花青素色素(Extended Data Fig. 8a–c)。通过比较两个蓝星睡莲品种中白色和蓝色花瓣中花青素生物合成途径中基因的表达谱,我们发现花青素合酶和翠雀素修饰酶基因的表达在蓝色花瓣中明显高于白色花瓣(Extended Data Fig. 8d, e)。这两种酶催化花青素生物合成的最后两个步骤,因此是蓝色素生物合成的关键酶。
睡莲在全球范围均有分布,包括寒冷地区(中国北部及加拿大北部),这与其他ANA被子植物类群不同,无油樟属仅在太平洋岛屿有分布,而八角茴香目仅在温带和热带地区有分布。与无油樟属及一些被子植物相比,我们发现蓝星睡莲中与免疫和应激反应相关的基因明显有扩张,包括编码核苷酸结合富亮氨酸重复(NLR)蛋白、蛋白激酶和WRKY转录因子基因(Extended Data Fig. 9, Supplementary Note 8)。这些基因数量的增加可能使睡莲适应了全球各种生态栖息地。
综上所述,蓝星睡莲(N. colorata)基因组为比较基因组学和解决被子植物间的系统发育关系提供了参考。它还揭示了睡莲科祖先发生的一次全基因组复制事件,并提供了关于被子植物早期发育及进化的重要见解,涉及诸如花的发育、花的气味和颜色等。
Zhang, L., Chen, F., Zhang, X. et al. The water lily genome and the early evolution of flowering plants. Nature 577, 79–84 (2020). https://doi.org/10.1038/s41586-019-1852-5
大多数的种系发生重建方法会产生无根树,但是观察树的拓扑结构无法识别树根应在哪一分支上。实际中,对于要证实哪一个分类单元的分支先于其他的分类单元,树根必须确定。
在无根树中设定一个根,最简单的方法是在数据集中增加一个外群(outgroup)。 外群是一种分类操作单元,且有外部信息表明外群在所有分类分类群之前就已分化。合适的外群与待分析的分类群关系不能相距太远,因为在比较关系较远的物种时,系统发生的信号会降低,这是核苷酸替换饱和的结果。使用一个以上的外群通常可以进一步改善推导的树状拓扑的准确度。
所谓的外类群就是与你研究的序列关系极为密切的序列,且外类群能很好的聚为一支(若外类群不止一条序列),若研究的是演化历史,一般应选择比目标序列具有较早进化历史的序列作为外类群。
另一种可选的引入外群的方法是,使用两套相同的、同时存在于待分析的所有分类操作单元中的并系同源基因。在这种方法中,第一个并系同源基因群中的基因可以成为第二个并系同源基因群中基因的外群。这种确定的系统已用于确定tree of life的第一层分支,树根可以置于通向生命树中细菌、古细菌以及真核细胞中任一分枝上。当使用单一外群时,根可以置于通向外群的分支上。另外,若使用多个外群,根必须置于连接外群和内群的分支上。
如果是鉴定物种,最好选一个外群。在缺少一个合适的外群时,根大约可以置于两个分类操作单元间最长支的中点上。这种确定根的方法叫做中点定根(midpoint rooting),当在树中所有分支的进化速度大致相同而且实际的外群与其它分类群间的支的长度不太短时,这种方法相当准确,但是中点生根这种方法慎用,它有一个假设前提:假设两个最不同的谱系以相同的速率进化。显然,这个假设现实中很可能不成立。
在进化过程中,新基因通常来自事先存在的基因,新基因的功能从先前基因的功能进化而来。新基因的原材料来自基因组区域的重复,这种重复可包括一个或多个基因。作为物种形成的伴随事件而被重复,并继续保持相同功能的基因,称为直系同源基因(orthologous gene)。新的基因功能可由在单个物种的基因组中发生的重复引起的。在一个基因组内部的重复导致旁系同源基因(paralogous gene)。
Orthology VS Paralogy
Relation of sequences
Orthologs: similar sequences that have arisen e to a speciation event.
Functionality Retained.
Orthologs: members of a gene (protein) family in various organisms.
Paralogs: Similar sequences that have arisen e to a gene plication event.
Paralogs are not necessarily to have the same or similar functions. Probably become pseudogenes.
Paralogs: members of a gene (protein) family within a species.
Xenologs: Similar sequences that have arisen out of horizontal transfer events.
Examples: Transformation; Conjugation; Transction; Transgene
所谓Bootstraping法 就是从整个序列的碱基(氨基酸)中任意选取一半,剩下的一半序列随机补齐组成一个新的序列。这样,一个序列就可以变成了许多序列,一个多序列组也就可以变 成许多个多序列组。根据某种算法(最大简约性法、最大可能性法、除权配对法或邻位相连法)每个多序列组都可以生成一个进化树。将生成的许多进化树进行比 较,按照多数规则(majority-rule)我们就会得到一个最“逼真”的进化树。
Jackknife则是另外一种随机选取序列的方法。它与Bootstrap法的区别是不将剩下的一半序列补齐,只生成一个缩短了一半的新序列。
在分支分类学中具有一个不为其他分类单元所共有的祖先的两个分类单元称为姐妹群。姐妹群是由一个祖种通过分裂产生的一对分支,是建立系统发育系统的基本结构,根据近裔共性加以识别。
D. “鸟枪法”
鸟枪法(Shotgun method):使用基因组中的随机产生的片段作为模板进行克隆的方法。
使用限制性内切酶将带有目的基因的DNA链切成若干小段
再使用DNA连接酶将其整合到载体的基因中,并使其表达
如果在某个细胞中得到了目的产物,就说明整合到该细胞中的DNA片断就是所需要的DNA片断
E. 基因组测序的具体过程
讲起来很复杂哈,简单来说,就是边合成边测序,把所要测序的序列通过各种手段来建库,然后加上5‘和3’都加上接头,放在一个小玻璃片一样的芯片上,上面一般是8个lane,把样品加入到lane上之后和上面长好的接头序列相结合(之前的建库时加的接头与之反向互补),然后再经过PCR扩增的过程形成一个DNA簇,然后再进行新一条链的合成,在此合成过程中边合成边测序,没加上去一个碱基,机器就读一次,根据其发到荧光来判断加的是那种碱基,这样就能得到所要的序列信息了。
测序会有误差的,就好比PCR扩增循环多了也会有误差一样,假如说在第一百个碱基合成时其误差是99.9%,那下一个的误差就是99.9%*99.9%.....依次以指数形式递增,到一定程度后测得的序列就不可靠了,所以会有一个极限,现在solexa的极限貌似比之前的150bp要长了。
当然以上都是指的Solexa的平台,454是焦磷酸测序原理,就又是一种方式了,这种方法得到的可靠序列长度之前说是500bp以上,现在估计都能测通1000bp了。
F. 矫正基因结构注释 - 做有良心的基因家族分析
半个月前,我推了一个《任何人都能掌握-基因家族分析》的腾讯课程(原本事实是开给课题组)。在热身课程(完全免费)中,将我个人对基因家族分析的认知和其意义均做了说明,感兴趣的可见 https://ke.qq.com/course/338062?tuin=72ed3eb
其中涉及到一点,即是,基因家族分析中一个常常被忽略甚至忽视的,对科研可能有所贡献的步骤(可能很多培训公司并不会涉及),那就是 矫正基因结构注释 。
基因组,尤其是植物基因组,从测序,到组装,到注释,每一个都不简单。甚至存在一种说法, 一篇基因组文章,一个组装算法 。而事实上,注释也是类似的。即使是拟南芥或者水稻这两个模式生物,都不能保证所有基因的结构都被注释出来。更何况刚发表的基因组?
在基因家族分析讲演中,有这么一个图
在半个月前的讲演中,其实也已经讲过了,用在线网页工具softberry就可以了,以上图的 Aco005453.1为例,基于motif pattern和domain info,可以明显地看出来,其缺少的是5端,GRAS结构域也被截断了( 注意,这个在几乎所有物种的基因组结构注释文件都会出现,因为软件永远不可能保证绝对的准确,至少目前这个事情上是的 ,而菠萝基因组已经做得很优秀了。)
1.首先,获得这个基因的位置信息
使用TBtools的gff3 gene info工具
恩。。。我又试了几个基因,基本无解。或者是正好基因与上下游有overlap,或者是基因结构预测后并没有太大的改善。不过似乎也OK,毕竟我们大概可以知道其中两个基因并不真实。
如果要证明他们是否真实, 可能还是需要RNAseq数据辅助,或者事实上,需要race实验 。得到序列之后,再按照上述操作,用TBtools重构gff3即可。
嗯,意料之外,情理之中。
据我个人了解,菠萝基因组应也是经过了基因结构注释的人工矫正。所以从某个角度来说,大部分基因结构应是正确。可能单纯从文本预测上,确实无法改善结构注释信息。而只能发现一些确实有问题的序列。对于这些序列,或者是事实并不存在的假基因,或者是需要race实验做进一步获取。
无论哪一种途径,得到序列之后可以直接用TBtools重构gff3信息,并修改原始gff3文件。做进一步分析。
G. 宏基因组shotgun入门笔记
目录
根据分析对象和实验目的,宏基因组的研究基本上可以分为
1. Pre-processing
2. Sequence analysis
包括两种分析策略: read-based (mapping) 和 assembly-based
简单来说,assembly-based approach 受到覆盖度的制约,因为组装时低覆盖度的区域是不会进行组装的,而是被丢弃,这样低丰度的细菌的信息就被丢弃了,反映在reads利用率上,就是往往reads利用率极低,往往低于50%
而 read-based (mapping) approach 则受到reference databases的制约,因为细菌的遗传多样性很高,即便是同一个菌种,它的不同菌株,其基因组的组成也是有相对比较大的差异的,那么在mapping的时候就会出现mapping不上的问题,使得mapping效率不够高;而且只能分析reference databases中有的物种,对于reference databases未收录的新物种,是无法进行分析的。
不过可用的微生物参考基因组正在迅速地增加,包括那些原先难以培养的细菌由于培养方法的改进,使得对其进行测序成为可能,再加上单细胞测序的途径和 metagenomic assembly的途径得到的基因组序列。现在一些类型的环境样品(如人肠道)的参考基因组的多样性已经可以满足 assembly-free taxonomic profiling 的要求。
随着测序成本的下降和测序深度的增加,其分析难度将会越来越大,制约效应也将会越来越明显
预计的单位测序成本将会以指数关系下降,但其中计算成本下降的幅度会远慢于测序成
在数据存储和数据处理的层面上,rDNA和扩增序列的分析难度较小,基本可以在个人电脑或者小型服务器上完成,但宏基因组全测序的分析却主要受限于计算技术的发展
即使在同一个环境中获取的不同样本,其微生物组成也会存在比较大的差异,这使得在样本集之间,寻找具有统计学显着性和生物学意义的差异变得很困难。因此如何做到,在即使其影响因素的作用程度很小的情况下,也能有效地检测出差异就显得十分重要。
一种策略是,构造 pilot data,即将不同浓度的绝对定量 control (spike-in) 加入到样本中,来评估实验与分析方法的稳健性(robust);
另一种策略:two-tiered approach,即挑取少部分样本,既做 16s rDNA 测序,又做 shotgun metagenomics 测序,对比这两个层次的结果来评估实验结果的稳健性。
两种研究策略:
由于在研究宏基因组过程中,比如研究人类的微生物群,影响其微生物群的因素众多,包括宿主基因型,年龄,饮食习惯等等,当进行两个环境微生物群横向比较时,很难做到控制变量,使得在进行比较分析时混入了许多干扰因素;此时如果进行单一环境微生物群多时间点采样的纵向比较,就可以从很大程度上消除这种影响。
1. 样本量与测序深度
当实验目的是检出显着性差异时,样本量与测序深度的选择取决于(1)不同样本间微生物组组成的一致性,(2)样本固有的微生物多样性,(3)影响因素的效应量(effect size)
建议:参考前人在类似环境中的研究。若没有可参照的类似研究,选择marker gene做预实验
2. Confounding variables and control groups
在进行宏基因组研究时,往往很难找到与目标样本集对应的没有其他干扰因素的对照组
建议:目前最佳的解决策略是,尽可能地搜集各个样本群体的元数据 (metadata),然后在随后的比较分析中将它们考虑进去。比如临床样本,包括性别、年龄、是否使用抗生素/药物、取样位置、饮食习惯等等。比如环境样本,包括地理位置、季节、pH、温度等等。
元数据的搜集可以参照MIMARKS (Minimum information about a marker gene sequence) 和 MIxS (minimum information about any (x) gene sequence) 标准
3. Sample collection/preservation
样本的处理和保存过程的差异会带来系统偏差,比如when samples are provided from a number of locations by different research groups,或者在纵向研究中,不同取样时间点的样本的保存时间长短不一。有时这些处理步骤的效应量可能比你感兴趣的生物学变量还大。
建议:尽可能按照相同的标准来进行取样和保存
4. Biomass/Contamination
当前采用的基于测序的方法具有很高的灵敏度 (highly sensitive),即使非常微量的DNA也能被检测出来。而实验室中使用到的常规仪器和试剂并不是无菌的,这样就很可能在实验操作过程中,人为地引入污染。由于检测方法的高灵敏度,当原样本的微生物量很少时,污染带来的信号很可能会盖过真实的信号。
建议:在上机测序前,做好微生物量的定量 (qPCR)。当样品中的微生物数量少于10 5 数量级时,其极有可能会受到背景污染的干扰。此时,可以参照以下的方法进行细胞/DNA的富集:
可以增设负对照实验 (Negative control),对其进行与实际样本相同的操作,使用相同的试剂,以此来找出污染的细菌类型,这样就可以在后续的生物信息学分析过程中将其过滤掉。
5. 选择合适的DNA提取方法
DNA提取的效果会直接对后续的实验和分析产生巨大的影响。DNA提取方法的选择依赖于样品中细胞类型的组成,然而即使是相同类型的样品其微生物组成也具有较大的差异(当人粪便中革兰氏阴性菌主导时,细胞很容易裂解,而当由相对顽强的革兰氏阳性菌主导时,则相反)。
因此不存在适用于所有样品的最佳的DNA提取方案。
若方案选择不当,则获得的DNA主要来自于那些易裂解的细菌
建议:
Illumina测序仪通量大 (up to 1.5 Tb per run),且准确率高 (with a typical error rate of 0.1–1%),通过在不同样本的序列上添加两重barcode,可以一次测序多个samples。
然而,Illumina测序仪存在carryover (between runs) 和 carry-between (within runs)的问题。最新的测序仪由于使用了新的扩增方法 (ExAmp),导致较高比例的‘index hopping’。
虽然没有一个明确的指导意见,告诉你在哪个特定的环境样品中应该测多大的覆盖度,但是一个基本的原则就是通量要尽可能地大,这样低丰度的细菌也能被测到。Illumina HiSeq 2500/4000, NextSeq 和 NovaSeq 的测序通量都很大,都适用于 metagenomics 的研究。
Metagenome de novo assembly 采用的策略与 whole-genome assembly 相同,均为 de Bruijn 图方法
用 de Bruijn 图方法进行宏基因组的从头组装时,面临着以下的挑战:
当进行单一基因组的组装时,其有一个前提假设:整个基因组的测序覆盖度是相对均匀的,这样就可以利用覆盖度信息来识别重复序列和鉴定测序错误和等位变异。
而metagenome中,各个组成基因组的覆盖度取决于它们的物种丰度,低丰度物种的基因组就会由于总体测序深度不够而使得最终组装出来的基因组是支离破碎的。使用更短的 k-mer 有助于低丰度基因组的组装,但是这会使得图中重复 k-mer 的频率大大增加,降低了组装的准确性。
这需要组装工具在考量低丰度物种与获得高丰度物种更长更准确的contig之间进行权衡,即选择合适的 k-mer :
同种细菌的不同菌株,它们的基因组组成很相近,常常就是一个碱基的变异或者整个基因/操纵子的丢失,当进行 de Bruijn 图组装时,就会在这些差异的位置出现分叉,组装工具在遇到这些分叉时,常常会停在这些位置,从而导致一个个不连续组装片段的产生。
Meta-IDBA:将图依据其拓扑结构拆分成各个元件,每个元件代表各个亚种的共有区域
解决计算能力与内存不足的策略:
Metagenome 组装完成后,我们得到的是成千上万的 contigs,我们需要知道哪些 contigs 来自哪一个基因组,或者都有哪些微生物的基因组。所以需要将 contigs 按照物种水平进行分组归类,称为 "bining"
一个很容易想到的策略就是,将组装得到的片段与已知物种的参考基因组进行比对,根据同源性进行归类。然而目前大多数的微生物的基因组还没有测序出来,因此限制了这种方法的可行性。
目前主流的 bining 策略利用的是 contigs 的序列组成特点。
依据:来自同一菌株的序列,其核酸组成是相似的
例如根 据核酸使用频率 (oligonucleotide frequency variations),通常是四核苷酸频率(tetranucleotide frequency), GC含量 和 必需的单拷贝基因 等
优势:即便只有一个样品的宏基因组数据也可以进行binning,这在原理上是可操作的
不足:由于很多微生物种内各基因型之间的基因组相似性很高,想利用1个样品的宏基因组数据通过核酸组成信息进行binning,效果往往并不理想或难度很大。利用核酸组成信息进行binning,基本上只适合那些群落中物种基因型有明显核酸组成差异的,例如低GC含量和一致的寡核苷酸使用频率
依据:来自同一个菌株的基因在不同的样品中 ( 不同时间或不同病理程度 ) 的丰度分布模式是相似的【PMID: 24997787】。
原因:比如,某一细菌中有两个基因,A和B,它们在该细菌基因组中的拷贝数比例为 A:B = 2:1,则不管在哪个样品中这种细菌的数量有多少,这两个基因的丰度比例总是为 2:1
优势:这种方法更有普适性,一般效果也比较好,能达到菌株的水平
不足:必须要大样本量,一般至少要50个样本以上,至少要有2个组能呈现丰度变化 ( 即不同的处理,不同的时间,疾病和健康,或者不同的采样地点等 ) ,每个组内的生物学重复也要尽量的多
对于像质粒这样的可移动遗传单元 (mobile genetic elements (MGEs)),由于其复制独立于细菌染色体,则同一种细菌的不同个体,该质粒的拷贝数可能存在差异,使得无法用丰度信息进行有效地bining
将核酸组成信息和丰度差异信息创建一个综合的距离矩阵,既能保证binning效果,也能相对节约计算资源,现在比较主流的binning软件多是同时依据核酸组成和丰度变化信息
依据:不同的细菌,其基因组甲基化模式不同,平均一种细菌有3种特意的甲基化 motif。MGEs (mobile genetic elements) 中含有 MTase 基因,其基因水平转移是细菌甲基化组多样性的驱动因素。虽然 MGEs 在不同个体的拷贝数不同,但是都存在,因此具有相同 MGEs 的细菌个体,其总遗传物质(包括染色体和 MGEs )都会受到相同的MTase的作用而得到相同的甲基化模式。
Q1:从哪些序列下手进行binning呢?
从原始的clean reads,还是从组装成的contig,还是从预测到的gene,都可以。根据基于聚类的序列类型的不同,暂且分为reads binning, contig binning和 genes binning
比较这三种binning的优劣:
总体来说应用最广泛的就是基于genes binning 和 contig binning
Genes binning的一般流程
在宏基因组做完组装和基因预测之后,把所有样品中预测到的基因混合在一起,去冗余得到unique genes集合,对这个unique genes集合进行binning,主要是根据gene在各个样品中的丰度变化模式,计算gene之间的相关性,利用这种相关性进行聚类
该图中的聚类过程类似于 K-means聚类 :随机选择几个seed genes作为诱饵,计算其他基因丰度分布模式与seed genes的相关性,按照固定的相关性值PCC>0.9,将它们归属于不同seed genes所代表的类,然后在聚好的类内重新选择seed genes,进行迭代,最终聚类得到一个个基因集合,较大的集合(超过700个基因)称为 metagenomic species (MGS),较小的集合称为 co-abundance gene group (CAG)
基于 bining 结果进行单菌组装:
比如对核酸组成信息的利用,开发得就不够充分,四碱基使用频率因简单而被广泛使用和接受,但现在已有研究表明k-mer丰度信息也是很好的种系特征,同时越长的k-mer含有越多的信息,还有基因和参考基因组间的同源关系也是有价值的种系信号,但这些都还没有被自动化的binning软件整合
想要获得高质量的bins经常需要手动调整
Taxonomic profiling: identifies which microbial species are present in a metagenome and estimates their abundance
优点:
当然它也有局限性:
对于与人类密切相关的样品,比如人肠道,可以使用该策略,而且已经有相关的成功实践
By looking at co-abundant markers from preassembled environment-specific gene catalogs
即前人研究 (MetaHIT consortium) 已经得出特定环境下的微生物的组成,这些微生物中有某些 co-abundant markers(这些 marker genes 的丰度与其物种的丰度成正比),这样就可以基于对这些 markers 的定量得到对应的物种丰度
选择 markers 的不同策略:
当样本量巨大,都进行组装是明显不切实际的,此时采用 marker-based approaches 是一个不错的选择;而且,如果该环境来源的样本其组成微生物是研究比较充分时,marker-based approaches 能得到比较准确的物种定量结果。
Gene identification
Characterization of the functional potential of the microbiome
局限性 : lack of annotations for accessory genes in most microbial species
因为在评估微生物群体的代谢潜能时,只对那些高度保守和 housekeeping 类型的功能进行了注释,这就解释了,为什么来自不同环境的不同样品,它们的功能特征常常是十分相似的,即使它们的物种组成有很大差异。
例如,鉴定出微生物群落中的抗生素抗性基因,该方法高度依赖特定功能相关基因集注释的质量。
参考资料:
(1) 魏子艳, 金德才, 邓晔. 环境微生物宏基因组学研究中的生物信息学方法[J]. 微生物学通报, 2015, 42(5):890-901.
(2) Quince C, Walker A W, Simpson J T, et al. Shotgun metagenomics, from sampling to analysis[J]. Nature Biotechnology, 2017, 35(9):833.
(3) 句句干货!一文读懂宏基因组binning
(4) Nielsen H B, Almeida M, Juncker A S, et al. Identification and assembly of genomes and genetic elements in complex metagenomic samples without using reference genomes[J]. Nature Biotechnology, 2014, 32(8):822-828.
(5) Sangwan N, Xia F, Gilbert J A. Recovering complete and draft population genomes from metagenome datasets[J]. Microbiome, 2016, 4(1):8.
(6) Abubucker, S. et al. Metabolic reconstruction for metagenomic data and its application to the human microbiome. PLoS Comput. Biol. 8, e1002358(2012).
(7) Beaulaurier J, Zhu S, Deikus G, et al. Metagenomic binning and association of plasmids with bacterial host genomes using DNA methylation.[J]. Nature Biotechnology, 2017, 36(1).
H. 简述完整基因组自上到下的测序过程必要的三个步骤
全基因组测序是指对全部基因组完整测序 是决定一套完整染色体基因组上核苷酸碱基准确顺序组成的过程 基因组测序是一项庞大的工程 其中三个必需关键技术是DNA大片段的克隆 测序的自动化和用生物信息学处理数据 当一个基因组完成测序之后 应该对其进行注释 全基因组测序主要应用于癌症
通常采用Fleischman等人在测定流感嗜血杆菌基因组中建立的 鸟枪随机测序法 而后再将测序结果进行整合
1.建立随机DNA文库
通过喷雾器进行机械剪切或使用超声波处理纯度高 完整性好的基因组DNA 制备随机片段 然后将随机片段插入到适宜的测序载体中 随机DNA文库建立后 对文库的随机质量和容量进行鉴定
2.高通量测序
最大限度地从文库中随机挑取克隆制备测序模板 并使用多台自动化测序仪进行高通量测序
3.随机片段的组装
将测序结果输入计算机 使用软件根据重叠序列将随机片段组装起来以还原整个基因组序列
4.缺口的补平
由于使用的是随机片段 因此在组装过程中可能出现物理缺口 对于这种情况 可以根据缺口两边的序列设计引物 以完整的DNA为模板进行PCR扩增 得到缺失部分的序列
I. 序列组装需不需要序列之间存在重叠部分
基因组组装(Genome assembly)是生物信息学领域的核心问题,基因组组装就是把序列测序产生的读取片段reads经过序列拼接组装,生成基因组的碱基序列。基因组组装软件可根据得到的所有读长组装成基因组。基因组组装这个步骤对于基因组分析是十分关键的,因为目前二代测序技术获得的测序序列一般都较短,需要组装拼接成较长的完整的序列用于进一步分析,例如长序列能提高物种注释分析的准确性。
宏观来说,基因组组装可以分为从头组装(De novo assembly) 和映射比对组装(mapping assembly), 从头组装是指不需要依靠任何已知的基因组信息,反过来,映射比对组装就是需要把测序序列和参考基因组来比对,找到序列的对应位置再进行组装,本文主要讲解的从头组装。 当然两种都有各个的用处,映射比对组装也有一些算法例如BWT算法。
由于目前组装技术的限制和实际情况的复杂性,最终组装得到的序列与真实基因组序列之间仍可能存在差异,甚至只能得到若干条无法进一步连接起来的序列。对组装效果的评价主要依由于据组装序列的连续性、完整性和准确性。连续性要求组装得到的(多条)序列长度尽可能长;完整性要求组装序列的总长度占基因组序列长度的比例尽可能大;准确性要求组装序列与真实序列尽可能符合。
目前基因组组装一般有基于OLC(Overlap-Layout-Consensus, 先重叠后扩展)和基于De Brujin Graph(DBG)两种组装算法,基于OLC的组装方法适合长序列组装,运行依赖的数据结构需要消耗大量的内存,且运行速度比较慢,错误率高,而DBG组装方法内存消耗相对较低,运算速度快,且准确率高。目前主流的基因组装算法都是基于后者改进设计的。
基本概念
在开始之前,有几个名词需要说明下:
reads:就是我们测序产生的短读序列,通常一代和三代的reads读长在几千到几万bp之间,二代的相对较短,平均是几十到几百bp。
contig:中文叫做重叠群,就是不同reads之间的overlap交叠区,拼接成的序列就是contig
scaffold:这是比contig还要长的序列,获得contig之后还需要构建paired-end或者mate-pair库,从而获得一定片段的两端序列,这些序列可以确定contig的顺序关系和位置关系,最后contig按照一定顺序和方向组成scaffold,其中形成scaffold过程中还需要填补contig之间的空缺。