A somatic variant analysis method based on individualized mhc haplotype reference

Abstract

The application discloses a somatic cell variation analysis method based on individualized MHC haplotype reference and belongs to the technical field of bioinformatics. In view of the problems that the prior art relies on a public reference genome to analyze an MHC polymorphism region, leading to serious alignment mismatch, and the existing analysis tool relies on known typing, leading to low resolution, the method obtains individualized haplotype sequences of subjects and extracts high-density heterozygous variation sites; sequencing data is mapped to the haplotype sequences, and the integrated copy number and the minor allele frequency are quantified based on the deduplication logic; after filtering the double-allele deletion region, the heterozygous deletion (LOH) state is finely determined based on a smoothing segmentation algorithm; in addition, competitive mapping filtering is performed by constructing a double-haplotype joint reference sequence to identify somatic mutations. The application eliminates the alignment bias caused by the public reference, breaks through the pre-typing limitation, and realizes the accurate analysis of the LOH and mutations in the MHC region with ultra-high resolution.

Description

Technical Field

[0001] This invention relates to a method for analyzing somatic variations based on individualized MHC haplotype references, and more particularly to a method for phasing, allele-specific copy number inference, and loss of heterozygosity (LOH) identification of somatic variations in the major histocompatibility complex (MHC) region of a subject's tumor tissue, belonging to the fields of bioinformatics and computational genomics. Background Technology

[0002] The Major Histocompatibility Complex (MHC) region, located on human chromosome 6p21.3 and approximately 5 Mb in length, is one of the most structurally complex functional regions in the human genome, exhibiting the highest allelic polymorphism and the highest density of paralogous genes. This region contains key human leukocyte antigen (HLA) genes that control antigen presentation and immune responses. The MHC region has long been considered one of the most difficult regions in the genome to annotate. Accurately resolving genetic variations in the MHC region, such as loss of heterozygosity (LOH), copy number variation (CNV), and somatic mutations, is of paramount clinical value for assessing tumor immune escape mechanisms, predicting the efficacy of immunotherapy, and guiding the design of neoantigen vaccines.

[0003] However, existing MHC region sequence analysis techniques face the following insurmountable technical bottlenecks: (1) Mapping bias due to reference genome dependence: Traditional somatic mutation detection and variant analysis strategies mainly rely on aligning short-read sequencing data (50–300 bp) to a public reference genome (such as GRCh38). However, due to the extremely high allelic polymorphism and extensive paralogous sequences in the MHC region, such alignment methods relying on a single reference genome face fundamental limitations. This alignment can lead to severe mismatches, resulting in an unacceptable false positive rate, causing the true mutation signal to be masked by noise. Multiple studies have confirmed that in short-read sequencing, the MHC region generates multi-mapping reads, some of which have the same alignment quality among these genes. Standard analysis procedures typically discard multi-mapping reads in order to achieve unique alignments. This approach not only significantly reduces the coverage of the target region but also introduces allele bias, underestimating the quantification of alleles that differ significantly from the reference genome: (I) Paralogous genes within the MHC region (such as HLA-A, HLA-B, and HLA-C) exhibit high sequence similarity, making it difficult to uniquely map some short reads to a single gene copy; (II) This region is rich in repetitive elements, structural variations, and insertion / deletion polymorphisms, and short read lengths are often insufficient to cross highly conserved paralogous regions or complex structural regions, making it difficult to establish accurate long-range connections; (III) Commonly used short-read alignment algorithms assume a difference rate of less than 1% between the sample and the reference genome, while the difference rate among alleles in the MHC region can be significantly higher than this assumption, leading to an increase in systematic alignment errors (including false negatives and false positives). Therefore, despite significant progress in long-read sequencing technology in recent years, the vast majority of sequencing data from large-scale clinical cohorts worldwide are still generated based on short-read sequencing platforms. Therefore, algorithms based on short read lengths must address the fundamental alignment bias caused by single reference genome dependence in order to maximize the scientific value of massive amounts of existing cohort data.

[0004] (2) LOH analysis tools that rely on known alleles have a very narrow scope of application: Existing allele-specific copy number inference tools (such as LOHHLA and MHC Hammer) require pre-identified HLA typing results before analysis. This design choice is a strategy adopted by researchers to cope with extreme MHC polymorphism: LOHHLA aligns reads to patient-specific HLA type information rather than a public reference genome, thus avoiding the alignment problem of a single reference genome. This dependence strictly limits the scope of analysis to known alleles already included in the database, and cannot cover the entire MHC region (especially complex non-classical genes), and usually can only provide a rough gene-level resolution, making it difficult to accurately locate the boundaries of copy number changes.

[0005] (3) LOH analysis tools based on heterozygous sites suffer from extremely low resolution: To circumvent the aforementioned limitation of "relying on pre-genotyping," existing analysis tools such as LILAC have emerged. These tools attempt to directly infer segment copy numbers using confirmed heterozygous variant sites. However, in highly polymorphic MHC regions, heterozygous sites that can be identified with high confidence using traditional alignment methods are extremely sparse. Although LILAC supports the detection of novel germline variations, including insertions and deletions, many short reads with low alignment quality or originating from gene transition regions in tumor samples cannot be effectively utilized, resulting in a sparse availability of heterozygous genetic markers. This sparsity forces existing algorithms to perform copy number inference only over a very large genomic span (tens to hundreds of kb). This large-span averaging process completely masks key variant events at the gene level (e.g., single gene) or subgenetic level (e.g., single exon), making the analysis results only usable for detecting LOH at the large-scale chromosome arm or region level, unable to finely locate event boundaries.

