Method for classifying adeno-associated virus data on multiple reference genomes based on next-generation sequencing

Abstract

The invention belongs to a second-generation sequencing data analysis technology in the field of bioinformatics, and particularly relates to a method for classifying adeno-associated virus data on multiple reference genomes based on second-generation sequencing. The method comprises the following steps: S1, constructing a position information file of a homologous sequence on a genome; s2, establishing an original BAM file; s3, constructing a comparison information summary list T; s4, a QNAME list is constructed; s5, filtering; and S6, counting. According to the method disclosed by the invention, the main components of AAV are identified and analyzed on the computer algorithm level through a short fragment sequencing method for the first time. Compared with an experiment-based method, the analysis result of the method has higher component quantitative accuracy and lower manual operation cost. The method provided by the invention is of great significance to clinical detection of AAV drugs.

Description

Technical Field

[0001] The present invention belongs to the field of bioinformatics and relates to a second-generation sequencing data analysis technology, and in particular to a classification method for adeno-associated virus data based on second-generation sequencing on multiple reference genomes. Background Art

[0002] Recombinant adeno-associated virus (rAAV) is a small virus that is often used as a vector for gene therapy. It is known as a powerful assistant for gene therapy because of its simple structure, high safety and low immunogenicity. AAV belongs to the genus Dependovirus of the family Parvoviridae and is a non-enveloped single-stranded DNA defective virus.

[0003] In the application of AAV as a gene vector, the AAV genome often needs to be modified to insert the target gene. The production process of rAAV vectors involves multiple complex steps, including cloning of viral genes, transfection of host cells, and collection and purification of viral particles. Therefore, the preparation process of rAAV vectors and the characteristics of gene therapy products also bring certain safety risks. In order to ensure the consistency of the quality and function of the rAAV vector finally obtained, it is necessary to identify the components of the cell fluid extracted from the experiment. Through component identification, we can identify pure rAAV vectors and exclude possible contaminants and non-target components.

[0004] At present, there are characterizations of impurities in adeno-associated virus (AAV) samples and formulation compositions for stabilizing AAV (e.g., patent application number CN202180073988.2). This method is mainly used to identify nucleic acid impurities in samples of rAAV vectors, specifically based on experimental methods of size exclusion chromatography (SEC) and spectrophotometry. The steps include: contacting a sample containing an AAV vector with a SEC column; then eluting the SEC column with a variety of eluents to separate molecules of different sizes; finally, using a spectrophotometer to detect nucleic acid impurities in the eluent, which is usually achieved by measuring the absorbance at a specific wavelength.

[0005] However, the current methods for identifying rAAV vector data components mainly rely on experimental methods, which have several significant shortcomings. For example, experimental methods lead to higher costs; the experimental process is also accompanied by more complicated manual operations; and when quantifying rAAV vector data components, experimental methods also have many limitations on the accuracy of quantitative data.

[0006] On the other hand, humans and AAV have a very long history of coexistence. AAV is a virus that is widely present in nature. It is also common in humans and other primates, and is distributed in large quantities in the liver, lungs, heart and other organs of the human body. There are 11 known serotypes of AAV. However, when using computer algorithms to identify the components of AAV samples, the sequence homology between various genomes has not yet been introduced as an analysis content. Directly comparing the sequencing sequence with each genome will result in a high false positive rate in the classification results. Due to the existence of homologous sequences between different genomes, a sequence may be compared to multiple genomic sources. This phenomenon is called a false positive in the comparison result.

[0007] In addition, host genome fragments exist as impurities in rAAV in two forms: one is the whole genome mixed in at a lower frequency, and the other is a smaller fragment integrated into each rAAV genome. The current method cannot effectively and intuitively distinguish the two.

[0008] In summary, at present, there is still a lack of a widely recognized bioinformatics algorithm to accurately determine the unique genomic source of a sequence, as well as to identify the components of AAV vector sequencing data and output a visual report. Therefore, it is necessary to improve it. Summary of the invention

[0009] In view of the above-mentioned problems in the actual analysis of the component identification of recombinant adeno-associated virus (rAAV) vector data, the present invention has developed AsCA (AAV sequence Classification Algorithm), i.e., an AAV sequence classification algorithm tool. The AsCA proposed a complete bioinformatics solution for the first time, introduced homology sequence judgment logic, effectively reduced the false positive of the comparison results, and improved the credibility of the classification and identification results. Therefore, the purpose of the present invention is to provide a classification method for adeno-associated virus data based on second-generation sequencing on multiple reference genomes.

[0010] In order to achieve the above object, the present invention adopts the following technical solution:

[0011] A classification method of adeno-associated virus data based on second-generation sequencing on multiple reference genomes, comprising the following steps;

[0012] S1. Construct the location information file of homologous sequences on the genome;

[0013] S2, constructing an original host BAM file, an original virus BAM file, and multiple original impurity BAM files;

[0014] S3, constructing a comparison information summary list T;

[0015] S4, build QNAME list;

[0016] S5, filtering: filtering the three original BAM files in S2 according to the QNAME list to obtain a host BAM file, a virus BAM file, and multiple impurity BAM files;

[0017] S6. Statistics.

[0018] According to the present invention, a classification method of adeno-associated virus data based on second-generation sequencing on multiple reference genomes includes the following steps:

[0019] S1. Construction of the location information file of homologous sequences on the genome: Perform homology analysis on the FASTA file of the host genome, the FASTA file of the viral genome, and the FASTA file of the impurity genome, respectively, to generate two files of the location information of homologous sequences on the genome, which record the location information of homologous sequences on each genome;

[0020] S2. Creation of original BAM files: Compare and analyze the double-ended Fastq files with the FASTA files of the host genome, the FASTA files of the virus genome, and the FASTA files of the impurity genome, respectively, to obtain three original BAM files, wherein the double-ended Fastq files are files storing sequence information of viruses, hosts, and impurities obtained after sequencing;

[0021] S3, construction of alignment information summary list T: read the alignment information in the three original BAM files in S2, and merge all alignment information of a sequence together according to QNAME to obtain the alignment information summary list T;

[0022] S4, construction of QNAME list: Classify each sequence in the comparison information summary list T in S3, and the classification judgment logic includes sequences aligned to a single genome and sequences aligned to multiple reference genome sources to determine the QNAME list of the virus, the QNAME list of the host, and one of the QNAME lists of multiple impurities, and obtain the QNAME list of the virus, the QNAME list of the host, and the QNAME list of multiple impurities; at this point, each sequence has and only has a unique reference genome source, and the sequence QNAME of each source is recorded separately in the QNAME list of the corresponding source;

[0023] S5, filtering: filtering the three original BAM files in S2 according to the QNAME list of the virus, the QNAME list of the host, and the QNAME list of multiple impurities to obtain a host BAM file, a virus BAM file, and multiple impurity BAM files, and each sequence appears in and only appears in one BAM file;

[0024] S6. Statistics: Perform statistics on the result data for the multiple BAM files obtained in S5.

[0025] According to the present invention, in S4, in the alignment information summary list T, for sequences that are aligned only to a single genome, they are directly classified into the genome source, and their QNAMEs are recorded in the corresponding QNAME list.

[0026] According to the present invention, in S4, in the comparison information summary list T, for sequences compared to multiple reference genome sources, it is determined whether they are in a homologous region based on the position information, and the sequences are divided into two parts, one part being "sequences not in a homologous region" and the other part being "sequences in a homologous region".

[0027] Furthermore, in the comparison information summary list T, for sequences that are not in the homologous region and are compared to multiple reference genome sources, the sequence matching lengths recorded by the CIGAR values ​​of the sequence on different reference genomes are compared, and the sequence is classified into the genome source with the longest matching length, and its QNAME is recorded in the corresponding QNAME list. For example, if the CIGAR value of a sequence in the viral genome is 150M, the CIGAR value in the host genome is 100M50D, and the CIGAR value in the impurity genome is 70M80D, it means that the matching lengths of the sequence on the three reference genomes are 150, 100, and 70, respectively, that is, the sequence has the longest matching length on the viral genome, so the sequence will be classified as a viral genome sequence, and its QNAME will also be recorded in the viral QNAME list.

[0028] Furthermore, in the comparison information summary list T, for sequences in the homologous region that are compared to multiple reference genome sources, they are preferentially considered not to be host sequences. The sequence matching lengths recorded by the CIGAR values ​​of the sequence on different reference genomes are compared, and the sequence is classified into the genome source with the longest matching length, and its QNAME is recorded in the corresponding QNAME list.

[0029] According to the present invention, the filtering step of S5 is: processing the three original BAM files obtained in S1 according to the QNAME list, if a sequence appears in the QNAME list of the virus, the corresponding sequences in the host original BAM file and the impurity original BAM file are filtered out; after filtering, a host BAM file, a virus BAM file and multiple impurity BAM files are obtained, wherein the number of impurity BAM files changes dynamically with the number of impurity types during actual analysis; at this point, each sequence appears and only appears in one BAM file.

[0030] According to the present invention, it also includes S7, mutation detection: for the virus BAM file obtained in S5, mutation detection is performed using a mutation detection tool, including SNV detection and InDel detection, to obtain a mutation detection result table in a VCF file format, and to visualize the SNV base ratio and the InDel distribution on the genome, respectively.

[0031] According to the present invention, the statistics in step S6 include respectively counting the following information:

[0032] 1) Count the number of sequences stored in each BAM file of a sample to obtain the quantitative results of the sample's components; 2) Count the sequence position information in the virus BAM file to obtain the sequence coverage on the virus genome; 3) Count the sequence depth of each chromosome in the host BAM file, calculate the depth ratio of the sample on each chromosome, use the same method to calculate the depth ratio of the standard on each chromosome, compare the depth ratio of the sample with the standard and perform log processing to obtain the sequence density information of the sample on each chromosome.

[0033] Furthermore, the statistical results in step S6 are collated and visualized, and the component identification results and visualization report of the AAV vector data are output.

[0034] According to the present invention, the double-ended Fastq file includes a double-ended Fastq file obtained after sequencing of an extracted sample containing AAV, impurities, and host cell fluid.

[0035] Compared with the prior art, the present invention has the following beneficial effects: for the first time, the main components of AAV are identified and analyzed at the computer algorithm level by short fragment sequencing. Compared with the experimental method, the analysis results of the present invention have higher component quantitative accuracy and lower human operation costs. This progress is of great significance to the clinical detection of AAV drugs. Specifically, it is embodied in:

[0036] 1) The present invention uses large-scale short-fragment sequencing for the first time, and uses information analysis algorithms to classify and process sequencing data for the first time, achieving accurate classification and quantification of rAAV data components, as well as mutation detection and visualization of rAAV data. Compared with the current experimental-based methods, the algorithm AsCA of the present invention greatly improves time efficiency and reduces the tediousness of manual operation in an automated process.

[0037] 2) When performing component identification, the present invention introduces a homology sequence comparison algorithm. First, the corresponding genomes of the host, AAV and impurities are used to generate corresponding inter-genomic homologous sequence files to assist in determining whether the sequencing sequence is a homologous sequence, and comprehensively consider its homology and the CIGAR value on each genome to determine the source. When a sequence is aligned to multiple genomes at the same time and is in a homologous region, it is determined based on the Bayesian theorem that it is more likely to belong to a viral genome. Based on this logic, each sequence is classified into a unique source, and the BAM file is extracted and filtered accordingly. This method improves the situation where a sequence is aligned to multiple genome sources due to homologous similar sequences, greatly improving the accuracy of the results.

[0038] 3) The present invention further subdivides the impurity sequences derived from the human genome according to the chromosome source, determines the chromosome from which the impurity sequences are derived, calculates, and compares the ratio of the sample to the standard in each chromosome, so as to more intuitively distinguish whether the host genome fragments are mixed in with a lower frequency in the whole genome as a whole, or are integrated in the rAAV sequence with a high frequency in smaller fragments. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 Schematic diagram of AsAC analysis steps; first, the reference genome FASTA comparison generates homologous sequence files, then the Fastq files are aligned to each reference genome respectively, and then the sequence source is determined based on the alignment information of each sequence on each reference genome. The original BAM file is filtered according to the classification results of the sequence, so that each sequence is stored and only stored in a BAM file from one source, mutation detection is performed on the BAM file of the virus source, and finally the corresponding mutation detection results and visualization reports are output;

[0040] Figure 2 This is a bar chart of the components of each sample. The horizontal axis is each sample, and the vertical axis is the sequence number statistics; the height of each column represents the total number of sequences that can be successfully matched for the sample, and each color represents a type of genome; standards and negative controls are used to confirm that the experimental results are in line with expectations, and further verify and improve the credibility of the sample results;

[0041] Figure 3 It is a scatter plot of the ratio of the depth of each sample to the standard on the host genome; the horizontal axis is the host chromosome, and the vertical axis is the ratio (after log processing);

[0042] Figure 4 It is a line graph of the coverage of each sample on the viral genome; the horizontal axis is the base position of the viral genome, and the vertical axis is the base coverage;

[0043] Figure 5It is a line graph of the frequency distribution of SNV mutation bases of samples on the viral genome; the horizontal axis is the viral genome, the vertical axis is the mutation frequency of each type of mutation base, and the 4 different broken lines represent 4 types of bases respectively;

[0044] Figure 6 The InDel mutation distribution of the sample on the viral genome; the horizontal axis is the viral genome, the vertical axis is the mutation frequency, and the size of the point represents the mutation length;

[0045] Figure 7 It is a comparison of the number of sequences before and after the removal of false positives determined by homologous sequences; the horizontal axis is the genome classification of the virus and the host, and the vertical axis is the number of sequences. DETAILED DESCRIPTION

[0046] In order to make the technical problems, technical solutions and advantages to be solved by the present invention more clear, a detailed description will be given below with reference to the accompanying drawings and specific embodiments.

[0047] The specific meanings of the terms in the present invention are as follows:

[0048] 1. AAV: Adeno-Associated Virus;

[0049] 2. AsCA: AAV sequence Classification Algorithm, AAV sequence classification algorithm tool;

[0050] 3. Fastq file: a file format for storing biological sequences (usually nucleic acid sequences) and their sequencing quality score information;

[0051] 4. FASTA file: a text-based file format used to represent nucleotide sequences or amino acid sequences;

[0052] 5. File of location information of homologous sequences on the genome: a file storing homologous sequence information, which is generated by performing homology comparison analysis on the sequences of two reference genomes. The stored content includes the location of homologous sequences on different genomes, homology scores, etc.

[0053] 6. BAM file: Binary Alignment / Map, a file format for storing sequence alignment results;

[0054] 7. CIGAR value: Compact Idiosyncratic Gapped Alignment Report, is one of the information recorded in the BAM file, which can reflect the alignment of the sequence, and is composed of English letters and numbers. The English letters represent the alignment status, M represents match, I represents insertion, D represents deletion, N represents variable shearing, S represents substitution, and H represents excision. The number before the English letter represents the number of bases in the corresponding alignment status. For example, 100M50D means that in a sequence, the first 100 bases are matched compared with the reference genome, and the last 50 bases are missing compared with the reference genome;

[0055] 8. QNAME: Query Name, the name of the sequence, is one of the information recorded in the BAM file. There is a one-to-one correspondence between QNAME and sequence. A sequence has only one unique and non-repeating QNAME.

[0056] 9. Alignment information table: a table that records all alignment information of sequences. Each row is a sequence, and stores the alignment position, alignment status, and other information of each sequence on all reference genomes.

[0057] 10. SNV mutation: Single nucleotide variant, refers to a mutation caused by the replacement of a single nucleotide in the genome;

[0058] 10. InDel mutation: insertion-deletion mutation, which refers to the mutation caused by the insertion and deletion of multiple nucleotides in the genome;

[0059] 11. VCF file: Variant Call Format, a file format used to store genomic variation data;

[0060] 12. REF base: the corresponding base on the reference genome at the variant site;

[0061] 13. ALT base: the base that mutates in the sequencing sequence.

[0062] Example 1

[0063] The main purpose of this example is to illustrate that the AsCA process, for the first time, uses a short-fragment sequencing method to identify and analyze the main components of AAV at the computer algorithm level, and its analysis results have higher component quantitative accuracy and lower human operation costs compared to experimental methods. And by comparing the ratio between human genome impurity chromosomes, it can be determined whether the human genome sequence is mixed with the whole genome at a lower frequency or integrated in each rAAV in smaller fragments.

[0064] The main analysis logic steps of AsCA are as follows: Figure 1 shown.

[0065] In this embodiment, the input data that needs to be prepared includes:

[0066] (1) Three AAV samples were used. The cell fluid containing AAV was extracted and sequenced to obtain three pairs of paired-end Fastq files of the samples.

[0067] (2) One standard sample was a mixture of sequences from various sources, including four types of impurity sequences, Helper, R2C, Ecoli, and E1A, as well as viral genome sequences. After sequencing, a pair of double-end Fastq files of the standard sample were obtained;

[0068] (3) One negative control, which contained only one impurity sequence and no viral or host sequences. After sequencing, a pair of impurity double-end Fastq files were obtained;

[0069] (4) Virus genome FASTA file, host genome FASTA file and impurity genome FASTA file. There are four types of impurity sequences in the impurity genome FASTA file, and the impurity names are Helper, R2C, Ecoli and E1A.

[0070] The main analysis logic steps of AsCA are as follows:

[0071] First, use homologous sequence comparison tools, such as NCBI BLAST, EBI hmmer, etc., to perform homology analysis on the viral genome and the host genome to obtain the location information file A of the homologous sequence on the genome; perform homology analysis on the impurity genome and the host genome to obtain the location information file B of the homologous sequence on the genome.

[0072] For each sample's Fastq file, perform the following operations.

[0073] Using alignment software such as BWA, Bowtie2, SOAP, etc., the double-ended Fastq files of the three samples were aligned to the viral genome to obtain the original viral BAM files; the double-ended Fastq files of the three samples were aligned to the host genome to obtain the original host BAM files; the double-ended Fastq files of the three samples were aligned to the impurity genome to obtain the original impurity BAM files. After the alignment analysis was completed, 3 original BAM files were obtained for each sample.

[0074] Read three original BAM files at the same time. When the impurity genome contains multiple impurities, for each type of impurity, AsCA will identify the types and names of all impurities, and merge the alignment information of the same sequence on the three genomes according to the sequence QNAME, obtain the alignment information summary list T, and classify each sequence in the table. The classification judgment logic is as follows:

[0075] (1) In the alignment information summary table T, if a sequence is aligned to only a single genome, the sequence is directly classified into the corresponding genome source, and the QNAME of the sequence is stored in the QNAME list of the corresponding genome source. For example, if a sequence is aligned to and only aligned to a viral genome, its QNAME is stored in the viral QNAME list.

[0076] (2) In the comparison information summary list T, when a sequence is compared to multiple reference genomes, homology information needs to be introduced to achieve accurate classification. Compare the position of the sequence on the reference genome with the position information in the homologous sequence file to determine whether the sequence is in the homologous region. At this point, the sequence will be divided into two parts, one part is "sequences not in the homologous region" and the other part is "sequences in the homologous region".

[0077] (3) In the comparison information summary list T, if a sequence is not in the homologous region and is simultaneously compared to multiple reference genomes, the sequence matching lengths recorded by the CIGAR values ​​of the sequence on different reference genomes are further compared, and the sequence will eventually be classified into the genome source with the longest matching length. For example, if a sequence is simultaneously compared to the viral genome and the host genome, and the CIGAR values ​​of the sequence on the two genomes are 150M and 100M50D, respectively, it means that the matching length of the sequence on the viral genome is 150 bases and the matching length on the host genome is 100 bases. Therefore, the sequence will be classified into the viral genome, and its QNAME will be stored in the viral QNAME list.

[0078] (4) In the comparison information summary list T, if a sequence is in a homologous region and is compared to multiple reference genomes at the same time, it is preferentially considered not to belong to the host sequence, and the judgment logic consistent with that in step (3) is used to classify the sequence.

[0079] After the above analysis is completed, each sample will get 7 QNAME lists, including 1 virus QNAME list, 1 host QNAME list, and 5 impurity QNAME lists (4 known impurities Helper, R2C, Ecoli, E1A and 1 other impurity). Each sequence will be classified into a unique genomic source, so there is no duplication of sequences in each QNAME list.

[0080] Use 7 QNAME lists to filter the 3 original BAM files so that the filtered BAM files contain and only contain sequence information that is accurately classified to the corresponding genome source. If a sequence is classified into the virus QNAME list, the host original BAM file and the impurity BAM file need to filter out the sequence information corresponding to the QNAME, and output the filtering results to the new BAM file, where the number of impurity BAMs changes dynamically with the number of impurity types in the actual input impurity reference genome FASTA file.

[0081] After filtering is completed, 1 virus BAM file, 1 host BAM file and 5 impurity BAM files (including 1 Helper impurity BAM file, 1 R2C impurity BAM file, 1 Ecoli impurity BAM file, 1 E1A impurity BAM file and 1 other impurity BAM file) will be obtained.

[0082] At this point, we start to count and visualize the result data. The specific operations are as follows:

[0083] (1) Accurate quantification of sample composition: For each sample, the seven BAM files obtained after filtering were counted, and the number of sequences stored in each BAM file was counted to obtain the quantitative results of the number of sequences from each reference genome, and a bar graph was drawn ( Figure 2 ), achieving accurate quantification of sample components;

[0084] (2) Sequence distribution information on host genome chromosomes: Count the sequence depth of each chromosome of each sample in the host BAM file, calculate the depth ratio of the sample on each chromosome, use the same method to calculate the depth ratio of the standard on each chromosome, compare the depth ratio of the sample with the standard and perform log processing to obtain the sequence density information of the sample on each chromosome, and draw a scatter plot ( Figure 3 ).from Figure 3 The sequence distribution of all samples on the host genome can be intuitively viewed on the PCR platform. The sequence quantity ratios of all samples and standards on each chromosome on the host genome are relatively consistent.

[0085] pass Figure 3 It can be seen that the host genome chromosome distribution is relatively consistent, and the sequence support number of each chromosome is within one time compared to the standard sample mixed with a certain whole genome, indicating that the host genome is distributed evenly on each chromosome. Therefore, it is determined that the mixing method of the host genome is low-frequency mixing of the whole genome, rather than high-frequency mixing of a certain chromosome fragment on rAAV.

[0086] (3) Base coverage of viral genome: For each sample's viral BAM file, use tools such as samtools to count the sequencing depth of each base position on the viral genome and draw a line graph ( Figure 4 ), which shows the distribution of sequencing depth on the viral genome. This method can cover the entire rAAV virus genome, proving that this method can intuitively and quickly observe and identify areas with abnormal coverage on the viral genome.

[0087] Example 2

[0088] This example aims to verify that this method (AsCA) can detect SNV mutations and InDel-type mutations in the rAAV genome.

[0089] After all samples were analyzed in Example 1, a virus BAM file was obtained for each sample, and the virus BAM file contained and only contained sequence information classified as a viral genome source.

[0090] For each sample's virus BAM file, perform the following operations:

[0091] (1) Use mutation detection tools, such as HaplotypeCaller, VarScan, VarDict, BCFtools, etc., to detect mutations in the sequences in the virus BAM file, and output the SNV mutation and InDel mutation results into a vcf file.

[0092] (2) The vcf file format is organized so that each line contains only one mutation site, and the specific mutation type is determined based on the REF base and ALT base of the mutation site. If the length of the REF base and the ALT base are both 1, the mutation is a SNP mutation; if the REF base length is longer than the ALT base, the mutation is a Deletion; if the REF base length is shorter than the ALT base, the mutation is an Insertion.

[0093] (3) For SNV mutations, count the ALT base types at each mutation site. There are only four types of bases, A, T, G, and C, and count the ALT base mutation frequency at each site. Draw a line graph of the mutation frequency of the four types of bases on the viral genome ( Figure 5 ), the line graph can intuitively reflect the mutation frequency of various bases on the viral genome, show whether the sample mutation has high-frequency SNV mutations on the viral gene, and observe the potential SNV mutation preference of the sample.

[0094] (4) For InDel mutations, the length of each mutation is counted and a mutation diagram of InDel mutations on the viral genome is drawn ( Figure 6 ), the size of each scatter point represents the length of the mutation, the color represents the type of mutation, the line under the scatter point indicates the location of the mutation site on the viral genome, and the height of the line represents the frequency of the mutation. This figure simply and clearly reflects the distribution of sample InDel mutations on the viral genome and the length of the mutation, which is conducive to quickly locating whether the InDel mutation is in an important coding region, and inferring the scope of the affected area from the length of the mutation, and further studying the subsequent protein compilation and gene expression effects.

[0095] Based on the above data, it can be seen that the present method AsCA can detect SNV mutations and InDel-type mutations in the rAAV genome, and additionally provides a visual report to quickly and intuitively reflect the status of mutation detection results, thereby improving the clinical usefulness of rAAV drug testing.

[0096] Embodiment 3:

[0097] Example 3 aims to verify that AsCA can effectively identify and filter multiple alignment sequences caused by homologous sequences in the original alignment result BAM file after the introduction of the homologous sequence classification algorithm.

[0098] In this embodiment, the input data that needs to be prepared includes:

[0099] (1) One sample was manually configured, with only viral sequences and a small amount of impurity sequences added, and no host sequences added. After sequencing, a pair of double-end Fastq files were obtained;

[0100] (2) Viral genome FASTA file, host genome FASTA file.

[0101] First, use alignment software, such as BWA, Bowtie2, SOAP, etc., to align the double-ended Fastq file to the viral genome to obtain the original viral BAM file A; align the double-ended Fastq file to the host genome to obtain the original host BAM file B.

[0102] The original virus BAM file A was counted to obtain the number of sequences that the sample was aligned to the virus genome; the original host BAM file B was counted to obtain the number of sequences that the sample was aligned to the host genome.

[0103] The original virus BAM file and the original host BAM file are processed separately using the homologous sequence classification algorithm to obtain two filtered BAM files, namely virus BAM file A' and host BAM file B'. This filtering step is called false positive removal.

[0104] The number of sequences in the virus BAM file A' is obtained by counting; the number of sequences in the host BAM file B' is obtained by counting.

[0105] Use the number of sequences contained in all BAM files before and after removing false positives to draw a bar graph ( Figure 7 ). As can be seen in the figure, before removing false positives, more comparison sequences can be observed in both the viral genome and the host genome of the sample; after removing false positives, the number of sequences in the viral genome of the sample increased, while the number of sequences in the host genome decreased greatly, which is in line with the sample design.

[0106] In summary, after introducing the homologous sequence classification algorithm, AsCA can effectively identify the source of sequences, improve the situation where one sequence is aligned to multiple genomes due to homologous sequences, and greatly improve the credibility and accuracy of the data.

[0107] The above is a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A classification method for adeno-associated virus data based on second-generation sequencing on multiple reference genomes, characterized in that: The method comprises the following steps: S1. Construct the location information file of homologous sequences on the genome; S2, constructing an original host BAM file, an original virus BAM file, and multiple original impurity BAM files; S3, constructing a comparison information summary list T; S4, build QNAME list; S5, filtering: filtering the three original BAM files in S2 according to the QNAME list to obtain a host BAM file, a virus BAM file, and multiple impurity BAM files; S6. Statistics.

2. The method for classifying adeno-associated virus data based on second-generation sequencing on multiple reference genomes according to claim 1, characterized in that: The method comprises the following steps: S1. Construction of the location information file of homologous sequences on the genome: Perform homology analysis on the FASTA file of the host genome, the FASTA file of the viral genome, and the FASTA file of the impurity genome, respectively, to generate two files of the location information of homologous sequences on the genome, which record the location information of homologous sequences on each genome; S2. Creation of original BAM files: Compare and analyze the double-ended Fastq files with the FASTA files of the host genome, the FASTA files of the virus genome, and the FASTA files of the impurity genome, respectively, to obtain three original BAM files, wherein the double-ended Fastq files are files storing sequence information of viruses, hosts, and impurities obtained after sequencing; S3, construction of alignment information summary list T: read the alignment information in the three original BAM files in S2, and merge all alignment information of a sequence together according to QNAME to obtain the alignment information summary list T; S4, construction of QNAME list: Classify each sequence in the comparison information summary list T in S3, and the classification judgment logic includes sequences aligned to a single genome and sequences aligned to multiple reference genome sources to determine the QNAME list of the virus, the QNAME list of the host, and one of the QNAME lists of multiple impurities, and obtain the QNAME list of the virus, the QNAME list of the host, and the QNAME list of multiple impurities; at this point, each sequence has and only has a unique reference genome source, and the sequence QNAME of each source is recorded separately in the QNAME list of the corresponding source; S5, filtering: filtering the three original BAM files in S2 according to the QNAME list of the virus, the QNAME list of the host, and the QNAME list of multiple impurities to obtain a host BAM file, a virus BAM file, and multiple impurity BAM files, and each sequence appears in and only appears in one BAM file; S6. Statistics: Perform statistics on the result data for the multiple BAM files obtained in S5.

3. The method for classifying adeno-associated virus data based on second-generation sequencing on multiple reference genomes according to claim 2, characterized in that: In S4, in the alignment information summary list T, for sequences that are aligned to only a single genome, they are directly classified into the genome source and their QNAME is recorded in the corresponding QNAME list.

4. The method for classifying adeno-associated virus data based on second-generation sequencing on multiple reference genomes according to claim 2, characterized in that: In S4, in the comparison information summary list T, for sequences that are compared to multiple reference genome sources, whether they are in the homologous region is determined based on the position information, and the sequence is divided into two parts, one part is "sequences that are not in the homologous region" and the other part is "sequences that are in the homologous region".

5. The method for classifying adeno-associated virus data based on second-generation sequencing on multiple reference genomes according to claim 4, characterized in that: In the alignment information summary list T, for sequences that are not in the homologous region and are aligned to multiple reference genome sources, the sequence matching lengths recorded by the CIGAR values ​​of the sequence on different reference genomes are compared, and the sequence is classified into the genome source with the longest matching length, and its QNAME is recorded in the corresponding QNAME list.

6. The method for classifying adeno-associated virus data based on second-generation sequencing on multiple reference genomes according to claim 4, characterized in that: In the comparison information summary list T, for sequences in the homologous region that are compared to multiple reference genome sources, they are preferentially considered not to be host sequences. The sequence matching lengths recorded by the CIGAR values ​​of the sequence on different reference genomes are compared, and the sequence is classified into the genome source with the longest matching length, and its QNAME is recorded in the corresponding QNAME list.

7. The method for classifying adeno-associated virus data based on second-generation sequencing on multiple reference genomes according to claim 4, characterized in that: The filtering steps of S5 are as follows: the three original BAM files obtained in S1 are processed according to the QNAME list. If a sequence appears in the QNAME list of the virus, the corresponding sequences in the host original BAM file and the impurity original BAM file are filtered out; after filtering, the host BAM file, the virus BAM file and multiple impurity BAM files are obtained, where the number of impurity BAM files changes dynamically with the number of impurity types during actual analysis; so far, each sequence appears and only appears in one BAM file.

8. The method for classifying adeno-associated virus data based on second-generation sequencing on multiple reference genomes according to claim 1 or 2, characterized in that: The statistics in step S6 include the following information: 1) Count the number of sequences stored in each BAM file of a sample to obtain the quantitative results of the sample's components; 2) Count the sequence position information in the virus BAM file to obtain the sequence coverage on the virus genome; 3) Count the sequence depth of each chromosome in the host BAM file, calculate the depth ratio of the sample on each chromosome, use the same method to calculate the depth ratio of the standard on each chromosome, compare the depth ratio of the sample with the standard and perform log processing to obtain the sequence density information of the sample on each chromosome.

9. The method for classifying adeno-associated virus data based on second-generation sequencing on multiple reference genomes according to claim 8, characterized in that: The statistical results in step S6 are collated and visualized, and the component identification results and visualization report of the AAV vector data are output.

10. The method for classifying adeno-associated virus data based on second-generation sequencing on multiple reference genomes according to any one of claims 1 to 7, characterized in that: It also includes S7, mutation detection: for the virus BAM file obtained in S5, mutation detection is performed using a mutation detection tool, including SNV detection and InDel detection, to obtain a mutation detection result table in a VCF file format, and to visualize the SNV base ratio and the InDel distribution on the genome.

Images