A macrovirus group analysis method based on second-generation and third-generation sequencing technology

Abstract

The application discloses a macrovirus group analysis method based on second-generation and third-generation sequencing technologies, which comprises the following steps: performing quality control and filtering on second-generation sequencing raw data to obtain clean reads, performing single-sample and multi-sample assembly on the quality-controlled data to obtain contigs sequences; performing correction and quality control on third-generation sequencing raw data to obtain clean long sequences, performing assembly on the quality-controlled third-generation data to obtain contigs sequences; performing mixed assembly on the quality-controlled data of the second-generation and third-generation sequencing to obtain contigs sequences; merging all contigs to construct a non-redundant contigs set; and finally performing virus identification and determination, virus species annotation and functional annotation. The application provides a reliable macrovirus group analysis method based on second-generation and third-generation sequencing technologies, and the implementation method is simple and the application range is wide.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of second-generation and third-generation sequencing technology and the field of macrovirology, in particular to a macrovirology analysis method based on second-generation and third-generation sequencing technology. BACKGROUND

[0002] Viruses are an important component of microorganisms, and their composition changes are related to many diseases, and are the most abundant species of biological circles on earth. Viruses that can invade bacteria and cause host lysis are also known as bacteriophages, which can transmit genetic material between bacteria as a mobile genetic element. In recent years, due to the widespread use of antibiotics, superbugs have emerged, and the threat of pathogenic bacteria to antibiotic resistance is growing, and the effectiveness of antibiotic treatment is severely challenged. Since bacteriophages can invade and change (usually reduce) the fitness of host bacteria, bacteriophage therapy has become a hot topic again. Therefore, viruses play an important role in regulating bacterial biomass, maintaining biodiversity, horizontal gene transfer, and biochemical cycles in the entire biosphere, and can also be used as an alternative to antibiotics, with strong practical value.

[0003] Macrovirology is a new branch of science that has emerged on the basis of the theory of metagenomics and existing molecular biology detection techniques for viruses. It is a general term for all viruses or virus-like genetic information carried by a certain sample. Macrovirology directly studies the genetic material of all viruses in the environment, and can quickly and accurately identify the composition of all viruses in the environment, playing an important role in virus discovery, virus tracing, and microbial early warning. Macrovirology research can be applied to the study of human or animal intestinal or blood samples, oceans, soil, etc., to uncover potential hazards to humans and the environment. The continuous progress of sequencing technology and the improvement of analysis methods have enabled people to glimpse the unknown scientific field. In addition, many characteristics of viruses are based on the mutual influence of the entire community environment and individuals. Therefore, macrovirology can better discover the characteristics of viruses than individual studies.

[0004] Due to the limitation of sequencing technology principle and the complexity of viral sequence, the original sequence read length of second-generation sequencing is short (usually only 50-250 bp), it is difficult to obtain low-abundance viruses, and the assembly fragment accuracy is affected by many factors; the read length of third-generation sequencing is longer (average 10-15 kb, maximum 80 kb), but because of the cost problem, it has not been widely used. How to balance the sequencing cost and the sequencing throughput, effectively process the original sequencing data and then obtain the complete and accurate viral sequence is a hot issue in the field of bioinformatics. With the development of metaviromics, analysis tools for viral sequencing data are emerging in an endless stream, but tools that can combine the advantages of second-generation and third-generation sequencing technologies to complete viral sequence identification, species analysis and functional analysis are not common. Therefore, it is an urgent need to integrate existing tools, complete the above analysis modules and make the analysis as efficient and accurate as possible. SUMMARY

[0005] The purpose of the present application is to provide a metaviromics analysis method based on second-generation and third-generation sequencing technology. The present application provides accurate analysis results for researchers by processing the data from second-generation and third-generation metaviromics sequencing to mixed assembly, virus identification and species, and functional annotation. It also allows non-bioinformatics professionals to independently complete viral high-throughput data analysis, optimizes the work efficiency of researchers, and reduces the cost of scientific research. The present application provides a reliable metaviromics analysis method based on second-generation and third-generation sequencing technology, which is simple to implement and widely applicable, and solves the technical problem of balancing sequencing cost and sequencing throughput in the prior art.

[0006] According to the purpose of the present application, a metaviromics analysis method based on second-generation and third-generation sequencing technology is provided, comprising the following steps:

[0007] (1) filtering the original high-throughput sequencing data obtained by second-generation genome sequencing to obtain filtered data; assembling the filtered data to obtain second-generation assembly contigs;

[0008] (2) filtering the original high-throughput sequencing data obtained by third-generation genome sequencing to obtain filtered data; assembling the filtered data to obtain third-generation assembly contigs; the assembly in this step is performed according to the size of the sample sequencing data, when the single sample sequencing data is more than 100G, the following method 1 or method 2 is used for single assembly; if it is less than or equal to 100G, the following two assembly methods are used simultaneously; wherein:

[0009] Method 1: The optimal K-mer is calculated by de Bruijn algorithm through reads length, sequencing depth and genome size, which is used as a fixed K-mer value to detect repetitive sequences and assemble contigs sequences;

[0010] Method 2: The reads are corrected according to the principle of majority over minority, and the regions with the same base are determined as high-quality regions according to the degree of overlap, and the low-quality regions that cannot completely overlap are trimmed; the sequence consistency of the overlapping region is recalculated, and the completely overlapping region is unitigs, and the overlapping region contains partially overlapping fragments, which is contigs; then the k-medoids clustering algorithm is used for iterative binning on the unitigs sequence file to obtain bins; and sequence determination: when all unitigs of a contig belong to the same bin, the contig is retained, otherwise, all unitigs are retained;

[0011] (3) Mixed assembly of second and third generation metaviromics sequencing data: the following two assembly methods are used at the same time to obtain contigs of second and third generation mixed assembly; wherein:

[0012] Method 1: Construct the assembly set of the filtered data in step (1) to obtain the second generation assembly set, and use the Bayesian model method to fill in the fragments that cannot be connected in the second generation assembly set with the filtered data in step (2) to obtain contigs;

[0013] Method 2: Use the HYBRIDSPADES algorithm to assemble the filtered data in step (1) and the filtered data in step (2) to obtain contigs;

[0014] (4) Construction of non-redundant contigs set: merge the second generation assembly contigs obtained in step (1), the third generation assembly contigs obtained in step (2), and the contigs of second and third generation mixed assembly obtained in step (3), cluster the sequences with more than 95% similarity into a cluster according to the similarity of the sequences, select the longest sequence as the representative sequence of the cluster, and all the representative sequences are the obtained non-redundant contigs set;

[0015] (5) Identification and determination of viral sequences: the non-redundant contigs set obtained in step (4) is subjected to viral identification and determination by no less than three methods, and the contigs identified as viruses by at least two methods are retained; the retained contigs are compared with the virus group database to determine the integrity, and the error matching more than 10% is filtered out, and the finally retained contigs are determined as viral sequences; the viral identification methods include:

[0016] Method 1: application of machine learning for virus identification using genomic features; the genomic features include structure, function, classification annotation and virus marker genes;

[0017] Method 2: identification of viral sequences based on K-mers word frequency of viral sequences and machine learning;

[0018] Method 3: based on deep learning algorithm, the network structure of double-channel convolutional neural network is designed, and the DNA fragments are effectively identified.

[0019] Preferably, the raw high-throughput sequencing data obtained by the second-generation genome sequencing in step (1) is filtered, specifically: the raw high-throughput data of the second-generation sequencing macrovirus group is filtered, the adapter sequence is removed, and the reads with a proportion of N greater than 10% are removed; after the data is filtered, the base quality analysis is performed by analyzing the composition and quality distribution of the bases, and the low-quality reads with a quality value Q≤30 are removed; the reads after quality control are removed based on the sequence alignment algorithm, the reads with a base number exceeding 85% aligned to the bacterial contaminant genome, and finally the high-quality clean reads are obtained.

[0020] Preferably, the assembly in step (1) is specifically: the assembly algorithm of de Bruijn graph based on K-mer iteration obtains the assembly result contig, and the short sequences with a length below 1.5 kb are filtered out; different data assembly methods are selected according to the number of samples and the sequencing amount; wherein:

[0021] Method 1: when the single sample sequencing amount is less than or equal to 20G and the number of samples is less than 100, single sample assembly and multi-sample assembly must be performed at the same time; single sample assembly means single sample assembly, and multi-sample assembly means mixed assembly of reads of multiple samples;

[0022] Method 2: when the single sample sequencing amount is greater than 20G, the single sample assembly method is selected.

[0023] Preferably, the raw high-throughput sequencing data obtained by the third-generation genome sequencing in step (2) is filtered, specifically: the raw high-throughput data of the third-generation sequencing macrovirus group is filtered, the adapter sequence is removed, and the subreads with less than 2 wrapping times are removed; based on the circular consistency sequencing CCS mode, the subreads derived from the same template chain are corrected for consistency to obtain high-precision single-molecule consistent reads; the reads after quality control are removed based on the sequence alignment algorithm, the reads with a base number exceeding 85% aligned to the bacterial contaminant genome, and finally the high-quality clean reads are obtained.

[0024] Preferably, step (5) is followed by species annotation of the viruses, using the final retained contigs obtained in step (5), calculating the average nucleotide identity ANI, retaining the viral group dataset with ANI > 95%, i.e. genus level; obtaining the annotation at the family level by aligning the genes at the amino acid level to a subset of viruses in the TrEMBL database; performing genus level species annotation of the viral genomes based on K-mer features by supervised machine learning methods; supplementing the annotation of the viral sequences at other taxonomic levels from the known taxonomy library according to the annotation results of the viruses at the family and genus levels.

[0025] Preferably, step (5) is followed by functional annotation of the viruses, detecting specific features of known benign and virulent viral sequences by a convolutional neural network, identifying whether the final retained contigs obtained in step (5) contain these features, and finally identifying the viral sequences as benign and virulent viruses.

[0026] Preferably, the retained contigs are aligned with a viral group database to determine the integrity, 0-5% error matching is considered as a high-quality contig, 5-10% error matching is a medium-quality contig, and more than 10% error matching is a low-quality contig that needs to be filtered.

[0027] Overall, compared with the prior art, the above technical solutions conceived by the present application mainly have the following technical advantages:

[0028] (1) The present application provides a complete process from processing of macroviral group data of second and third generation sequencing to mixed assembly of second and third generation, virus determination and identification, and species and functional annotation, which provides accurate analysis results for researchers and enables non-biological information professionals to independently complete viral high-throughput data analysis, thereby optimizing the work efficiency of researchers and reducing the cost of scientific research. The present application proposes a reliable analysis method for macroviral group based on second and third generation sequencing technology, which is simple to implement and widely applicable.

[0029] (2) The analysis process of the present application is more reasonable and comprehensive, and the analysis result is more accurate, which solves the problem of non-standard process of viral group assembly and identification by combining multiple generations of sequencing technology, and provides convenience and technical support for researchers. BRIEF DESCRIPTION OF DRAWINGS

[0030] Figure 1 is a flowchart of the present application. DETAILED DESCRIPTION

[0031] In order to make the objectives, technical solutions and advantages of the present application clearer, the present application will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and should not be used to limit the present application. In addition, the technical features involved in the various embodiments of the present application described below can be combined with each other as long as they do not conflict with each other.

[0032] The specific process of some examples of the present application is as follows:

[0033] Filtering of second-generation sequencing data: filtering the original high-throughput data of the second-generation sequencing macrovirus group, removing adapter sequences, removing reads with a proportion of N greater than 10%; performing base quality analysis based on the composition and quality distribution of the bases, removing low-quality reads with a quality value Q≤30; based on the short sequence alignment algorithm, removing reads with more than 85% of the bases aligned to contaminant genomes such as bacteria, and finally obtaining high-quality clean reads;

[0034] Specifically, the raw sequencing data is first evaluated for quality using fastqc, and based on the evaluation results, the software Trimmomatic and the software bowtie2 are used for quality filtering and DNA filtering.

