A method for analyzing macroviral group data

Abstract

The application discloses a macrovirus group data analysis method, and belongs to the technical field of macrovirology. The method comprises the following steps: sequence quality control, sequence assembly, sequence clustering, virus sequence identification, virus sequence checking, virus abundance calculation, species annotation, virus lifestyle judgment, virus host prediction and virus auxiliary metabolism gene analysis. The application uses Trimmomatic software, BWA-MEN, Megahit software, CD-HIT and other tools to execute the analysis process of macrovirus group data. Practice proves that the application can accurately identify and annotate virus species, comprehensively and systematically deeply analyze and mine macrovirus group data, the steps are simple and clear, the analysis time is short, and the effect of macrovirus identification research is greatly improved.

Description

TECHNICAL FIELD

[0001] The present application relates to a macrovirus group data analysis method, belonging to the technical field of macrovirology. BACKGROUND

[0002] Viruses are the most numerous life-like organisms on Earth, and they have an important contribution to the biochemical cycle of the Earth. Viral metagenomics, also known as virome, is a new branch of science that has emerged on the basis of metagenomics theory and in combination with existing viral molecular biology detection techniques.

[0003] The macrovirus group directly takes all the genetic material (DNA or RNA) of viruses in the sample as the research object, enriches the viral particles first, then obtains the genomic sequence information, and identifies all the viral components and relative abundance, which is a powerful means for discovering new viruses, virus infection early warning and control. It has important significance in the research fields of virus origin and evolution mode, genetic diversity and geographical distribution, and the relationship between viruses and hosts.

[0004] Macro-virus group research can be applied to human or animal blood, tissue, feces and other samples, plant tissue samples, and various environmental samples such as water and soil, to analyze the viral community therein. However, due to the generally small size of viral genomes, the relative content of viral nucleic acid in the sample is very low, the host's genomic sequence is severely interfered, and the number of known viruses and genomic information is limited, which makes it difficult to analyze macro-virus group research data.

[0005] The current macro-virus group data analysis method is not accurate, not systematic, and relatively cumbersome, and the analysis process takes a long time, which hinders the comprehensive understanding of the viral community and is not conducive to the in-depth mining of sequencing data. SUMMARY

[0006] In order to solve the problems of complicated macro-virus group data analysis process, poor accuracy, and lack of systematization, the present application provides a macro-virus group data analysis method, comprising:

[0007] Step 1: sequence quality control;

[0008] Use Trimmomatic software to remove low-quality sequences, and use BWA-MEN to remove sequences that may contain eukaryotic genomes;

[0009] Step 2: sequence assembly and clustering;

[0010] Use Megahit software to assemble the high-quality sequences obtained in step 1 to obtain contigs sequences, and use CD-HIT to remove redundancy and cluster the contigs sequences;

[0011] Step three: Virus sequence identification and inspection;

[0012] For the contigs sequences obtained after clustering, use VirSorter2 to identify viral sequences, and use Check-V to check the integrity and degree of contamination of the viral sequences;

[0013] Step four: Virus species annotation and abundance calculation;

[0014] Use Hmmer v3.1b2 to align the hmm model of the reference protein, and perform species annotation according to the protein alignment result, and use BWA software to calculate the abundance of the viral sequence;

[0015] Step five: Virus lifestyle judgment: use PhaTYP software to judge the lifestyle of the virus;

[0016] Step six: Virus host prediction: use CRISPRone to search for CRISPR intervals, tRNA matches, and genomic homology matches in bacterial sequences to predict viral hosts;

[0017] Step seven: Virus auxiliary metabolism gene analysis: use Vibrant and DRAM-V to determine the auxiliary metabolism genes encoded by the virus.

[0018] Optionally, the step six comprises:

[0019] Use CRISPRone to search for CRISPR intervals in bacterial sequences, use Blastn to align the CRISPR interval sequences with the viral contigs, select the alignment results with homology greater than or equal to 95% and single nucleotide polymorphism less than or equal to 2, and identify the bacterial sequence corresponding to the CRISPR interval as a putative viral host;

[0020] Use tRNAscan-SE software to predict tRNA genes in the viral contigs, use Blastn to align the tRNA sequences with the bacterial scaffolds, and retain the alignment results with sequence identity greater than or equal to 95% as viral hosts;

[0021] Directly align the bacterial sequences with the viral contigs by Blastn, and filter the predicted hosts of the virus according to the parameters: bitscore≥50, e-value≥0.003, identity≥70%, and matching length≥2500bp.

[0022] Optionally, the step one uses Trimmomatic software to trim adapters and other Illumina-specific sequences from reads; performs sliding window trimming, window size is 4bp, moving from the first base position at 5' end, the average quality of bases in window is greater than or equal to Q15, and the sequence is truncated at the 3' end of the first window with an average quality less than Q15; after the above quality screening, sequences with a length of less than 36bp are removed, and low-quality sequences are removed.

[0023] Optionally, the step one uses BWA-MEN software to establish an index of the host genome, then aligns the host genome with sample reads, and finally uses samtools to remove host sequences.

[0024] Optionally, in the step two, Megahit software is used, the minimum kmer is set to 35, the maximum kmer is set to 95, the kmer interval is set to 20, and the minimum value of the output contig is set to 500.

[0025] Optionally, in the step two, the est mode of CD-HIT is used, the similarity is 0.9, the word size is 5 when comparing each pair of sequences, 16GB RAM is used, and the number of threads is 4.

[0026] Optionally, in the step three, VirSorter2 is used, the minimum sequence length is selected to be 1500bp, and the number of threads is 4.

[0027] Optionally, in the step three, the end_to_end mode of Check-V software is used, and the number of threads is 40.

[0028] Optionally, in the step four, BWA software is used to compare the virus contig file with the clean read, the mem mode of BWA is run, the kmer length is 30, the number of threads is 40, and the output file is calculated using pileup.sh.

[0029] Optionally, in the step four, Hmmer v3.1b2 software is used for virus species annotation, including:

[0030] First, download the VOGDB database and the species corresponding classification table, then use the program Hmmscan under Hummer to compare the virus contig with the hmm model of the database, set E-value to 0.0001, and the number of cpus to 40;

[0031] It is required that at least 20% of the genes of each contig have annotations, and then at least 60% of the 20% genes come from the same genus, otherwise it is considered to be unannotated.

[0032] The present application has the following advantages:

[0033] The virus group data analysis method of the present application removes low-quality sequences using Trimmomatic software and removes sequences that may exist in eukaryotic genomes using BWA-MEN, thereby improving the quality of virus sequences and reducing the interference of host genomes; high-quality sequences obtained in step one are assembled using Megahit software to obtain contigs sequences, and the contigs sequences are de-redundant and clustered using CD-HIT; for the contigs sequences obtained after clustering, VirSorter2 is used to identify virus sequences, Check-V is used to check the integrity and degree of contamination of the virus sequences; Hmmer v3.1b2 is used to align the hmm model of the reference protein, and species annotation is performed according to the protein alignment result; BWA software is used to calculate the abundance of virus sequences; PhaTYP software is used to determine the lifestyle of viruses; CRISPRone is used to search for CRISPR intervals, tRNA matches, and genomic homology matches in bacterial sequences to predict virus hosts; Vibrant and DRAM-V are used to determine the auxiliary metabolic genes encoded by viruses.

[0034] The method of the present application can accurately identify and annotate virus species, comprehensively and systematically analyze and mine macrovirus group data, has simple and clear steps, short analysis time, and greatly optimizes the effect of macrovirus identification research. BRIEF DESCRIPTION OF DRAWINGS

[0035] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the drawings needed in the embodiment description will be briefly introduced. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor.

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

[0037] In order to make the purpose, technical solutions and advantages of the present application clearer, the embodiments of the present application will be further described in detail below with reference to the drawings.

[0038] Embodiment one:

[0039] The present embodiment provides a macrovirus group data analysis method, comprising:

[0040] Step one: sequence quality control;

[0041] Using Trimmomatic software, low-quality sequences are removed, and BWA-MEN is used to remove sequences that may exist in eukaryotic genomes.

[0042] Step two: sequence assembly and clustering

[0043] Using Megahit software to assemble high-quality sequences obtained in the above step one, contigs sequences are obtained, and CD-HIT is used to remove redundancy and cluster the contigs sequences;

[0044] Step three: virus sequence identification and inspection

[0045] For the contigs sequences obtained after clustering, VirSorter2 is used to identify viral sequences, and Check-V is used to check the integrity and degree of contamination of the viral sequences;

[0046] Step four: virus species annotation and abundance calculation

[0047] Using Hmmer v3.1b2 to align the hmm model of the reference protein, the species annotation is performed according to the protein alignment result, and the abundance of the viral sequence is calculated using BWA software;

[0048] Step five: virus lifestyle judgment: using PhaTYP software to judge the lifestyle of the virus

[0049] Step six: virus host prediction: using CRISPRone to search for CRISPR interval, tRNA matching, and genomic homology matching in bacterial sequences to predict the virus host;

[0050] Step seven: virus auxiliary metabolism gene analysis: using Vibrant and DRAM-V to determine the auxiliary metabolism genes encoded by the virus.

[0051] Example two

[0052] The present embodiment provides a method for analyzing macroviral group data, as shown in Figure 1 , comprising:

[0053] Step one: sequence quality control

[0054] Using Trimmomatic software, low-quality sequences are removed, and BWA-MEN is used to remove sequences that may exist in eukaryotic genomes, and the specific operation is as follows:

[0055] The adapter and other Illumina-specific sequences, i.e. index sequences, are trimmed from the reads using the Trimmomatic software; sliding window trimming is performed, with a window size of 4 bp, starting from the first base position at the 5' end, and moving the window, with an average base quality of Q15 in the window, and the sequence is truncated at the 3' end of the first window with an average quality lower than Q15; after the above quality screening, sequences with a length of less than 36 bp are removed, and low-quality sequences are removed, i.e. paired reads with adapters are removed, duplicate reads generated by PCR amplification are removed, and ployX sequences are removed.

[0056] The host genome index is established using the BWA-MEN software, and then the host genome is aligned with the sample reads, and finally the host sequences are removed using the samtools.

[0057] Step two: sequence assembly and clustering

[0058] The high-quality sequences obtained in the above step one are assembled using the Megahit software to obtain contigs sequences, and the contigs sequences are de-redundant and clustered using the CD-HIT, with the following specific parameter settings:

[0059] The Megahit software is used, with the minimum kmer set to 35, the maximum kmer set to 95, the kmer interval set to 20, and the minimum contig output value set to 500;

[0060] The est mode of the CD-HIT is used, with a similarity of 0.9, a word size of 5 when comparing each pair of sequences, 16 GB RAM, and 4 threads.

[0061] Step three: virus sequence identification and inspection

[0062] The contigs sequences obtained after clustering are used for virus sequence identification using the VirSorter2, and the integrity and contamination degree of the virus sequences are inspected using the Check-V, with the following specific parameter settings:

[0063] The VirSorter2 selects a minimum sequence length of 1500 bp and a thread number of 4; the Check-V software is set to the end_to_end mode and a thread number of 40.

[0064] Step four: virus species annotation and abundance calculation

[0065] The Hmmer v3.1b2 is used to align the reference protein hmm model, the species annotation is performed according to the protein alignment result, and the BWA software is used to calculate the abundance of the virus sequences.

[0066] Using BWA software to compare the virus contig file with clean read, running the mem mode of BWA, kmer length is 30, thread number is 40, output file uses pileup.sh to calculate coverage.

[0067] Using Hmmer v3.1b2 software to annotate virus species, including:

[0068] First, download VOGDB database and species corresponding classification table, then use the program Hmmscan under Hummer to compare virus contig and database hmm model, set E-value to 0.0001, cpu number to 40;

[0069] It is required that at least 20% of the contigs have annotations, and then at least 60% of the 20% genes come from the same genus, otherwise it is considered to be unannotated.

[0070] Step five: Virus lifestyle judgment: use PhaTYP software to judge the lifestyle of virus;

[0071] Step six: Virus host prediction: use CRISPRone to search for CRISPR interval, tRNA matching, and genomic homology matching in bacterial sequences to predict the virus host, as follows:

[0072] Using CRISPRone to search for CRISPR interval in bacterial sequences, using Blastn to compare CRISPR interval sequence with virus contig, selecting comparison results with homology greater than or equal to 95% and single nucleotide polymorphism less than or equal to 2, corresponding to the bacterial sequence where the CRISPR interval is located, is identified as a putative virus host;

[0073] Using tRNAscan-SE software to predict tRNA genes in virus contig, using Blastn to compare tRNA sequence with bacterial scaffold, retaining sequence identity greater than or equal to 95% comparison results identified as virus host;

[0074] By Blastn, bacterial sequences are directly compared with virus contig, and according to the parameters: bitscore≥50, e-value≥0.003, identity≥70%, matching length≥2500bp, the predicted host of virus is screened.

[0075] At the same time, three methods are used to predict the host of virus contig, and the results of the three methods are integrated, and the results are taken as the union.

[0076] Step seven: Virus-assisted metabolic gene analysis: determine virus-encoded auxiliary metabolic genes using Vibrant and DRAM-V.

[0077] The present embodiment improves the quality of viral sequences and reduces the interference of host genomes by sequence quality control, removing low-quality sequences using Trimmomatic software, and removing sequences that may exist in eukaryotic genomes using BWA-MEN; the high-quality sequences obtained in step one are assembled using Megahit software to obtain contigs sequences, and the contigs sequences are de-redundant and clustered using CD-HIT; for the contigs sequences obtained after clustering, VirSorter2 is used to identify viral sequences, Check-V is used to check the integrity and degree of contamination of the viral sequences; Hmmer v3.1b2 is used to align the reference protein hmm model, and species annotation is performed according to the protein alignment result, BWA software is used to calculate the abundance of viral sequences; PhaTYP software is used to determine the lifestyle of viruses; CRISPRone is used to search for CRISPR intervals, tRNA matches, and genomic homology matches in bacterial sequences to predict viral hosts; Vibrant and DRAM-V are used to determine virus-encoded auxiliary metabolic genes.

[0078] The method of the present embodiment can accurately identify and annotate viral species, comprehensively and systematically analyze and mine macroviral group data, and has simple and clear steps, short analysis time, and greatly improved macrovirus identification research effect.

[0079] Some steps in the embodiments of the present application can be implemented using software, and the corresponding software programs can be stored in a readable storage medium, such as an optical disc or a hard disk.

[0080] The above only describes the preferred embodiments of the present application and does not limit the present application, and any modification, equivalent replacement, improvement, etc. made within the spirit and principles of the present application shall be included in the protection scope of the present application.

Claims

1. A method of analyzing macroviral group data, characterized by, The method comprises: Step 1: sequence quality control; Using Trimmomatic software, low-quality sequences are removed, and BWA-MEN is used to remove sequences that may exist in eukaryotic genomes; Step 2: sequence assembly and clustering; Using Megahit software, the sequences obtained in step 1 are assembled to obtain contigs sequences, and CD-HIT is used to remove redundancy and cluster the contigs sequences; Step 3: virus sequence identification and inspection; Using VirSorter2, the contigs sequences obtained after clustering are used to identify virus sequences, and Check-V is used to check the integrity and degree of contamination of the virus sequences; Step 4: virus species annotation and abundance calculation; Using Hmmer v3.1b2 to align the hmm model of the reference protein, the species annotation is performed according to the protein alignment result, and the abundance of the virus sequence is calculated using BWA software; Step 5: virus lifestyle judgment: using PhaTYP software to judge the virus lifestyle; Step 6: virus host prediction: using CRISPRone to search for CRISPR intervals, tRNA matches, and genomic homology matches in bacterial sequences to predict virus hosts; Step 7: virus auxiliary metabolism gene analysis: using Vibrant and DRAM-V to determine virus-encoded auxiliary metabolism genes.

2. The method of analyzing macroviral group data according to claim 1, wherein, The step 6 comprises: Using CRISPRone to search for CRISPR intervals in bacterial sequences, using Blastn to align CRISPR interval sequences with virus contigs, selecting alignment results with homology greater than or equal to 95% and single nucleotide polymorphism less than or equal to 2, and identifying the bacterial sequence corresponding to the CRISPR interval as a putative virus host; Using tRNAscan-SE software to predict tRNA genes in virus contigs, using Blastn to compare tRNA sequences with bacterial scaffolds, and retaining alignment results with sequence identity greater than or equal to 95% to identify virus hosts; By Blastn, bacterial sequences are directly aligned with virus contigs, and according to the parameters: bitscore≥50, e-value≥0.003, identity≥70%, and matching length≥2500bp, the predicted hosts of the virus are screened.

3. The method of claim 1, wherein the macrovirus group data is analyzed by, In step 1, Trimmomatic software is used to trim adapters and other Illumina-specific sequences from sample reads; sliding window trimming is performed, with a window size of 4bp, starting from the first base position at the 5' end, and moving the window. The average quality of the bases in the window is greater than or equal to Q15, and the sequence is truncated at the 3' end of the first window with an average quality less than Q15. After the above quality screening, sequences with a length of less than 36bp are removed, and low-quality sequences are removed.

4. The method of claim 1, wherein the macrovirus group data is analyzed by, In step 1, BWA-MEN software is used to establish an index of the host genome, then the host genome is aligned with the sample reads, and finally samtools is used to remove host sequences.

5. The method of claim 1, wherein the macrovirus group data is analyzed by, In the second step, the minimum kmer is set as 35, the maximum kmer is set as 95, the kmer interval is set as 20, and the minimum contig sequence is set as 500.

6. The method of claim 1, wherein the macrovirus group data is analyzed by, In the second step, the est mode of CD-HIT is used, the similarity is set as 0.9, the word size is set as 5 when comparing two sequences, 16 GB RAM is used, and the thread number is set as 4.

7. The method of claim 1, wherein the macrovirus group data is analyzed by a method comprising: In the third step, VirSorter2 is used, the minimum sequence length is set as 1500 bp, and the thread number is set as 4.

8. The method of claim 1, wherein the macrovirus group data is analyzed by, In the third step, the end_to_end mode of Check-V software is used, and the thread number is set as 40.

9. The method of claim 1, wherein the macrovirus group data is analyzed by, In the fourth step, the BWA software is used to compare the virus contig file with the clean read, the mem mode of BWA is run, the kmer length is set as 30, the thread number is set as 40, and the output file is calculated by using pileup.sh.

10. The method of claim 1, wherein the macrovirus group data is analyzed by, In the fourth step, the Hmmer v3.1b2 software is used for virus species annotation, including: First, download the VOGDB database and the species corresponding classification table, then use the program Hmmscan under Hummer to compare the virus contig with the hmm model of the database, set the E-value as 0.0001, and the cpu number as 40; Require at least 20% of the genes of each virus contig to have annotation, then at least 60% of the 20% genes come from the same genus, otherwise it is considered as unannotated.

Images