HLA typing method based on high-throughput sequencing

Abstract

The invention belongs to the technical field of bioinformatics, and particularly relates to an HLA typing method based on high-throughput sequencing. According to the method, the aligned HLA sequence is innovatively processed to form an algorithm database with the same length as the sequencing read length, and each basic group of each allele is placed in the absolute coordinate of the whole gene locus through the aligned HLA, so that the allelic gene is obtained. Therefore, the sequence which is obtained after processing and has the same length as the sequencing read length directly has absolute coordinates of one or more HLA gene loci. Meanwhile, high-throughput sequencing data is directly compared with an algorithm database, the base coverage of each absolute coordinate can be directly obtained on the locus, and then the loci are subjected to the steps of locus coincident fragment disassembly, position coincident base disassembly, haploid phasing, candidate sequence construction, database comparison and the like. And obtaining one or more pairs of clear HLA typing sequences according to the clearness of the original data and the mutation position density of the gene locus.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of bioinformatics, and particularly relates to an HLA typing method based on high-throughput sequencing. BACKGROUND

[0002] Human Leukocyte Antigen (HLA) is the expression product of human major histocompatibility complex (MHC), a group of closely linked genes located on the short arm of chromosome 6. It plays a crucial role in the immune system. The HLA system is very complex and highly polymorphic, which makes the HLA of different individuals vary greatly, and determines the immune response ability of the body to pathogens and other foreign antigens and the important characteristics of HLA molecules in medical scenarios such as organ transplantation.

[0003] HLA is usually divided into three categories, I, II and III genes. HLA-I: almost distributed on the surface of all nucleated cells, the encoded product is mainly responsible for activating CD8 + T cells, including classic antigen genes such as HLA-A, HLA-B and HLA-C. HLA-II: mainly distributed on the surface of professional antigen presenting cells (APC), including HLA-DR, HLA-DQ and HLA-DP genes; HLA-III: mainly the complement of antigens, involved in immune regulation and inflammatory response processes.

[0004] In the current HLA typing method based on high-throughput sequencing technology, some methods only use a simple scoring system, lack of verification mechanism, resulting in sometimes incorrect typing results. For the type not recorded in the database, no prompt is given, but an incorrect typing result is given; in addition, although the same principle of aligning HLA is applied, the database is processed for the entire HLA sequence, and when compared, an algorithm allowing mismatch such as Bowtie2 is used for comparison, and the sequence needs to be compared twice, the running speed is slow, and the cross-locus fragments are not carefully processed, which can easily produce errors and affect the accuracy of typing. SUMMARY

[0005] The purpose of the present application is to provide an HLA typing method based on high-throughput sequencing, which can construct an algorithm database with the same sequence length as the sequencing read length, based on which the present application can place each base of each allele at the absolute coordinates of the entire locus, and then obtain one or several clear HLA typing sequences, with very high accuracy.

[0006] The present application provides a method for constructing an HLA typing algorithm database, comprising the following steps: 1) Download the IMGT / HLA alignment database from the International Immunogenetic Database (IMGT), where each locus in the IMGT / HLA alignment database includes alleles with known full sequences and alleles with known partial exon sequences; 2) Process the alignment database of the IMGT / HLA to correct alignment errors; 3) After completing the processing described in step 2), for each locus, complete the partial exon sequence with known alleles using the known alleles of the full sequence. 4) Cut the alleles in each locus according to the high-throughput sequencing read length with a step size of 1, and remove duplicates to obtain the database for each locus; 5) Compare all locus databases. When duplicate base sequences are found in two or more locus databases, delete the duplicate base sequences and classify them into a joint database to obtain the HLA typing algorithm database. The HLA typing algorithm database includes a database for each locus, as shown in one of the following two categories: I: Independent database of gene loci; II: Independent databases of loci and joint databases containing the loci.

[0007] Preferably, the method for completing step 3) includes: using the intron sequences of the alleles with known full-sequence sequences to complete both ends of the alleles with known partial exon sequences.

[0008] Preferably, the method for selecting the allele with a known full sequence includes: selecting the allele with a known full sequence that is the same as or closest to the exon sequence of the partial exon sequence.

[0009] Preferably, after step 3), the method further includes deleting the absolute positions of the unsequencing and deleted portions shown in the alleles of each locus, retaining only the sequence of bases and the absolute positions of the bases.

[0010] The present invention also provides the application of the construction method described in the above technical solution or the HLA typing algorithm database constructed using the construction method described in the above technical solution in HLA typing.