[0006] In summary, there is currently a lack of a comprehensive solution in the field that can completely eliminate the limitations of public reference genomes, eliminate the need for pre-provided HLA typing, and obtain sufficiently dense heterozygous sites to achieve high-resolution MHC full-region analysis. Summary of the Invention

[0007] [Technical Issues] The technical problem to be solved by this invention is to overcome the mismatch caused by paralogous sequences in the MHC region of short read sequencing data, eliminate the dual dependence on public reference genome and pre-typed HLA alleles, obtain a high-density set of heterozygous variant sites specific to the subject, and support high-precision identification of subgenetic LOH and somatic mutations.

[0008] [Technical Solution] The purpose of this invention is to overcome the shortcomings of existing technologies and provide a variation analysis method based on individualized MHC haplotypes. This method introduces subject-specific assembled haplotypes as an absolute benchmark, eliminating mismatches and obtaining a high density of true heterozygous sites. This solves the problem of low LOH resolution caused by the sparse heterozygous sites in existing technologies (such as LILAC), achieving accurate allele-specific copy number quantification and somatic mutation detection.

[0009] This invention provides a method for LOH variant analysis based on individualized MHC haplotypes, comprising the following steps: S1. Obtain two individualized assembled pseudo-haplotype sequences of the target region of the subject; perform double sequence alignment on the two pseudo-haplotype sequences to identify and extract subject-specific high-density heterozygous variant sites.

[0010] S2. Obtain short-read sequencing data (sequence reads) of the target sample and control sample from the subject; independently align the sequencing data (sequence reads) of the tumor cells with the sequencing data (sequence reads) of the control sample to the two pseudo-haplotype sequences described in step 1; based on the haplotype-specific alignment results, calculate the comprehensive copy number (CLOSC) at each genomic position of the target sample. CN combined ) and calculate the minor allele frequency (BAF).

[0011] S3, based on the above CN combined The statistical distribution of biallelic deletions is used to calculate the lower limit threshold LB for biallelic deletion detection, thereby identifying and filtering biallelic deletion regions in the target sample; wherein, the formula for calculating the lower limit threshold LB is: , in, μ MHC region CN combined The mean, σ Its standard deviation; S4. In the remaining region after excluding heterozygous variant sites within the biallelic deletion region, combined with the... CN combined BAF and tumor purity ρ Calculate the specific copy number of allele 1 ( ) and the specific copy number of allele 2 ( CN 2) Clustering algorithms are used to smoothly segment heterozygous variant sites, and loss-of-heterozygosity (LOH) regions are determined based on allele copy number differences; wherein, the following equations are used to calculate and CN 2: , .

[0012] In one implementation, in S1, a de novo assembly method based on whole-genome sequencing data or targeted MHC sequencing data is used to obtain two haplotype sequences of the subject's MHC region.

[0013] Preferably, two haplotype sequences from the subject's MHC region are obtained with reference to the patent with publication number CN110777192A.

[0014] In one implementation, the nucme / dnadiffr module of MUMmer4 is used for double sequence alignment.

[0015] In one implementation, in S1, heterozygous variant sites are identified by performing double sequence alignment on two haplotype sequences and applying the following filtering conditions: (I) alignment length ≥ 1500 bp; (II) in the multiple alignment region, only the alignment result with the highest consistency score is retained.

[0016] In one implementation, in S1, for complex regions where two haplotypes cannot be directly compared due to extremely high heterogeneity, the standard reference genome is used as a computational bridge to infer and supplement additional heterozygous variant sites by comparing the two haplotypes to the standard reference genome respectively. Specifically, each haplotype assembly sequence is mapped to a standard reference genome, and the physical coordinates based on the assembly sequence are converted to standard coordinates of the standard reference genome, thereby generating a sample-specific variant call file; In the process of generating the variant call file, regions where only one haplotype can be aligned to the standard reference genome are excluded.

[0017] In one implementation, the reference genome version includes, but is not limited to, hg19 (GRCh37), hg38 (GRCh38), or CHM13.

[0018] In one embodiment, the target sample refers to sequencing data of a tumor tissue sample from the same individual; the control sample refers to sequencing data of a peripheral blood leukocyte sample, saliva sample, oral swab sample, cell-free DNA sample from a fingernail, or adjacent normal tissue sample from the same individual.

[0019] In one embodiment, in S2, the target sample is a tumor tissue sample, and the control sample is a peripheral blood leukocyte sample.

[0020] In one implementation, in S2, the sequencing data of the target sample and the control sample are independently aligned to the two pseudo-haplotype sequences using the Razers3 alignment tool, with the parameters set as follows: the number of mismatches is set to 0~1, and uniquely mapped read pairs are retained.