[0035] Second-generation data assembly: the assembly algorithm based on the K-mer iterative de Bruijn graph obtains the assembly results contig, and filters out short sequences with a length of less than 1.5 kb; according to the number of samples and the sequencing amount, different data assembly strategies are selected. Among them:

[0036] Strategy 1: when the single-sample sequencing amount is less than 20G and the number of samples is less than 100, single-sample assembly and multi-sample assembly must be performed at the same time. Single-sample assembly means assembly of a single sample, and multi-sample assembly means mixed assembly of multiple sample reads;

[0037] Strategy 2: when the single-sample sequencing amount is greater than 20G, the single-sample assembly strategy can be selected;

[0038] Specifically, the software IDBA_UD is used to complete single-sample assembly, and the software megahit is used to complete multi-sample assembly.

[0039] Third-generation sequencing data filtering: filter the original high-throughput data of the third-generation sequencing macrovirus group, remove the adapter sequence, remove the subreads with less than 2 times of the number of loops; based on the circular consistency sequencing CCS mode, the subreads derived from the same template chain are corrected by consistency, and high-precision single-molecule consistent reads are obtained; the quality-controlled reads are based on sequence alignment algorithm, and the reads with more than 85% of the number of bases aligned to the contaminant genomes such as bacteria are removed, and finally the high-quality clean reads are obtained;

[0040] Specifically, the CCS software is used to complete quality filtering, and the bowtie2 software is used for DNA filtering.

[0041] Third-generation sequencing assembly: according to the size of the sample sequencing data, when the single sample sequencing data is more than 100G, a single assembly strategy can be used; if less than 100G, two assembly strategies are used at the same time; wherein:

[0042] Strategy 1: calculate the optimal K-mer based on the length of reads, sequencing depth and estimated genome size, and use it as a fixed K-mer value to detect repeated sequences and assemble contigs sequences, which is realized by using the software flye;

[0043] Strategy 2: adopt the principle of minority obeying majority to correct and modify reads, determine the region with base consistency as high-quality region according to the degree of overlap, and modify the low-quality region that cannot completely overlap; recalculate the sequence consistency of the overlapping region, and assemble to obtain unitigs with shorter fragments but more accurate bases and contigs with longer fragments but relatively poor quality, which is realized by using the software CANU; use k-medoids clustering algorithm for iterative binning of unitigs sequence file to obtain bins, which is realized by using MetaBAT; sequence determination: when all unitigs of contig belong to the same bin, the contig is retained, otherwise, all unitigs are retained, which is realized by using R, python or shell language.

[0044] Second-generation and third-generation macrovirus group sequencing data assembly: two assembly strategies are adopted. Among them:

[0045] Strategy 1: construct a splicing set of short reads (reads after filtering of second-generation sequencing data), and use Bayesian model method to improve the fragments that cannot be connected in the second-generation splicing set by long reads (reads after filtering of third-generation sequencing data) with low coverage, to obtain contigs, which is realized by using the software OPERA-MS;

[0046] Strategy 2: Assemble contigs using HYBRIDSPADES algorithm with both short reads and long reads, specifically using software metaSPAdes.

[0047] Constructing non-redundant contigs set: merge the contigs assembled by second generation, the contigs assembled by third generation and the contigs assembled by hybrid of second and third generation, cluster the sequences with more than 95% similarity into one cluster, preferably, select the longest sequence as the representative sequence of the cluster, all the representative sequences are the obtained non-redundant contigs set, specifically using software CD-hit to achieve.

[0048] Virus identification and determination: use no less than three methods to identify and determine the viruses of the non-redundant contigs, the contigs identified as viruses by at least two methods are retained, specifically using R, python or shell language to achieve; compare the retained contigs with the public complete virus group database to estimate the integrity, 0-5% error matching is considered as high-quality contig, 5-10% error matching is medium-quality contig, more than 10% error matching is low-quality contig and needs to be filtered, the finally retained contigs are determined as virus sequences, specifically using software CheckV to achieve. The virus set is used for subsequent species annotation and functional annotation. The virus identification method includes:

