CN113035269B - Genome metabolism model construction, optimization and visualization method based on high-throughput sequencing technology - Google Patents

Abstract

The invention discloses a method for constructing, optimizing and visualizing a genome metabolism model based on a high-throughput sequencing technology, which combines the I l umina second-generation sequencing technology and the PacB I o third-generation sequencing technology in the high-throughput sequencing field, carries out genome splicing and annotation with complementary advantages, compares a biological metabolism pathway database KEGG and MetaCyc through annotated protein information on the basis, and further obtains a whole genome scale metabolism network model according to the correlation of gene-enzyme-reaction according to genes corresponding to the proteins; and further searching a metabolic network model of the existing genetic species in the KEGG database, constructing relatively accurate genome-scale metabolic network models (GEMs) by an integration strategy of gene homology comparison genetic models, realizing flow automation and generating related result report files. The construction of the genome-scale metabolic network models GEMs based on genome sequence annotation and metabolic biochemical information integration is helpful for accelerating the reconstruction of microorganisms.

Description

Genome metabolic model construction, optimization and visualization based on high-throughput sequencing technology Methods

technical field

The invention belongs to the fields of metabolic genomics and computational simulation biology in bioinformatics, and relates to a method for constructing, optimizing and visualizing a genome metabolic model, in particular to a method for constructing, optimizing and visualizing a genome metabolic model based on high-throughput sequencing technology Methods.

Background technique

In order to promote the development of industrial biotechnology, it is urgent to further improve the metabolic ability of designing and modifying microbial cells. With the rapid development of whole-genome high-throughput sequencing and the continuous emergence of bioinformatics analysis strategies, it is possible to design and regulate microbial metabolic functions in a global and systematic manner. Among them, the construction of genome-scale metabolic network model GEMs based on the integration of genome sequence annotation and metabolic and biochemical information will help to better understand and accelerate the transformation of microorganisms. Therefore, how to use high-throughput sequencing technology and large-scale computing simulation to quickly construct a high-precision GEM and effectively visualize and analyze it, so that the biological simulation transformation process can be quickly grasped has become an urgent problem in this field.

Contents of the invention

The purpose of the present invention is to integrate the respective technical advantages of the second- and third-generation sequencing, improve the quality of genome assembly, and the accuracy of genome annotation, and calculate the gene model expression to obtain the comprehensive scoring value of the gene prediction results, and determine the highest score. The gene prediction model is the preliminary gene function annotation result; combined with the transcriptome data of similar species, the preliminary genome splicing annotation result is matched, and the genome sequence annotation result is optimized; the genome-scale metabolic network model GEMs is quickly constructed to visualize the metabolic network , which is convenient for users to quickly grasp the simulation transformation process of microorganisms and scientific research, analysis and search.

The present invention is a method for constructing, optimizing and visualizing a genome metabolic model based on high-throughput sequencing technology, comprising the following steps:

Step 1. Use the second-generation Illumina and third-generation PacBio high-throughput sequencing technologies to sequence the species, and obtain the original sequence files with quality control;

Step 2. Perform quality statistics and evaluation on the original sequence files respectively. If the evaluation is qualified, execute step 3; otherwise, return to execute step 1;

Step 3. Remove primers, linker sequences, and read sequence reads whose quality/length is lower than the set requirement for the original sequence file respectively, and obtain the sequence data after screening;

Step 4. Splicing the second-generation Illumina high-throughput screening data;

Splice the screened second-generation Illumina sequencing data with different k-mer size parameters, and select the splicing result with the largest parameter N50 as the contig file of the second-generation Illumina sequencing data;

Step 5. Perform the first splicing of the third-generation PacBio high-throughput screening data;

Use the splicing software to splice the screened third-generation PacBio sequencing data for the first time, and obtain the contig file for the first splicing. The splicing steps include:

Step 51, correction: stacking the screened reads together for base correction to improve the accuracy of the bases in the reads;

Step 52, pruning: determine the high-quality region and the low-quality region of the reads, and trim the low-quality region of the reads, retain a single highest-quality sequence block, and delete the low-quality sequence block;

Step 53, assembly: sort the reads sequence into a contigs file, generate corresponding consensus sequences consensus sequences, and create a path connecting the consensus sequences;

Step 6. Fusion and embed the contig file of the second-generation Illumina into the third-generation PacBio high-throughput screening data for the second splicing to obtain high-quality genome splicing results. The specific steps are as follows:

Step 61, performing sequence error correction on the first splicing data of PacBio, and obtaining the corrected reads sequence based on PacBio;

Step 62, using the contig file spliced by the second-generation Illumina as a supplement to the reads sequence corrected by PacBio, filling the gap in the third-generation PacBio sequencing, performing second splicing of the genome, and obtaining the contig file of PacBio;

Step 63: For the contig file of PacBio, use Illumina original sequencing data to further compare, filter the reads of high-quality comparison, use the bam file generated by the align part through the base error correction software Pilon, and make the second splicing result base errors were corrected;

Step 64, using the extension software to further extend the corrected splicing result of Illumina original reads data PacBio, and finally obtain a high-quality genome splicing result;

Step 7. Predict and match the genome assembly results, and complete the genome sequence annotation;

According to the splicing results, use the gene prediction software Augustus, Maker and Braker respectively to obtain all the gene prediction results, compare the predicted results with the NCBI nt database, BUSCO database and Pfam database to obtain the corresponding scores, and pass the gene model The expression is calculated to obtain the comprehensive scoring value of the gene prediction model, and the gene prediction model with the highest score is determined to be the preliminary genome assembly annotation result. Combined with the transcriptome data of similar species, the preliminary genome assembly annotation result is matched. If the sequencing read sequence If the matching degree between reads and the genome is >80%, the preliminary genome assembly annotation result is determined as the final result, and the genome sequence annotation is completed; otherwise, change the gene prediction software and re-predict and select;

Step 8, constructing the genome metabolic network model GEMs;

Step 9, build a visual network platform of metabolic network model GEMs through MetExploreViz, and perform visual display of metabolic model GEMs;

Step 10, using optFlux software to perform interface visualization FBA flow balance analysis on the metabolic network model GEMs.

