CN113035269A - 基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法 - Google Patents

基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法 Download PDF

Info

Publication number
CN113035269A
CN113035269A CN202110422896.2A CN202110422896A CN113035269A CN 113035269 A CN113035269 A CN 113035269A CN 202110422896 A CN202110422896 A CN 202110422896A CN 113035269 A CN113035269 A CN 113035269A
Authority
CN
China
Prior art keywords
genome
splicing
sequence
quality
reads
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.)
Granted
Application number
CN202110422896.2A
Other languages
English (en)
Other versions
CN113035269B (zh
Inventor
汤雷翰
申志勇
周豪魁
叶羽
罗明晶
谢秀秀
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Computational Science Research Center
Shenzhen Institute of Advanced Technology of CAS
Original Assignee
Beijing Computational Science Research Center
Shenzhen Institute of Advanced Technology of CAS
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Beijing Computational Science Research Center, Shenzhen Institute of Advanced Technology of CAS filed Critical Beijing Computational Science Research Center
Priority to CN202110422896.2A priority Critical patent/CN113035269B/zh
Publication of CN113035269A publication Critical patent/CN113035269A/zh
Application granted granted Critical
Publication of CN113035269B publication Critical patent/CN113035269B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B5/00—ICT specially adapted for modelling or simulations in systems biology, e.g. gene-regulatory networks, protein interaction networks or metabolic networks
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G16B30/10—Sequence alignment; Homology search
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G16B30/20—Sequence assembly
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B45/00—ICT specially adapted for bioinformatics-related data visualisation, e.g. displaying of maps or networks
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B50/00—ICT programming tools or database systems specially adapted for bioinformatics
    • G16B50/10—Ontologies; Annotations

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Biophysics (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Molecular Biology (AREA)
  • Bioethics (AREA)
  • Databases & Information Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Genetics & Genomics (AREA)
  • Physiology (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
  • Information Retrieval, Db Structures And Fs Structures Therefor (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法，将高通量测序领域中I l l umi na二代测序和PacB i o三代测序技术相结合，优势互补地进行基因组拼接和注释，并在此基础上通过注释的蛋白质信息比对生物代谢通路数据库KEGG和MetaCyc,依据蛋白质对应的基因,进而根据基因‑酶‑反应的关联获取到全基因组规模代谢网络模型；并进一步寻找KEGG数据库中已有亲缘物种的代谢网络模型,通过基因同源比对亲缘模型的整合策略构建了相对精准的基因组规模代谢网络模型GEMs，实现流程自动化并生成相关结果报告文件。本发明基于基因组序列注释和代谢生化信息整合的基因组规模代谢网络模型GEMs的构建将有助于加快改造微生物。

Description

基于高通量测序技术的基因组代谢模型构建、优化及可视化 的方法
技术领域
本发明属于生物信息学中代谢基因组学和计算模拟生物学领域，涉及一种基因组代谢模型构建、优化及可视化的方法，特别涉及一种基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法。
背景技术
为推动工业生物技术的发展，迫切需要进一步提高设计和改造微生物细胞的代谢能力。伴随着全基因组高通量测序的快速发展和生物信息学分析策略的不断涌现，使全局性、系统化设计和调控微生物代谢功能成为可能。其中基于基因组序列注释和代谢生化信息整合的基因组规模代谢网络模型GEMs的构建将有助于更好理解和加快改造微生物。因此如何利用高通量测序技术和大型计算模拟快速构建一个高精确度的GEM并能有效可视化分析，让生物的模拟改造过程能被迅速掌握微成为该领域目前急需解决的问题。
发明内容
本发明的目的在于整合二三代测序的各自技术优势，提高基因组组装质量，以及基因组注释的准确性，并通过基因模型表达式进行计算，得到基因预测结果的综合打分数值，确定分值最高的基因预测模型为初步基因功能注释结果；结合相近物种的转录组数据,对初步基因组拼接注释结果进行匹配，完成基因组序列注释结果的优化；快速构建基因组规模的代谢网络模型GEMs，对代谢网络进行可视化，方便用户迅速掌握微生物的模拟改造过程及科研、分析查找。
本发明为一种基于高通量测序技术的基因组代谢模型的构建、优化及可视化的方法，包括以下步骤：
步骤1、分别利用二代Illumina和三代PacBio高通量测序技术对物种进行序列测定，得到带有质量控制的原始序列文件；
步骤2、分别对原始序列文件进行质量统计及评价，评价合格则执行步骤3，否则返回执行步骤1；
步骤3、分别对原始序列文件去除引物、接头序列和质量/长度低于设定要求的读序reads序列，获得经过筛选的序列数据；
步骤4、对二代Illumina高通量筛选数据进行拼接；
以不同的k-mer size参数对经过筛选的二代Illumina测序数据进行拼接，选择参数N50最大的拼接结果作为二代Illumia测序数据拼接的contig文件；
步骤5、对三代PacBio高通量筛选数据进行第一次拼接；
通过拼接软件对经过筛选的三代PacBio测序数据进行第一次拼接，获取到第一次拼接的contig文件，其拼接步骤包括：
步骤51、修正：将经过筛选的reads堆垛到一起进行碱基修正，提高reads中碱基的准确性；
步骤52、修剪：确定reads的高质量区域和低质量区域，并对reads的低质量区域进行修整，保留单个最高质量的序列块，删除低质量序列块；
步骤53、组装：将reads序列排序为contigs文件，生成对应的共有的一致序列consensus sequences，并创建共有序列互相相连的路径；
步骤6、将所述二代Illumina的contig文件融合嵌入三代PacBio高通量筛选数据进行第二次拼接，获取高质量的基因组拼接结果，具体步骤如下：
步骤61、对PacBio第一次拼接数据进行序列纠错，得到基于PacBio纠正后的reads序列；
步骤62、将二代Illumina拼接的contig文件作为PacBio纠正后的reads序列的补充，填补三代PacBio测序中gap，进行基因组第二次拼接，得到PacBio的contig文件；
步骤63、对PacBio的contig文件，利用Illumina原始测序数据进一步比对，过滤高质量比对的reads，通过碱基错误修订软件Pilon以align部分产生的bam文件为依据，对第二次拼接结果中的碱基错误进行更正；
步骤64、通过延伸软件利用Illumina原始reads数据PacBio的纠正后拼接结果进一步延伸，最终获取到高质量的基因组拼接结果；
步骤7、对基因组拼接结果进行预测、匹配，完成基因组序列注释；
依据所述拼接结果，分别使用基因预测软件Augustus,Maker和Braker，获取所有基因预测结果，将预测出的结果分别比对NCBI的nt数据库，BUSCO数据库以及Pfam数据库获取对应分值，并通过基因模型表达式进行计算，得到基因预测模型的综合打分数值，确定分值最高的基因预测模型为初步基因组拼接注释结果，结合相近物种的转录组数据,对初步基因组拼接注释结果进行匹配，如果测序读序reads与基因组genome的匹配度>80％，则确定初步基因组拼接注释结果为最终结果，完成基因组序列注释，否则改变基因预测软件，重新进行预测和选取；
步骤8、构建基因组代谢网络模型GEMs；
步骤9、通过MetExploreViz搭建代谢网络模型GEMs的可视化网络平台，进行代谢模型GEMs的可视化展示；
步骤10、利用optFlux软件对代谢网络模型GEMs进行界面可视化的FBA流平衡分析。
进一步地，所述步骤3中质量低于设定要求指质量低于20的碱基，长度低于设定要求指长度小于100bp的序列。
进一步地，所述步骤4中不同的k-mer size参数分别为21、33、55、77、99。
进一步地，所述步骤7中基因模型表达式为：
Evidence score(gene model)＝BLASTp_score*cov(query)*cov(target)+BUSCO_score+Pfam_scores+BLASTn_score*cov(query)*cov(target)
其中，BLASTp_score为来源于BLASTp蛋白序列比对程序的分值；
BLASTp_score为来源于BLASTn序列碱基比对程序的分值；BUSCO_score为来源于BUSCO的Benchmarking Universal Single-Copy Orthologs数据库比对得到的分值；Pfam_scores为来源于Pfam数据库中比对得到的分值；cov(query)为查询序列的全长比对覆盖度；cov(target)为目标序列的全长比对覆盖度；Evidence score为最后基因模型的综合打分数值。
进一步地，所述步骤8中具体步骤包括：
步骤81、利用RAVEN对基因组注释出的基因组成的蛋白序列分别比对生物通路数据库KEGG和MetaCyc,获取到对应的代谢网络模型涉及的反应Rxns，代谢物Mets以及基因数目，将相同的反应，代谢物和涉及基因合并，进行去重融合，最终获取到基于KEGG和MetaCyc的代谢网络模型GEMs；
步骤82、在KEGG数据库中寻找与所述物种最接近的物种的代谢网络模型GEMs，通过基因同源比对策略，获取到牵扯的反应Rxns和代谢物Mets及基因，补充到所述获取的基于KEGG和MetaCyc代谢网络模型GEMs中；
步骤83、通过GapFind和GapFill,完善整个代谢网络模型GEMs使模型获得更多的基因和反应。
进一步地，所述步骤2中二代Illumina高通量测序技术对物种进行序列质量测定评价标准为Clean Data Q20>85％；三代PacBio高通量测序技术对物种进行序列质量测定评价标准为Mean Concordance>85％；
进一步地，所述步骤5中的高质量区域是指reads的碱基质量大于等于20的区域，低质量区域是指reads的碱基质量小于20的区域。
进一步地，所述物种为酵母。
进一步地，所述酵母为解脂耶氏酵母。
本发明的效果如下：
本发明取长补短将Illumina二代测序和PacBio三代测序数据有效整合，获取到高质量的基因组染色体水平组装结果；整合多个真核基因预测软件结果并结合物种真实的转录组数据得到最可靠有效的基因组序列注释结果；基于前面的基因组序列注释结果快速构建基因组规模的代谢网络模型GEMs；对代谢网络进行可视化，方便科研人员或者用户分析查找；对已有GEMs模型进行快速有效的流平衡分析，完成微生物的模拟改造预测；显著提高微生物改造工程工作效率，降低科研成本，可广泛应用于高通量测序中的各种微生物研究工作，特别适用于在实际微生物应用领域中各种工程菌的分析和改良研究上。
附图说明
图1为本发明基于高通量测序技术的基因组代谢模型的构建、优化及可视化的方法的基因组拼接组装流程图；
图2为本发明基因组序列注释流程图；
图3为本发明基因组注释结果示意图；
图4为本发明代谢网络模型构建流程图；
图5为本发明459个独特基因在染色体上的分布情况图；
图6为本发明分别对三个组装结果的评估图；
图7为本发明序列拼接过程图解；
图8为参考染色体基因组与推测的对应染色体基因组的序列比对图；
图9为本发明一个高度匹配的scaffold contig-Scaffold86_1结果图；
图10为本发明数据结果存放说明图。
具体实施方式
以下，参照附图1-10对本发明的实施方式进行说明。
为具体说明本技术方案，提供方案中涉及的中英文对照表如下：
Figure BDA0003024821100000051
Figure BDA0003024821100000061
本发明公开了一种基于高通量测序领域中，广泛应用的以Illumina为代表的二代测序和PacBio的三代测序技术相结合，优势互补的基因组拼接和注释方法，并在此基础上构建精准的基因组规模代谢网络模型GEMs。本方法流程首先由系统生成自定义的参数文件，在依据用户修改设定参数文件后和高通量数据处理流程模块生成批处理可执行文件，由系统执行批处理可执行文件，实现整个流程自动化并生成相关结果报告文件。从而能够高效的让生物信息科研人员，甚至是不太懂使用编程脚本进行数据分析的科研人员，完成高通量测序数据分析，其中基于基因组序列注释和代谢生化信息整合的基因组规模代谢网络模型(GEMs)的构建将有助于科研人员全局理解和加快改造微生物。本发明流程不仅仅可以用于酵母这类真核微生物的研究，并可推广应用到其他真核生物，甚至细菌等原核生物如大肠杆菌的GEMs重构，应用范围较为广泛。
为分别解决上述问题，本发明提供了一种基于高通量测序技术的基因组代谢模型的构建、优化及可视化的方法，包括如下步骤：
1、获取高质量染色体水平组装结果步骤如附图1：用户根据所分析的物种，输入设定的拼接参数配置文件，分别导入二代测序和三代测序的高通量测序原始序列数据，经过各自专有处理流程软件，进行筛选和拼接得到有效的基因组拼接结果，即对应的contig文件。并进一步分析整合结果，拼接得到染色体水平上的基因组全长序列，在此基础上进行后续的生物信息学分析和代谢网络模型的构建，具体包括如下步骤：
(a)分别导入物种二代和三代测序的两种高通量测序原始序列文件；
(b)分别用各自基本处理流程对原始序列文件进行质量控制与统计，并剔除低质量的序列数据，获得经过筛选的序列数据；
(c)分别将经过筛选的数据进行各自拼接、组装成基因组的contig文件；
Illumina原始测序数据在去掉接头序列，以及尾部低质量(q<20)碱基序列后，通过SPAdes拼接软件，计算不同k-mer size参数结果，选择有效contig序列长度大于500bp,并选取N50最大的拼接结果作为最终Illumina组装结果。PacBio原始测序数据通过开源的软件canu标准流程，获取到PacBio数据的拼接组装contigs结果，接下来将二者contigs文件合并，再次通过碱基错误修订软件Pilon以及contig进一步延伸软件SSPACE,再次利用Illumina原始reads数据和组装contig结果，以PacBio的拼接结果为骨架，进一步提高组装结果的质量。步骤(C)中，使用SPAdes(v3.14.1)拼接软件，以不同的k-mer长度(21,33,55,77,99)对Illumina的高通量测序数据进行拼接，对于一些GC含量高，覆盖率低或参差不齐的数据集合可进一步选择单细胞single-cell模式(设定--sc参数),最后选择N50最大的拼接结果作为Illumia测序数据拼接的contig文件进行下一步分析。使用Canu(v2.1.1)拼接软件对PacBio的高通量测序数据进行拼接，其拼接过程中分为三步，第一步是修正(Correct),即将read"堆垛"(Pileup)到一起进行修正，提高reads中碱基的准确性，一般三代测序的错误以多一个碱基或者少一个碱基为主，因为单分子测序有时候可能测两个碱基的时候信号连到一起了或者有时候对同一个碱基测了多次导致这种错误；第二步是修剪(Trim),即使用overlap确定read哪些区域是高质量区域，哪些区域质量较低需要修整。最后保留单个最高质量的序列块，并删除可疑的区域；第三步是组装(assemble),即将reads排序为contigs，生成对应的一致序列(consensus sequences)并创建可能的共有序列互相相连的路径。
(d)然后整合两种测序拼接结果，以测序较长的三代结果为骨架，利用测序错误率更低的二代测序reads进行矫正和延长，最终拿到最佳的基因组组装结果；步骤(d)中，先通过各自标准拼接流程拿到初步contig文件，然后在以PacBio的长contig文件为骨架，以Illumina的contig文件和高质量的reads测序数据为补充辅助，在碱基错误纠正软件Pilon(v1.23)进行校正和序列进一步延伸软件SSPACE(v2.1.1)的延伸协助下，拿到基因组拼接组装的最终版本，进行下一步基因组注释和代谢模型构建。
(1)获取可靠有效的基因组序列注释步骤如附图2：
依据上一步拼接结果，用户根据所分析的物种，选取物种转录组参考数据，依据真核生物中多种基因预测软件结果，对所有预测出的结果进行验证，并依据综合打分结果，去除低分数的基因预测模型，最终确定最可靠的基因预测结果完成注释；
(a)分别导入参考物种的转录组高通量测序原始序列文件和步骤(1)中最终的基因组拼接结果文件；
(b)使用RepeatModeler从头对基因组的重复序列家族进行建模注释，并且使用RepeatMasker来检测和屏蔽基因组中的重复序列，因为基因组中存在重复序列结构区，特别是高等真核生物，重复序列占了相当大的比例，会影响基因预测的质量，也会带来不必要的资源消耗。
然后通过三种常见的基因预测软件Augustus,Maker和Braker来完成基因组拼接后的基因注释结果。对于参考的转录组数据，通过Hisat/Trinity流程，拿到RNA-seq的拼接组装结果，在利用Maker和Braker完成对该结果的基因预测。所有的预测结果，将全部进行评估打分，通过BLASTp,BUSCO,InterProScan和BLASTn获取到对应的分值，最后通过基因模型的Evidence Score表达式如下，获取到Evidence Score分值最高的预测结果作为最终的基因预测模型。
Evidence score(gene model)＝BLASTp_score*cov(query)*cov(target)+BUSCO_score+Pfam_scores+BLASTn_score*cov(query)*cov(target)
其中，BLASTp_score为来源于BLASTp蛋白序列比对程序的分值；BLASTp_score为来源于BLASTn序列碱基比对程序的分值；BUSCO_score为来源于Benchmarking UniversalSingle-Copy Orthologs(BUSCO)数据库比对得到的分值；Pfam_scores为来源于Pfam数据库中比对得到的分值；cov(query)为查询序列的全长比对覆盖度；cov(target)为目标序列的全长比对覆盖度；Evidence score为最后基因模型的综合打分数值。
(2)构建代谢网络模型GEMs步骤如附图4：根据参数配置文件和基因组序列注释结果，选取开源的RAVEN工具集，生成对应的自动化GEMs构建流程，具体包括如下步骤：
(a)利用RAVEN2.0工具集平台，首先对基因组注释出的蛋白序列分别通过比对KEGG和MetaCyc两大主要的生物通路数据库,分别获取到对应的代谢网络模型所涉及的反应Rxns，代谢物Mets以及牵扯到的基因数目。将相同的反应，代谢物和涉及基因合并，最终获取到一个基于以上两大主流平台的代谢网络模型GEMs。
(b)利用KEGG数据库中已经存在的代谢网络模型GEMs，寻找最接近的物种，通过基因的同源比对策略，获取到牵扯的反应Rxns和代谢物Mets及基因，作为进一步补充，整合到前面获取到的代谢网络模型GEMs中。
(c)最后对整合得到的代谢网络模型GEMs，通过平台流程中的GapFind和GapFill,完善整个代谢网络模型GEMs，
(3)执行及输出步骤：执行所描述的自动化GEMs构建流程，获得基于序列注释结果的基因组规模的代谢网络模型GEMs报告并完成可视化展示，具体包括如下步骤:
使用但不限于利用ModelExplore2.0完成GEM网络可视化评估分析，通过评估所涉及的基因数目，网络的连通性以及基因敲除后所预测的生物学功能的准确性来确定代谢网络模型的好坏。并通过MetExploreViz进一步搭建代谢模型GEMs的可视化网络平台，方便用户具体查看某些特定代谢通路和反应。
(4)代谢网络模型GEMs的界面可视化的流平衡分析主要但不限于利用开源的optFlux软件进行相关FBA流平衡分析。
利用本发明，在拼接模块流程上，将二代Illumina为主高通量测序数据和三代PacBio的高通量测序数据进行有效结合，优势互补，使得物种基因组拼接能够接近或达到染色体组装水平，在基因组注释模块流程上，结合了真实的相近物种转录组数据，并使用多个基因预测软件，获取到所有可能预测结果，将预测出的结果分别去比对NCBI的nt数据库，BUSCO数据库以及Pfam数据库获取对应分值，通过Evidence score公式计算，得到基因预测模型的综合打分数值，并最终确定综合打分值最高的基因预测模型为最可靠的基因组拼接注释结果。最后在构建基因组代谢网络模型流程上，利用基因组注释的蛋白序列，通过比对KEGG和MetaCyc数据库内同源基因所涉及的代谢网络，以及提取已有相近物种GEMs中同源蛋白所涉及的代谢网络，最终获取到一个较完整的基因组规模代谢网络模型GEMs.整个模块流程能够帮助科研人员和检测人员迅速完成一整套从基因组测序到重建基因组规模代谢网络模型GEMs的分析，并模拟物种代谢途径从而为物种的优化改造提供方向。本发明分析流程思路清晰，着重于整合二三代测序的各自技术优势，提高基因组组装质量，以及基因组注释的准确性，最终显著提高微生物改造工程工作效率，降低科研成本，可广泛应用于高通量测序中的各种微生物研究工作，特别适用于在实际微生物应用领域中各种工程菌的分析和改良研究上。
本发明的方法，首先由系统流程提供自定义参数配置文件，在根据用户修改设定参数后的自定义参数文件和高通量数据处理流程模块生成与数据流程对应的批处理可执行文件，由流程系统执行批处理可执行文件，实现数据流程分析，生成结果报告文件，从而能够高效的帮助生物信息分析人员完成一套标准化的高通量数据整合分析流程，甚至可以让不太懂高通量数据分析和代谢网络建模的科研人员自己完成从测序数据拼接组装，注释分析，到建立基因组规模代谢网络模型。本发明的拼接技巧和注释方法可以进一步推广和应用到后续更多的高通量测序领域中。
以下结合具体实施例子对上述方案做进一步说明，该实施例子将被用于说明本发明而不是限制本发明的应用范围，实施例子中采用的相关条件和参数设定可以根据具体应用要求的条件做进一步调整，未注明的实施条件通常为常规实验中的条件。
第一，对原始数据进行过滤处理，即去除接头序列和低质量的reads序列，先分别利用各自标准流程完成对Illumina和PacBio前期reads的处理和拼接。
Illumina测序数据质量的统计见下表：
Figure BDA0003024821100000111
上表中：Sample ID:测序样品的名称；Insert Size(bp):建库时候插入的DNA片段长度；Clean Reads Length(bp)：reads测序的长度；Raw Data(Mb)：原始数据测序碱基总量；Clean Q20(％):测序错误率小于1％的碱基数量百分比；Clean Q30(％):测序错误率小于0.1％的碱基数量百分比；GC(％):C+G的碱基占所有碱基的数量百分比。
PacBio测序数据质量的统计见下表：
Figure BDA0003024821100000112
上表中：Sample ID:测序样品的名称；Mean Concordance：平均一致度得分；Number of Reads：reads个数；Number of Bases：数据量；Mean Read Length：平均测序读长：N50 ReadLength(bp)：测序读长N50长度，即将测序数据从大到小依次排列，在总体50％处的测序读长。
测序数据质量优化：高通量测序中通常会出现一些点突变等测序错误，而且序列末端的质量一般比较低，而且含有接头序列，因此为了提高拼接质量，拿到更准确的生物信息分析结果，需要对原始数据进行优化处理。
具体步骤及参数设定：使用Trimmomatic(v0.32)，对Illumina测序原始数据去除引物和接头序列，去除两端质量低于20的碱基，去除长度小于100bp的序列；使用SPAdes(v3.14.1)对Illumina测序数据进行拼接，选择paired reads输入形式，默认参数,完成contig的组装拼接。使用Canu(v2.1.1)对PacBio测序数据进行拼接，设定genomeSize为20m(酵母参考基因组大小),correctedErrorRate默认设定(0.045for PacBio reads,and0.144for Nanopore reads)，将依据测序深度，覆盖度低于30X时候，可以增加1％的corrected ErrorRate，当测序覆盖度高于60X时候，减少1％的corrected ErrorRate。其他参数设定，均以默认参数为准，如minOverlapLength<integer＝500>，minReadLength<integer＝1000>。使用Pilon(v1.23)，以FASTA和BAM文件作为输入，根据比对结果对输入的参考基因组进行提高，包括查找单碱基差异，小的插入缺失(indels)，较大的插入缺失或者block，填充参考序列中的N以及找到局部的错误组装。其contig的fasta文件，来自Illumina和PacBio各自标准流程组装结果，合并去除冗余，最终获取到经过Pilon处理后，拼接质量进一步提高的基因组拼接组装序列。使用SSPACE(v2.1.1)，进一步对拼接后的contig利用原始reads序列paired-end信息，进行可能的延伸，最终获取到延伸后的基因组拼接版本。
分别使用QUAST(v5.0.2)，BUSCO(v4.1.4)对延伸后的基因组拼接版本进行质量评估。其中用QUAST计算组装结果里的mismatch和InDel比例来评估碱基的准确度。通过查询BUSCO数据库比对结果来对组装结果进行完整性评估。其中C代表能匹配上的数据库中完整的单拷贝直系同源数目，F代表匹配上但不完整，M代表没能匹配上的单拷贝直系同源基因。n代表对应数据库中包含的单拷贝直系同源基因总数。
考虑到实际拼接的基因组来自酵母，这里选择比较适合真菌的专属注释流程Fungap(v1.1.1)来完成，并选取了已有转录组数据的参考邻近酵母基因组Y.lipolyticaW29的转录组RNAseq序列作为辅助，流程中对以下三个不同基因预测软件Augustus(v3.3.3),Braker(v2.1.5),Maker(v2.31.10)的结果进行评估,通过构建evidence score分值公式，确定最佳基因模型预测结果。
使用RAVEN Toolbox(v2.0)，基于预测到的基因组注释的蛋白序列文件，通过对KEGG数据库和MetaCyc数据库中比对查找，获取到蛋白对应的同源基因，依据基因-酶-反应的关联，分别重头构建基因组代谢规模网络模型，包含各自所涉及的反应Rxns，代谢物Mets和基因Genes.通过对相同反应，代谢物和基因的整合，获取到整合的基因组代谢规模网络模型。
通过ModelExplore(v2.0)完成GMEs的评估,并通过MetExploreViz进一步搭建代谢模型GEMs的可视化网络平台展示。使用OptFlux(v3.4.0)导入代谢模型，实现FBA流平衡的分析,帮助科研人员进行微生物的代谢模拟和改造。
整套流程所涉及使用到的分析软件：Trimmomatic(v0.32),SPAdes(v3.14.1),Canu(v2.1.1),Pilon(v1.23),SSPACE(v2.1.1)，QUAST(v5.0.2)，BUSCO(v4.1.4)，Fungap(v1.1.1),Augustus(v3.3.3),Braker(v2.1.5),Maker(v2.31.10),Samtools(v1.10),Bamtools(v2.5.1),Pfam_scan(v1.6),GeneMark-ES/ET(v4.59),RepeatModeler(v2.0.1),Hisat2(v2.2.0)，RAVEN Toolbox(v2.0),ModelExplore(v2.0),MetExploreViz,OptFlux(v3.4.0)。
下面以解脂耶氏酵母为例说明本发明技术方案：
第一部分：解脂耶氏酵母基因组拼接
序列拼接过程如附图7所示，序列组装拼接、评估流程图如附图1所示。
基因组序列拼接将分别以Illumina二代测序和PacBio三代测序完成相应平台的基因组组装结果，并将以PacBio三代测序结果为骨架，进一步使用Illumina二代测序的高质量reads通过pilon和SSPACE程序进行序列纠错和衍生，从而得到我们的第一个基因组分析组装版本Poly_version 1.0。最后分别利用N50参数、QUAST和BUSCO数据库完成整个组装评估工作。
(1)Illumina平台建库测序结果：
Illumina测序深度评估
7,567,4601502/20,000,000＝113X(酵母基因组按照20MB预估)
Illumina NovaSeq PE150测序结果见下表：
Figure BDA0003024821100000141
Sample ID:测序样品的名称；Insert Size(bp):建库时候插入的DNA片段长度；Clean Reads Length(bp)：reads测序的长度；Raw Data(Mb)：原始数据测序碱基总量；Clean Q20(％):测序错误率小于1％的碱基数量百分比；Clean Q30(％):测序错误率小于0.1％的碱基数量百分比；GC(％):C+G的碱基占所有碱基的数量百分比
PacBio平台建库测序结果见下表：
Figure BDA0003024821100000142
PacBio平台建库测序结果文件解压缩后对应文件为data/data_jie/jiequ05.subreads.bamPacBio测序深度评估
3,000,081,562/20,000,000＝150X(酵母基因组按照20MB预估)
(2)基因组拼接结果：
Illumina测序初步拼接结果见下表：
序号 项目 Scaffold片段 Contig片段
1 总数目(>500bp) 117 147
2 总长度(bp) 20,295,846 20,295,546
3 N50长度(bp) 372,953 316,305
4 N90长度(bp) 139,075 87,734
5 最大长度(bp) 949,861 949,861
6 最小长度(bp) 574 547
7 序列GC含量(％) 49 49
过程：通过SPAdes等二代测序常用软件进行拼接组装。
PacBio测序拼接结果见下表：
Figure BDA0003024821100000151
过程：通过PacBio官方提供的SMRT Link v5.0.1软件对reads进行组装，利用变异检测软件polish完成矫正后，进一步使用resequencing程序将没用到的reads与组装结果polished_assembly.fasta进行mapping,最后得到6个contigs的PacBio拼装结果。
Illumina和PacBio数据结合拼接结果见下表：
Figure BDA0003024821100000152
Figure BDA0003024821100000161
过程：Pilon和bowtie程序对组装结果错误程序进行修订，详细结果解释举例：Contig1:1146572Contig1_pilon:1146573-1146574.AT
在PacBio拼接出的最初Contig1序列文件，经过pilon和bowtie程序的核查，发现在其1146572个核甘酸的位置，遗漏了两个核甘酸AT，即1146573-1146574位置。
纠正后结果进一步衍生
(3)基因组拼接结果评估：
Quast评估结果见下表：
Figure BDA0003024821100000162
依据BUSCO数据库评估结果：基因组拼接和注释的完整性评估通过使用BUSCO(通用单拷贝直系同源基因为基准)数据库。构成每个主要谱系的BUSCO集的基因是选自直系同源组，其基因在至少90％的物种中以单拷贝直系同源物形式存在。
过程：针对我们拼接的基因组Yarrowia lipolytica属于ascomycota phyla.选取fungi_odb9和ascomycota_odb9分别对三个组装结果进行评估。结果如(附图6)。三个组装结果，在依赖BUSCO数据库评估上，相差不大，在ascomycota_odb9中都能找到97.5％的BUSCO集合.我们选取二代和三代测序相结合的组装结果Poly_version 1.0进行后续分析。
(4)染色体编号的确定：由实验信息所知拼接的酵母总共有6个染色体，我们的组装拼接结果对应的contig数目也正好是6，可以确定拼接结果达到了染色体水平，利用网上KEGG数据库中的一个Yarrowia lipolytica的yli菌株信息所知，信息如下，我们首先依据大小对应关系，初步确定了我们染色体的编号。
序号 Ref Chromosome ID Size(bp) Scaffold ID Scaffold size(bp)
1 Chromosome E 4,224,103 scaffold1 4,202,664
2 Chromosome F 4,003,362 scaffold2 4,015,094
3 Chromosome D 3,633,272 scaffold3 3,641,580
4 Chromosome C 3,272,609 scaffold4 3,323,792
5 Chromosome B 3,066,374 scaffold5 3,073,838
6 Chromosome A 2,303,261 scaffold6 2,275,489
并通过flak-1.0软件对参考染色体基因组和我们推测的对应染色体基因组进行了序列比对，结果如附图8：由以上比对结果可知，我们组装的染色体基因组序列和对应的参考染色体基因组序列高度匹配，除了scaffold5和ChrB在150kb附近出现了个倒转reverse的匹配。我们针对该区段进一步查看Illumina reads和组装的scaffold5 contig的mapping情况，对所得到的bam文件利用可视化软件IGV查看alignment的情况，我们能够找到大部分高度匹配的灰色reads覆盖了该区域，这里reads会用不同的颜色表示(插入大小小于期望值；插入大小大于期望值；倒置、重复、易位事件)，由此推测该区段确实是个与参考yli菌株不同的地方。对PacBio的数据也和scaffold5 contig进行了alignmentmapping,结果可视化相同区段around 160kb,同样能找到长的PacBio reads顺利跨过匹配。如果组装有问题，PacBio全部不能顺利跨过去。
(5)线立体染色体序列的查找：由于参考基因组有yli还存在一个线粒体染色体序列，鉴于他们6条对应染色体序列的高度匹配，我们有理由相信拼接的基因组中是否漏掉了线粒体序列，首先通过reads和参考线粒体序列进行mapping查找，结果如下：
我们找到有7.2％的reads可以匹配，接下来，我们将这些匹配的reads抓取出来，通过SPAdes软件进行denovo的再次拼接。得到size大小是47,912bp的线粒体染色体基因组文件mt_chr.fasta.为了验证是否为环状结构序列，我们截取首尾一段DNA序列，创建test序列的fasta文件。利用test文件与前面Illumina拼接得到的scaffold contigs文件进行blat查询，找到了一个高度匹配的scaffold contig-Scaffold86_1，见附图9，Test测试序列与scaffolds86_1的reverse序列可以找到很好的跨度匹配,从而证明了为一个完整的环状mt_chr.fasta，该序列与参考基因组线粒体染色体序列比对，高度相似，证明线粒体染色体基因组序列高度保守，基本没有发生变化，最后的基因组组装序列结果为(我们将线粒体环状比对中多出的8个碱基去掉了)
序号 Ref Chromosome ID Size(bp) Scaffold ID Scaffold size(bp)
1 Chromosone MT 47,916 scaffold_mt 47,904
2 Chromosome E 4,224,103 scaffold1 4,202,664
3 Chromosome F 4,003,362 scaffold2 4,015,094
4 Chromosome D 3,633,272 scaffold3 3,641,580
5 Chromosome C 3,272,609 scaffold4 3,323,792
6 Chromosome B 3,066,374 scaffold5 3,073,838
7 Chromosome A 2,303,261 scaffold6 2,275,489
以上整个第一部分流程参数选择说明如下：
对Illumina二代测序数据，首先通过Trimmomatic软件对Illumina原始测序数据在去掉接头序列，以及尾部低质量(q<20)碱基序列，去除长度小于36pb序列，通过SPAdes拼接软件，计算不同k-mer size参数结果后默认选择有效contig序列长度大于500bp,并选取N50最大的拼接结果作为最终Illumina组装contig结果.对于PacBio三代测序数据，则是通过开源的canu标准流程中先进行序列纠错，对我们的酵母基因组数据进行处理，并依据网上已有参考酵母基因组大小设定参考基因组大小为20M。在拿到纠正后的PacBio数据后，我们发现利用Circlator软件，整合Illumina拼接的contig文件作为有益补充，基于PacBio纠正后的reads序列来组装，其效果会好于canu默认自带拼接组装步骤结果，特别是对于环状的原核生物基因组，甚至可以拿到整个完整的环状基因组序列，而无需要在用PCR实验去补基因组上的gap。纠正部分，这里我们对拿到的PacBio的contig组装结果，利用Illumina原始测序数据，经过bowtie和samtools软件进一步完成比对，过滤高质量比对的read，通过碱基错误修订软件Pilon以align部分产生的bam文件为依据，对拼接结果中一些可能的碱基错误进行更正，并通过进一步的延伸软件SSPACE,利用Illumina原始reads数据PacBio的纠正后拼接结果进一延伸，最终获取到了一个高质量的基因组组装结果。基因组组装评估部分，我们通过quast和BUSCO数据库信息完成了对组装结果的评估工作。基因组注释部分，我们通过fungap调用RepeatModeler从头对基因组的重复序列家族进行建模注释，并且使用RepeatMasker来检测和屏蔽基因组中的重复序列，因为基因组中存在重复序列结构区，特别是高等真核生物，重复序列占了相当大的比例，会影响基因预测的质量，也会带来不必要的资源消耗。然后，通过三种常见的基因预测软件Augustus,Maker和Braker来完成基因组拼接后的基因注释结果。对于参考的转录组数据，通过Hisat/Trinity流程，拿到RNA-seq的拼接组装结果，在利用Maker和Braker完成对该结果的基因预测。所有的预测结果，将全部进行评估打分，通过BLASTp,BUSCO,InterProScan和BLASTn获取到对应的分值，最后通过基因模型的Evidence Score公式如下，获取到Evidence Score分值最高的预测结果作为最终的基因预测模型。在这一步参数选择上，我们还加入了sister_organisms,即下载了NCBI数据库中所有9个解脂耶氏酵母基因组做注释参考，以提高注释的准确性。以我们的酵母基因组拼接为例子，获取到的拼接，评估和注释比对结果详见研究报告文档。这里我们可以看到与NCBI数据库中已有的完整染色体水平的解脂耶氏酵母基因组参考的Yarrowialipolytica CLIB122比较我们最后的拼接组装结果也是得到了6个染色体和一个质粒序列。
第二部分：解脂耶氏酵母基因组注释
(1)基因组注释过程
整套流程利用了转录组数据做辅助，针对目前我们没有自己的转录组序列，选取了高度相似的同种酵母Yarrowia lipolytica Strain W29一套转录组序列为辅助。
依据文章提供的信息，当选取相同species或genus水平，对应平均DNA相似度达到96％或92.7％，则预测结果与真实转录组数据辅助结果非常接近。
(2)基因组注释具体过程步骤如下
Step1：输入数据的预处理。
重复遮蔽；基因组区域如转座子重复，经常会使错误的排列，干扰基因预测。
使用RepeatModeler构建的基因组特异性重复库。
mRNA读数的组装；
用户提供的mRNA读数由Trinity程序组装。
Step2：基因预测
FunGAP使用的制造者和默认参数；用迭代SNAP基因模型训练运行四次。由于真菌基因组的基因密度较高，在mRNA组装中纠正相邻转录本的融合。Augustus和FunGAP使用的默认参数；与用户指定的augustus_species参数。允许重叠的CDS预测。FunGAP使用的刹车器和默认参数；使用GeneMark-ET和Augustus进行基于无监督RNA测序的基因组注释。
Step3：基因模型评估
BLASTp；BUSCO；InterProScan(Pfam域预测)；BLASTn；证据得分(基因模型)＝BLASTp_score*cov(query)*cov(target)+BUSCO_score+Pfam_scores+BLASTn_score*cov(query)*cov(target)
Step4：基因模型过滤
提供的信息，注释流程中,基因预测程序augustus选取的参考模型为相同的种yarrowia_lipolytica进行预测，并下载9个yarrowia proteome数据进行同源查找。
下载姐妹蛋白质组(Download sister proteome)
(3)基因组注释结果:FunGAP报告：创建于2019-09-17，该基因组中预测了6970个基因。
1.基因结构Gene structure见下表：
序号 Attributes Values
1 Total protein-coding genes 6,970
2 Transcript length(avg/med) 1,463.9/1,201.0
3 CDS length(avg/med) 1,362.7/1,131.0
4 Protein length(avg/med) 454.2/377.0
5 Exon length(avg/med) 1,035.0/810.0
6 Intron length(avg/med) 319.5/168.0
7 Spliced genes 1,833(26.3％)
8 Gene density(genes/Mb) 338.67
9 Number of introns 2,207
10 Number of introns per gene(med) 1.0
11 Number of exons 9,177
12 Number of exons per gene(med) 1.0
2.转录组读长序列组装(Transcriptome reads assembly)见下表：
序号 Attributes Values
1 Number of mapped reads 20,207,847
2 Number of assembled contigs 10,448
3 Number of contigs>1kbp 5,295
4 Total transcript size(Mbp) 13,320,738
3.转录本长度分布(Transcript length distribution)见附图3
4.蛋白质长度分布(Protein length distribution)见附图3总阅读量：22,385,814次；映射读数：20,207,847次；对齐率：90.27％。Total reads:22,385,814；Mappedreads:20,207,847；Alignment rate:90.27％
依据现有技术可知，当选取的transcriptome reads alignment rate不低于80％时候，基因预测效果就比较可靠，这里由于我们选择了species水平的酵母转录组数据，reads与genome的匹配度高达90.27％.说明预测的结果是对的。
第三部分：解脂耶氏酵母基因组注释与参考基因组比较
(1)将注释完成的基因组通过BUSCO数据库评估组装结果5个Missing BUSCOs(M):EOG092D12CU、EOG092D2L06、EOG092D3PHH、EOG092D4JRM、EOG092D4SCY对于KEGG数据库中存在的完成参考基因组Yarrowia lipolytica Strain CLIB122的进行同样的BUSCO数据库核查，结果如下以下9个Missing BUSCOs(M)(EOG092D1EU1、EOG092D1Q04、EOG092D1RDH、EOG092D2L06、EOG092D2RP6、EOG092D37V8、EOG092D3C2F、EOG092D4H30、EOG092D4SCY)与前面的5个Missing BUSCOs存在2个重叠BUSCOs，可见他们都缺乏以下2个BUSCOs.
第1个：EOG092D2L06"Zinc finger,RING-like"156 155 1 329 1.3703 NAIPR011513；IPR014857
第2个：EOG092D4SCY"EF-hand domain"164 162 2 205 1.0865NA IPR002048；IPR011992；IPR018247
(2)选取的参考基因组为Yarrowia lipolytica Strain W29，下载其对应的蛋白序列注释文件进行在线比对，找到的共有的6511pair genes。
接下来，我们以注释的6970基因文件fungap_out_prot.faa为基准，进行VENN图分析，其中参考基因组YarliW29与其共有的基因map数目为6511.另外我们也比较了前面BUSCO中ascomycota_odb9数据库中的1315个BUSCO group匹配情况，总共是匹配到了1310个BUSCO group(包含完整和片段匹配),对应6970个基因中的1318个基因。其中2个为0的数值，是由于reference bias我们只选择我们拼接注释的6970个基因为共同参考比较，因此YarliW29与BUSCO(1315groups)的完全比较存在遗漏可能。并没有见到比较明显的基因在染色体上的聚类情况，每条染色体上与参考基因组的比较，存在相当数目的独特基因分布，附图5为459个独特基因在染色体上的分布情况图。
第四部分：解脂耶氏酵母基因组代谢模型重构
(1)基因组规模的代谢模型(GEMs)
基因组尺度代谢模型(GEMs)通过计算描述生物体内整个代谢基因的基因-蛋白-反应关联，可以模拟预测各种系统级代谢研究的代谢通量。最新的GEM，iMK1208，是一个高质量的模型，包括1208个基因和1643个反应，并成功地用于预测增加肌动蛋白生产的代谢工程目标。天蓝链霉菌(Streptomyces coelicolor)基因组尺度模型(GEMs)用于系统生物学研究，创业板自动重建草案Automated draft GEM reconstruction
基于两个数据库的重构
(a)基于MetaCyc的重建模块：MetaCyc是一个精心策划的数据库，收录了来自所有生命领域的实验阐明的代谢途径。MetaCyc包含来自3009种不同生物体的2722条途径。
(b)基于KEGG的重建模块：(1)基于KEGG注释的基因组信息的GEM草案。(2)用KEGG训练的HMMs从同源性中重建GEM草案；
ScoKEGGAnnotation＝getKEGGModelForOrganism('yli',”,”,”,0,0)；
ScoKEGGAnnotation.id＝'ScoKEGGAnnotation'；
(c)基于MetaCyc的模型和基于KEGG的模型的合并。KEGG依赖于基于序列的注释，而MetaCyc只收集经过实验验证的途径。.
(2)参考Yarrowia lipolytica CLIB122基因组信息
(3)JGI数据库目前已有的全部8个Yarrowia lipolytica信息
Figure BDA0003024821100000241
(4)半自动创建代谢模型迭代过程
(5)解脂耶氏酵母基因组代谢模型具体重构流程
基于以上流程，我们目前的代谢模型第一个版本信息如下
Figure BDA0003024821100000242
对于后期进一步的迭代优化，弥补gap需要进一步查看相关文献以及后期合作探讨才能深入优化进行。
第五部分：数据结果存放说明见附图10
以上显示描述了本发明的基本原理，特征和主要优点。本发明不受上述实例的限制，上述实例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都将落入要求保护的本发明范围内。

Claims (9)

1.一种基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法，其特征在于，其包括以下步骤：
步骤1、分别利用二代Illumina和三代PacBio高通量测序技术对物种进行序列测定，得到带有质量控制的原始序列文件；
步骤2、分别对原始序列文件进行质量统计及评价，评价合格则执行步骤3，否则返回执行步骤1；
步骤3、分别对原始序列文件去除引物、接头序列和质量/长度低于设定要求的读序reads序列，获得经过筛选的序列数据；
步骤4、对二代Illumina高通量筛选数据进行拼接；
以不同的k-mer size参数对经过筛选的二代Illumina测序数据进行拼接，选择参数N50最大的拼接结果作为二代Illumia测序数据拼接的contig文件；
步骤5、对三代PacBio高通量筛选数据进行第一次拼接；
通过拼接软件对经过筛选的三代PacBio测序数据进行第一次拼接，获取到第一次拼接的contig文件，其拼接步骤包括：
步骤51、修正：将经过筛选的reads堆垛到一起进行碱基修正，提高reads中碱基的准确性；
步骤52、修剪：确定reads的高质量区域和低质量区域，并对reads的低质量区域进行修整，保留单个最高质量的序列块，删除低质量序列块；
步骤53、组装：将reads序列排序为contigs文件，生成对应的共有的一致序列consensus sequences，并创建共有序列互相相连的路径；
步骤6、将所述二代Illumina的contig文件融合嵌入三代PacBio高通量筛选数据进行第二次拼接，获取高质量的基因组拼接结果，具体步骤如下：
步骤61、对PacBio第一次拼接数据进行序列纠错，得到基于PacBio纠正后的reads序列；
步骤62、将二代Illumina拼接的contig文件作为PacBio纠正后的reads序列的补充，填补三代PacBio测序中gap，进行基因组第二次拼接，得到PacBio的contig文件；
步骤63、对PacBio的contig文件，利用Illumina原始测序数据进一步比对，过滤高质量比对的reads，通过碱基错误修订软件Pilon以align部分产生的bam文件为依据，对第二次拼接结果中的碱基错误进行更正；
步骤64、通过延伸软件利用Illumina原始reads数据PacBio的纠正后拼接结果进一步延伸，最终获取到高质量的基因组拼接结果；
步骤7、对基因组拼接结果进行预测、匹配，完成基因组序列注释；
依据所述拼接结果，分别使用基因预测软件Augustus,Maker和Braker，获取所有基因预测结果，将预测出的结果分别比对NCBI的nt数据库，BUSCO数据库以及Pfam数据库获取对应分值，并通过基因模型表达式进行计算，得到基因预测模型的综合打分数值，确定分值最高的基因预测模型为初步基因组拼接注释结果，结合相近物种的转录组数据,对初步基因组拼接注释结果进行匹配，如果测序读序reads与基因组genome的匹配度>80％，则确定初步基因组拼接注释结果为最终结果，完成基因组序列注释，否则改变基因预测软件，重新进行预测和选取；
步骤8、构建基因组代谢网络模型GEMs；
步骤9、通过MetExploreViz搭建代谢网络模型GEMs的可视化网络平台，进行代谢模型GEMs的可视化展示；
步骤10、利用optFlux软件对代谢网络模型GEMs进行界面可视化的FBA流平衡分析。
2.根据权利要求1所述的基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法，其特征在于，所述步骤3中质量低于设定要求指质量低于20的碱基，长度低于设定要求指长度小于100bp的序列。
3.根据权利要求1所述的基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法，其特征在于，所述步骤4中不同的k-mer size参数分别为21、33、55、77、99。
4.根据权利要求1所述的基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法，其特征在于，所述步骤7中基因模型表达式为：
Evidence score(gene model)＝BLASTp_score*cov(query)*cov(target)+BUSCO_score+Pfam_scores+BLASTn_score*cov(query)*cov(target)
其中，BLASTp_score为来源于BLASTp蛋白序列比对程序的分值；
BLASTp_score为来源于BLASTn序列碱基比对程序的分值；BUSCO_score为来源于BUSCO的Benchmarking Universal Single-Copy Orthologs数据库比对得到的分值；Pfam_scores为来源于Pfam数据库中比对得到的分值；
cov(query)为查询序列的全长比对覆盖度；cov(target)为目标序列的全长比对覆盖度；Evidence score为最后基因模型的综合打分数值。
5.根据权利要求1所述的基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法，其特征在于，所述步骤8中具体步骤包括：
步骤81、利用RAVEN对基因组注释出的基因组成的蛋白序列分别比对生物通路数据库KEGG和MetaCyc,获取到对应的代谢网络模型涉及的反应Rxns，代谢物Mets以及基因数目，将相同的反应，代谢物和涉及基因合并，进行去重融合，最终获取到基于KEGG和MetaCyc的代谢网络模型GEMs；
步骤82、在KEGG数据库中寻找与所述物种最接近的物种的代谢网络模型GEMs，通过基因同源比对策略，获取到牵扯的反应Rxns和代谢物Mets及基因，补充到所述获取的基于KEGG和MetaCyc代谢网络模型GEMs中；
步骤83、通过GapFind和GapFill,完善整个代谢网络模型GEMs使模型获得更多的基因和反应。
6.根据权利要求1所述的基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法，其特征在于，所述步骤2中二代Illumina高通量测序技术对物种进行序列质量测定评价标准为Clean Data Q20>85％；三代PacBio高通量测序技术对物种进行序列质量测定评价标准为Mean Concordance>85％。
7.根据权利要求1所述的基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法，其特征在于，所述步骤5中的高质量区域是指reads的碱基质量大于等于20的区域，低质量区域是指reads的碱基质量小于20的区域。
8.根据权利要求1所述的基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法，其特征在于，所述物种为酵母。
9.根据权利要求8所述的基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法，其特征在于，所述酵母为解脂耶氏酵母。
CN202110422896.2A 2021-04-16 2021-04-16 基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法 Active CN113035269B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110422896.2A CN113035269B (zh) 2021-04-16 2021-04-16 基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110422896.2A CN113035269B (zh) 2021-04-16 2021-04-16 基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法

Publications (2)

Publication Number Publication Date
CN113035269A true CN113035269A (zh) 2021-06-25
CN113035269B CN113035269B (zh) 2022-11-01

Family

ID=76456928

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110422896.2A Active CN113035269B (zh) 2021-04-16 2021-04-16 基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法

Country Status (1)

Country Link
CN (1) CN113035269B (zh)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111161797A (zh) * 2019-12-31 2020-05-15 北京百迈客生物科技有限公司 一种基于三代测序检测多样本量比较转录组分析方法
CN113903411A (zh) * 2021-08-11 2022-01-07 东北林业大学 一种基于后缀数组与单调栈的基因组组装预处理方法
CN114168708A (zh) * 2021-11-15 2022-03-11 哈尔滨工业大学 一种基于多域特征的个性化生物通路检索方法
CN114333987A (zh) * 2021-12-30 2022-04-12 天津金匙医学科技有限公司 一种基于宏基因组测序的预测耐药表型的数据分析方法
CN114822700A (zh) * 2022-04-25 2022-07-29 至本医疗科技（上海）有限公司 用于呈现重排或融合结构亚型的方法、设备和介质
CN114875118A (zh) * 2022-06-30 2022-08-09 北京百图智检科技服务有限公司 确定细胞谱系的方法、试剂盒和装置
CN115691673A (zh) * 2022-10-25 2023-02-03 广东省农业科学院蔬菜研究所 一种端粒到端粒的基因组组装方法
CN116072222A (zh) * 2023-02-16 2023-05-05 湖南大学 病毒基因组鉴定和拼接的方法及应用
CN117275590A (zh) * 2023-11-10 2023-12-22 华东师范大学 有机固废体系中大分子的降解功能基因数据库与分析平台

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106778076A (zh) * 2016-11-15 2017-05-31 上海派森诺生物科技股份有限公司 一种高效的针对于放线菌基因组拼接的方法
CN107609347A (zh) * 2017-08-21 2018-01-19 上海派森诺生物科技股份有限公司 一种基于高通量测序技术的宏转录组数据分析方法
US20180260521A1 (en) * 2017-03-09 2018-09-13 Shine Biopharma Inc. Method and apparatus for multiple dot plot analysis
CN108804875A (zh) * 2018-06-21 2018-11-13 中国科学院北京基因组研究所 一种利用宏基因组数据分析微生物群体功能的方法
US20190144852A1 (en) * 2017-11-13 2019-05-16 The Board Of Trustees Of The University Of Illinois Combinatorial Metabolic Engineering Using a CRISPR System
CN111192630A (zh) * 2019-12-24 2020-05-22 中国科学院生态环境研究中心 一种宏基因组数据挖掘方法
CN112133368A (zh) * 2020-10-13 2020-12-25 南开大学 一种基于三代测序技术的宏基因组测序数据的自动化分析方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106778076A (zh) * 2016-11-15 2017-05-31 上海派森诺生物科技股份有限公司 一种高效的针对于放线菌基因组拼接的方法
US20180260521A1 (en) * 2017-03-09 2018-09-13 Shine Biopharma Inc. Method and apparatus for multiple dot plot analysis
CN107609347A (zh) * 2017-08-21 2018-01-19 上海派森诺生物科技股份有限公司 一种基于高通量测序技术的宏转录组数据分析方法
US20190144852A1 (en) * 2017-11-13 2019-05-16 The Board Of Trustees Of The University Of Illinois Combinatorial Metabolic Engineering Using a CRISPR System
CN108804875A (zh) * 2018-06-21 2018-11-13 中国科学院北京基因组研究所 一种利用宏基因组数据分析微生物群体功能的方法
CN111192630A (zh) * 2019-12-24 2020-05-22 中国科学院生态环境研究中心 一种宏基因组数据挖掘方法
CN112133368A (zh) * 2020-10-13 2020-12-25 南开大学 一种基于三代测序技术的宏基因组测序数据的自动化分析方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
赵刚等: "基于Illumina高通量测序的岩原鲤转录组分析", 《西南农业学报》 *

Cited By (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111161797B (zh) * 2019-12-31 2023-06-06 北京百迈客生物科技有限公司 一种基于三代测序检测多样本量比较转录组分析方法
CN111161797A (zh) * 2019-12-31 2020-05-15 北京百迈客生物科技有限公司 一种基于三代测序检测多样本量比较转录组分析方法
CN113903411A (zh) * 2021-08-11 2022-01-07 东北林业大学 一种基于后缀数组与单调栈的基因组组装预处理方法
CN114168708A (zh) * 2021-11-15 2022-03-11 哈尔滨工业大学 一种基于多域特征的个性化生物通路检索方法
CN114168708B (zh) * 2021-11-15 2022-06-14 哈尔滨工业大学 一种基于多域特征的个性化生物通路检索方法
CN114333987A (zh) * 2021-12-30 2022-04-12 天津金匙医学科技有限公司 一种基于宏基因组测序的预测耐药表型的数据分析方法
CN114822700A (zh) * 2022-04-25 2022-07-29 至本医疗科技（上海）有限公司 用于呈现重排或融合结构亚型的方法、设备和介质
CN114875118A (zh) * 2022-06-30 2022-08-09 北京百图智检科技服务有限公司 确定细胞谱系的方法、试剂盒和装置
CN114875118B (zh) * 2022-06-30 2022-10-11 北京百图智检科技服务有限公司 确定细胞谱系的方法、试剂盒和装置
CN115691673A (zh) * 2022-10-25 2023-02-03 广东省农业科学院蔬菜研究所 一种端粒到端粒的基因组组装方法
CN115691673B (zh) * 2022-10-25 2023-08-15 广东省农业科学院蔬菜研究所 一种端粒到端粒的基因组组装方法
CN116072222A (zh) * 2023-02-16 2023-05-05 湖南大学 病毒基因组鉴定和拼接的方法及应用
CN116072222B (zh) * 2023-02-16 2024-02-06 湖南大学 病毒基因组鉴定和拼接的方法及应用
CN117275590A (zh) * 2023-11-10 2023-12-22 华东师范大学 有机固废体系中大分子的降解功能基因数据库与分析平台
CN117275590B (zh) * 2023-11-10 2024-03-26 华东师范大学 有机固废体系中大分子的降解功能基因数据库与分析平台

Also Published As

Publication number Publication date
CN113035269B (zh) 2022-11-01

Similar Documents

Publication Publication Date Title
CN113035269B (zh) 基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法
Venturini et al. Leveraging multiple transcriptome assembly methods for improved gene structure annotation
Lischer et al. Reference-guided de novo assembly approach improves genome reconstruction for related species
Scalzitti et al. A benchmark study of ab initio gene prediction methods in diverse eukaryotic organisms
Hara et al. Optimizing and benchmarking de novo transcriptome sequencing: from library preparation to assembly evaluation
Sun et al. LTR retrotransposons contribute to genomic gigantism in plethodontid salamanders
Rudd Expressed sequence tags: alternative or complement to whole genome sequences?
Henson et al. Next-generation sequencing and large genome assemblies
Wang et al. The draft nuclear genome assembly of Eucalyptus pauciflora: a pipeline for comparing de novo assemblies
Rupp et al. Construction of a public CHO cell line transcript database using versatile bioinformatics analysis pipelines
Rödelsperger et al. CYNTENATOR: progressive gene order alignment of 17 vertebrate genomes
Vendrell-Mir et al. A benchmark of transposon insertion detection tools using real data
Singh Fundamentals of bioinformatics and computational biology
Bowman et al. A modified GC-specific MAKER gene annotation method reveals improved and novel gene predictions of high and low GC content in Oryza sativa
Lozano-Fernandez A practical guide to design and assess a phylogenomic study
EP4018450A1 (en) Methods and systems for assessing genetic variants
Beier et al. Multiplex sequencing of bacterial artificial chromosomes for assembling complex plant genomes
Sandhya et al. Methods and tools for plant organelle genome sequencing, assembly, and downstream analysis
Wang et al. BAUM: improving genome assembly by adaptive unique mapping and local overlap-layout-consensus approach
Haiminen et al. Assessing pooled BAC and whole genome shotgun strategies for assembly of complex genomes
Gao et al. Functional and structural divergence of an unusual LTR retrotransposon family in plants
Williams et al. Methods and strategies for gene structure curation in WormBase
Dubchak et al. Multi-species sequence comparison: the next frontier in genome annotation
Lu et al. Predicting functional alternative splicing by measuring RNA selection pressure from multigenome alignments
Parisi et al. The structurally constrained protein evolution model accounts for sequence patterns of the LβH superfamily

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