[0011] This invention also provides an HLA typing method, characterized by comprising the following steps: S1) The HLA typing algorithm database is constructed using the construction method described in the above technical solution; S2) The high-throughput sequencing data after quality inspection is compared with the HLA typing algorithm database for typing; S21) When performing genotyping for a specific locus, the high-throughput sequencing data after quality control is compared with both independent databases and joint databases containing the locus, yielding the following comparison results: S211) When both sequences in a fragment pair are aligned to independent databases, the fragment pair is confirmed to originate from the locus, and the base coverage of the locus is determined based on the alignment information. S212) When one sequence in a fragment pair aligns to an independent database and the other sequence aligns to a joint database, the fragment pair is confirmed to originate from the locus. S213) When both fragments in a fragment pair are aligned to a joint database, the locus of the fragment pair cannot be identified. S22) After aligning the S211) and S212) fragments to specific loci, the coverage and base types at each position are statistically analyzed. S221) If the coverage at any location is less than 10% of the average coverage depth of the exon or intron region it is located in, mark the location as a low-coverage location; if the low-coverage location is located in an exon, report an error and stop typing. S222) If the abundance of the third most abundant base in the coverage data at any location exceeds 5%, microsatellite detection is performed at that location; if no microsatellite exists at that location, an error is reported and typing is stopped; if a microsatellite exists at that location, the two most abundant bases at that location are recorded. S223) If no error is triggered at any position as in S241) and S242), count the dibase positions at those positions; S224) Bases that can be placed in multiple positions are temporarily omitted and recorded separately; S23) The coverage data pairs consisting of the fragment pairs in S213) and the fragment pairs with known loci are compared position by position. If the fragment pair is different from the known coverage data of the locus at any position and the coverage is less than 10%, it is confirmed that the fragment pair does not come from the locus. When there is only one candidate locus left for the fragment pair, it is confirmed that the fragment pair comes from the locus. After confirmation, the fragment pair should be added to the coverage data of the correct locus. S24) For the bases (wobbly bases) that can be placed in multiple positions encountered in S224), compare them with the coverage data of these positions: if it is confirmed that there are bases in the coverage data that do not match the bases that can be placed in multiple positions, delete the position from the candidate positions of the bases that can be placed in multiple positions; when there is only one candidate position left for the bases that can be placed in multiple positions, confirm that this base is located in the remaining position; S25) Haploid phasing and polyploid phenomena were analyzed for each pair of adjacent two dibase positions; S26) Based on the haploid phasing analysis described in S25), one or more pairs of complete sequences are generated for each locus; S27) Based on the complete sequence described in S26), the typing result is confirmed by comparing it with the HLA typing algorithm database.

[0012] Preferably, step S25) haploid phasing includes the following steps: for each pair of adjacent two dibase positions, extract all fragment pairs that simultaneously cover the two dibase positions; The number of fragment pairs supporting cis and trans is counted separately. When the number of fragment pairs supporting cis or trans exceeds 90% of the total, the two positional relationships are marked as cis or trans. When the number of fragment pairs supporting cis or trans is less than 90% of the total, microsatellite detection is performed on the fragment pairs. If microsatellites exist, the dibase position is recorded as an unknown haploid phase determination result. If microsatellites do not exist, an error is reported and typing is stopped.

[0013] Preferably, the method for generating the complete sequence in step S26) includes the following steps: If no dibase position is detected in step S22), a complete sequence is generated; If a double base position is detected, a pair of sequences is generated. When a double base position is encountered during the generation process, the two bases are assigned to the two sequences according to the haploid phasing result. If the haploid phasing result of the double base position encountered during the generation process is unknown compared to the previous pair of double base positions, the sequence pair is doubled, and the double base position is assigned to the doubled sequence pair in cis and trans configurations, respectively.

[0014] Preferably, the comparison method in step S27) includes the following steps: comparing the complete sequence with the HLA typing algorithm database in the order of HLA type I genes, HLA type II genes, and HLA type III genes.

[0015] The present invention also provides a storage medium or processor, the storage medium including a stored program, the stored program executing the construction method as described in the above technical solution or the HLA typing method as described in the above technical solution; The processor is used to run a program, which executes the construction method or the HLA typing method described in the above technical solution.

[0016] Beneficial effects: This invention innovatively processes aligned HLA sequences to form an algorithmic database of the same length as the sequencing reads. After alignment, each base of each allele is placed at the absolute coordinates of the entire locus. Therefore, the processed sequences, identical in length to the sequencing reads, directly carry the absolute coordinates of one or more HLA loci. Simultaneously, this invention directly compares high-throughput sequencing data with the algorithmic database, directly determining the base coverage at each absolute coordinate at the locus. Following steps such as locus overlap fragment decomposition, positional overlap base decomposition, haploid phasing, candidate sequence construction, and database comparison, depending on the clarity of the original data and the mutation density at the locus, one or more pairs of clear HLA genotyping sequences are obtained. It has the following advantages: 1. It can clearly identify the vast majority of overlapping segments at the loci, providing clear typing and ensuring accuracy; 2. Able to identify new subtypes and provide suggestions, and able to identify the closest allele between the new subtype and the previously reported subtype, as well as the differences between the new allele and the new allele; 3. It can clearly indicate the presence of contamination, chimerism, or copy number variation, i.e., polyploidy, and has a high sensitivity in indicating this. 4. When the two alleles at the sample locus do not contain microsatellites, the results can be accurate to the quartile, i.e., allele precision. 5. The method uses direct comparison in the process of comparing data and algorithm database, which simplifies the comparison difficulty, increases the running speed, and significantly reduces the program running time; 6. Partially resolves the issue of incomplete allele sequences in the IMGT database. Attached Figure Description