Further, in the step 3, the quality lower than the set requirement refers to the bases whose quality is lower than 20, and the length lower than the set requirement refers to the sequence whose length is less than 100 bp.

Further, the different k-mer size parameters in step 4 are 21, 33, 55, 77, and 99, respectively.

Further, the gene model expression in the step 7 is:

Evidence score(gene model)＝BLASTp_score*cov(query)*cov(target)+BUSCO_score+Pfam_scores+BLASTn_score*cov(query)*cov(target)

Among them, BLASTp_score is the score derived from the BLASTp protein sequence alignment program;

BLASTp_score is the score derived from the BLASTn sequence base alignment program; BUSCO_score is the score derived from the comparison of the Benchmarking Universal Single-Copy Orthologs database of BUSCO; Pfam_scores is the score derived from the comparison of the Pfam database; cov (query) is the full-length alignment coverage of the query sequence; cov(target) is the full-length alignment coverage of the target sequence; Evidence score is the comprehensive scoring value of the final gene model.

Further, the specific steps in the step 8 include:

Step 81. Use RAVEN to compare the protein sequences of the genes annotated by the genome to the biological pathway databases KEGG and MetaCyc, and obtain the reactions Rxns, metabolite Mets and gene numbers involved in the corresponding metabolic network model, and convert the same reactions, metabolic Combining substances and related genes, performing deduplication fusion, and finally obtaining metabolic network model GEMs based on KEGG and MetaCyc;

Step 82, search for the metabolic network model GEMs of the species closest to the species in the KEGG database, and obtain the involved reaction Rxns, metabolite Mets and genes through the gene homology comparison strategy, and supplement the obtained based on KEGG and MetaCyc metabolic network model GEMs;

Step 83, through GapFind and GapFill, improve the entire metabolic network model GEMs so that the model can obtain more genes and reactions.

Further, in the step 2, the second-generation Illumina high-throughput sequencing technology evaluates the sequence quality of the species as Clean Data Q20>85%; the third-generation PacBio high-throughput sequencing technology evaluates the sequence quality of the species as Mean Concordance >85%;

Further, the high-quality region in the step 5 refers to the region where the base quality of the reads is greater than or equal to 20, and the low-quality region refers to the region where the base quality of the reads is less than 20.

Further, the species is yeast.

Further, the yeast is Yarrowia lipolytica.

Effect of the present invention is as follows:

The present invention learns from each other and effectively integrates Illumina second-generation sequencing and PacBio third-generation sequencing data to obtain high-quality genome chromosome-level assembly results; integrates the results of multiple eukaryotic gene prediction software and combines the real transcriptome data of species to obtain the most reliable and effective genome Sequence annotation results; quickly construct genome-scale metabolic network model GEMs based on the previous genome sequence annotation results; visualize the metabolic network, which is convenient for researchers or users to analyze and search; perform rapid and effective flow balance analysis on existing GEMs models, and complete microbial It can significantly improve the efficiency of microbial transformation engineering and reduce the cost of scientific research. It can be widely used in various microbial research work in high-throughput sequencing, especially suitable for the analysis and improvement of various engineering bacteria in the field of actual microbial applications. research.

Description of drawings

Fig. 1 is the genome splicing and assembling flowchart of the method for constructing, optimizing and visualizing the genome metabolic model based on high-throughput sequencing technology of the present invention;

Fig. 2 is the flowchart of genome sequence annotation of the present invention;

Figure 3 is a schematic diagram of the genome annotation results of the present invention;

Fig. 4 is the construction flow chart of metabolic network model of the present invention;

Fig. 5 is the distribution situation map on chromosome of 459 unique genes of the present invention;

Fig. 6 is the evaluation figure of the present invention respectively to three assembly results;

Figure 7 is a schematic illustration of the sequence splicing process of the present invention;

Figure 8 is a sequence comparison diagram of the reference chromosomal genome and the inferred corresponding chromosomal genome;

Fig. 9 is a highly matched scaffold contig-Scaffold86_1 result diagram of the present invention;

Fig. 10 is an explanatory diagram for storing data results in the present invention.

Detailed ways

Hereinafter, embodiments of the present invention will be described with reference to FIGS. 1-10 .

In order to specifically illustrate this technical solution, the Chinese-English comparison table involved in the solution is provided as follows:

The invention discloses a genome splicing and annotation method based on the combination of the second-generation sequencing represented by Illumina and the third-generation sequencing technology of PacBio, which are widely used in the field of high-throughput sequencing, with complementary advantages. Genome-scale metabolic network models of GEMs. The process of this method first generates a custom parameter file by the system, and generates a batch executable file with the high-throughput data processing process module after modifying the parameter file according to the user, and executes the batch executable file by the system to realize the automation of the entire process And generate relevant result report files. In this way, bioinformatics researchers, even those who do not know how to use programming scripts for data analysis, can efficiently complete high-throughput sequencing data analysis. The genome-scale metabolic network model based on genome sequence annotation and metabolic and biochemical information integration ( The construction of GEMs) will help researchers understand and accelerate the transformation of microorganisms globally. The process of the present invention can not only be used in the research of eukaryotic microorganisms such as yeast, but also can be extended and applied to other eukaryotic organisms, even bacteria and other prokaryotic organisms such as Escherichia coli for reconstruction of GEMs, and the application range is relatively wide.

In order to solve the above problems respectively, the present invention provides a method for constructing, optimizing and visualizing a genome metabolic model based on high-throughput sequencing technology, including the following steps:

1. The steps to obtain high-quality chromosome level assembly results are shown in Figure 1: According to the analyzed species, the user enters the set splicing parameter configuration file, and imports the high-throughput sequencing raw sequence data of the second-generation sequencing and the third-generation sequencing respectively. Proprietary processing software for screening and splicing to obtain effective genome splicing results, that is, the corresponding contig files. The integration results were further analyzed, and the full-length genome sequence at the chromosome level was obtained by splicing. On this basis, the subsequent bioinformatics analysis and the construction of the metabolic network model were carried out, which specifically included the following steps:

(a) Import two high-throughput sequencing original sequence files of the second-generation and third-generation sequencing of the species;