[0049] Strategy 1: apply machine learning to use genomic features (including structural / functional / classification annotation and virus marker genes) to identify viruses, specifically using software VirSorter to achieve;

[0050] Strategy 2: identify virus sequences based on virus sequence K-mers word frequency and machine learning, specifically using software VirFinder to achieve;

[0051] Strategy 3: based on deep learning algorithm, design the network structure of "double-channel convolutional neural network", effectively identify viruses from DNA fragments, specifically using software PPP-META to achieve.

[0052] Virus species annotation: Calculate the average nucleotide identity (ANI) of all virus contigs, retain the virus group dataset with ANI>95%, that is genus level, which is realized by using software FastANI; Get the annotation of family level by aligning the amino acid level genes with the virus subset in the TrEMBL database, which is realized by using software demovir; Based on K-mer characteristics, the genus level species annotation of virus genome is realized by using supervised machine learning method, which is realized by using software VirusTaxo; According to the annotation results of virus at family and genus levels, the annotation of other classification levels of virus sequences is perfected from the known taxonomy library, which is realized by using R or python language.

[0053] Virus function annotation: The valuable local features of virus sequence characteristics are detected by using convolutional neural network, and the virus sequence is recognized as mild and toxic virus, which is realized by using DeePhage.

[0054] In summary, the application constructs a more comprehensive and more accurate virus identification and recognition method. The application develops a metavirome assembly method based on second-generation sequencing, third-generation sequencing and second-third-generation sequencing, thereby solving the problem of non-standardization of the macro-virus assembly identification process of multi-generation sequencing, and integrating various assembly and identification tools, so that the analysis result is more accurate.

[0055] Those skilled in the art can easily understand that the above description is only a preferred embodiment of the application, and is not intended to limit the application, and any modification, equivalent replacement and improvement made within the spirit and principle of the application shall be included in the protection scope of the application.

Claims

1. A method for macroviral panel analysis based on second and third generation sequencing technology, characterized in that, Comprise the following steps: (1) The raw high-throughput sequencing data obtained by second-generation genome sequencing is filtered to obtain filtered data; the filtered data is assembled to obtain second-generation assembled contigs; (2) The raw high-throughput sequencing data obtained by third-generation genome sequencing is filtered to obtain filtered data; the filtered data is assembled to obtain third-generation assembled contigs; the assembly in this step is performed according to the size of the sample sequencing data, when the single sample sequencing data is more than 100G, the following method 1 or method 2 is used for single assembly method; if less than or equal to 100G, the following two assembly methods are used at the same time; wherein: Method 1: Calculate the optimal K-mer by de Bruijn algorithm through reads length, sequencing depth and genome size, and use it as fixed K-mer value to detect repetitive sequences and assemble contigs sequences; Method 2: Adopt the principle of minority to majority to correct reads, determine the base consistent area as high quality area according to the degree of overlap, and modify the low quality area that cannot completely overlap; recalculate the sequence consistency of the overlapping area, and the completely overlapping area is unitigs, and the overlapping area containing partially overlapping fragments is contigs; then use k-medoids clustering algorithm for iteration binning on unitigs sequence file to obtain bins; then sequence determination: when all unitigs of contig belong to the same bin, the contig is retained, otherwise, all unitigs are retained; (3) Mixed assembly of second-generation and third-generation macrovirus sequencing data: the following two assembly methods are used at the same time to obtain contigs assembled by second-generation and third-generation; wherein: Method 1: Construct the splicing set of the filtered data in step (1) to obtain the second-generation splicing set, use the Bayesian model method to fill the fragments in the second-generation splicing set that cannot be connected to obtain contigs; Method 2: Use HYBRIDSPADES algorithm to assemble the filtered data in step (1) and the filtered data in step (2) to obtain contigs; (4) Construct a non-redundant contigs set: merge the second-generation assembled contigs obtained in step (1), the third-generation assembled contigs obtained in step (2) and the contigs obtained by mixed assembly of second-generation and third-generation in step (3), cluster the sequences with more than 95% similarity into a cluster according to the similarity of the sequences, select the longest sequence as the representative sequence of the cluster, and all the representative sequences are the obtained non-redundant contigs set; (5) Virus sequence identification and determination: the non-redundant contigs set obtained in step (4) is subjected to virus identification and determination by not less than three methods, and the contigs identified as viruses by at least two methods are reserved; the reserved contigs are compared with the virus group database to determine the integrity, and the reads with more than 10% of the error matching are filtered out, and the finally reserved contigs are determined as virus sequences; the virus identification methods include: Method 1: machine learning is applied to identify viruses using genomic features; the genomic features include structure, function, classification annotation and virus marker genes; Method 2: virus sequence K-mers word frequency and machine learning are used to identify virus sequences; Method 3: based on deep learning algorithm, a network structure of double-channel convolutional neural network is designed to effectively identify DNA fragments as viruses.