[0017] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the accompanying drawings used in the embodiments will be briefly described below.

[0018] Figure 1 This is the aligned database structure of IMGT / HLA downloaded from the IMGT database in Example 1; Figure 2 This is a diagram illustrating alignment errors in the IMGT / HLA alignment database in Example 1; Figure 3 This is a diagram showing the sequence after completion during the algorithm database construction process in Example 1; Figure 4 This is a diagram illustrating the joint database in Example 1; Figure 5 This is a diagram showing the results of comparing the quality-inspected data with the algorithm database in Example 1; Figure 6This is a diagram illustrating the matching of databases containing more than one type of absolute coordinate in Example 1. Figure 7 This is a schematic diagram of the initial absolute coverage of HLA-B loci after statistical analysis in Example 1; Figure 8 This is a diagram showing the location area coverage data in Example 1; Figure 9 This is a diagram showing the phase determination results of haploids in Example 1; Figure 10 This is a diagram illustrating the haploid sequence to be compared generated based on the coverage data and the dibase positions of the haploid phase in Example 1. Detailed Implementation

[0019] This invention provides a method for constructing an HLA typing algorithm database, comprising the following steps: 1) Download the IMGT / HLA alignment database from the International Immunogenetic Database (IMGT), where each locus in the IMGT / HLA alignment database includes alleles with known full sequences and alleles with known partial exon sequences; 2) Process the alignment database of the IMGT / HLA to correct alignment errors; 3) After completing the processing described in step 2), for each locus, complete the partial exon sequence with known alleles using the known alleles of the full sequence. 4) Cut the alleles in each locus according to the high-throughput sequencing read length with a step size of 1, and remove duplicates to obtain the database for each locus; 5) Compare all locus databases. When duplicate base sequences are found in two or more locus databases, delete the duplicate base sequences and classify them into a joint database to obtain the HLA typing algorithm database. The HLA typing algorithm database includes a database for each locus, as shown in one of the following two categories: I: Independent database of gene loci; II: Independent databases of loci and joint databases containing the loci.

[0020] As one implementation, in step 2), the correction of comparison errors includes processing unnecessary insertions and deletions. The processing method is not particularly limited and can be performed according to conventional methods in the art.

[0021] In one implementation, the method for completing step 3) includes: filling in the gaps between the ends of the allele with a known partial exon sequence using the intron sequence of the allele with a known full sequence. In another implementation, the method for selecting the allele with a known full sequence includes: selecting the allele with a known full sequence that is identical to or closest to the exon sequence of the partial exon sequence.

[0022] In one implementation, after step 3), the absolute positions of the unsequencing and deleted portions shown in the alleles of each locus are deleted, leaving only the sequence of bases and the absolute positions of the bases; in the IMGT / HLA alignment database, the unsequencing portion is represented by "*" and the deleted portion is represented by "·".

[0023] As one implementation, the high-throughput sequencing read length in step 4) can be 150 nt, 250 nt, or other common sequencing read lengths in the field of high-throughput sequencing. In the embodiments of the present invention, a sequencing method with a sequencing strategy of PE150 and a sequencing read length of 150 nt is used as an example for illustration, but it should not be limited to the entire scope of protection of the present invention. As one implementation, the deduplication method in step 4) is as follows: when two or more sequencing reads have the same base sequence length but different absolute coordinate positions, one entry is retained, and the two or more absolute coordinate positions are recorded in the entry.

[0024] The present invention also provides the application of the construction method described in the above technical solution or the HLA typing algorithm database constructed using the construction method described in the above technical solution in HLA typing.

[0025] This invention also provides an HLA typing method, characterized by comprising the following steps: S1) The HLA typing algorithm database is constructed using the construction method described in the above technical solution; S2) The high-throughput sequencing data after quality inspection is compared with the HLA typing algorithm database for typing; S21) When performing genotyping for a specific locus, the high-throughput sequencing data after quality control is compared with both independent databases and joint databases containing the locus, yielding the following comparison results: S211) When both sequences in a fragment pair are aligned to independent databases, the fragment pair is confirmed to originate from the locus. S212) When one sequence in a fragment pair aligns to an independent database and the other sequence aligns to a joint database, the fragment pair is confirmed to originate from the locus. S213) When both fragments in a fragment pair are aligned to a joint database, the locus of the fragment pair cannot be identified. S22) After aligning the S211) and S212) fragments to specific loci, the coverage and base types at each position are statistically analyzed. S221) If the coverage at any location is less than 10% of the average coverage depth of the exon or intron region it is located in, mark the location as a low-coverage location; if the low-coverage location is located in an exon, report an error and stop typing. S222) If the abundance of the third most abundant base in the coverage data at any location exceeds 5%, microsatellite detection is performed at that location; if no microsatellite exists at that location, an error is reported and typing is stopped; if a microsatellite exists at that location, the two most abundant bases at that location are recorded. S223) If no error is triggered at any position as in S241) and S242), count the dibase positions at those positions; S224) Bases that can be placed in multiple positions are temporarily omitted and recorded separately; S23) The coverage data pairs consisting of the fragment pairs in S213) and the fragment pairs with known loci are compared position by position. If the fragment pair is different from the known coverage data of the locus at any position and the coverage is less than 10%, it is confirmed that the fragment pair does not come from the locus. When there is only one candidate locus left for the fragment pair, it is confirmed that the fragment pair comes from the locus. After confirmation, the fragment pair should be added to the coverage data of the correct locus. S24) For the bases (wobbly bases) that can be placed in multiple positions encountered in S224), compare them with the coverage data of these positions: if it is confirmed that there are bases in the coverage data that do not match the bases that can be placed in multiple positions, delete the position from the candidate positions of the bases that can be placed in multiple positions; when there is only one candidate position left for the bases that can be placed in multiple positions, confirm that this base is located in the remaining position; S25) Haploid phasing and polyploid phenomena were analyzed for each pair of adjacent two dibase positions; S26) Based on the haploid phasing analysis described in S25), one or more pairs of complete sequences are generated for each locus; S27) Based on the complete sequence described in S26), the typing result is confirmed by comparing it with the HLA typing algorithm database.

[0026] In one implementation, the high-throughput sequencing data in step S2) is next-generation sequencing high-throughput data. In one implementation, the high-throughput sequencing library is an HLA gene-enriched library; in another implementation, the enrichment method for the HLA gene-enriched library can be PCR amplification or probe fishing, but this should not be limited to the entire scope of protection of this invention. In one implementation, the quality control method in step 2) includes removing all low-quality fragment pairs containing any bases with a confidence level below 95%. In one implementation, the alignment method in step S2) is direct alignment, and no mismatches of any kind are allowed.

[0027] As one implementation method, the standard for microsatellites in step S222) is that a single base appears 7 or more times, or a double base appears 5 or more times. As one implementation method, after step S222), if the ratio of the most abundant base to the second most abundant base in the data coverage at any position is less than 6, the position is recorded as a double base position; at the same time, this position is recorded as an unbalanced position, and the library construction scheme is manually checked and adjusted based on this.

[0028] As one implementation method, step S25) of haploid phasing includes the following steps: for each pair of adjacent dibasic positions, extract all fragment pairs that simultaneously cover the two dibasic positions; count the number of fragment pairs supporting cis and trans respectively; when the number of fragment pairs supporting cis or trans exceeds 90% of the total, mark the relationship between the two positions as cis or trans; when the number of fragment pairs supporting cis or trans does not reach 90% of the total, perform microsatellite detection on the fragment pairs; if microsatellites exist, record the dibasic position as an unknown haploid phasing result; if microsatellites do not exist, report an error and stop typing.

[0029] As one implementation, the method for generating the complete sequence in step S26) includes the following steps: if no dibase position is detected in step S22), a complete sequence is generated; if a dibase position is detected, a pair of sequences is generated. When a dibase position is encountered during the generation process, the two bases are assigned to the two sequences according to the haploid phasing result; if the haploid phasing result of the dibase position encountered during the generation process is unknown compared to the previous pair of dibase positions, the sequence pair is doubled, and the dibase position is assigned to the doubled sequence pair in cis and trans configurations, respectively; as another implementation, when any base (not just a dibase position) is assigned to the sequence, "deletion" is treated the same as regular bases.

[0030] In one implementation, the alignment method in step S27) includes the following steps: aligning the complete sequence with the HLA genotyping algorithm database sequentially according to the order of HLA genes of type I, type II, and type III. In another implementation, the alignment order for type I HLA genes is as follows: antigen recognition domain (ARD), exon 4, exon 5, exon 1, exon 6, exon 7, exon 8 (skipped if absent), intron 2, intron 3, intron 4, intron 1, intron 5, intron 6, and intron 7 (skipped if absent); the alignment order for type II HLA genes is: ARD, exon 3, exon 4, exon 1, exon 5 (skipped if absent), exon 6 (skipped if absent), intron 2, intron 3, intron 1, intron 4 (skipped if absent), and intron 5 (skipped if absent).