(b) Carry out quality control and statistics on the original sequence files with their respective basic processing procedures, and eliminate low-quality sequence data to obtain screened sequence data;

(c) respectively splicing and assembling the screened data into contig files of the genome;

After the Illumina original sequencing data removes the linker sequence and the low-quality (q<20) base sequence at the tail, use the SPAdes splicing software to calculate the results of different k-mer size parameters, select the effective contig sequence length greater than 500bp, and select the largest N50 The splicing results were used as the final Illumina assembly results. The PacBio original sequencing data obtained the splicing and assembly contigs results of the PacBio data through the open source software canu standard process, and then merged the two contigs files, and then further extended the software SSPACE through the base error revision software Pilon and contig, and again used the original Illumina The reads data and assembly contig results use PacBio's splicing results as the backbone to further improve the quality of the assembly results. In step (C), use the SPAdes (v3.14.1) assembly software to assemble the Illumina high-throughput sequencing data with different k-mer lengths (21, 33, 55, 77, 99), for some GC content For data sets with low coverage or uneven coverage, you can further select the single-cell mode (set --sc parameter), and finally select the splicing result with the largest N50 as the contig file of the Illumia sequencing data splicing for the next step of analysis. Use Canu (v2.1.1) splicing software to splice PacBio's high-throughput sequencing data. The splicing process is divided into three steps. The first step is to correct (Correct), which is to read "stack" (Pileup) together. Correction to improve the accuracy of bases in reads. Generally, the error of third-generation sequencing is mainly one base more or one less, because single-molecule sequencing may sometimes measure two bases when the signals are connected together or have When the same base is measured multiple times, this error is caused; the second step is trimming (Trim), that is, using overlap to determine which regions of the read are high-quality regions and which regions are of low quality and need to be trimmed. Finally, a single sequence block with the highest quality is retained, and suspicious regions are deleted; the third step is assemble, which is to sort reads into contigs, generate corresponding consensus sequences, and create paths where possible consensus sequences are connected to each other.

(d) Then integrate the two sequencing results, take the longer third-generation sequencing results as the skeleton, use the second-generation sequencing reads with lower sequencing error rates to correct and extend, and finally get the best genome assembly results; step (d) ), the preliminary contig files are first obtained through their respective standard splicing processes, and then the base error correction software Pilon ( v1.23) for correction and sequence extension with the assistance of the extension software SSPACE (v2.1.1), the final version of genome assembly was obtained, and the next step of genome annotation and metabolic model construction was carried out.

(1) The steps to obtain reliable and effective genome sequence annotation are shown in Figure 2:

According to the splicing results of the previous step, the user selects the species transcriptome reference data according to the analyzed species, verifies all the predicted results according to the results of multiple gene prediction software in eukaryotes, and removes low scores based on the comprehensive scoring results Gene prediction model, and finally determine the most reliable gene prediction results to complete the annotation;

(a) respectively importing the original sequence file of the high-throughput sequencing of the transcriptome of the reference species and the final genome assembly result file in step (1);

(b) Use RepeatModeler to model and annotate the repeat sequence family of the genome from scratch, and use RepeatMasker to detect and mask the repeat sequence in the genome, because there are repeat sequence structural regions in the genome, especially in higher eukaryotes, the repeat sequence accounts for A considerable proportion will affect the quality of gene prediction and cause unnecessary resource consumption.

Then three common gene prediction software Augustus, Maker and Braker are used to complete the gene annotation results after genome splicing. For the reference transcriptome data, through the Hisat/Trinity process, the RNA-seq assembly results are obtained, and Maker and Braker are used to complete the gene prediction of the results. All prediction results will be evaluated and scored, and the corresponding scores are obtained through BLASTp, BUSCO, InterProScan and BLASTn. Finally, the Evidence Score expression of the gene model is as follows, and the prediction result with the highest Evidence Score is obtained as the final result Gene prediction models.

Evidence score(gene model)＝BLASTp_score*cov(query)*cov(target)+BUSCO_score+Pfam_scores+BLASTn_score*cov(query)*cov(target)

Among them, BLASTp_score is the score derived from the BLASTp protein sequence alignment program; BLASTp_score is the score derived from the BLASTn sequence base alignment program; BUSCO_score is the score derived from the Benchmarking Universal Single-Copy Orthologs (BUSCO) database comparison ; Pfam_scores is the score derived from the comparison in the Pfam database; cov (query) is the full-length alignment coverage of the query sequence; cov (target) is the full-length alignment coverage of the target sequence; Evidence score is the last gene The comprehensive scoring value of the model.

(2) The steps for constructing metabolic network model GEMs are shown in Figure 4: According to the parameter configuration file and genome sequence annotation results, the open source RAVEN tool set is selected to generate the corresponding automated GEMs construction process, which specifically includes the following steps:

(a) Using the RAVEN2.0 tool set platform, firstly, the protein sequences annotated by the genome are compared with the two main biological pathway databases KEGG and MetaCyc, and the corresponding metabolic network models involved in the reaction Rxns and metabolites are respectively obtained Mets and the number of genes involved. Merge the same reactions, metabolites and genes involved, and finally obtain a metabolic network model GEMs based on the above two mainstream platforms.

(b) Use the existing metabolic network model GEMs in the KEGG database to find the closest species, and obtain the involved reaction Rxns, metabolite Mets and genes through the homology comparison strategy of the genes, and integrate them into the front as a further supplement The obtained metabolic network model GEMs.

(c) Finally, for the integrated metabolic network model GEMs, through GapFind and GapFill in the platform process, improve the entire metabolic network model GEMs,

(3) Execution and output steps: Execute the described automated GEMs construction process, obtain the genome-scale metabolic network model GEMs report based on the sequence annotation results and complete the visual display, specifically including the following steps:

Use but not limited to ModelExplore2.0 to complete the GEM network visualization evaluation analysis, and determine the quality of the metabolic network model by evaluating the number of genes involved, the connectivity of the network, and the accuracy of the predicted biological function after gene knockout. And through MetExploreViz to further build a visualization network platform of metabolic model GEMs, it is convenient for users to view certain specific metabolic pathways and reactions.