[0021] In one implementation, the sequence alignment tool includes razers3.

[0022] In one implementation, in S2, the CN combined Calculate using the following formula: , in, Cov tumor This refers to the overall coverage of genomic locations in tumor samples. Cov WBCThis refers to the overall coverage of genomic locations in the control sample, where M is the ratio of the size of the control sequencing library to the size of the tumor sequencing library. ρ For tumor purity.

[0023] In one implementation, the overall coverage refers to the sum of coverage mapped to each haplotype, minus the shared coverage of equivalent mappings on two pseudo-haplotypes.

[0024] In one embodiment, the ρ Obtained through pathological evaluation, whole-exome sequencing evaluation, etc.

[0025] In one implementation, simultaneously, within a window with a radius of 75 bp centered on the high-density heterozygous variant site, the minor allele frequency (BAF) is calculated based on the allele-specific read that is uniquely aligned to a single pseudohaplotype.

[0026] In one implementation, in S3, the confirmation of the biallelic deletion region must satisfy: (I) comprehensive copy number CN combined <LB; (II) Continuous region length > 25 kb; (III) No assembly gap exists within the region. The assembly gap is defined as a length > 12.5 kb within the continuous fragment where no sequencing read can be aligned to any pseudohaplotype.

[0027] In one implementation, in S4, genomic regions that simultaneously meet the following conditions are selected as regions lacking heterozygosity: (a)| - CN 2|≥Set value, where the set value is used to measure the degree of difference between two haplotype copy numbers; (b) and CN The smaller of the two values ​​is below the preset lower limit for the number of copies.

[0028] In one embodiment, the set value is 0.5; preferably 0.7.

[0029] In one implementation, the lower limit of the copy number is 0.5.

[0030] In one implementation, a deletion LOH is defined as a major allele copy number ≤ 1.5, and a copy-neutral LOH is defined as a major allele copy number > 1.5.

[0031] This invention also provides a method for somatic mutation variation analysis based on individualized MHC haplotypes, comprising the following steps: S1. Obtain two individualized assembled pseudohaplotype sequences from the target region of the subject; S2. Obtain short-read sequencing data of the target tumor sample, align the sequencing data to a standard reference genome, identify and remove PCR duplicate reads based on the alignment results, and obtain deduplicated sequencing reads. S3. Merge the two pseudo-haplotype sequences obtained in S1 to construct a double haplotype joint reference sequence; map the deduplication sequencing reads obtained in S2 to the double haplotype joint reference sequence, and screen reads that meet the following conditions: the sequence identity between the read and the double haplotype joint reference sequence is not less than 95%, and the read has only one optimal matching position in the double haplotype joint reference sequence; S4. Based on the reads obtained after screening, apply the somatic mutation detection algorithm to output the set of somatic mutations in the MHC region of the target tumor sample.

[0032] In one implementation, a de novo assembly method based on whole-genome sequencing data or targeted MHC sequencing data is used to obtain two haplotype sequences of the subject's MHC region.

[0033] Preferably, two haplotype sequences from the subject's MHC region are obtained with reference to the patent with publication number CN110777192A.

[0034] In one implementation, the sequencing data is aligned to a standard reference genome using alignment tools such as Razers3 or BWA-MEM.

[0035] In one implementation, the standard reference genome includes, but is not limited to, hg19 (GRCh37), hg38 (GRCh38), or CHM13.

[0036] In one implementation, the read segment having only one optimal alignment position in the double haplotype joint reference sequence means that: In the joint reference sequence composed of two haplotype sequences, there is one and only one optimal alignment result for this read, excluding multi-mapped reads that can be aligned to two haplotypes simultaneously and reads that cannot be aligned.

[0037] In one implementation, the somatic mutation detection algorithm is the Strelka2 mutation invocation tool.

[0038] In one embodiment, the control sample includes a peripheral blood leukocyte sample, saliva sample, oral swab sample, cell-free DNA sample from a fingernail, or adjacent normal tissue sample from the same individual.

[0039] In one embodiment, the MHC region is the 6p21.3 region of human chromosomes.

[0040] Beneficial effects (1) Solving the problem of sparse loci and achieving ultra-high resolution: Unlike existing technologies such as LILAC that rely on standard comparisons to obtain sparse variants, this invention directly performs high-precision cross-comparison between two individualized pseudo-haplotypes, which can exhaustively extract all true heterozygous loci of the individual. The dense locus distribution means that copy number inference is no longer forced to cross huge intervals, realizing true gene-level or even sub-gene-level LOH breakpoint detection.

[0041] (2) Completely eliminate mapping bias and reduce false positives: Completely abandon public reference genomes and directly align to individual-specific sequences. Through a strict "unique alignment" filtering strategy, homologous sequence mismatches are fundamentally eliminated, which greatly improves the detection accuracy of somatic mutations in the MHC region.

[0042] (3) Global analysis that breaks free from genotyping limitations: No prior HLA genotyping input is required; comprehensive genetic analysis of the entire MHC region, including classical and non-classical HLA genes and any alleles not included in the database, can be directly achieved. Attached Figure Description