[0031] As one implementation method, for some alleles that lack the complete gene sequence but only have partial exon sequences, if there are sequence missing parts during alignment, they should be considered as a perfect match.

[0032] As one implementation method, after the comparison, if multiple candidate sequence pairs are not eliminated, all results should be reported, prioritizing them according to the proportion of each allele in the population. However, all results may be the genotyping results for this instance, differing only in probability. If only one candidate sequence pair is not eliminated after comparison, the comparison result is the genotyping result for this instance. If all candidate sequence pairs are ultimately eliminated, all eliminated sequence pairs in the last round should be reported, along with any new allele sequence discovery. Simultaneously, the eliminated candidate sequence pairs in the last round should be compared with their remaining candidate alleles in the last round, reporting the differences from all candidate alleles. As another implementation method, when reporting new allele sequence discoveries, other genotyping methods in the field can be used for detection.

[0033] The present invention also provides a storage medium or processor, the storage medium including a stored program, the stored program executing the construction method as described in the above technical solution or the HLA typing method as described in the above technical solution; The processor is used to run a program, which executes the construction method or the HLA typing method described in the above technical solution.

[0034] To further illustrate the present invention, the technical solutions provided by the present invention will be described in detail below with reference to the accompanying drawings and embodiments, but these should not be construed as limiting the scope of protection of the present invention.

[0035] Example 1 An HLA typing method 1. Construct the algorithm database, the steps are as follows: 1) Download the database; Download the latest version of the IMGT / HLA alignment database from the official website of the International Immunogenetics Database (IMGT). The download path is https: / / github.com / ANHIG / IMGTHLA / tree / Latest / alignments. The downloaded database structure is as follows: Figure 1 As shown; when downloading the database, each locus requires two files, namely the allele with known full sequence (the corresponding file has the suffix .gen) and the allele with known partial exon sequence (the corresponding file has the suffix .nuc).

[0036] 2) Align, verify, and repair the downloaded database; For example, in the IMGT / HLA database of IMGT version 3.6, there is an alignment error in the allele B*51:01:01:124 (e.g.) Figure 2 (as shown in the display position), which leads to a large number of unnecessary insertions and deletions. Therefore, such errors are checked and repaired, and the algorithm database is constructed based on the checked and repaired database.

[0037] 3) Confirm the absolute coordinates of each locus sequence and generate the algorithm database; 3.1) Since the alignment database was downloaded in step 1), each allele is of equal length, and the start and end positions of its exons and introns are also the same. The alleles with known full-sequence sequences and those with known partial exon sequences at each locus are compared. The intron sequences of the alleles with known full-sequence sequences are used to fill in the gaps between the ends of the alleles with known partial exon sequences, with a fill length of 149 nt. In this step, the alleles with known full-sequence sequences are selected when their exon sequences are identical to or closest to the exon sequences of the alleles with known partial exon sequences to be filled. The filled sequence is as follows: Figure 3 As shown; 3.2) After the sequence is completed, since all the bases of the alleles correspond to statistical absolute coordinates, the absolute positions of the unsequenced parts (represented by *) and the deleted parts (represented by ·) are removed, leaving only the sequence of bases and the absolute positions corresponding to these bases. That is, all the bases of the alleles correspond to statistical absolute coordinates.

[0038] 3.3) After processing in step 3.2), all bases of all alleles have corresponding absolute coordinates. Take a 150nt string of these sequences (it should be noted that the selection of 150nt string in this embodiment is for data with a common read length of 150nt in current high-throughput sequencing (NGS), and should not be limited to the entire scope of protection of this patent). With a step size of 1, remove duplicates from the base sequences. For example, if two or more 150nt base sequences are the same but have different absolute coordinate positions, only one entry is retained, and both or more absolute coordinate positions are recorded in this entry.

[0039] 3.4) Compare databases generated from different loci. If databases containing 150nt repeating bases are found at two or more loci, these repeating databases are removed from their respective databases and reclassified into a combined database. For example... Figure 4 The image shown is a partial list of entries from the HLA-B / HLA-C repeat database. The first entry indicates that this 150nt base sequence appears simultaneously at both HLA-B and HLA-C loci, and its absolute position at the HLA-B locus has two possible variations.

[0040] After completing the above processing, the algorithm database is built.

[0041] 2. Analyze the high-throughput sequencing data; The high-throughput sequencing data were second-generation sequencing data, and the libraries were HLA gene-enriched libraries constructed using the probe fishing method, covering HLA-A, HLA-B, HLA-C, HLA-E, HLA-F, HLA-G, HLA-H, HLA-DRB1, HLA-DRB3, HLA-DRB4, HLA-DRB5, HLA-DQB1, HLA-DPB1, HLA-DQA1, HLA-DPA1, and MICA and MICB; the sequencing strategy was PE150.

[0042] 1) Perform quality control on the raw data: discard all low-quality PE150 fragment pairs containing any bases with a confidence level below 95%.

[0043] 2) Compare the quality-checked data with the algorithm database constructed in step 1. When genotyping a specific locus, it should be compared with both the independent database for that locus and the combined database containing that locus. For example, when genotyping the HLA-B locus, it should be compared with the HLA-B database and all duplicate databases containing HLA-B. In this embodiment, the HLA-B, HLA-A / HLA-B, and HLA-B / HLA-C databases need to be compared.

[0044] The comparison method is direct comparison, disallowing any form of mismatch. Therefore, this invention avoids using traditional bioinformatics algorithms that allow mismatches, such as Bowtie2, BWA, and HISAT2, significantly accelerating the software comparison speed.

[0045] 3) In this embodiment, after alignment, all fragment pairs that can be perfectly aligned to the above three databases (HLA-B, HLA-A / HLA-B, and HLA-B / HLA-C) are recorded, and the results are as follows: Figure 5 As shown.

[0046] The comparison results can be divided into the following three categories of fragment pairs: (1) The success ratio of the two fragment pairs in the HLA-B database; (2) One fragment is compared to the HLA-B database, and the other is compared to the duplicate database; (3) The two pairs in the fragment pair are compared to the duplicate database; For the three types of fragment pairs above, the first two can be identified as originating from the HLA-B locus, while the locus of the third cannot be determined.

[0047] It should also be noted that for partially matched segments, the database contains more than one type of absolute coordinate, requiring the extraction of the common portions from these multiple absolute coordinates, such as... Figure 6 As shown: For this fragment pair, since one fragment matches the HLA-B database, it can be known that the fragment pair originates from the HLA-B locus. Fragment 1 corresponds to one set of absolute coordinates. Fragment 2 corresponds to two sets of absolute coordinates. Comparing the two sets of absolute coordinates, 149 out of 150 absolute positions are the same, which are called definite positions. One different position is called a wobble position. Fragment 1 of this fragment pair can increase the coverage of 150 HLA-B absolute positions, and fragment 2 can increase the coverage of 149 HLA-B absolute positions. The remaining 1 base is recorded in the wobble base library (1 guanine originates from position 1536 or 1561, denoted as 1G, 1536, 1561).

[0048] 4) A schematic diagram of the initial absolute coverage of HLA-B loci after statistical analysis is shown below. Figure 7 At this point, further allocation is needed for fragment pairs with uncertain locus origins and wobble bases.

[0049] Taking the HLA-B locus genotyping of this sample as an example, all fragment pairs with uncertain locus origins were aligned to the HLA-B / HLA-C database. These uncertain fragment pairs should be temporarily considered as originating from HLA-B and compared with existing HLA-B coverage data, while simultaneously being temporarily considered as originating from HLA-C and compared with existing HLA-C coverage data. If, at any absolute position, an uncertain fragment pair differs from existing coverage data and has a coverage of less than 10%, it should be confirmed that it does not originate from that locus. If, after the above comparisons, only one candidate locus remains for the fragment pair, it can be considered to originate from that locus.

[0050] For positional rocking bases, a similar method is needed for allocation. If a rocking base does not appear in the confirmed coverage data at any absolute position, that position can be removed from the candidate positions for rocking bases. If only one rocking base remains after processing, it can be considered to have originated from that absolute position.

[0051] The resulting coverage data is the final coverage data for that locus.

[0052] 5) Classify the regions (e.g., exon 2, intron 3, etc.) and states at each location; If the coverage depth at any location is less than 10% of the average coverage depth in that region, that location is marked as a low-coverage location. If any low-coverage location occurs in an exon, an error should be reported and typing should be stopped.

[0053] If, in any location's coverage data, the third most common base (in the coverage data, deletion is considered a single base) accounts for more than 5%, that location must be marked as contaminated / chimerous. Locations marked as contaminated / chimerous should be checked for the presence of microsatellites nearby. Microsatellites are defined as seven or more repeated bases, or five or more repeated dibases. If microsatellites are found at a contaminated / chimerous location, it should be relabeled as bimodal. If no microsatellites are found at a contaminated / chimerous location, an error should be reported and typing should be stopped. The final coverage data is as follows: Figure 8 As shown.

[0054] If, in any location covered by data, the ratio of the first and second most abundant bases is less than 6, that location is recorded as a dibase and marked as imbalanced. All imbalanced locations should be manually checked to adjust library construction strategies, such as probe design.

[0055] 6) Haploid phasing is performed at all dibase positions, following these steps: List all adjacent dibase pairs and determine whether the relationship between these dibase positions and the previous adjacent dibase position in the table is cis (Cis) or trans (Trans) by detecting the fragment pairs that cover these positions.

[0056] like Figure 9As shown, for a pair of dibase positions, the dibases at the preceding position are recorded as F1 and F2, and the dibases at the following position are recorded as R1 and R2. Fragment pairs covering this pair of dibase positions can support the F1-R1 cophase, F2-R2 cophase, F1-R2 cophase, F2-R1 cophase, and others. If the number of fragment pairs covering this pair of dibase positions is less than 50, the two positions are marked as distant. Otherwise, if more than 90% of the fragment pairs supporting the F1-R1 and F2-R2 cophases are marked as cis, and if more than 90% of the fragment pairs supporting the F2-R1 and F1-R2 cophases are marked as trans, the two positions are marked as contradiction. Otherwise, it is marked as contradiction. If a contradiction is found, it should be checked whether there are microsatellites near the contradiction. If not, an error should be reported and typing should be stopped.

[0057] 7) Generate the haploid sequence to be aligned based on the coverage data and the dibase positions of the haploid phase. The specific steps are as follows: Determine if a dibase position exists. If not, generate a haploid sequence. If so, initially generate a pair of haploid sequences, adding the bases from the overlay data one by one to the two haploid sequences. When the first dibase position marked as too far away is encountered, duplicate the existing sequence pair, adding the new dibase in cis and trans arrangements to the two duplicated sequences respectively. If more dibase positions are encountered, double the existing sequence pairs. The generated sequence is as follows. Figure 10 As shown.

[0058] 8) Compare the generated haploid sequences with the algorithm database: (1) Although theoretically all generated sequence pairs could be correctly genotyped, considering that different types account for too low a proportion in the population or have never been recorded, these sequence pairs need to be removed. In this embodiment, a total of 8 sequence pairs were generated for genotyping of the HLA-B locus. First, the ARD region was compared, and 4 pairs could be successfully matched with alleles in the database. These 4 pairs were compared on exon 4, and the range of matched alleles was narrowed down to the alleles that each sequence successfully matched in the ARD region. A total of 2 pairs of sequences were matched with alleles on exon 4. A total of 2 pairs of sequences were matched with alleles on exon 5. A total of 2 pairs of sequences were matched with alleles on exon 6. A total of 2 pairs of sequences were matched with alleles on exon 7.

[0059] (2) Then, the non-coded regions are compared: Two sequence pairs aligned to alleles in the 5' untranslated region. Two sequence pairs aligned to alleles in intron 1. Two sequence pairs aligned to alleles in intron 2. Two sequence pairs aligned to alleles in intron 3. Two sequence pairs aligned to alleles in intron 4. Two sequence pairs aligned to alleles in intron 5. Two sequence pairs aligned to alleles in intron 6. One sequence pair aligned to alleles in the 3' untranslated region.

[0060] The sequence verification results are HLA-B*37:01:01:01 and HLA-B*40:02:01:01. The above genotyping results are the genotyping results of this sample.

[0061] Example 2 To verify the accuracy of the HLA genotyping method of this invention, 350 samples were genotyped using the method described herein, and these samples were analyzed using the industry-recognized open-source analysis software HLA-HD. Due to limitations of the probes used at the time, rather than the algorithm itself, genotyping was performed on 11 loci at these sites: HLA-A, HLA-B, HLA-C, HLA-DRB1, HLA-DRB3, HLA-DRB4, HLA-DRB5, HLA-DQB1, HLA-DPB1, HLA-DQA1, and HLA-DPA1, totaling 3850 sample loci. Among these, 61 sample loci showed different results from the two algorithms, accounting for 1.58%. After investigation, the genotyping accuracy of the algorithm of this invention was 100% for these sample loci.

[0062] Although the above embodiments have provided a detailed description of the present invention, they are only some embodiments of the present invention, and not all embodiments. People can obtain other embodiments based on these embodiments without creative effort, and these embodiments all fall within the protection scope of the present invention.

Claims

1. A method of constructing a database of HLA typing algorithms, characterized in that, The method comprises the following steps: 1) downloading the IMGT / HLA alignment database in the international immunogenetics database, each locus of the IMGT / HLA alignment database comprising full sequence known alleles and partial exon sequence known alleles; 2) processing the IMGT / HLA alignment database to repair alignment errors; 3) after completing the processing in step 2), for each locus, filling in the partial exon sequence known alleles with the full sequence known alleles; 4) cutting each locus of alleles according to a high-throughput sequencing read length with a step length of 1 and removing duplicates to obtain a database of each locus; 5) comparing all locus databases, and when a repeated base sequence is found in two or more locus databases, deleting the repeated base sequence and classifying it into a joint database to obtain the HLA typing algorithm database; In the HLA typing algorithm database, each locus comprises a database shown in one of the following two items: I: a locus independent database; II: a locus independent database and a joint database comprising the locus.

2. The construction method according to claim 1, characterized in that, The method for filling in in step 3) comprises filling in both ends of the partial exon sequence known alleles with intron sequences of the full sequence known alleles.