(4) The flow balance analysis of the interface visualization of the metabolic network model GEMs is mainly but not limited to using the open source optFlux software for related FBA flow balance analysis.

Utilizing the present invention, in the splicing module process, the second-generation Illumina-based high-throughput sequencing data and the third-generation PacBio high-throughput sequencing data are effectively combined, and the advantages are complementary, so that the genome splicing of species can approach or reach the level of chromosome assembly. In the process of the genome annotation module, the real transcriptome data of similar species are combined, and multiple gene prediction software are used to obtain all possible prediction results, and the predicted results are compared with NCBI's nt database, BUSCO database and Pfam database The corresponding score was obtained, calculated by the Evidence score formula, and the comprehensive score value of the gene prediction model was obtained, and finally the gene prediction model with the highest comprehensive score value was determined as the most reliable genome splicing annotation result. Finally, in the process of constructing the genome metabolic network model, using the annotated protein sequence of the genome, comparing the metabolic network involved in the homologous genes in the KEGG and MetaCyc databases, and extracting the metabolic network involved in the homologous proteins in GEMs of similar species , and finally obtain a relatively complete genome-scale metabolic network model GEMs. The entire module process can help researchers and testing personnel quickly complete a set of analysis from genome sequencing to reconstruction of genome-scale metabolic network model GEMs, and simulate species metabolic pathways to provide The optimal transformation of species provides direction. The analysis process of the present invention has a clear idea, focuses on integrating the respective technical advantages of the second and third generation sequencing, improves the quality of genome assembly, and the accuracy of genome annotation, and finally significantly improves the efficiency of microbial transformation projects and reduces scientific research costs. It can be widely used in Qualcomm Various microbial research work in quantitative sequencing, especially suitable for the analysis and improvement research of various engineering bacteria in the field of practical microbial applications.

In the method of the present invention, firstly, the system flow provides a custom parameter configuration file, and the user-defined parameter file and the high-throughput data processing flow module generate a batch executable file corresponding to the data flow according to the user-defined parameter file after modifying the set parameters. The process system executes batch executable files, realizes data process analysis, and generates result report files, so that it can efficiently help bioinformatics analysts complete a set of standardized high-throughput data integration analysis processes, and even help people who do not understand high-throughput Data analysis and metabolic network modeling researchers complete the sequence data splicing and assembly, annotation analysis, and the establishment of genome-scale metabolic network models by themselves. The splicing technique and annotation method of the present invention can be further extended and applied to more subsequent fields of high-throughput sequencing.

Below in conjunction with specific implementation example, above-mentioned scheme is described further, and this implementation example will be used for illustrating the present invention rather than limiting the scope of application of the present invention, and the relevant conditions and parameter setting that adopt in the implementation example can be done according to the condition of concrete application requirement For further adjustments, the unspecified implementation conditions are usually the conditions in routine experiments.

First, filter the raw data, that is, remove adapter sequences and low-quality reads sequences, and use their respective standard procedures to complete the processing and splicing of Illumina and PacBio's previous reads.

The statistics of Illumina sequencing data quality are shown in the table below:

In the above table: Sample ID: the name of the sequencing sample; Insert Size (bp): the length of the DNA fragment inserted when building the library; Clean Reads Length (bp): the length of the reads sequencing; Raw Data (Mb): the original data sequencing base Total; Clean Q20(%): the percentage of bases whose sequencing error rate is less than 1%; Clean Q30(%): the percentage of bases whose sequencing error rate is less than 0.1%; GC(%): the bases of C+G % of all bases.

The statistics of PacBio sequencing data quality are shown in the table below:

In the above table: Sample ID: the name of the sequencing sample; Mean Concordance: the average consistency score; Number of Reads: the number of reads; Number of Bases: the amount of data; Mean Read Length: the average sequencing read length: N50 ReadLength (bp): Sequencing read length N50 length, that is, the sequencing data are arranged in descending order, and the sequencing read length at 50% of the total.

Sequencing data quality optimization: Sequencing errors such as point mutations usually occur in high-throughput sequencing, and the quality of the end of the sequence is generally low and contains adapter sequences. Therefore, in order to improve the quality of splicing and obtain more accurate biological information analysis results, Raw data needs to be optimized.

Specific steps and parameter settings: Use Trimmomatic (v0.32) to remove primers and adapter sequences from Illumina sequencing raw data, remove bases with a quality lower than 20 at both ends, and remove sequences with a length less than 100 bp; use SPAdes (v3.14.1 ) to splice the Illumina sequencing data, select the paired reads input form, default parameters, and complete the assembly and splicing of the contig. Use Canu (v2.1.1) to splice PacBio sequencing data, set genomeSize to 20m (yeast reference genome size), correctedErrorRate default setting (0.045for PacBio reads, and 0.144for Nanopore reads), will be based on sequencing depth, coverage When the sequencing coverage is lower than 30X, the corrected ErrorRate can be increased by 1%, and when the sequencing coverage is higher than 60X, the corrected ErrorRate can be reduced by 1%. Other parameter settings are subject to the default parameters, such as minOverlapLength<integer=500>, minReadLength<integer=1000>. Use Pilon (v1.23), take FASTA and BAM files as input, improve the input reference genome according to the comparison results, including finding single base differences, small indels, large indels or blocks , filling N in the reference sequence and finding local misassemblies. The fasta file of its contig comes from the assembly results of Illumina and PacBio's respective standard processes, which are merged to remove redundancy, and finally obtain the genome assembly sequence with further improved splicing quality after Pilon processing. Using SSPACE (v2.1.1), further use the paired-end information of the original reads sequence for the spliced contig to perform possible extensions, and finally obtain the spliced version of the extended genome.

QUAST (v5.0.2) and BUSCO (v4.1.4) were used to evaluate the quality of the extended genome assembly version respectively. Among them, QUAST is used to calculate the ratio of mismatch and InDel in the assembly result to evaluate the accuracy of the base. The integrity of the assembly results was evaluated by querying the comparison results of the BUSCO database. Among them, C represents the number of complete single-copy orthologous genes in the database that can be matched, F represents the number of single-copy orthologous genes that are matched but incomplete, and M represents the single-copy orthologous genes that cannot be matched. n represents the total number of single-copy orthologous genes contained in the corresponding database.

Considering that the actual spliced genome comes from yeast, Fungap (v1.1.1), a dedicated annotation process suitable for fungi, is selected here to complete it, and the transcriptome RNAseq sequence of the reference adjacent yeast genome Y. lipolyticaW29 with existing transcriptome data is selected as an auxiliary In the process, the results of the following three different gene prediction software Augustus (v3.3.3), Braker (v2.1.5), Maker (v2.31.10) are evaluated, and the best gene model prediction is determined by constructing the evidence score formula result.

Using RAVEN Toolbox (v2.0), based on the predicted protein sequence files of genome annotations, by comparing and searching the KEGG database and the MetaCyc database, the homologous genes corresponding to the proteins were obtained, and based on the gene-enzyme-reaction association, Re-construct the genome metabolic scale network model from scratch, including the respective reaction Rxns, metabolite Mets and gene Genes. By integrating the same reactions, metabolites and genes, an integrated genome metabolic scale network model is obtained.

The evaluation of GMEs was completed through ModelExplore (v2.0), and the visual network platform display of metabolic model GEMs was further built through MetExploreViz. Use OptFlux (v3.4.0) to import the metabolic model to realize the analysis of FBA flow balance and help researchers to simulate and transform the metabolism of microorganisms.

The analysis software involved in the whole process: 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).

Below take Yarrowia lipolytica as an example to illustrate the technical scheme of the present invention:

Part I: Yarrowia lipolytica genome assembly

The process of sequence assembly is shown in Figure 7, and the flowchart of sequence assembly and evaluation is shown in Figure 1.

Genome sequence assembly will use Illumina next-generation sequencing and PacBio third-generation sequencing to complete the genome assembly results of the corresponding platforms, and will use the PacBio third-generation sequencing results as the backbone to further use the high-quality reads of Illumina next-generation sequencing to perform sequence correction through pilon and SSPACE programs Wrong and derived, resulting in our first genome analysis assembly version Poly_version 1.0. Finally, the whole assembly evaluation work is completed by using N50 parameters, QUAST and BUSCO databases respectively.

(1) Illumina platform library construction and sequencing results:

Illumina Sequencing Depth Assessment

7,567,4601502/20,000,000=113X (the yeast genome is estimated at 20MB)

Illumina NovaSeq PE150 sequencing results are shown in the table below:

Sample ID: the name of the sequencing sample; Insert Size (bp): the length of the DNA fragment inserted when building the library; Clean Reads Length (bp): the length of the reads sequencing; Raw Data (Mb): the total number of bases sequenced in the raw data; Clean Q20(%): The percentage of bases whose sequencing error rate is less than 1%; Clean Q30(%): The percentage of bases whose sequencing error rate is less than 0.1%; GC(%): C+G bases account for all bases % of quantity

The results of library construction and sequencing on the PacBio platform are shown in the table below:

The PacBio platform library construction sequencing result file is decompressed and the corresponding file is data/data_jie/jiequ05.subreads.bamPacBio sequencing depth evaluation

3,000,081,562/20,000,000=150X (the yeast genome is estimated at 20MB)

(2) Genome splicing results:

The preliminary splicing results of Illumina sequencing are shown in the table below:

serial number project Scaffold Fragment Contig Fragment 1 Total number (>500bp) 117 147 2 total length (bp) 20,295,846 20,295,546 3 N50 length (bp) 372,953 316,305 4 N90 length (bp) 139,075 87,734 5 Maximum length (bp) 949,861 949,861 6 Minimum length (bp) 574 547 7 Sequence GC content (%) 49 49

Process: Splicing and assembly is performed by commonly used next-generation sequencing software such as SPAdes.

The splicing results of PacBio sequencing are shown in the table below:

Process: The reads were assembled through the SMRT Link v5.0.1 software officially provided by PacBio. After the correction was completed using the variation detection software polish, the resequencing program was further used to map the unused reads and the assembly result polished_assembly.fasta, and finally 6 were obtained. PacBio assembly results of contigs.

The combined splicing results of Illumina and PacBio data are shown in the table below:

Process: Pilon and bowtie programs revised the assembly result error program, detailed result explanation example: Contig1:1146572Contig1_pilon:1146573-1146574.AT

In the original Contig1 sequence file spliced by PacBio, after verification by the pilon and bowtie programs, it was found that two nucleotide ATs were missing at the position of 1146572 nucleotides, that is, positions 1146573-1146574.

The corrected result is further derived

(3) Evaluation of genome splicing results:

The Quast evaluation results are shown in the table below:

Results were assessed against the BUSCO database: The completeness of genome assembly and annotation was assessed by using the BUSCO (universal single-copy orthologous genes as benchmark) database. The genes comprising the BUSCO set of each major lineage were selected from ortholog groups whose genes were present as single-copy orthologs in at least 90% of the species.

Process: For our spliced genome Yarrowia lipolytica belongs to ascomycota phyla. Select fungi_odb9 and ascomycota_odb9 to evaluate the three assembly results respectively. The result is as (accompanying drawing 6). The three assembly results are not much different in terms of relying on the evaluation of the BUSCO database, and 97.5% of the BUSCO collection can be found in ascomycota_odb9. We selected the assembly result Poly_version 1.0 of the combination of second-generation and third-generation sequencing for subsequent analysis.

(4) Determination of chromosome number: According to the experimental information, the spliced yeast has a total of 6 chromosomes, and the number of contigs corresponding to our assembly and splicing results is exactly 6. It can be determined that the splicing results have reached the chromosome level. Using the online KEGG database The information of a yli strain of Yarrowia lipolytica is known, and the information is as follows. We first determined the number of our chromosomes based on the size correspondence.

serial number 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

And through the flak-1.0 software, the reference chromosomal genome and our speculated corresponding chromosomal genome were compared, and the results are shown in Figure 8: From the above comparison results, we can see that the chromosomal genome sequence we assembled and the corresponding reference chromosomal genome sequence are highly Matching, except for scaffold5 and ChrB, there is a reverse match near 150kb. We further checked the mapping situation of Illumina reads and assembled scaffold5 contig for this section, and used the visualization software IGV to check the alignment of the obtained bam file. We were able to find that most of the highly matched gray reads covered this area, here reads It will be indicated by different colors (insertion size smaller than expected; insertion size larger than expected; inversion, duplication, translocation event), so it is speculated that this segment is indeed different from the reference yli strain. Alignment mapping was also performed on the PacBio data and the scaffold5 contig, and the results visualized the same segment around 160kb, and the long PacBio reads could also be found to cross the match smoothly. If there is a problem with the assembly, PacBio will not be able to pass through smoothly.

(5) Search for linear three-dimensional chromosome sequences: Since the reference genome has yli and there is also a mitochondrial chromosome sequence, in view of the high matching of their six corresponding chromosome sequences, we have reason to believe whether the mitochondrial sequence is missing in the spliced genome. First, through reads and reference mitochondrial sequence for mapping search, the results are as follows:

We found that 7.2% of the reads could be matched. Next, we grabbed these matched reads and used the SPAdes software to perform denovo splicing again. The mitochondrial chromosome genome file mt_chr.fasta with a size of 47,912bp was obtained. In order to verify whether it is a circular structure sequence, we intercepted a DNA sequence from the beginning and the end to create a fasta file of the test sequence. Using the blat query between the test file and the scaffold contigs file spliced by Illumina, a highly matching scaffold contig-Scaffold86_1 was found, as shown in Figure 9. A good span match can be found between the Test test sequence and the reverse sequence of scaffolds86_1, thus proving In order to create a complete circular mt_chr.fasta, the sequence is highly similar to the mitochondrial chromosome sequence of the reference genome, which proves that the mitochondrial chromosome genome sequence is highly conserved and basically unchanged. The final genome assembly sequence result is (we put the mitochondrial ring The extra 8 bases in the comparison were removed)

serial number 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

The selection of process parameters in the first part of the above is explained as follows:

For Illumina next-generation sequencing data, first use Trimmomatic software to remove adapter sequences and tail low-quality (q<20) base sequences from Illumina original sequencing data, remove sequences less than 36pb in length, and use SPAdes splicing software to calculate different k-mers After the size parameter results, the effective contig sequence length is selected by default to be greater than 500bp, and the splicing result with the largest N50 is selected as the final Illumina assembly contig result. For the PacBio third-generation sequencing data, the sequence error correction is performed first through the open source canu standard process. The yeast genome data were processed, and the reference genome size was set to 20M based on the existing reference yeast genome size on the Internet. After getting the corrected PacBio data, we found that using the Circlator software to integrate the contig files spliced by Illumina as a useful supplement, based on the corrected reads sequence of PacBio to assemble, the effect will be better than the default splicing and assembly steps of canu. Especially for circular prokaryotic genomes, even the entire complete circular genome sequence can be obtained without the need to use PCR experiments to fill in the gaps on the genome. For the correction part, here we use the Illumina original sequencing data to further complete the alignment through the bowtie and samtools software for the obtained PacBio contig assembly results, filter the high-quality aligned reads, and use the base error correction software Pilon to generate the align part Based on the bam file of Illumina, some possible base errors in the splicing results were corrected, and through the further extension software SSPACE, the corrected splicing results were further extended using Illumina original reads data PacBio, and finally a high-quality genome was obtained Assembly result. In the genome assembly evaluation part, we completed the evaluation of the assembly results through the quast and BUSCO database information. In the genome annotation part, we use Fungap to call RepeatModeler to model and annotate the repetitive sequence family of the genome from scratch, and use RepeatMasker to detect and mask the repetitive sequences in the genome, because there are repetitive sequence structural regions in the genome, especially in higher eukaryotes, Repeated sequences account for a considerable proportion, which will affect the quality of gene prediction and cause unnecessary resource consumption. Then, three common gene prediction software Augustus, Maker and Braker are used to complete the gene annotation results after genome splicing. For the reference transcriptome data, through the Hisat/Trinity process, the RNA-seq splicing and assembly results are obtained, and Maker and Braker are used to complete the gene prediction of the results. All prediction results will be evaluated and scored, and the corresponding scores will be obtained through BLASTp, BUSCO, InterProScan and BLASTn. Finally, the Evidence Score formula of the gene model is as follows, and the prediction result with the highest Evidence Score is obtained as the final gene predictive model. In this step of parameter selection, we also added sister_organisms, that is, downloaded all 9 Yarrowia lipolytica genomes in the NCBI database for annotation reference to improve the accuracy of annotation. Taking our yeast genome assembly as an example, the obtained assembly, evaluation and annotation comparison results are detailed in the research report document. Here we can see that compared with Yarrowia lipolytica CLIB122, which is a complete chromosome-level Yarrowia lipolytica genome reference in the NCBI database, our final splicing and assembly results also obtained 6 chromosomes and a plasmid sequence.

Part II: Yarrowia lipolytica genome annotation

(1) Genome annotation process

The entire process uses transcriptome data as an aid. Since we do not have our own transcriptome sequence at present, we selected a set of transcriptome sequence from the highly similar yeast Yarrowia lipolytica Strain W29 as an assistant.

According to the information provided in the article, when the same species or genus level is selected, and the corresponding average DNA similarity reaches 96% or 92.7%, the predicted result is very close to the auxiliary result of the real transcriptome data.

(2) The specific process steps of genome annotation are as follows

Step1: Preprocessing of input data.

Repeat cloaking; genomic regions such as transposon repeats, which often result in erroneous alignments, interfere with gene predictions.

Genome-specific repeat library built using RepeatModeler.

Assembly of mRNA reads;

User-supplied mRNA reads were assembled by the Trinity program.

Step2: Gene prediction

Maker and default parameters used by FunGAP; four runs with iterative SNAP gene model training. Correction of fusions of adjacent transcripts in mRNA assemblies due to the high gene density of fungal genomes The default parameters used by Augustus and FunGAP; with the user-specified augustus_species parameter. Overlapping CDS predictions are allowed. Brakes and default parameters used by FunGAP; Unsupervised RNA-Seq-based genome annotation using GeneMark-ET and Augustus.

Step3: Gene model evaluation

BLASTp; BUSCO; InterProScan(Pfam domain prediction); BLASTn; Evidence score(gene model) = BLASTp_score*cov(query)*cov(target)+BUSCO_score+Pfam_scores+BLASTn_score*cov(query)*cov(target)

Step4: Gene model filtering

Provided information, in the annotation process, the reference model selected by the gene prediction program augustus is the same species yarrowia_lipolytica for prediction, and 9 yarrowia proteome data are downloaded for homology search.

Download sister proteome

(3) Genome annotation results: FunGAP report: Created on 2019-09-17, 6970 genes were predicted in this genome.

1. Gene structure See the table below:

serial number 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 is shown in the table below:

serial number 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 (Transcript length distribution) see Figure 3

4. Protein length distribution (Protein length distribution) is shown in Figure 3. Total reads: 22,385,814; mapped reads: 20,207,847; alignment rate: 90.27%. Total reads: 22,385,814; Mapped reads: 20,207,847; Alignment rate: 90.27%

According to the existing technology, when the selected transcriptome reads alignment rate is not less than 80%, the gene prediction effect is more reliable. Here, because we have selected the yeast transcriptome data at the species level, the matching degree between reads and genome is as high as 90.27%. It shows that the predicted result is correct.

Part III: Yarrowia lipolytica genome annotation and reference genome comparison

(1) Evaluate the assembly results of the annotated genome through the BUSCO database 5 Missing BUSCOs (M): EOG092D12CU, EOG092D2L06, EOG092D3PHH, EOG092D4JRM, EOG092D4SCY For the completed reference genome Yarrowia lipolytica Strain CLIB122 existing in the KEGG database, perform the same BUSCO database核查，结果如下以下9个Missing BUSCOs(M)(EOG092D1EU1、EOG092D1Q04、EOG092D1RDH、EOG092D2L06、EOG092D2RP6、EOG092D37V8、EOG092D3C2F、EOG092D4H30、EOG092D4SCY)与前面的5个Missing BUSCOs存在2个重叠BUSCOs，可见他们都缺乏以下2 BUSCOs.

1st: EOG092D2L06 "Zinc finger, RING-like" 156 155 1 329 1.3703 NAIPR011513; IPR014857

2nd: EOG092D4SCY "EF-hand domain" 164 162 2 205 1.0865NA IPR002048; IPR011992; IPR018247

(2) The reference genome selected was Yarrowia lipolytica Strain W29, and its corresponding protein sequence annotation file was downloaded for online comparison, and a total of 6511 pair genes were found.

Next, we use the annotated 6970 gene file fungap_out_prot.faa as a benchmark to perform VENN map analysis, in which the reference genome YarliW29 shares 6511 gene maps with it. In addition, we also compared the 1315 BUSCO groups in the ascomycota_odb9 database in the previous BUSCO In terms of matching, a total of 1310 BUSCO groups (including complete and fragment matches) were matched, corresponding to 1318 genes in 6970 genes. Two of the values are 0, because we only select the 6970 genes that we spliced and annotated as the common reference comparison for reference bias, so there may be omissions in the complete comparison between YarliW29 and BUSCO (1315groups). There is no obvious clustering of genes on chromosomes. Compared with the reference genome, there are a considerable number of unique gene distributions on each chromosome. Attached Figure 5 is a map of the distribution of 459 unique genes on chromosomes.

Part IV: Reconstruction of Yarrowia lipolytica genome metabolic model

(1) Genome-scale metabolic models (GEMs)

Genome-scale metabolic models (GEMs) can simulate and predict the metabolic fluxes of various system-level metabolic studies by computationally describing the gene-protein-response associations of the entire metabolic gene in an organism. The latest GEM, iMK1208, is a high-quality model comprising 1208 genes and 1643 reactions and was successfully used to predict metabolic engineering targets that increase actin production. Streptomyces coelicolor Genome Scale Models (GEMs) for Systems Biology Research, Automated draft GEM reconstruction

Refactoring based on two databases

(a) MetaCyc-based reconstruction module: MetaCyc is a curated database of experimentally elucidated metabolic pathways from all domains of life. MetaCyc contains 2722 pathways from 3009 different organisms.

(b) KEGG-based reconstruction module: (1) GEM draft based on KEGG-annotated genome information. (2) Reconstruct the GEM draft from the homology with KEGG-trained HMMs;

ScoKEGGAnnotation = getKEGGModelForOrganism('yli',",",",",0,0);

ScoKEGGAnnotation.id = 'ScoKEGGAnnotation';

(c) Merge of MetaCyc-based and KEGG-based models. KEGG relies on sequence-based annotations, while MetaCyc only collects experimentally validated pathways. .

(2) Refer to the genome information of Yarrowia lipolytica CLIB122

(3) All 8 Yarrowia lipolytica information currently available in the JGI database

(4) Iterative process of creating metabolic model semi-automatically

(5) Specific reconstruction process of Yarrowia lipolytica genome metabolic model

Based on the above process, the information of the first version of our current metabolic model is as follows

For further iterative optimization in the later stage, making up for the gap requires further review of relevant literature and later cooperation and discussion before in-depth optimization can be carried out.

Part V: For the storage instructions of data results, please refer to Attachment 10

The above presentation describes the basic principles, features and main advantages of the invention. The present invention is not limited by the above-mentioned examples, and what described in the above-mentioned examples and the description only illustrates the principle of the present invention, without departing from the spirit and scope of the present invention, the present invention also has various changes and improvements, these changes and improvements All will fall within the scope of the claimed invention.

Claims (9)