[0043] Figure 1 A schematic diagram of the overall process of a high-resolution MHC region variation analysis method provided in an embodiment of the present invention; Figure 2 The diagrams are schematic diagrams of specific examples of biallelic deletion and loss of heterozygosity (LOH) identified in the embodiments of the present invention; (A, B) are schematic diagrams of biallelic deletion; (C) is a schematic diagram of loss of heterozygosity (LOH).

[0044] Figure 3 This is a schematic diagram of the LOH detection results for the entire MHC region of a single sample; Figure 4 A schematic diagram for PacBio HiFi long read sequencing validation of LOH in the MHC region; (A, B) are views of the control sample (WBC), target sample (tumor) HLA-C and HLA-DPB1 locus IGV; (C) paired rank sum test of coverage ratio. Figure 5 A schematic diagram illustrating the accuracy of LOH region definition in this invention for PacBio HiFi long-read sequencing verification. Detailed Implementation

[0045] Many details are described in order to enable this application to be better understood; some features may be omitted in different circumstances, or implemented in other ways, or used in combination with other features; those skilled in the art can make various modifications and variations without departing from the spirit and scope of this invention.

[0046] As used herein, “comprising” means “including but not limited to”; “for example” means “for example but not limited to”; “preferred” means an embodiment that has better technical effects within the scope of the present invention; “more preferred” means an embodiment that further improves performance based on the preferred embodiment; and “most preferred” means that the best performance is achieved under the current technical conditions based on the more preferred embodiment.

[0047] As used in this article, the “MHC region” covers chromosome 6:29942530-33109100 (GRCh38), including 23 classical and non-classical HLA genes and their complete exon and intron sequences: HLA-A, HLA-B, HLA-C, HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DQB1, HLA-DRA, HLA-DRB1, HLA-DRB3 / 4 / 5, HLA-E, HLA-F, HLA-G, HLA-H, HLA-J, HLA-K, HLA-L, HLA-P, HLA-V, HLA-W, HLA-Y, and HLA-Z.

[0048] As used herein, the term "heterozygous variant site" refers to a genomic location where there is a base difference between two MHC haplotype sequences. Since the MHC region is one of the most polymorphic regions in the human genome, its SNP density is much higher than the average level of the whole genome. Therefore, in this invention, densely distributed heterozygous variant sites within the MHC region are referred to as "high-density heterozygous variant sites".

[0049] As used herein, the term "haploidy" refers to a nucleotide sequence in the diploid genome of a subject that originates from a single chromosome of the same parent. In the context of this invention, haploidy specifically refers to a single chromosome sequence of an MHC region obtained through assembly.

[0050] As used in this article, "biallelic deletion" refers to a situation in which both haplotypes in a certain genomic region of the target sample have lost their copy number.

[0051] As used herein, the terms "major allele" and "minor allele" do not refer to frequency in a population, but rather to the allele with a relatively higher copy number and the allele with a relatively lower copy number on two homologous chromosomes of the same individual. Specifically, a "major allele" is an allele or haplotype with a higher copy number in a genomic region; a "minor allele" is an allele or haplotype with a lower copy number in the same region.

[0052] As used in this article, “coverage” refers to the number or density of sequencing reads aligned to a target genome region.

[0053] As used in this article, “assembly gap” refers to a genomic region in which a continuous assembly sequence could not be obtained during individualized MHC haplotype assembly due to genomic sequence complexity (such as excessively long repetitive sequences or GC bias); in sequencing data, it appears as a continuous blank region without any read alignment support.

[0054] As used in this article, “tumor” refers to a malignant neoplasm originating from epithelial tissue, including but not limited to non-small cell lung cancer, melanoma, colorectal cancer, gastric cancer, and esophageal cancer.

[0055] As used in this article, "tumor purity" ρ "Tumor cells" refers to the proportion of tumor cells in the target sample (tumor sample) tissue; its value ranges from 0 to 1; the measurement methods include pathological evaluation or whole exome sequencing evaluation.

[0056] Example 1: Application of a high-resolution MHC region variation analysis method in the detection of loss of heterozygosity This embodiment provides a high-resolution method for MHC region variation analysis, incorporating the analytical workflow. Figure 1 Taking a tumor sample as the target sample and a white blood cell (WBC) sample as the control sample as an example, the specific implementation steps are as follows: Step S1: Identification of high-density heterozygous variant sites and cross-sample standardization S1.1: The method of the patent with publication number CN110777192A was used to obtain two haplotype sequences of the MHC region pre-assembled by the subject.

[0057] S1.2: The two false haplotype sequences obtained in S1 are submitted to the nucmer / dnadiff module of MUMmer4 for cross-alignment. To filter out false positive matches, sequence results shorter than 1500 bp are strictly excluded; in regions with multiple alignments, only the alignment result with the highest consistency is retained, and heterozygous variant sites are extracted accordingly.

[0058] S1.3: To enable cross-sample comparisons, each haplotype assembly sequence is mapped back to a standard reference genome (e.g., GRCh38). The physical coordinates based on the assembly sequence are converted to GRCh38 standard coordinates, generating a sample-specific variant call file. During this process, regions with only a single haplotype that can be aligned to the reference genome are excluded. Specifically, for complex regions where direct alignment between two haplotypes fails due to extremely high heterogeneity, this method introduces GRCh38 as a computational bridge, inferring and supplementing additional heterozygous variant sites by aligning two haplotypes to GRCh38 separately.