3. The construction method of claim 2, wherein, The selection method of the full sequence known alleles comprises selecting the closest full sequence known alleles with the same exon sequence as the partial exon sequence.

4. The construction method of claim 1, wherein, After step 3), the absolute positions of unsequenced parts and deleted parts shown in each locus of alleles are deleted, and only the sequence with base composition and the absolute position of the base are retained.

5. The application of the HLA typing algorithm database constructed by the construction method in any one of claims 1-4 or constructed by using the construction method in any one of claims 1-4 in HLA typing.

6. A method of HLA typing, characterized by, The method comprises the following steps: S1) constructing an HLA typing algorithm database by using the construction method in any one of claims 1-4; S2) comparing and typing the quality inspected high-throughput sequencing data with the HLA typing algorithm database; S21) when comparing and typing a specific locus, comparing the quality inspected high-throughput sequencing data with both the independent database of the locus and the joint database comprising the locus to obtain the following comparison results: S211) when both sequences in a fragment pair are compared to the independent database, confirming that the fragment pair is from the locus; S212) when one sequence in a fragment pair is compared to the independent database and the other sequence is compared to the joint database, confirming that the fragment pair is from the locus; S213) when both sequences in a fragment pair are compared to the joint database, the locus of the fragment pair cannot be confirmed; S22) after comparing and typing the fragment pairs in S211) and S212) to a specific locus, counting the coverage and base types at each position. S221) If the coverage of any position is less than 10% of the average coverage depth of the exon or intron region where the position is located, mark the position as a low-coverage position; if the low-coverage position is located in an exon, report an error and stop the typing; S222) If the abundance of the third most abundant base in the coverage data of any position is more than 5%, perform microsatellite detection on the position; if the position does not contain a microsatellite, report an error and stop the typing; if the position contains a microsatellite, record the second and third most abundant bases at the position; S223) If the report of errors in S241) and S242) is not triggered at any position, count the double-base positions of the position; S224) Temporarily do not place the base that can be placed in multiple positions, and record it separately; S23) Compare the coverage data of the fragment pairs in S213) and the fragment pairs of the already determined loci one by one; if the fragment pairs are different from the determined coverage data of the locus at any position and the coverage is less than 10%, confirm that the fragment pairs do not come from the locus; when there is only one candidate locus for the fragment pairs, confirm that the fragment pairs come from the locus, and after confirmation, add the fragment pairs to the coverage data of the correct locus; S24) Compare the coverage data of the positions where the base that can be placed in multiple positions is encountered in S224) with the determined coverage data: if the base that is confirmed to exist in the determined coverage data is inconsistent with the base that can be placed in multiple positions, delete the position from the candidate positions of the base that can be placed in multiple positions; when there is only one candidate position for the base that can be placed in multiple positions, confirm that the base is located in the remaining position; S25) Perform haplotype phasing and polyploidy phenomenon analysis on each pair of adjacent double-base positions; S26) Based on the haplotype phasing analysis in S25), generate one or more complete sequences for each locus; S27) Based on the complete sequences in S26), compare the complete sequences with the HLA typing algorithm database to confirm the typing results.

7. The HLA typing method of claim 6, wherein, The haplotype phasing in step S25) includes the following steps: for each pair of adjacent double-base positions, extract all fragment pairs that cover the two double-base positions simultaneously; Respectively count the number of fragment pairs supporting cis and trans, when the number of fragment pairs supporting cis or trans is more than 90% of the total number, mark the relationship between the two positions as cis or trans; when the number of fragment pairs supporting cis or trans does not reach 90% of the total number, perform microsatellite detection on the fragment pairs; if there is a microsatellite, record the double-base positions as unknown haplotype phasing results; if there is no microsatellite, report an error and stop the typing.

8. The HLA typing method of claim 6, wherein, The method for generating complete sequences in step S26) includes the following steps: If no double-base position is detected in step S22), generate one complete sequence; If a double base position is detected, a pair of sequences is generated, in the process of generating, when a double base position is encountered, two bases are assigned to two sequences according to the haplotype phasing result; if the haplotype phasing result of the double base position encountered in the generating process is unknown from the previous pair of double base positions, the sequence pair is doubled, and the double base position is respectively assigned to the doubled sequence pair in cis and trans.

9. The HLA typing method of claim 6, wherein, The method of the alignment in step S27) comprises the following steps: sequentially aligning the complete sequence with the HLA typing algorithm database according to the order of class I HLA genes, class II HLA genes and class III HLA genes.

10. A storage medium or processor, characterized by, The storage medium comprises a storage program, and the storage program executes the construction method according to any one of claims 1-4 or the HLA typing method according to any one of claims 6-9. The processor is configured to run a program, and the program is configured to execute the construction method according to any one of claims 1-4 or the HLA typing method according to any one of claims 6-9.

Images