1. A method for constructing, optimizing and visualizing a genome metabolic model based on high-throughput sequencing technology, characterized in that it comprises the following steps: Step 1. Use the second-generation Illumina and third-generation PacBio high-throughput sequencing technologies to sequence the species, and obtain the original sequence files with quality control; Step 2. Perform quality statistics and evaluation on the original sequence files respectively. If the evaluation is qualified, execute step 3; otherwise, return to execute step 1; Step 3. Remove primers, linker sequences, and read sequence reads whose quality/length is lower than the set requirement for the original sequence file respectively, and obtain the sequence data after screening; Step 4. Splicing the second-generation Illumina high-throughput screening data; Splice the screened second-generation Illumina sequencing data with different k-mer size parameters, and select the splicing result with the largest parameter N50 as the contig file of the second-generation Illumina sequencing data; Step 5. Perform the first splicing of the third-generation PacBio high-throughput screening data; Use the splicing software to splice the screened third-generation PacBio sequencing data for the first time, and obtain the contig file for the first splicing. The splicing steps include: Step 51, correction: stacking the screened reads together for base correction to improve the accuracy of the bases in the reads; Step 52, pruning: determine the high-quality region and the low-quality region of the reads, and trim the low-quality region of the reads, retain a single highest-quality sequence block, and delete the low-quality sequence block; Step 53, assembly: sort the reads sequence into a contigs file, generate corresponding consensus sequences consensus sequences, and create a path connecting the consensus sequences; Step 6. Fusion and embed the contig file of the second-generation Illumina into the third-generation PacBio high-throughput screening data for the second splicing to obtain high-quality genome splicing results. The specific steps are as follows: Step 61, performing sequence error correction on the first splicing data of PacBio, and obtaining the corrected reads sequence based on PacBio; Step 62, using the contig file spliced by the second-generation Illumina as a supplement to the reads sequence corrected by PacBio, filling the gap in the third-generation PacBio sequencing, performing second splicing of the genome, and obtaining the contig file of PacBio; Step 63: For the contig file of PacBio, use Illumina original sequencing data to further compare, filter the reads of high-quality comparison, use the bam file generated by the align part through the base error correction software Pilon, and make the second splicing result base errors were corrected; Step 64, using the extension software to further extend the corrected splicing result of Illumina original reads data PacBio, and finally obtain a high-quality genome splicing result; Step 7. Predict and match the genome assembly results, and complete the genome sequence annotation; According to the splicing results, use the gene prediction software Augustus, Maker and Braker respectively to obtain all the gene prediction results, compare the predicted results with the NCBI nt database, BUSCO database and Pfam database to obtain the corresponding scores, and pass the gene model The expression is calculated to obtain the comprehensive scoring value of the gene prediction model, and the gene prediction model with the highest score is determined to be the preliminary genome assembly annotation result. Combined with the transcriptome data of similar species, the preliminary genome assembly annotation result is matched. If the sequencing read sequence If the matching degree between reads and the genome is >80%, the preliminary genome assembly annotation result is determined as the final result, and the genome sequence annotation is completed; otherwise, change the gene prediction software and re-predict and select; Step 8, constructing the genome metabolic network model GEMs; Step 9, build a visual network platform of metabolic network model GEMs through MetExploreViz, and perform visual display of metabolic model GEMs; Step 10, using optFlux software to perform interface visualization FBA flow balance analysis on the metabolic network model GEMs. 2. The method for constructing, optimizing and visualizing a genome metabolism model based on high-throughput sequencing technology according to claim 1, characterized in that, in said step 3, the quality is lower than the set requirement refers to bases whose quality is lower than 20 , the length is lower than the set requirement refers to the sequence whose length is less than 100bp. 3. The method for constructing, optimizing and visualizing a genome metabolic model based on high-throughput sequencing technology according to claim 1, wherein the different k-mer size parameters in step 4 are 21, 33, 55 respectively. , 77, 99. 4. The method for constructing, optimizing and visualizing a genome metabolic model based on high-throughput sequencing technology according to claim 1, wherein the expression of the gene model in the step 7 is: Evidence score(gene model)＝BLASTp_score*cov(query)*cov(target)+BUSCO_score+Pfam_scores+BLASTn_score*cov(query)*cov(target) Among them, BLASTp_score is the score derived from the BLASTp protein sequence alignment program; BLASTp_score is the score derived from the BLASTn sequence base alignment program; BUSCO_score is the score derived from the comparison of the Benchmarking Universal Single-Copy Orthologs database of BUSCO; Pfam_scores is the score derived from the comparison of the Pfam database; cov(query) is the full-length alignment coverage of the query sequence; cov(target) is the full-length alignment coverage of the target sequence; Evidence score is the comprehensive scoring value of the final gene model. 5. The method for constructing, optimizing and visualizing a genome metabolism model based on high-throughput sequencing technology according to claim 1, wherein the specific steps in the step 8 include: Step 81. Use RAVEN to compare the protein sequences of the genes annotated by the genome to the biological pathway databases KEGG and MetaCyc, and obtain the reactions Rxns, metabolite Mets and gene numbers involved in the corresponding metabolic network model, and convert the same reactions, metabolic Combining substances and related genes, performing deduplication fusion, and finally obtaining metabolic network model GEMs based on KEGG and MetaCyc; Step 82, search for the metabolic network model GEMs of the species closest to the species in the KEGG database, and obtain the involved reaction Rxns, metabolite Mets and genes through the gene homology comparison strategy, and supplement the obtained based on KEGG and MetaCyc metabolic network model GEMs; Step 83, through GapFind and GapFill, improve the entire metabolic network model GEMs so that the model can obtain more genes and reactions. 6. The method for constructing, optimizing and visualizing the genome metabolism model based on high-throughput sequencing technology according to claim 1, characterized in that, in the step 2, the second-generation Illumina high-throughput sequencing technology measures the sequence quality of the species The evaluation standard is Clean Data Q20>85%; the third-generation PacBio high-throughput sequencing technology determines the sequence quality of the species and the evaluation standard is Mean Concordance>85%. 7. The method for constructing, optimizing and visualizing a genome metabolism model based on high-throughput sequencing technology according to claim 1, wherein the high-quality region in the step 5 means that the base quality of the reads is greater than or equal to 20 The low-quality region refers to the region where the base quality of the reads is less than 20. 8. The method for constructing, optimizing and visualizing a genome metabolic model based on high-throughput sequencing technology according to claim 1, wherein the species is yeast. 9. The method for constructing, optimizing and visualizing a genome metabolic model based on high-throughput sequencing technology according to claim 8, wherein the yeast is Yarrowia lipolytica.

Images