CN115810395B - 一种基于高通量测序动植物基因组t2t组装方法 - Google Patents
一种基于高通量测序动植物基因组t2t组装方法 Download PDFInfo
- Publication number
- CN115810395B CN115810395B CN202211549592.3A CN202211549592A CN115810395B CN 115810395 B CN115810395 B CN 115810395B CN 202211549592 A CN202211549592 A CN 202211549592A CN 115810395 B CN115810395 B CN 115810395B
- Authority
- CN
- China
- Prior art keywords
- genome
- data
- assembly
- sequencing
- ont
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Landscapes
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本发明公开了一种基于高通量测序动植物基因组T2T组装方法，所述方法包括：S1.基于样本测序，分别进行ONT超长测序数据组装及HiFi测序数据组装；S2.基于基因组大小评估结果对两条路线的组装结果去冗余；S3.分别对两条路线组装结果进行HiC挂载，得到染色体级别参考基因组其gap；S4.对ONT线路基因组着丝粒区域进行三代纠错，基于二代数据对ONT线路基因组进行二代纠错；S5.基于原始数据比对末端片段，根据端粒重复基序列的丰度大小进行重组装，再经比对进行末端端粒序列替换，对两条线路的端粒延伸；S6.分别对基因组序列补gap、HiFi纠错；S7.结合已知物种基因组排序对染色体进行调整得到组装后的基因组；所述T2T基因组端粒完整、无空缺，推动了动植物的研究进展。
Description
技术领域
本发明涉及基因组装技术领域，尤其涉及一种基于高通量测序动植物基因组T2T组装方法。
背景技术
高通量测序技术(High-throughputsequencing)又称下一代测序技术(Next-generationsequencingtechnology)。高通量测序技术目前已经来到***，纳米孔测序技术(又称***测序技术)是最近几年兴起的新一代测序技术。目前测序长度可以达到4Mb。这项技术开始于90年代，经历了三个主要的技术革新：一、单分子DNA从纳米孔通过；二、纳米孔上的酶对于测序分子在单核苷酸精度的控制；三、单核苷酸的测序精度控制。随着物种基因组研究的不断发展，测序技术的不断提升，物种基因组的连续性和完整性也都得到大幅提升。物种基因组从普通二代测序组装的draftgenome(草图基因组，基因组1.0)，到使用PacBio(PacificBiosciences公司)或ONT(Oxfordnanopore公司)(R9.4)测序技术结合Hi-C(High-throughput/resolutionchromosomeconformationcapture，高通量染色体构象捕获技术)组装的high-qualitygenome(高质量基因组，基因组2.0)，再到由ONT(N50>50Kb或R10)测序技术组装的nearcompletegenome(近完全基因组，基因组3.0)，以及结合ONTultra-long(牛津纳米孔超长测序技术)和PacBioHiFi(HighFidelityreads,高准确度读长测序技术)技术组装的T2Tgenome(端粒到端粒基因组，基因组4.0)。动植物T2T基因组，指通过ONTultra-longN50>100Kb(测序读长N50大于1000000，N50：N50是基因组拼接之后一个评价指标，将拼接得到的所有的序列，根据序列大小从大到小进行排序，然后逐步开始累加，当加和长度超过总长一半时，加入的序列长度即为N50长度。)结合HiFi和二代数据(高通量测序，如Illumina HiSeqTM/MiseqTM)进行混合组装，得到的有一条或者多条染色体达到端粒到端粒(Telomere-to-Telomere，T2T)的水平基因组，T2T基因组完成图是基因组组装的终极目标。
T2T基因组组装的要求包括：染色体组装无空缺，QV值(基因组组装准确度评估标准)大于40，BUSCO(基因组组装完整度评估标准)评估大于95％，端粒完整。以往单一测序手段ONT普通测序或者PacBioHiFi都难以达到T2T水平。普通的ONT测序组装精度偏低、难以达到QV大于40的质量标准，且难以达到染色体组装无空缺；HiFi测序组装虽然能带来很高的组装精度，但是无法达到染色体组装无空缺，且着丝粒区域往往组装空缺。
因此，有必要针对现有技术的缺陷，提出一种基于高通量测序动植物基因组T2T组装方法。
发明内容
本发明的目的在于提供一种基于高通量测序动植物基因组T2T组装方法，通过ONT超长测序技术结合HiFi测序技术和高标准、高质量的组装技术能克服以往的缺点，基因组组装达到T2T水平，推动了整个物种的研究进展。
有鉴于此，本发明的方案如下：
一种基于高通量测序动植物基因组T2T组装方法，包括如下步骤：
S1.基于样本测序，分别进行ONT超长测序数据组装及HiFi测序数据组装；
S2.基于基因组大小评估结果对两条路线的组装结果去冗余；
S3.分别对两条路线组装结果进行HiC挂载，得到染色体级别参考基因组其gap；
S4.对ONT线路基因组着丝粒区域进行三代纠错，基于二代数据对ONT线路基因组进行二代纠错；
S5.基于原始数据比对末端片段，根据端粒重复基序列的丰度大小进行重组装，再经比对进行末端端粒序列替换，对两条线路的端粒延伸；
S6.分别对基因组序列补gap、HiFi纠错；
S7.结合已知物种基因组排序对染色体进行调整得到组装后的基因组。
在本发明的一个实施例中，所述步骤S4中三代纠错过程为：对ONT基因组着丝粒区域使用k-mer锚定方法进行迭代多轮三代纠错。
作为优选的实施例，所述ONT基因组着丝粒区域基于二代数据和基因组中k-mer的频次，标记在组装中出现一次且在二代数据中出现14到46次的21-kmer。
在本发明的一个实施例中，所述步骤S4中二代纠错过程为：将二代原始数据进行切分并比对到基因组上，通过深度神经网络进行变异检测得到变异信息，变异信息经过滤对齐后导出一致性序列从而得到二代纠错后的基因组。
在本发明的一个实施例中，所述步骤S5端粒延伸过程为：
S51.将每条染色体分别与原始数据进行比对，收集在染色体末端定长内比对一次的所有片段，计算每条片段中出现端粒重复基序列的次数，定义出现最多次数的片段为ref，其他为query，将ref和query重新组装，得到一致性序列；
S52.将该一致性序列分别比对到每条染色体上，取最佳比对结果对末端端粒序列进行替换。
作为优选的实施例，所述步骤S52中，如果identity低于80阈值或比对上区域不在染色体末端20Kbp，则不进行替换。
在本发明的一个实施例中，步骤S6中所述基因组序列补gap步骤为：基于补洞程序对基因组进行补洞，将补洞数据与基因组数据进行比对，按其他纠错后基因组版本、HiFi数据、ONT数据的先后顺序分别进行补洞。
作为优选的实施例，所述补gap步骤中，在数据量不够的情况下，用subreads/ont原始片段对基因组进行补洞；如果比对上的位置能刚好跨过gap两端，则选取比对上的最长长度区域的最佳比对区域用补gap数据对基因组上包含gap区域的序列进行替换。
在本发明的一个实施例中，步骤S6中所述基因组HiFi纠错步骤为：先过滤掉低于10kbp的HiFi数据，将过滤后的数据与补完洞的基因组进行比对得到合并排序后的文件；再分别过滤掉二次比对、嵌合比对的片段后进行三代纠错。
在本发明的一个实施例中，所述T2T组装方法还包括共线性分析步骤：对步骤S7中所述调整后的基因组与已发表的物种的基因组做共线性比对分析，分析指标包括完整度和准确度。
相比现有技术，本发明的有益效果包括但不限于：
1.本发明提出的基于高通量测序动植物基因组T2T组装方法通过ONT超长测序技术结合HiFi测序技术和高标准、高质量的组装技术能克服以往的缺点，基因组组装达到T2T水平，得到端粒完整、无空缺的基因组，推动了整个物种的研究进展。
2.本发明提出的基于高通量测序动植物基因组T2T组装方法基于多种测序手段并辅以复杂的开源软件与自研软件达到超出单一测序手段简单组装的效果，为基因组T2T组装提供新的选择。
附图说明
图1为本发明所述动植物基因组T2T组装方法流程图。
图2为本发明所述HiC具体具体步骤流程图。
图3为本发明实施例中基于contig互作强度与位置关系得到的互作热图。
图4为本发明实施例中最终T2T版本Hi-C互作热图。
图5为本发明实施例中对T2T组装方法得到的基因组共线性分析结果示意图。
具体实施方式
为了使本发明的目的、技术方案和有益技术效果更加清晰明白，以下结合具体实施方式，对本发明进行进一步详细说明。应当理解的是，本说明书中描述的具体实施方式仅仅是为了解释本发明，并不是为了限定本发明。
在一个实施例中，提出一种基于高通量测序动植物基因组T2T组装方法，流程如图1所示，步骤如下：
1.测序：同一份动植物生物样本需要进行二代测序(illumina)或者BGI(华大基因)、ONT超长测序(OxfordNanoporeTechnology平台仪器测序)、HiC测序(Hi-C技术源于基因组捕获技术)、PacBioHiFi测序(PacificBiosciences)。
2.评估：通过二代测序进survey基因组分析，估计出基因组大小。
3.初组装：分两条线进行。
3.1以ONT超长测序数据组装结果为骨架(后简称ONT骨架)，分别使用Flye、Nextdenovo及Necat三款基因组组装软件进行组装；
Flye：https://github.com/fenderglass/Flye；
Nextdenovo：https://github.com/Nextomics/NextDenovo；
Necat：https://github.com/xiaochuanle/NECAT。
3.2以HIFI测序数据组装结果为骨架(后简称HIFI骨架)，使用hifiasm(https://github.com/chhylp123/hifiasm)进行组装。
4.去冗余与评估
组装后基因组大小有超过survey评估大小，为把基因组减少到评估大小，需对两条路线的组装结果进行去冗余；使用merqury(https://github.com/search？q＝merqury)进行QV值计算(组装后基因组准确度评估值)；使用BUSCO(https://busco.ezlab.org/)进行组装后基因组完整度评估；使用reads回比到组装后基因组；去除低质量contig(组装后的序列)；在NCBI上下载本物种线粒体、叶绿体，并使用minimap2(https://github.com/lh3/minimap2)比对组装后基因组，去除比对上的contig。
5.HiC挂载
5.1ONT骨架路线：选取三个组装结果中最好的数据进行HiC挂载(挂载就是利用3维情况下HiC互做,指导二维基因组组装的一种应用)；二、HiFi骨架路线：HiC挂载。经过此步得出染色体数目与基因组Gap(未组装出来的序列，在基因组染色体层面体现为序列空白以100个N填充)数量。HiC挂载具体具体步骤如图2所示，为本领域常用的路线。
6.本步骤应用于ONT骨架路线
三代纠错：先将去线粒体叶绿体低reads支持contig流程中生成的中间结果按40Kbp长度过滤后的三代filterreads比对到hic挂载后的基因组上。通过meryl(https://github.com/marbl/meryl)来统计二代数据的kmer出现次数和基因组中kmer的出现次数，找到在组装中出现一次且在二代数据中出现14到46次的21-kmer，并认为是需要标记的着丝粒区域。用T2T_polish(https://github.com/arangrhie/T2T-Polish)下的子程序filter_by_marker_nosplit.sh对bam文件进行标记，并用medaka consensus(https://github.com/nanoporetech/medaka)进行纠错，medaka stitch得到一致性序列，以达到避免着丝粒过度纠错的目的。
二代纠错：将二代原始数据进行切分，用bwa(https://github.com/lh3/bwa)分别比对到基因组上，然后sort，merge得到总的sorted.merged.bam，并建立索引文件。用deepvariant(https://github.com/***/deepvariant)进行变异检测得到vcf文件，经过过滤对齐后的vcf文件被bcftoolsconsensus(https://github.com/samtools/bcftools)导出一致性序列从而得到二代纠错后的基因组。
7.端粒延伸应用于两条路线
将每条染色体分别与原始数据用winnowmap(https://github.com/marbl/Winnowmap)比对，收集在染色体末端50bp内比对一次的所有reads，计算每条read中出现端粒重复基序(‘CCCTAAA’/‘TTTAGGG’)的次数，定义出现最多次数的read为ref，其他为query，使用medaka_consensu将ref和query重新组装，得到一致性序列。将该一致性序列分别比对到每条染色体上，取最佳比对结果对末端端粒序列进行替换。如果identity低于80阈值或比对上区域不在染色体末端20Kbp，则不进行替换。
8.补gap
Gap_closer流程用于对基因组上gap的填补，用winnowmap将补洞数据(不含N)与基因组数据比对(含N)，流程分为三个水平对gap进行填补，其优先级为：其他纠错后基因组版本>HiFi数据>ONT数据。其他版本纠错后的基因组可以是nextdenovo、necat、hifiasm、canu…等任何版本基因组；HiFi数据最好为ccs.fasta，既去环化后的HiFi数据；ONT数据最好为cns.fa，即ont数据经过nextdenovo自纠错后生成的一致性序列consensus.fa。在数据量不够的情况下，可以考虑选用subreads/ont原始reads对基因组进行补洞。如果比对上的位置能刚好跨过gap两端，则选取比对上的最长长度区域的最佳比对区域用补gap数据对基因组上包含gap区域的序列进行替换。从理论上说，只要所提供补洞用的HiFi/ONT数据量足够大，该流程就能把gap给补上。但是流程也是有一定局限性的，对于很短一片区域存在多个gap时不太适用，因此gap向两端延申的长度不建议超过两个gap间的最短距离。
9.HIFI纠错
上述步骤中做了填补gap的操作，为了使得我们填补上的序列可靠，以及保证基因组质量的一致性，这一步使用HiFi数据对补洞后的基因组进行纠错。过滤掉HiFi数据10Kbp以下的原始数据，用winnowmap将过滤后的数据与补完洞的基因组进行比对，得到比对完成的merge.sorted.bam文件后。用samtools view-F256(https://github.com/samtools/samtools)对bam文件进行过滤，并用falconcbam-filter-clipped(https://github.com/PacificBiosciences/FALCON)过滤掉嵌合比对片段，最后用racon(https://github.com/isovic/racon)的liftover分支进行纠错，得到经过HiFi纠错后的基因组序列。
10.染色体排序与共线性
将上述组装结果mapping到已发表的物种基因组上，并将染色体编号和方向按照已发表的基因组染色体编号和方向进行对应，然后将调整之后的基因组与已发表的物种的基因组做共线性。同样，调整位置后的基因组和ref也可以输入到共线性流程中做mummer共线性图或jcvi共线性图，这样就能得到一条直线或者是瀑布式的共线性图，可以更好的看出两版本之间共线性的情况。
11.基因组最终确认
两条线路中选择完整度高、准确度高的基因组为最终基因组。
实施例
1.使用太平洋生物公司(PacBio)高保真测序数据(HiFi)进行基因组组装，生成重叠群(contig)。使用从美国国立生物技术信息中心(NCBI)下载了大约20Gb水稻太平洋生物公司(PacBio)高保真测序(HiFi)测序数据。使用hifiasm(https://github.com/chhylp123/hifiasm)进行组装。组装结果如表1所示。
表1：hifiasm初步组装结果表
| num_seqs | sum_len | min_len | avg_len | max_len | N50 |
| 141 | 401950498 | 9114 | 2850712.8 | 45002514 | 31388662 |
2.使用从美国国立生物技术信息中心(NCBI)下载了大约50Gb水稻ONT超长数据，分别使用三款基因组组装软件进行组装Flye(https://github.com/fenderglass/Flye)、Nextdenovo(https://github.com/Nextomics/NextDenovo)、Necat(https://github.com/xiaochuanle/NECAT)进行组装。组装结果如表2所示。
表2：ONT数据初步组装结果表
| num_seqs | sum_len | min_len | avg_len | max_len | N50 | |
| Nextdenovo | 23 | 398481786 | 110804 | 17325295.0 | 44949911 | 31426128 |
| flye | 77 | 395339674 | 1472 | 5134281.5 | 52039294 | 24103585 |
| necat | 61 | 401818451 | 23811 | 6587187.7 | 44884577 | 24129713 |
3.对于高杂合物种，需要使用软件Purge_haplotigs(v1.0.4；https://github.com/skingan/purge_haplotigs_multiBAM)/Purge_dups(v1.2.5；https://github.com/df guan/purge_dups)对基因组进行去杂合。
4.使用minimap2(2.17-r941)(H.Li 2018)软件，比对线粒体、叶绿体，去除碱基比对超过50％的序列；通过blast refseq库去除细菌污染；去除低reads支持的contig(Minimap2将ONT reads(>＝40kbp)与contig比对，如果一条contig上超过50％的位点的深度低于15，就进行移除)，该阈值会根据数据量适当调整。
5.利用互作关系进行辅助组装：
1)Contig聚类
使用Hi-C互作关系，确定有效数据中不同Contig间关联的紧密程度，对Contig进行聚类。对于核型为2n的基因组草图，利用ALLHiC(v0.9.8)(Zhang etal.2019)软件通过agglomerativehierarchicalclustering(自下而上的层次聚类算法)，将草图的Contig序列聚类为n个染色体群。
2)Contig定序和定向
利用ALLHIC(v0.9.8)(Zhangetal.2019)软件对n个染色体群内部的Contig进行定序和定向，再通过软件3D-DNA(v180419)(Dudchenkoetal.2017)和juicer(v1.6)(NevaC.etal.2016b)将Contig两两之间的互作关系转化为指定的二进制文件(即.hic文件)。再通过windows软件Juicebox(v1.11.08)(NevaC.etal.2016a)对已经定序和定向Contig(.assembly)进行手动定序与定向(.review.assembly)。
3)利用互作关系进行Contig去冗余
若某段序列与一相同大小的序列完全没有互作，但与其它序列互作正常，那么它极大可能为杂合序列，需要手动去除。
将定序、定向、去冗余的Contig序列，使用100个N补Gap得到最终的染色体级别基因组序列。详细每条染色体及其长度统计结果如表3所示，没有互作或者有互作噪音的Contig为未挂载的片段(记为chrUnn)。
4)染色体挂载结果统计：经过Hi-C组装和人工调整后，共有397,785,342bp的序列长度的基因组序列被定位到12条染色体上，占比99.25％。每条染色体及其长度的详细统计结果如表3所示。
表3：染色体长度统计表
| Chromosome | Length(bp) | contignumber |
| chr1 | 25,815,008 | 6 |
| chr2 | 31,882,477 | 1 |
| chr3 | 38,842,104 | 2 |
| chr4 | 32,745,327 | 13 |
| chr5 | 32,075,978 | 5 |
| chr6 | 24,932,023 | 1 |
| chr7 | 37,439,256 | 3 |
| chr8 | 30,588,819 | 2 |
| chr9 | 27,101,476 | 2 |
| chr10 | 31,388,662 | 1 |
| chr11 | 45,002,514 | 1 |
| chr12 | 39,971,698 | 1 |
| chrUnn | 3,021,772 | 73 |
5)使用软件HiCExplorer(v3.6)(Joachim et al.2020)基于contig互作强度与位置关系进行绘图，如图3所示。
6.组装纠错
6.1三代纠错：
针对着丝粒区域使用k-mer锚定方法进行三代纠错(注意此纠错手段仅针对T2T基因组ont骨架版)；
1)使用minimap2(v2.17-r941)(H.Li 2018)将ONT reads(>＝40kbp)比对到伪染色体上。
2)为了避免着丝粒序列的过度纠错，找到在组装中出现一次且在Illumina reads中出现14到46次的21-kmer，并使用T2T-Polish
(https://github.com/malonge/T2T-Polish)的子程序filter_by_marker_nosplit.sh进行对比对好的bam进行标记。
3)使用medakaconsensus(v1.5.0；参数：–modelr941_prom_hac_g507
–batch_size200；https://github.com/nanoporetech/medaka)进行纠错。
4)使用medakastitch(v1.5.0；https://github.com/nanoporetech/medaka)得到三代纠错后的一致性序列。
重复上述步骤，即T2T基因组ONT骨架版需要迭代两轮三代纠错。
6.2二代纠错：
使用deepVariant(v1.3.0)(Poplinetal.2018)进行二代纠错(注意此纠错手段仅针对T2T基因组ONT骨架版)；
1)使用bwamem(v0.7.17-r1188)(HengLi2013)，samtools(v1.9)
(Daneceketal.2021)得到三代纠错后的基因组和二代数据比对好的bam文件。
2)使用deepvariant(v1.3.0；参数：–model_type＝WGS)(Poplinetal.20
18)得到vcf文件。
3)使用bcftools(v1.15；参数：view-e‘type＝“ref”’，view-i‘QUAL>1
&&(GT＝“AA”||GT＝“Aa”)’)(Daneceketal.2021)对vcf文件进行过滤，并用bcftoolsnorm(Daneceketal.2021)将vcf文件左对齐并且规范化。
4)将vcf文件压缩并生成索引文件。
5)bcftoolsconsensus(Daneceketal.2021)导出二代纠错后的一致性序列。
7.端粒延伸
1)使用Winnowmap(v1.11,参数：k＝15,–MD)(Chiragetal.2020)将所有的ontreads比对到ref,收集在染色体末端50bp内比对一次的所有r eads。
2)计算每条read中端粒重复基序(‘CCCTAAA’/‘TTTAGGG’)的出现次数(端粒数据库：http://telomerase.asu.edu/sequences_telomere.html)，并定义出现最多次数的read为ref，其他为query。
3)medaka_consensu(v1.2.1，参数：-mr941_min_high_g360；https://github.
com/nanoporetech/medaka)，将ref端粒read和query端粒read进行重新组装，得到一致性序列。
4)用nucmer(v3.1)(Kurtzetal.2004)将上述端粒的一致性序列分别比对到每条染色体上，取最佳比对结果对末端端粒序列进行替换。如果identity低于80阈值或比对上区域不在染色体末端20Kbp，则不进行替换。
8.gap填补
用于对基因组上gap的填补，用winnowmap(v1.11,参数：k＝15,–MD)
(Chiragetal.2020)将补洞数据(不含N)与基因组gap区间比对(含N)，该步骤分为三个水平对gap进行填补，其优先级为：其他纠错后基因组版本>hifi数据>ont数据。其他版本纠错后的基因组可以是nextdenovo、necat、hifiasm、canu…等任何版本基因组；hifi数据最好为ccs.fasta，既去环化后的hifi数据；ont数据最好为cns.fa，即ont数据经过nextdenovo自纠错后生成的一致性序列consensus.fa。再数据量不够的情况下，可以考虑选用subreads/ontreads对基因组进行补洞。如果比对上的位置能刚好跨过gap两端，则选取比对上的最长长度区域的最佳比对区域用补gap数据对基因组上包含gap区域的序列进行替换。从理论上说，只要所提供补洞用的hifi/ont数据量足够大，用于补洞的基因组其他组装版本足够多，就能得到0gap基因组。
9.HiFi reads纠错
1)使用Winnowmap2(Chirag et al.2020)将>＝10kbp的HiFi reads与填补后gap的版本基因组比对(参数：k＝15greater-than distinct＝0.9998
–MD-ax map-pb)；
2)使用samtools“view”(v1.10,参数：-F 256)(Danecek et al.2021)对比对片段进行过滤；
3)使用“falconc bam-filter-clipped”来删除嵌合比对片段(-t-F 0x104)；
4)使用这些过滤后的比对信息，使用racon的特殊分支进行纠错，(v1.6.0,
-L-u，https://github.com/isovic/racon/tree/liftover)。
10.T2T基因组评估
1)通过基因组上gap的位置和个数来评估基因组组装的连续性，基因组gap情况如表4所示。
表4：0gap基因组
2)利用软件bwa(version：0.7.17-r1188)将二代高通量测序(如Illumina测序平台)得到的质控后reads(clean reads)比对(mapping)到参考基因组上。统计二代比对率和覆盖度如表5所示。
表5：
注：Sample name：样本名称；Mapping rate(％)：比对率；Paired mapping rate(％)：双端比对率；Average sequencing depth：平均测序深度；Coverage(％)：覆盖度；Coverage at least 4X(％)：4X深度以上的覆盖度；Coverage at least 10X(％)：10X深度以上的覆盖度；Coverage at least 20X(％)：20X深度以上的覆盖度。
3)使用软件BUSCO(version：4.1.4；参数：–evalue 1e-05)来评估基因组组装的完整性，结果如表6所示。在基因组完整性较好时，T2T基因组中B USCO评估中CompleteBUSCOs(C)的占比一般大于95％。
表6：基因组组装的完整性评估结果
注：Complete BUSCOs为完整的BUSCO；Complete and single-copy BUSCOs为完整且单拷贝的BUSCO；Complete and duplicated BUSCOs为完整但多拷贝的BUSCO；Fragmented BUSCOs为不完整的BUSCO；Missing BUSCOs为缺失的BUSCO；Total BUSCOgroups searched为本次使用BUSCO库中收录的保守蛋白基因的总数目。
4)通过二代数据的kmer统计结果和基因组的kmer统计结果，结果如表7所示。比较其kmer种类数频数以及分布来评判基因组组装的质量和准确度。
表7：基因组组装QV值统计表
| qv | completeness(％) | length | N50 | ReadsNum |
| 45.1739 | 99.1838 | 397,733,901 | 32,070,994 | 12 |
注：qv为用于评估基因组质量的qv值；completeness用于评估基因组完整度；length为基因组长度；N50为contig级别的N50，即将所有的contigs按照从长到短排序后依次累加，当相加的长度达到基因组总长度的一半时，最后一个加上的contig长度即为N50；Reads Num为contig级别基因组reads条数。
5)根据端粒区域碱基重复序列(5’端为CCCTAAA，3’端为TTTAGGG)，对组装结果进行端粒序列鉴定，结果如表8所示。
表8：端粒重复序列统计结果
/>
6)着丝粒序列鉴定
(I)ChIP-Seq
建立基于CenH3抗体的ChIP-seq的着丝粒分离分析方法，在获得着丝粒ChIP-seq数据后，通过与基因组序列比对筛选，获得在基因组上具有单一拷贝的数据，并将其锚定到对应物种的基因组上，即可获得ChIP-seq富集的峰图，该峰图则指示着丝粒区。这一方法目前已被用于着丝粒的基因组精确定位。
(II)FISH验证
结合细胞学的荧光原位杂交技术(FISH)验证，最终可鉴定出着丝粒区域特异的重复序列。
11.最终T2T版本Hi-C互作热图
在经过端粒延申，补gap，纠错后，以最终T2T版本基因组，再一次通过Hi-C数据来得到新的染色体级别的Hi-C互作热图，如图4所示，通过热图中互作信号的异常与否来鉴定我们前面的操作是否正确。
12.共线性分析
本次组装的基因组与目标物种历史基因组版本进行比较；该分析仅限有历史基因组版本的物种，共线性分析结果如图5所示。该历史基因组版本为另一种T2T组装方法得到的基因组。
本发明并不仅仅限于说明书和实施方式中所描述，因此对于熟悉领域的人员而言可容易地实现另外的优点和修改，故在不背离权利要求及等同范围所限定的一般概念的精神和范围的情况下，本发明并不限于特定的细节、代表性的方案和这里描述的示例。
Claims (9)
1.一种基于高通量测序动植物基因组T2T组装方法，其特征在于，包括如下步骤：
S1.基于样本测序，分别进行ONT超长测序数据组装及HiFi测序数据组装；
S2.基于基因组大小评估结果对两条路线的组装结果去冗余；
S3.分别对两条路线组装结果进行HiC挂载，得到染色体级别参考基因组其gap；
S4.对ONT线路基因组着丝粒区域进行三代纠错，基于二代数据对ONT线路基因组进行二代纠错；
S5.基于原始数据比对末端片段，根据端粒重复基序列的丰度大小进行重组装，再经比对进行末端端粒序列替换，对两条线路的端粒延伸；
S6. 分别对基因组序列补gap、HiFi纠错；
S7. 结合已知物种基因组排序对染色体进行调整得到组装后的基因组；
所述步骤S5 端粒延伸过程为：
S51.将每条染色体分别与原始数据进行比对，收集在染色体末端定长内比对一次的所有片段，计算每条片段中出现端粒重复基序列的次数，定义出现最多次数的片段为ref，其他为query，将 ref 和 query 重新组装，得到一致性序列；
S52.将该一致性序列分别比对到每条染色体上，取最佳比对结果对末端端粒序列进行替换。
2.根据权利要求1所述的动植物基因组T2T组装方法，其特征在于，所述步骤S4中三代纠错过程为：对ONT基因组着丝粒区域使用 k-mer 锚定方法进行迭代多轮三代纠错。
3.根据权利要求2所述的动植物基因组T2T组装方法，其特征在于，所述ONT基因组着丝粒区域基于二代数据和基因组中k-mer的频次，标记在组装中出现一次且在二代数据中出现14到46次的21-kmer。
4.根据权利要求1所述的动植物基因组T2T组装方法，其特征在于，所述步骤S4中二代纠错过程为：将二代原始数据进行切分并比对到基因组上，通过深度神经网络进行变异检测得到变异信息，变异信息经过滤对齐后导出一致性序列从而得到二代纠错后的基因组。
5.根据权利要求1所述的动植物基因组T2T组装方法，其特征在于，所述步骤S52中，如果identity低于80阈值或比对上区域不在染色体末端20Kbp，则不进行替换。
6.根据权利要求1所述的动植物基因组T2T组装方法，其特征在于，步骤S6中所述基因组序列补gap步骤为：基于补洞程序对基因组进行补洞，将补洞数据与基因组数据进行比对，按其他纠错后基因组版本、HiFi数据、ONT数据的先后顺序分别进行补洞。
7.根据权利要求6所述的动植物基因组T2T组装方法，其特征在于，所述补gap步骤中，在数据量不够的情况下，用subreads/ont原始片段对基因组进行补洞；如果比对上的位置能刚好跨过gap两端，则选取比对上的最长长度区域的最佳比对区域用补gap数据对基因组上包含gap区域的序列进行替换。
8.根据权利要求1所述的动植物基因组T2T组装方法，其特征在于，步骤S6中所述基因组HiFi纠错步骤为：先过滤掉低于10kbp的HiFi数据，将过滤后的数据与补完洞的基因组进行比对得到合并排序后的文件；再分别过滤掉二次比对、嵌合比对的片段后进行三代纠错。
9.根据权利要求1所述的动植物基因组T2T组装方法，其特征在于，还包括共线性分析步骤：对步骤S7中所述调整后的基因组与已发表的物种的基因组做共线性比对分析，分析指标包括完整度和准确度。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211549592.3A CN115810395B (zh) | 2022-12-05 | 2022-12-05 | 一种基于高通量测序动植物基因组t2t组装方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211549592.3A CN115810395B (zh) | 2022-12-05 | 2022-12-05 | 一种基于高通量测序动植物基因组t2t组装方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115810395A CN115810395A (zh) | 2023-03-17 |
| CN115810395B true CN115810395B (zh) | 2023-09-26 |
Family
ID=85484938
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211549592.3A Active CN115810395B (zh) | 2022-12-05 | 2022-12-05 | 一种基于高通量测序动植物基因组t2t组装方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115810395B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116343919B (zh) * | 2023-04-11 | 2023-12-08 | 天津大学四川创新研究院 | 一种全基因组图谱绘制测序方法 |
| CN117672354B (zh) * | 2023-12-21 | 2024-05-28 | 北京诺禾致源科技股份有限公司 | 比较哺乳动物近源物种完整基因组组装质量的方法和装置 |
Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105121661A (zh) * | 2013-02-01 | 2015-12-02 | 加利福尼亚大学董事会 | 用于基因组组装及单体型定相的方法 |
| CN107250356A (zh) * | 2014-12-16 | 2017-10-13 | 加文医学研究所 | 测序对照 |
| CN109326323A (zh) * | 2018-09-13 | 2019-02-12 | 北京百迈客生物科技有限公司 | 一种基因组的组装方法及装置 |
| CN111584004A (zh) * | 2020-05-12 | 2020-08-25 | 西藏自治区农牧科学院水产科学研究所 | 一种基于三维组学数据的西藏特色鱼类基因组组装方法 |
| CN112996927A (zh) * | 2018-10-31 | 2021-06-18 | 罗格斯新泽西州立大学 | Gramc：顺式调节模块的基因组规模报道子测定方法 |
| CN113151426A (zh) * | 2021-04-16 | 2021-07-23 | 中国农业科学院兰州畜牧与兽药研究所 | 一种基于三代PacBio和Hi-C技术组装和注释霍巴藏绵羊基因组的方法 |
| CN113767438A (zh) * | 2019-02-28 | 2021-12-07 | 加利福尼亚太平洋生物科学股份有限公司 | 使用均聚物折叠测序读段改进对齐 |
| CN113808668A (zh) * | 2021-11-18 | 2021-12-17 | 北京诺禾致源科技股份有限公司 | 提升基因组组装完整性的方法、装置及其应用 |
| CN115148290A (zh) * | 2022-08-12 | 2022-10-04 | 武汉希望组生物科技有限公司 | 一种基于三代测序数据的补洞方法 |
-
2022
- 2022-12-05 CN CN202211549592.3A patent/CN115810395B/zh active Active
Patent Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105121661A (zh) * | 2013-02-01 | 2015-12-02 | 加利福尼亚大学董事会 | 用于基因组组装及单体型定相的方法 |
| CN107250356A (zh) * | 2014-12-16 | 2017-10-13 | 加文医学研究所 | 测序对照 |
| CN109326323A (zh) * | 2018-09-13 | 2019-02-12 | 北京百迈客生物科技有限公司 | 一种基因组的组装方法及装置 |
| CN112996927A (zh) * | 2018-10-31 | 2021-06-18 | 罗格斯新泽西州立大学 | Gramc：顺式调节模块的基因组规模报道子测定方法 |
| CN113767438A (zh) * | 2019-02-28 | 2021-12-07 | 加利福尼亚太平洋生物科学股份有限公司 | 使用均聚物折叠测序读段改进对齐 |
| CN111584004A (zh) * | 2020-05-12 | 2020-08-25 | 西藏自治区农牧科学院水产科学研究所 | 一种基于三维组学数据的西藏特色鱼类基因组组装方法 |
| CN113151426A (zh) * | 2021-04-16 | 2021-07-23 | 中国农业科学院兰州畜牧与兽药研究所 | 一种基于三代PacBio和Hi-C技术组装和注释霍巴藏绵羊基因组的方法 |
| CN113808668A (zh) * | 2021-11-18 | 2021-12-17 | 北京诺禾致源科技股份有限公司 | 提升基因组组装完整性的方法、装置及其应用 |
| CN115148290A (zh) * | 2022-08-12 | 2022-10-04 | 武汉希望组生物科技有限公司 | 一种基于三代测序数据的补洞方法 |
Non-Patent Citations (2)
| Title |
|---|
| 单倍体基因组序列组装方法研究;官登峰;《全国优秀博硕士学位论文全文库(博士) 基础科学辑》;全文 * |
| 芥菜型油菜A9染色体物理图谱构建及重要基因研究;陆赢;《全国优秀博硕士学位论文全文库(博士) 农业科技辑》;全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115810395A (zh) | 2023-03-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN115810395B (zh) | 一种基于高通量测序动植物基因组t2t组装方法 | |
| Qiu et al. | The genome of Populus alba x Populus tremula var. glandulosa clone 84K | |
| Li et al. | Genome sequence of cultivated Upland cotton (Gossypium hirsutum TM-1) provides insights into genome evolution | |
| Shi et al. | The complete reference genome for grapevine (Vitis vinifera L.) genetics and breeding | |
| JP5960917B1 (ja) | 水稲全ゲノム育種チップ及びその応用 | |
| Zheng et al. | QTL mapping combined with bulked segregant analysis identify SNP markers linked to leaf shape traits in Pisum sativum using SLAF sequencing | |
| CN106939342B (zh) | 一种与谷子米色连锁的snp标记、引物及应用 | |
| CN108642201B (zh) | 与谷子株高性状相关的snp标记及其检测引物和应用 | |
| Wang et al. | Variation burst during dedifferentiation and increased CHH-type DNA methylation after 30 years of in vitro culture of sweet orange | |
| KR101539737B1 (ko) | 유전체 정보와 분자마커를 이용한 여교잡 선발의 효율성 증진 기술 | |
| Ercolano et al. | Complex migration history is revealed by genetic diversity of tomato samples collected in Italy during the eighteenth and nineteenth centuries | |
| CN111088389B (zh) | 玉米叶宽紧密连锁的ssr分子标记及其扩增引物和应用 | |
| Zhong et al. | Genome-wide identification of sequence variations and SSR marker development in the Munake grape cultivar | |
| CN109055593B (zh) | 提高棉花衣分的snp标记以及高产棉的鉴定和育种方法 | |
| KR20200055666A (ko) | 배추, 무 및 이들의 속간교배체의 유전형을 구분할 수 있는 dna 마커 및 이의 용도 | |
| CN108707685B (zh) | 与谷子分蘖数性状相关的snp标记及其检测引物和应用 | |
| CN114566214B (zh) | 检测基因组缺失***变异的方法及检测装置和计算机可读存储介质与应用 | |
| CN116004898A (zh) | 一种花生40K液相SNP芯片PeanutGBTS40K及其应用 | |
| CN108642199B (zh) | 与谷子旗叶长性状相关的snp标记及其检测引物和应用 | |
| CN108715901B (zh) | 一种与谷子株高性状相关的snp标记及其检测引物和应用 | |
| CN113718342A (zh) | 一种重组自交系群体高密度遗传图谱的构建方法 | |
| CN108642203B (zh) | 与谷子茎粗性状相关的snp标记及其检测引物和应用 | |
| CN108707684B (zh) | 一种与谷子旗叶长性状相关的snp标记及其检测引物和应用 | |
| KR101911307B1 (ko) | 유전자 단위에서 단상형을 구분하는 태그-snp 선발 및 활용 기술 | |
| CN109988769B (zh) | 基因及其应用 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |