Sequence typing method based on BLAST alignment

Abstract

The invention relates to a sequence typing method based on BLAST alignment, and the method specifically comprises the following steps: step 1, preparing an FASTA format file containing a plurality of known HIV subtype reference sequences, each sequence corresponding to a known subtype, the HIV subtypes including but not limited to CRF07BC, CRF08BC, B and CRF01AE; step 2, establishing a local nucleic acid database for the reference sequence by using a makeblastdb tool in BLAST, wherein the database is used for subsequent sequence alignment; 3, target sequences to be classified are input and subjected to BLAST comparison with the sequences in the reference database, a large number of HIV sequences can be automatically classified, the workload of manual analysis is omitted, and the method is suitable for epidemiological monitoring and gene sequence analysis.

Description

Technical Field

[0001] The present invention belongs to the technical field of bioinformatics, and in particular, relates to a sequence typing method based on BLAST comparison. Background Art

[0002] With the rapid development of high-throughput sequencing technology, a large amount of genomic and transcriptomic data has been accumulated in the field of bioinformatics. In infectious disease research, viral gene sequence information is crucial for monitoring, tracking and understanding viral transmission. For example, the genetic diversity of human immunodeficiency virus (HIV) makes studying its subtype distribution and evolutionary path a complex but necessary task. Different HIV subtypes are distributed differently around the world, and they may be closely related to disease progression, treatment efficacy and drug resistance. Therefore, accurate HIV subtype classification is crucial for epidemiological research and disease prevention and control.

[0003] Traditional HIV subtype classification methods mostly rely on manual analysis or the use of specific software tools, which are inefficient when processing large-scale data and may have difficulty maintaining high accuracy for different subtypes. BLAST (Basic Local Alignment Search Tool), as a commonly used sequence alignment tool, has the characteristics of fast speed and high accuracy and has been widely used in the alignment and analysis of genome and protein sequences. However, when using BLAST for sequence classification in the prior art, there is often a lack of a systematic automated process, especially when processing large-scale sequence data, and it is difficult to effectively classify sequences with low similarity or partial alignment coverage. Summary of the invention

[0004] The purpose of the present invention is to provide a sequence typing method based on BLAST comparison, which solves the problems existing in the prior art.

[0005] The purpose of the present invention can be achieved through the following technical solutions: A sequence typing method based on BLAST comparison, characterized in that it specifically comprises the following steps: Step 1, prepare a FASTA format file containing multiple known HIV subtype reference sequences, each sequence corresponding to a known subtype, and the HIV subtypes include but are not limited to CRF07_BC, CRF08_BC, B, CRF01_AE; Step 2: Use the makeblastdb tool in BLAST to establish a local nucleic acid database for the reference sequence, and the database is used for subsequent sequence alignment; Step 3: Input the target sequence to be classified and perform BLAST comparison with the sequences in the reference database. The comparison uses the following output format: sequence ID, comparison sequence ID, similarity, comparison length, number of mismatches, number of gaps, query sequence start and end positions, reference sequence start and end positions, expected value and score; Step 4: Calculate the similarity and coverage of each target sequence with the reference sequence based on the comparison results. The similarity is the percentage of the matched base pairs in the comparison to the total length, and the coverage is the percentage of the comparison length to the total length of the reference sequence. Step 5: classify the target sequence according to the preset similarity threshold and coverage threshold; Step 6. When the similarity and coverage of the target sequence are greater than or equal to the preset similarity threshold and coverage threshold, the target sequence is classified as the subtype corresponding to the reference sequence with the highest similarity; otherwise, the target sequence is classified as the "other" type.

[0006] Preferably, the FASTA format file in step 1 is as follows: >CRF07_BC ATGGAAGAGAAAT... >CRF08_BC CTGAGATGGGATG... >B GCTAGCTTGACTG... >CRF01_AE TACGATGACTCGT... FASTA contains reference sequences for multiple HIV subtypes, and the sequence of each subtype is consistent with its subtype name.

[0007] Preferably, in step 2, the command for using the makeblastdb tool in BLAST to establish a local nucleic acid database for the reference sequence is as follows: makeblastdb -in sequences.fasta -dbtype nucl -out hiv_subtype_db2 -in: input reference sequence file sequences.fasta; -dbtype: specifies the nucleic acid database nucl; -out: The database file prefix is ​​hiv_subtype_db2; This step generates a series of files prefixed with hiv_subtype_db2, including .nhr, .nin, and .nsq, which together constitute the BLAST database of HIV subtype reference sequences.

[0008] Preferably, the command for performing BLAST comparison with the sequence in the reference database in step 3 is as follows: blastn -query temp_fasta_file -db hiv_subtype_db2 -out result_file -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstartsend evalue bitscore" -query: target sequence file temp_fasta_file; -db: reference database hiv_subtype_db2; -out: output comparison result file result_file; -outfmt: Set the output format to a tab-delimited table.

[0009] Preferably, the similarity calculation in step 5 is as follows: obtaining the similarity from the pident field of the BLAST result, which indicates the percentage of bases matching between the target sequence and the reference sequence; Coverage calculation: The coverage is calculated based on the alignment length in the BLAST result and the total length of the reference sequence. The calculation formula is: Coverage = (total length of reference sequence / alignment length) × 100 Among them, the total length of the reference sequence can be calculated from the send and sstart fields of the BLAST result; The target sequence classification criteria are as follows: the similarity threshold is set to 95% and the coverage threshold is set to 95%.

[0010] Beneficial effects of the present invention: Automation and efficiency: This method can automatically classify a large number of HIV sequences, eliminating the workload of manual analysis, and is suitable for epidemiological monitoring and gene sequence analysis. High accuracy: By setting similarity and coverage thresholds, the misclassification problem caused by incomplete alignment of some sequences is effectively avoided. Scalability and flexibility: This method is not limited to HIV sequence classification, but can also be applied to the genetic typing of other viruses or organisms, and the database can be dynamically updated according to newly discovered subtypes. Improved processing speed: Through multi-threaded parallel processing, the data processing speed is greatly improved, which is suitable for the comparison and analysis of large-scale high-throughput data. DETAILED DESCRIPTION

[0011] The following will be combined with the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0012] The present invention is a sequence typing method based on BLAST comparison, which is characterized by comprising the following steps: Step 1, prepare a FASTA format file containing multiple known HIV subtype reference sequences, each sequence corresponding to a known subtype, and the HIV subtypes include but are not limited to CRF07_BC, CRF08_BC, B, CRF01_AE; The FASTA format file is as follows: >CRF07_BC ATGGAAGAGAAAT... >CRF08_BC CTGAGATGGGATG... >B GCTAGCTTGACTG... >CRF01_AE TACGATGACTCGT... The FASTA file is named sequences.fasta. FASTA contains reference sequences of multiple HIV subtypes. The sequence of each subtype is consistent with its subtype name, which facilitates subsequent comparison and analysis.

[0013] Step 2: Use the makeblastdb tool in BLAST to establish a local nucleic acid database for the reference sequence, and the database is used for subsequent sequence alignment; The command to establish a local nucleic acid database for the reference sequence using the makeblastdb tool in BLAST is as follows: makeblastdb -in sequences.fasta -dbtype nucl -out hiv_subtype_db2 -in: input reference sequence file sequences.fasta; -dbtype: specifies the nucleic acid database nucl; -out: The database file prefix is ​​hiv_subtype_db2; This step generates a series of files prefixed with hiv_subtype_db2, including .nhr, .nin, and .nsq, which together constitute the BLAST database of HIV subtype reference sequences.

[0014] Step 3: Input the target sequence to be classified and perform BLAST comparison with the sequences in the reference database. The comparison uses the following output format: sequence ID, comparison sequence ID, similarity, comparison length, number of mismatches, number of gaps, query sequence start and end positions, reference sequence start and end positions, expected value and score; The command for BLAST comparison with the sequence in the reference database is as follows: blastn -query temp_fasta_file -db hiv_subtype_db2 -out result_file -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstartsend evalue bitscore" -query: target sequence file temp_fasta_file; -db: reference database hiv_subtype_db2; -out: output comparison result file result_file; -outfmt: Set the output format to a tab-delimited table. Includes the following fields: qseqid: query sequence ID (i.e., the sequence ID to be classified).

[0015] sseqid: the reference sequence ID for the alignment (i.e., the subtype category of the reference sequence).

[0016] pident: Alignment similarity, indicating the proportion of matching bases.

[0017] length: comparison length.

[0018] mismatch: number of mismatched bases.

[0019] gapopen: number of gaps.

[0020] qstart and qend: the start and end positions of the query sequence alignment.

[0021] sstart and send: the start and end positions of the reference sequence.

[0022] evalue: expected value, used to measure the significance of the comparison.

[0023] bitscore: comparison score, reflecting the matching quality of the comparison.

[0024] Step 4: Calculate the similarity and coverage of each target sequence with the reference sequence based on the comparison results. The similarity is the percentage of the matched base pairs in the comparison to the total length, and the coverage is the percentage of the comparison length to the total length of the reference sequence. Step 5: classify the target sequence according to the preset similarity threshold and coverage threshold; Read the BLAST comparison results. The generated BLAST comparison results are saved in the result_file file, and each line represents a comparison record. The format is as follows: Sample_001 CRF07_BC 97.5 900 22 1 1 900 1 900 0.0 1600 Sample_002 B 95.2 850 40 3 50 900 1 850 0.0 1500 The first column of each row indicates the ID of the target sequence, the second column indicates the reference sequence (i.e., isoform) on which the alignment is made, the third column indicates the similarity, and the fourth column indicates the length of the alignment.

[0025] The similarity was calculated as follows: The similarity was obtained from the pident field of the BLAST result, which represents the percentage of bases that matched the target sequence with the reference sequence; Coverage calculation: The coverage is calculated based on the alignment length in the BLAST result and the total length of the reference sequence. The calculation formula is: Coverage = (total length of reference sequence / alignment length) × 100 Among them, the total length of the reference sequence can be calculated from the send and sstart fields of the BLAST result; The target sequence classification criteria are as follows: the similarity threshold is set to 95% and the coverage threshold is set to 95%.

[0026] Step 6. When the similarity and coverage of the target sequence are greater than or equal to the preset similarity threshold and coverage threshold, the target sequence is classified as the subtype corresponding to the reference sequence with the highest similarity; otherwise, the target sequence is classified as the "other" type.

[0027] The generation of classification results uses the following algorithm to extract the classification results and generate the final output: import pandas as pd def get_subtype_list(results_file): SIMILARITY_THRESHOLD = 95.0 COVERAGE_THRESHOLD = 95.0 data = [] with open(results_file, 'r') as file: for line in file: parts = line.strip().split('\t') sample_id = parts[0] subtype = parts[1] similarity = float(parts[2]) alignment_length = int(parts[3]) total_length = int(parts[8]) - int(parts[7]) + 1 coverage = (alignment_length / total_length) * 100 data.append([sample_id, subtype, similarity, coverage]) df = pd.DataFrame(data, columns=['SampleID', 'Subtype', 'Similarity','Coverage']) result_list = [] for sample_id in df['SampleID'].unique(): sample_data = df[df['SampleID'] == sample_id] max_similarity = sample_data['Similarity'].max() max_coverage = sample_data['Coverage'].max() if max_similarity <SIMILARITY_THRESHOLD or max_coverage<COVERAGE_THRESHOLD: result_list.append({"id": sample_id, "result": "Other"}) else: dominant_subtype = sample_data.loc[sample_data['Similarity'].idxmax(), 'Subtype'] result_list.append({"id": sample_id, "result": dominant_subtype}) return result_list Based on the BLAST alignment results, the algorithm classifies each target sequence according to the preset similarity and coverage thresholds to generate the final classification results.

[0028] In order to improve the efficiency of processing large-scale sequence data, the present invention introduces multi-threaded parallel processing in the BLAST comparison process. By adding the -num_threads parameter in the blastn command, multiple CPU cores can be specified for parallel computing. For example, the command for comparison using 4 threads is as follows: blastn -query temp_fasta_file -db hiv_subtype_db2 -out result_file -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstartsend evalue bitscore" -num_threads 4 After using multithreading, the alignment time is significantly reduced, especially when dealing with thousands of target sequences.

[0029] The database of the present invention can be updated according to the newly discovered HIV subtype. By adding new subtype reference sequences and re-running the makeblastdb command, the reference database can be dynamically expanded. For example, when a new subtype CRF09_CX is found, the sequence of the subtype can be added to the FASTA file and the database can be re-established.

[0030] makeblastdb -in updated_sequences.fasta -dbtype nucl -out updated_hiv_subtype_db.

[0031] The updated database will contain the latest HIV subtype information for subsequent comparison and classification.

[0032] In the description of the specification, the description with reference to the terms "one embodiment", "example", "specific example", etc. means that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representation of the above terms does not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described can be combined in any one or more embodiments or examples in a suitable manner.

[0033] The above contents are merely examples and explanations of the concept of the present invention. The technicians in this technical field may make various modifications or additions to the specific embodiments described or replace them in a similar manner. As long as they do not deviate from the concept of the invention or exceed the scope defined by the claims, they should all fall within the protection scope of the present invention.

Claims

1. A sequence typing method based on BLAST comparison, characterized in that: The specific steps include: Step 1, prepare a FASTA format file containing multiple known HIV subtype reference sequences, each sequence corresponding to a known subtype, and the HIV subtypes include but are not limited to CRF07_BC, CRF08_BC, B, CRF01_AE; Step 2: Use the makeblastdb tool in BLAST to establish a local nucleic acid database for the reference sequence, and the database is used for subsequent sequence alignment; Step 3: Input the target sequence to be classified and perform BLAST comparison with the sequences in the database. The comparison uses the following output format: sequence ID, comparison sequence ID, similarity, comparison length, number of mismatches, number of gaps, query sequence start and end positions, reference sequence start and end positions, expected value and score; Step 4: Calculate the similarity and coverage of each target sequence with the reference sequence based on the comparison results. The similarity is the percentage of the matched base pairs in the comparison to the total length, and the coverage is the percentage of the comparison length to the total length of the reference sequence. Step 5: classify the target sequence according to the preset similarity threshold and coverage threshold; Step 6: When the similarity and coverage of the target sequence are greater than or equal to the preset similarity threshold and coverage threshold, the target sequence is classified as the subtype corresponding to the reference sequence with the highest similarity; otherwise, the target sequence is classified as "other" type; The FASTA format file in step 1 is as follows: >CRF07_BC ATGGAAGAGAAAT... >CRF08_BC CTGAGATGGGATG... >B GCTAGCTTGACTG... >CRF01_AE TACGATGACTCGT... FASTA contains reference sequences for multiple HIV subtypes, and the sequence of each subtype is consistent with its subtype name; In step 2, the command for using the makeblastdb tool in BLAST to establish a local nucleic acid database for the reference sequence is as follows: makeblastdb -in sequences.fasta -dbtype nucl -out hiv_subtype_db2 -in: input reference sequence file sequences.fasta; -dbtype: specifies the nucleic acid database nucl; -out: The database file prefix is ​​hiv_subtype_db2; This step generates a series of files prefixed with hiv_subtype_db2, including .nhr, .nin, and .nsq, which together constitute the BLAST database of HIV subtype reference sequences; The command for BLAST comparison with the sequence in the reference database in step 3 is as follows: blastn -query temp_fasta_file -db hiv_subtype_db2 -out result_file -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstartsend evalue bitscore" -query: target sequence file temp_fasta_file; -db: reference database hiv_subtype_db2; -out: output comparison result file result_file; -outfmt: Set the output format to a tab-delimited table.

2. A sequence typing method based on BLAST comparison according to claim 1, characterized in that: The similarity calculation in step 5 is as follows: The similarity is obtained from the pident field of the BLAST result, which represents the percentage of bases that match the target sequence with the reference sequence; Coverage calculation: The coverage is calculated based on the alignment length in the BLAST result and the total length of the reference sequence. The calculation formula is: Coverage = (total length of reference sequence / alignment length) × 100 Among them, the total length of the reference sequence can be calculated from the send and sstart fields of the BLAST result; The target sequence classification criteria are as follows: the similarity threshold is set to 95% and the coverage threshold is set to 95%.