2. The macrovirus panel analysis method based on second and third generation sequencing technology according to claim 1, wherein, The filtering of the raw high-throughput sequencing data obtained by the second-generation genome sequencing in step (1) is as follows: the raw high-throughput data of the second-generation sequencing macrovirus group is filtered to remove adapter sequences and reads with a proportion of N greater than 10%; after filtering, the data is analyzed for base composition and quality distribution, and low-quality reads with a quality value Q≤30 are removed; The quality-controlled reads are based on sequence alignment algorithm, and the reads with more than 85% of the base number aligned to the bacterial contaminant genome are removed, and finally the clean reads of high quality are obtained.

3. The macrovirus panel analysis method based on second and third generation sequencing technology according to claim 1 or 2, characterized in that, The assembly in step (1) is as follows: the assembly algorithm based on K-mer iteration de Bruijn graph is used to obtain the assembly results contig, and short sequences with a length of less than 1.5 kb are filtered out; different data assembly methods are selected according to the number of samples and the sequencing amount; wherein: Method 1: when the single sample sequencing amount is less than or equal to 20G and the number of samples is less than 100, single sample assembly and multi-sample assembly must be performed simultaneously; single sample assembly means single sample assembly, and multi-sample assembly means mixed assembly of multiple sample reads; Method 2: when the single sample sequencing amount is greater than 20G, the single sample assembly method is selected.

4. The method of macroviral panel analysis based on second and third generation sequencing technology according to claim 1, wherein, The filtering of the raw high-throughput sequencing data obtained by the third-generation genome sequencing in step (2) is as follows: the raw high-throughput data of the third-generation sequencing macrovirus group is filtered to remove adapter sequences and subreads with less than 2 loop numbers; based on the circular consistency sequencing CCS mode, the subreads derived from the same template chain are corrected for consistency to obtain high-precision single-molecule consistent reads; The quality-controlled reads are based on sequence alignment algorithm, and the reads with more than 85% of the base number aligned to the bacterial contaminant genome are removed, and finally the clean reads of high quality are obtained.

5. The method of macrovirus panel analysis based on second and third generation sequencing technology according to claim 1, wherein, Step (5) is followed by species annotation of viruses, using the final retained contigs from step (5), average nucleotide identity (ANI) is calculated, and viral group datasets with ANI > 95% are retained, i.e. genus level; annotations at family level are obtained by aligning genes at amino acid level against a subset of viruses in the TrEMBL database; genus level species annotation of viral genomes is performed based on K-mer features by supervised machine learning methods; annotations at other taxonomic levels of viral sequences are supplemented from known taxonomy databases based on the annotations at family and genus levels of viruses.

6. The macrovirus panel analysis method based on second and third generation sequencing technology according to claim 1, wherein, Step (5) is followed by functional annotation of viruses, specific features of known mild pathogenic viral sequences are detected by convolutional neural networks, and whether the final retained contigs from step (5) contain these features are identified, and finally the viral sequences are identified as mild pathogenic viruses.

7. The macrovirus panel analysis method based on second and third generation sequencing technology according to claim 1, wherein, The retained contigs are aligned with a viral group database to determine integrity, 0-5% error matching is considered as high-quality contigs, 5-10% error matching is medium-quality contigs, and more than 10% error matching is low-quality contigs which need to be filtered.

Images