[0059] Step S2: Independent Mapping and Synthetic Copy Quantization S2.1: Acquire short-read sequencing data from tumor samples and peripheral blood leukocyte (WBC) samples, and use an alignment tool (razers3) to independently map the sequencing data to the two pseudo-haplotype sequences in S1. Extremely strict error tolerance (e.g., a maximum of one mismatch) is allowed, and uniquely mapped read pairs are preserved, meaning that reads can only be aligned to a single haplotype with no alternative positions.

[0060] S2.2: Use Samtools to calculate read depth. To avoid duplicate calculations of reads that do not cover heterozygous sites, the overall coverage at a location is defined as: the sum of coverage mapped to each haplotype, minus the shared coverage of equivalent mappings on two haplotypes.

[0061] S2.3: In tumor samples and WBC samples, calculate the overall coverage centered on each genomic location (e.g., a smooth window of 10kb). Cov tumor and Cov normal Calculate the global normalization factor M for the target analysis region (M is the ratio of the WBC sequencing library size to the tumor sequencing library size). Then calculate the [missing information - likely a specific value or feature] at that location. CN combined : , in, ρ Tumor purity can be obtained through pathological evaluation, whole-exome sequencing evaluation, etc.

[0062] S2.4: BAF calculations involve only allele-specific unique alignments. Within a window (e.g., 150 bp) centered on each heterozygous variant site, the coverage of reads strictly and uniquely aligned to haplotype 1 is extracted. Cov hap1_uniq ) and the coverage of reads that are strictly and uniquely matched to haplotype 2 ( Cov hap2_uniq ). Calculate the unique alignment coverage of minor alleles relative to the total specific unique alignment coverage of that locus (i.e. Cov hap1_uniq + Cov hap2_uniq The proportion of ) is used to obtain an accurate BAF.

[0063] Step S3: Detection and filtering of biallelic deletions S3.1: Calculate the total copy number of tumor samples across the entire MHC analysis region. CN combined mean ( μ ) and standard deviation ( σ The lower limit threshold (LB) for biallelic deletion is set as follows:

[0064] S3.2: Identify continuous regions with a length greater than a set value (e.g., 25 kb) and a combined copy number. CN combined Regions shorter than LB are considered candidate biallelic deletion regions. Further, if a candidate region contains a fragment exceeding a set length (e.g., 12.5 kb) lacking any evidence of false haplotype alignment, it is identified as an assembly gap and excluded. The remaining regions are confirmed as true biallelic deletion regions.

[0065] Step S4: Fine determination of allele-specific copy number and LOH status S4.1: Heterozygous sites falling within the biallelic deletion region identified in step S3 are removed and not included in subsequent LOH analysis. For the remaining heterozygous sites, combined with... ρ , CN combined and BAF, calculate the specific copy number of allele 1 based on the following set of equations ( CN 1) and the specific copy number of allele 2 ( CN 2):

[0066]

[0067] S4.2: Introduce and optimize a location-based clustering algorithm (bumphunter algorithm logic) to identify LOH regions: First, cluster heterozygous sites based on genomic physical location (using ClusterMaker, with the maximum gap set to 3,000 bp). Second, cluster based on the overall copy number ( CN combined The numerical values ​​divide the loci into three categories (e.g.) CN combined <0.7, 0.7≤ CN combined <1.3, CN combined ≥1.3). Within each classification segment, the copy number is smoothed using a local multinomial regression fitting and other smoothing algorithms (loessByCluster, span bpSpan = 200), and the smoothed values ​​are normalized to ensure... + CN 2=2. Identifying the maximum absolute difference in allele copy number (| - CN The region with a copy number of 2|>0.7 and a minor allele copy number of very low (<0.5) is considered as a candidate LOH region.

[0068] S4.3: Short LOH regions containing fewer than 3 heterozygous sites are removed. Finally, LOH events are classified according to the copy number of the major allele: when the copy number of the major allele is ≤1.5, it is classified as a deletion LOH; when the copy number of the major allele is >1.5, it is classified as a copy-neutral LOH (usually reflecting mitotic recombination or gene conversion). In terms of clinical significance, if more than half of the heterozygous sites in a specific gene (such as a specific HLA gene) fall into an LOH region, then the gene as a whole is considered to have undergone a loss of heterozygosity.

[0069] Example 2: Allele-specific copy number analysis and loss of heterozygosity detection 1. Identification and verification of copy number missing regions To verify the ability of this method to detect copy number deletions in the MHC region, a full-region copy number analysis was performed on a representative sample. The results are as follows: Figure 2 As shown.

[0070] Figure 2 The copy number distribution and coverage verification results of the MHC region of this sample are presented. Figure 2 A and Figure 2 B is calculated using this method. CN combined Distribution across the entire region: Blue scatter dots represent the original genomic loci. CN combined The calculated values, with the purple curve representing the smoothed copy number trend line. In the normal control region, the smoothed copy number fluctuates around the diploid level (copy number = 2), consistent with the baseline copy number characteristics of normal samples.

[0071] Further analysis revealed that within a specific interval of the MHC region (highlighted in orange in the image), the smoothed... CN combined The value is significantly lower than the preset missing threshold, and the abnormal region has a large continuous span, which is consistent with the typical characteristics of copy number loss. Therefore, it is determined that a copy number loss event has occurred in this region.

[0072] To rule out the possibility that the abnormal copy number signal originated from assembly or alignment errors, the coverage distribution of tumor samples and paired normal samples in this region was further compared. Figure 2A shows the sequencing coverage comparison of the region: the blue curve represents the coverage of normal samples, whose coverage level in the target region is consistent with that of the surrounding non-deleted regions, indicating that the genome assembly and alignment results in this region are reliable and there is no systematic technical bias; the red curve represents the coverage of tumor samples, whose coverage in the same target region is significantly lower than that of the surrounding regions and normal samples, which is consistent with the expected performance of copy number deletion.

[0073] Based on the copy number analysis and coverage verification results, this method can accurately identify copy number loss events in the MHC region, and the detection results have high reliability.

[0074] Figure 2 C shows the allele-specific copy number distribution at polymorphic sites in the MHC region of the target sample. Normal hap1 / Normal hap2 represents the allele-specific copy number calculation results for the normal control sample (peripheral blood leukocytes): in a normal diploid genome, the sum of the copy numbers of the two homologous chromosomes is 2, therefore the theoretical copy number of a single haplotype at each polymorphic site is 1. As shown in the figure, the estimated copy numbers of the two haplotypes (hap1 and hap2) in the normal sample fluctuate around 1, consistent with the expected distribution pattern of normal diploid tissue.

[0075] In the target samples (Tumor hap1 / Tumor hap2), this method identifies a continuous region by scanning a sliding window. Figure 2 The area highlighted in orange (bump region) shows a significant deviation in allele-specific copy number distribution. Within the orange-highlighted target region, the copy number estimates for hap1 and hap2 exhibit a clear imbalance: one haplotype shows a significantly reduced copy number, with its lowest value below the 0.5 threshold; while the other haplotype shows a correspondingly increased copy number or remains at a high level. This significant difference in copy number between the two haplotypes is consistent with the molecular characteristics of a LOH event, thus indicating that a loss of heterozygosity has occurred in this region.

[0076] The above results demonstrate that this method can achieve copy number resolution at the allele level within the MHC region, effectively identifying and accurately locating the occurrence interval of LOH events.

[0077] 2. Results of LOH detection at each locus in the entire MHC region of a single sample Taking sample 1652776 as an example, the detection was carried out by performing allele-specific copy number analysis on polymorphic sites at each locus in the MHC region. The number of heterozygous sites that can be used for LOH determination and the number of sites that are actually determined to be LOH in each locus were counted. The LOH detection results of the whole region were presented on a locus-by-locus basis.

[0078] The results showed that in sample 1, at both the HLA-C and HLA-DPB1 loci, and CN There was a significant imbalance between the two haplotypes, with the copy number of one haplotype significantly lower than that of the other and below the threshold, both meeting the criteria for LOH. Specifically, at the HLA-C locus (HLA-I class), 50 heterozygous loci were detected, of which 43 were suitable for LOH analysis, and 35 loci met the LOH criteria (81.4% LOH positive rate), indicating that LOH occurred at this locus. At the HLA-DPB1 locus (HLA-II class), 226 heterozygous loci were detected, of which 213 were suitable for LOH analysis, and all 213 loci met the LOH criteria (100% LOH positive rate), indicating that LOH occurred at this locus (Table 1). Figure 3 ).

[0079] The above results demonstrate that the method of the present invention can perform single-gene resolution LOH detection on the entire MHC region (covering HLA-I and HLA-II genes).

[0080] Table 1. LOH Detection Results of Samples

[0081] Note: total_het: total polymorphic sites; available_het: polymorphic sites used for analysis (normal sites with significant bias are removed); LOH_het: polymorphic sites where LOH occurred (within the bump, min( , CN 2) < 0.5, max( , CN 2>0.5) 3. Independent validation of LOH detection results by long-read sequencing To eliminate the potential impact of short-read sequencing alignment bias on LOH detection results, this method further employed PacBio HiFi long-read sequencing technology to independently validate the above samples, with the results as follows: Figure 4 As shown.

[0082] Figure 4 The visualization results in the IGV genome browser are shown, after aligning long-read sequencing data using two assembled haplotype sequences (hap1 and hap2) as reference sequences. The LOH detection results were orthogonally validated by comparing the read distribution and coverage of the two haplotypes aligned to the target regions (HLA-C and HLA-DPB1 loci).

[0083] Figure 4 A and Figure 4 The left subplot of B shows the alignment results of the tumor sample: In the HLA-C and HLA-DPB1 locus regions, the number of sequencing reads aligned to hap1 was significantly greater than the number aligned to hap2, and the mean coverage of hap1 was significantly higher than that of hap2. This uneven distribution of coverage indicates that in this tumor sample, the allele corresponding to hap2 experienced copy number loss, while the allele corresponding to hap1 was preserved, consistent with the aforementioned LOH detection conclusion based on short read data.

[0084] Figure 4 A and Figure 4 The right subplot of B shows the alignment results of paired normal samples (peripheral blood leukocytes): in the same HLA-C and HLA-DPB1 locus regions, the number of sequencing reads and coverage levels aligned to hap1 and hap2 are basically the same, and no significant imbalance was observed, which is consistent with the expected distribution pattern of normal diploid tissue.

[0085] To further quantify and validate the above observations, this method calculated the coverage ratio (hap1 coverage / hap2 coverage) of each heterozygous site in normal and tumor samples, and performed a paired rank-sum test. The results showed that ( Figure 4 C), the ratio in tumor samples deviated significantly from 1 (P<0.05), while the ratio in normal samples was close to 1 and there was no statistically significant difference. This further confirmed from a statistical perspective that tumor samples had allele-specific loss of heterozygosity at the HLA-C and HLA-DPB1 loci.

[0086] The above-mentioned long-read sequencing validation results not only confirm the accuracy of the LOH detection results based on short-read sequencing data, but also show that the method can effectively cover HLA-I genes (represented by HLA-C) and HLA-II genes (represented by HLA-DPB1) in the MHC region, demonstrating the versatility of high-resolution LOH analysis in the entire MHC region.

[0087] 4. Long-read sequencing validates the high-resolution LOH region delimitation capability of this method. To further verify the accuracy of this method in defining the boundary of LOH events, the resolution of the aforementioned detection results was validated using PacBio HiFi long-read sequencing data. The results are as follows: Figure 5 As shown.

[0088] Figure 5The figure shows the coverage distribution curves of normal samples (top) and tumor samples (bottom) after PacBio HiFi long-read sequencing, with the sequencing reads aligned to hap1 and Hap2, respectively, at various genomic locations. The orange highlighted areas in the figure represent the LOH regions (bump regions) detected by this method in short-read data.

[0089] In normal samples, the coverage curves of hap1 and hap2 are basically consistent throughout the MHC region, with no obvious coverage deviation, which is consistent with the expectation that two homologous chromosomes exist in equal amounts in normal diploid tissue.

[0090] In the tumor sample, within the orange-highlighted LOH region, the coverage curves of hap1 and hap2 showed a significant separation: the coverage of hap1 was significantly higher than that of hap2, indicating that the allele corresponding to hap2 had been lost within this region, consistent with the results of the short-read data. Of particular note, in the non-target regions on either side of this orange LOH region (white background areas in the figure), the coverage curves of hap1 and hap2 reverted to the same pattern, with no significant deviation observed.

[0091] The above results demonstrate that the LOH region boundaries identified by this method closely match the coverage difference region boundaries independently validated by long-read sequencing. Due to its read length advantage, long-read sequencing technology can completely traverse the complex structural regions of the MHC, avoiding the alignment ambiguity problems of short-read sequencing, and its coverage analysis results have higher reliability. The LOH region boundaries detected by this method on short-read data are precisely consistent with the boundaries of the independently validated results from third-generation long-read sequencing, fully demonstrating that this method has the ability to define high-resolution boundaries of LOH events in the MHC region.

[0092] Example 3: Application of a high-resolution MHC region variation analysis method in highly accurate identification of somatic mutations. S1: Two haplotype sequences from the subject's MHC region were obtained using the method published by patent CN110777192A.

[0093] S2: Preprocessing and deduplication based on standard reference: Obtain short-read sequencing data of target tumor samples; First, use alignment tools (razers3 or BWA-MEM) to align the sequencing data to a standard reference genome (such as GRCh38). Based on the initial alignment results, identify and remove PCR duplicates to obtain high-quality deduplicated sequencing reads.

[0094] S3: Construction of the Joint Reference Sequence and Competitive Alignment Filtering: To eliminate alignment bias in the public reference genome, the two haplotype sequences generated in step S1 were merged to construct a double haplotype joint reference sequence. Subsequently, the high-quality sequencing reads after deduplication in S2 were remapped to the double haplotype joint reference sequence, and extremely strict sequence specificity thresholds were set: (I) the identity between the read and the double haplotype joint reference sequence must be ≥95%; (II) the read must have only one optimal matching position in the joint reference sequence composed of the two haplotypes (i.e., excluding multi-mapped reads that can be aligned to both haplotypes simultaneously and reads that cannot be aligned). After this screening, only high-specificity reads that can be clearly attributed to a certain haplotype were retained for subsequent analysis.

[0095] S4: Somatic mutation detection: Based on the sequencing data that has undergone double filtering (duplicate removal and extremely high specificity mapping), a somatic mutation detection algorithm (such as Strelka mutation retrieval tool) is applied to calculate and output a set of high-confidence somatic mutations in the MHC region of the target sample.

[0096] Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Anyone skilled in the art can make various modifications and alterations without departing from the spirit and scope of the present invention. Therefore, the scope of protection of the present invention should be determined by the claims.

Claims

1. A method for LOH variation analysis based on individualized MHC haplotypes, characterized in that, Includes the following steps: S1. Obtain two individualized assembled pseudo-haplotype sequences of the target region of the subject; perform double sequence alignment on the two pseudo-haplotype sequences to identify and extract subject-specific high-density heterozygous variant sites; S2. Obtain short-read sequencing data of the target sample and control sample from the subject; independently align the sequencing data of the tumor cells and the sequencing data of the control sample to the two pseudo-haplotype sequences described in step 1; based on the haplotype-specific alignment results, calculate the comprehensive copy number at each genomic position of the target sample (S2). CN combined ) and calculation of second allele frequency (BAF); S3, based on the above CN combined The statistical distribution of biallelic deletions is used to calculate the lower limit threshold (LB) for biallelic deletion detection, thereby identifying and filtering biallelic deletion regions in the target sample; wherein, the formula for calculating the lower limit threshold LB is: ， in, MHC region CN combined The mean of σ is its standard deviation; S4. In the remaining region after excluding heterozygous variant sites within the biallelic deletion region, combined with the... CN combined BAF and tumor purity ρ Calculate the specific copy number of allele 1 ( ) and the specific copy number of allele 2 ( CN 2) Clustering algorithms are used to smoothly segment heterozygous variant sites, and loss-of-heterozygosity (LOH) regions are determined based on allele copy number differences; wherein, the following equations are used to calculate and CN 2: ， 。 2. The LOH variation analysis method according to claim 1, characterized in that, In S1, heterozygous variant sites are identified by performing double sequence alignment on two haplotype sequences and applying the following filtering conditions: (I) alignment length ≥ 1500 bp; (II) in multiple alignment regions, only the alignment result with the highest consistency score is retained.

3. The LOH variation analysis method according to claim 1, characterized in that, The target sample refers to sequencing data of a tumor tissue sample from the same individual; the control sample refers to sequencing data of a peripheral blood leukocyte sample or adjacent normal tissue sample from the same individual.

4. The LOH variation analysis method according to claim 1, characterized in that, In S2, the sequencing data of the target sample and the control sample are independently aligned to the two pseudo-haplotype sequences using a sequence alignment tool. The parameters are set as follows: the number of mismatches is set to 0~1, and uniquely mapped read pairs are retained.

5. The LOH variation analysis method according to claim 1, characterized in that, In S2, the CN combined Calculate using the following formula: ， in, Cov tumor This refers to the overall coverage of genomic locations in tumor samples. Cov normal This refers to the overall coverage of genomic locations in the control sample, where M is the ratio of the size of the control sequencing library to the size of the tumor sequencing library. ρ For tumor purity; The overall coverage refers to the sum of coverage mapped to each haplotype minus the shared coverage of equivalent mappings on two pseudo-haplotypes; Meanwhile, within a window of 75 bp centered on the high-density heterozygous variant site, the minor allele frequency (BAF) is calculated based on the allele-specific read that is uniquely aligned to a single pseudohaplotype.

6. The LOH variation analysis method according to claim 1, characterized in that, In S3, the confirmation of the biallelic deletion region must meet the following requirements: (I) CN combined <LB; (II) Continuous region length > 25 kb; (III) No assembly gap exists within the region; the assembly gap is defined as a length > 12.5 kb within the continuous fragment and no sequencing read can be aligned to any pseudohaplotype.

7. The LOH variation analysis method according to claim 1, characterized in that, In S4, genomic regions that simultaneously meet the following conditions are selected as regions lacking heterozygosity: (a)| - CN 2|≥Set value, where the set value is used to measure the degree of difference between two haplotype copy numbers; (b) and CN The smaller of the two values ​​is below the preset lower limit for the number of copies; Preferably, the set value is 0.5; more preferably, it is 0.

7. Preferably, the lower limit of the copy number is 0.

5.

8. The LOH variation analysis method according to claim 1, characterized in that, When the major allele copy number is ≤1.5, it is determined to be a deletion-type LOH; when the major allele copy number is >1.5, it is determined to be a copy-neutral LOH.

9. A method for analyzing somatic mutation variations based on individualized MHC haplotypes, characterized in that, Includes the following steps: S1. Obtain two individualized assembled pseudohaplotype sequences from the target region of the subject; S2. Obtain short-read sequencing data of the target tumor sample, align the sequencing data to a standard reference genome, identify and remove PCR duplicate reads based on the alignment results, and obtain deduplicated sequencing reads. S3. Merge the two pseudo-haplotype sequences obtained in S1 to construct a double haplotype joint reference sequence; map the deduplication sequencing reads obtained in S2 to the double haplotype joint reference sequence, and screen reads that meet the following conditions: the sequence identity between the read and the double haplotype joint reference sequence is not less than 95%, and the read has only one optimal matching position in the double haplotype joint reference sequence; S4. Based on the reads obtained after screening, apply the somatic mutation detection algorithm to output the set of somatic mutations in the MHC region of the target tumor sample.

10. The somatic cell mutation and variation analysis method according to claim 9, characterized in that, In step S3, the read segment has only one optimal alignment position in the double haplotype joint reference sequence, which means that the read segment has one and only one optimal alignment result in the joint reference sequence composed of two haplotype sequences, excluding multi-mapped read segments that can be aligned to two haplotypes at the same time and read segments that cannot be aligned.

Images