A method and system for accurate identification of loss of heterozygosity in prenatal diagnosis

Abstract

The application provides a calculation method and system for accurately identifying loss of heterozygosity in prenatal diagnosis, and is applied to the technical field of medical data processing. Bam files and variation information are obtained through whole genome library construction and sequencing analysis through data processing. SNP sites are screened, allele mutation frequency is calculated, and effective sites are screened. The genome window is divided, the homozygosity rate HR value is calculated, and the nonlinear weight reference is constructed based on normal samples. Finally, the improved CBS algorithm is used to identify significant AOH segments combined with the weight reference, generate a table containing position and HR average information, complete the identification of loss of heterozygosity, and do not need additional CMA detection, which is accurate and efficient.

Description

Technical Field

[0001] This invention relates to the field of medical data processing technology, and in particular to a computational method and system for accurately identifying loss of heterozygosity in prenatal diagnosis. Background Technology

[0002] Absence of heterozygosity (AOH) in embryos can lead to miscarriage and a range of birth defects, many of which are undetectable by routine ultrasound examinations. More importantly, AOH significantly increases the risk of genetic diseases caused by homozygous recessive genes. In regions associated with genetic imprinting, AOH can cause uniparental diploidy (UPD), leading to imprinted diseases such as UPD on chromosomes 6, 11, 14, 15, and 20. Typical examples include Prader–Willi syndrome (maternally derived chromosome 15 UPD), Angelman syndrome (paternally derived chromosome 15 UPD), Silver–Russell syndrome (maternally derived chromosome 7 UPD), and Beckwith–Wiedemann syndrome (paternally derived chromosome 11 UPD). Furthermore, female X-chromosome-related latent hemorrhage (LOH) may be closely related to the occurrence of X-linked recessive genetic diseases. If the mother is a carrier of the disease-causing gene and the fetus has LOH, the disease phenotype may manifest. Therefore, early screening for LOH in prenatal diagnosis is crucial for reducing the risk of birth defects.

[0003] Currently, first-line methods for prenatal diagnosis include karyotyping, chromosomal microarray analysis (CMA), and low-depth copy number variation equating (CNV-seq) based on next-generation sequencing (NGS). Karyotyping is the traditional "gold standard" for diagnosing chromosomal aberrations. It mainly detects most numerical abnormalities, large deletions / duplications, and other chromosomal variations by observing chromosome morphology and structure. However, karyotyping is complex to perform, has a long testing cycle, and struggles to detect CNVs smaller than 5Mb. CMA is the most commonly used diagnostic method for AOH individuals. It can detect CNVs and UPDs and has high clinical application value for AOH screening. However, CMA has limited resolution when detecting chromosomal copy number abnormalities, making it difficult to detect low proportions (<30%) of mosaicism. Furthermore, its high cost limits its application in large-scale screening.

[0004] In contrast, CNV-seq offers higher throughput, a wider detection range, and lower cost, effectively detecting low-proportion mosaicism. It can identify >10% of chromosomal aneuploid mosaicism in clinical samples, and under ideal conditions, the detected mosaicism rate can be as low as 5%. Furthermore, CNV-seq has lower sample requirements, accurately detecting as little as 10-50 ng of DNA, making it suitable for prenatal diagnosis or microsample analysis. However, while CNV-seq excels in copy number variation detection, currently used clinically cannot identify aneuploidy of chromosomes (AOH), including undesirable copy number variations (UPDs). Patients often require additional chromosomal microarray analysis (CMA) or STR analysis, increasing both financial burden and time costs.

[0005] It should be noted that the information disclosed in the background section above is only used to enhance the understanding of the background of this disclosure, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0006] The purpose of this application is to provide a computational method and system for accurately identifying loss of heterozygosity in prenatal diagnosis, overcoming at least to some extent the problems of existing technologies. It eliminates the need for additional CMA testing and can accurately identify AOH. The process involves extracting DNA from prenatal samples, constructing a library, sequencing, processing the data to obtain a BAM file, generating sequencing results, obtaining genotype and VCF datasets, screening SNP sites and calculating mutation frequencies, dividing the genomic window, calculating HR values, establishing a weighted reference using normal samples, and using an improved CBS algorithm to divide AOH segments, generating a table to complete the identification, thus improving accuracy and efficiency.

[0007] Other features and advantages of this application will become apparent from the following detailed description, or may be learned in part by practice of the invention.

[0008] According to one aspect of this application, a computational method for accurately identifying loss of heterozygosity in prenatal diagnosis is provided, comprising: extracting genomic DNA from prenatal samples and performing fragmentation, end repair, A-tailing, adapter ligation, PCR amplification, and magnetic bead purification before sequencing to complete whole-genome library construction; removing impurities from the sequencing data using trim_galore to obtain a pure FASTQ file, and obtaining a BAM file after genome alignment, duplicate removal, Indel re-alignment, and base correction; calculating depth and coverage; and using GATK to perform variant identification, filtering, and annotation to generate sequencing data analysis results; re-aligning the sample BAM file to insert the deleted region to obtain genotypes and generating a VCF dataset; using VCFtools to screen for high-quality SNP sites; calculating the proportion of mutant alleles at the SNP sites; and screening for sites with a total read count ≥3. The algorithm generates allele mutation frequency calculation and screening results. The entire genome is divided into 0.1Mb ​​windows. The number of valid SNP sites with a total read count ≥3 within each window is counted, and low-quality windows are excluded. The HR value of the retained windows is calculated, and the homozygosity rate of each window is calculated. CNV-seq data from normal samples validated by CMA are processed to obtain the mean HR value for each window. A non-linear weighted reference is constructed, and the median is used to fill in missing value windows. This is integrated to obtain a weighted reference list, generating non-linear weighted reference results. The HR data is input into the target CBS algorithm and introduced into the weighted reference list. By calculating the significant changes in HR values ​​between windows and recursively segmenting significant AOH regions, a table containing the chromosome location, start and end positions, mean HR value, and number of valid bins within the region for AOH fragments is generated, completing the identification of heterozygosity loss.

[0009] Another aspect of this application discloses a computational device for accurately identifying loss of heterozygosity in prenatal diagnosis, comprising: an acquisition module for extracting genomic DNA from prenatal samples and performing fragmentation, end repair, A-tailing, adapter ligation, PCR amplification, and magnetic bead purification before sequencing to complete whole-genome library construction; a processing module for removing impurities from the sample data using trim_galore to obtain a clean FASTQ file, performing genome alignment, duplicate removal, Indel re-alignment, and base correction to obtain a BAM file, calculating depth and coverage, and then using GATK to perform variant identification, filtering, and annotation to generate sequencing data analysis results; re-aligning the sample BAM file to insert the missing region to obtain genotypes and generate a VCF dataset, using VCFtools to screen for high-quality SNP sites, calculating the proportion of mutated alleles at the SNPs, and screening the total reads. For loci with ≥3 reads, allele mutation frequency calculation and screening results are generated. The whole genome is divided into 0.1Mb ​​windows. The number of effective SNP loci with ≥3 reads in each window is counted and low-quality windows are excluded. The HR value of the retained windows is calculated, and the homozygosity of the windows is calculated. CNV-seq data of normal samples verified by CMA are processed to obtain the mean HR value of each window. A nonlinear weighted reference is constructed and the median is used to fill the missing value windows. The weighted reference list is integrated to generate the nonlinear weighted reference result. The HR data is input into the target CBS algorithm and introduced into the weighted reference list. By calculating the significant change of HR value between windows and recursively segmenting significant AOH segments, a table of information containing the chromosome position, start and end positions, mean HR value and number of effective bins in the region of AOH segments is generated to complete the identification of heterozygosity loss.

[0010] According to another aspect of this application, an electronic device includes: a first processor; and a memory for storing executable instructions of the first processor; wherein the first processor is configured to execute the executable instructions to perform the above-described computational method for accurately identifying loss of heterozygosity in prenatal diagnosis.

[0011] According to another aspect of this application, a computer-readable storage medium is provided having a computer program stored thereon, which, when executed by a second processor, implements the above-described computational method for accurately identifying loss of heterozygosity in prenatal diagnosis.

[0012] This application provides a computational method and system for accurately identifying loss of heterozygosity (LOH) in prenatal diagnosis. The method involves extracting genomic DNA from prenatal samples, processing it, and then sequencing it to build a library. After sequencing, the data is processed to obtain a BAM file, and sequencing analysis results are generated using GATK. Genotypes are obtained to generate a VCF dataset, SNP loci are screened, and allele mutation frequencies are calculated. The genome is divided into 0.1Mb ​​windows, and the retention window (HR) value is calculated. A non-linear weighted reference is constructed using normal samples validated by CMA. The HR data is input into an improved CBS algorithm, combined with the weighted reference, to divide AOH regions and generate a table for identification. Multi-step processing and algorithm optimization improve accuracy and efficiency.

[0013] It should be understood that the above general description and the following detailed description are exemplary and explanatory only, and are not intended to limit this disclosure. Attached Figure Description

[0014] Figure 1 A flowchart illustrating a calculation method for accurately identifying loss of heterozygosity in prenatal diagnosis provided in an embodiment of this application is shown.

[0015] Figure 2 This is a schematic diagram showing the weight values ​​for linear and nonlinear mappings.

[0016] Figure 3 The whole-genome HR value map of positive control sample C1105 after visualization in this application;

[0017] Figure 4 The whole-genome HR value map of negative control sample C20230189 after visualization in this application;

[0018] Figure 5 A visualization of chromosome 6 HR values ​​for sample C749 in this application;

[0019] Figure 6 This illustration shows a schematic diagram of a computational device for accurately identifying loss of heterozygosity in prenatal diagnosis, provided in an embodiment of this application. Detailed Implementation

[0020] The preferred embodiments of the present invention will be described below with reference to the accompanying drawings. It should be understood that the preferred embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.

[0021] The following is combined with Figure 1 This application describes a computational method for accurately identifying loss of heterozygosity in prenatal diagnosis according to an exemplary embodiment of the present application. It should be noted that the following application scenarios are shown only to facilitate understanding of the spirit and principles of the present application, and the embodiments of the present application are not limited in any way. Rather, the embodiments of the present application are applicable to any suitable scenario.

[0022] In one embodiment, this application also proposes a computational method and system for accurately identifying loss of heterozygosity in prenatal diagnosis. Figure 1 A schematic flowchart illustrating a calculation method for accurately identifying loss of heterozygosity in prenatal diagnosis according to an embodiment of this application is shown.

[0023] S101 involves extracting genomic DNA from prenatal samples and then performing fragmentation, end repair, A-tailing, adapter ligation, PCR amplification, and magnetic bead purification before sequencing to complete the whole genome library construction.

[0024] In one embodiment, genomic DNA is extracted from prenatal samples (such as amniotic fluid, umbilical cord blood, and aborted chorionic villus tissue). The specific process is as follows: a conventional genomic DNA extraction kit is used to isolate pure genomic DNA from the collected prenatal samples, ensuring DNA purity (OD260 / 280≈1.8-2.0) and integrity (agarose gel electrophoresis shows no obvious degradation).

[0025] Genomic DNA is fragmented using nonspecific endonucleases (such as Tn5) to cut it into fragments approximately 300-500 bp in length to fit sequencing library requirements. The sticky ends of the fragmented DNA are then repaired using DNA polymerase, transforming them into blunt ends to prepare for subsequent tailing and adapter ligation. A single adenosine nucleotide (A) is added to the 3' end of the repaired DNA fragment, forming a 3' protruding A-tail structure, which facilitates ligation to sequencing adapters with T tails.

[0026] Sequencing adapters with specific indices are ligated to A-tailed DNA fragments. The adapters contain sequencing primer binding sites and sample identification sequences for sample differentiation during subsequent sequencing. The ligation product is then amplified by PCR using primers adapted to the adapter sequences. The amplification cycle number is controlled at 10-12 to reduce amplification bias and obtain a sufficiently concentrated library product.

[0027] PCR amplification products were purified using magnetic beads to remove unligated adapters, primer dimers, and short fragment impurities, ultimately yielding sequencing libraries with a concentration ≥20 ng / μL and a fragment distribution consistent with expectations. The purified libraries were then mixed in a specific ratio and loaded onto a sequencer (such as an Illumina platform) for whole-genome sequencing, completing the whole-genome library construction process.

[0028] S102 uses trim_galore to remove impurities from the sequencing data to obtain a pure FASTQ file. After genome alignment, deduplication, Indel re-alignment, and base correction, a bam file is obtained. Depth and coverage are statistically analyzed, and GATK is used to complete variant identification, filtering, and annotation to generate sequencing data analysis results.

[0029] In one implementation, the offline sequencing data is processed using `trim_galore` to remove end sequencing adapters, amplification primers, and low-quality bases, generating a clean FASTQ file. The `trim_galore` tool (with parameters set to `--quality20--phred33--stringency3`) processes the raw offline sequencing data, automatically identifying and deleting end sequencing adapter sequences, residual PCR amplification primer fragments, and bases with a quality value below 20, ultimately generating a clean FASTQ file free of impurities, ensuring the accuracy of subsequent data analysis.

[0030] The clean FASTQ file is aligned with a reference genome. After steps including read removal, indel re-alignment, and base quality correction, a BAM file is generated. The clean FASTQ file is then aligned with a human reference genome (e.g., hg38) using an alignment tool (e.g., BWA) to obtain preliminary alignment results. Subsequently, Picard is used to remove duplicate reads generated by PCR amplification. GATK is used for indel re-alignment to optimize the alignment accuracy of inserted and deleted regions. Finally, GATK's Base Recalibrator is used for base quality correction, resulting in a multi-optimized BAM file.

[0031] The sequencing depth and genome coverage of the BAM files were statistically analyzed to obtain depth and coverage information. The samtools tool was used to perform statistical analysis on the BAM files, calculating the sequencing depth of the whole genome (the average sequencing depth of the samples in this method is approximately 0.67-fold) and genome coverage.

[0032] The BAM file was processed using GATK software for variant identification, filtering, and annotation to generate a dataset containing variant site information. The HaplotypeCaller module of GATK was used for variant calling, yielding a raw variant dataset containing SNPs and indels. The VariantFiltration module of GATK was used for filtering (with filtering parameters set to DP<5, QD<2.0, FS>60.0, MQ<40.0) to retain high-quality variant sites. The ANNOVAR tool was then used to annotate the filtered variant sites, adding information such as gene function and population frequency, generating a dataset with detailed variant site information. Finally, the bcftools tool was used to extract the location information, wild-type bases, mutant bases, and their corresponding read counts (i.e., AD values) of all filtered SNP sites, which served as the input matrix for subsequent calculation modules.

[0033] The above processing results are integrated to generate complete sequencing data analysis results. First, the quality report of the clean FASTQ files is summarized. This report includes information such as data integrity and base quality distribution (e.g., Q20 and Q30 ratios) after trim_galore processing, reflecting the basic quality of the original sequencing data after impurity removal. Second, the sequencing depth and genome coverage statistics of the bam files are included, including the average sequencing depth of the whole genome (e.g., 0.67-fold), the depth distribution of each chromosome, and the proportion of the genome that is effectively covered, to assess the uniformity and coverage integrity of the sequencing data. Then, the high-quality variant site dataset and its annotation information obtained through GATK screening are integrated. The variant site dataset contains basic information such as the location and genotype of the SNPs and Indels retained after filtering, while the annotation information covers the genes and functional regions (e.g., coding and non-coding regions) where the sites are located.

[0034] Through the above integration, the resulting complete sequencing data analysis results not only include a comprehensive assessment of data quality (such as whether the sequencing depth meets the requirements of subsequent analysis and whether the coverage meets the standards), but also cover the accuracy information of sequence alignment (such as the alignment quality after Indel re-alignment and base correction) and the detailed characteristics of variant sites. This information together provides reliable basic data support for the screening of SNP sites (such as screening sites with a total number of reads ≥3) and the identification of loss of heterozygosity (AOH) in subsequent steps, ensuring that subsequent analysis is carried out on the basis of high-quality data.

[0035] S103, re-align the sample BAM files, insert missing regions to obtain genotypes and generate VCF datasets. After screening high-quality SNP sites using VCFtools, calculate the proportion of mutated alleles at the SNPs, screen sites with a total read count ≥3, and generate allele mutation frequency calculation and screening results.

[0036] In one implementation, GATK software is used to re-align the insertion and missing regions in the sample BAM file to obtain the sample genotypes and generate a VCF format dataset, where each variant site has corresponding genotype information. GATK software (with parameters set to -stand_call_conf30-stand_emit_conf10) is used to re-align the insertion and missing regions in the sample BAM file to optimize the accuracy of genotype identification. This step obtains the genotype information of the sample, including the homozygous or heterozygous status of each variant site, and generates a VCF format dataset, where each variant site (e.g., SNP, Indel) corresponds to a specific genotype label (e.g., AA, AB, BB).

[0037] The VCF dataset was filtered using vcftools to retain high-quality SNP loci, generating a filtered SNP locus dataset. By setting quality thresholds (e.g., QUAL≥30, DP≥3), high-quality SNP loci were retained, while low-quality variants (e.g., loci with insufficient sequencing depth or low genotype confidence) were excluded, thus generating a filtered SNP locus dataset and ensuring the reliability of subsequent analyses.

[0038] The selected SNP loci were calculated using R software, based on the formula. The allele mutation frequency of each SNP locus was obtained, generating a dataset of SNP locus allele mutation frequency values. This represents the proportion of mutated alleles at this SNP. This is a mutation-type read. This represents the wild-type read. The AF value reflects the proportion of the mutant allele at that SNP (ranging from 0 to 1). For example, if a certain SNP site has 2 wild-type reads and 1 mutant read, then AF = 1 / (2+1) ≈ 0.33. This step generates a dataset containing the AF values ​​for each SNP site.

[0039] The SNP allele mutation frequency dataset was filtered to retain sites with at least 3 wild-type reads plus at least 3 mutant reads, generating a dataset of SNP sites and their allele mutation frequencies that met the criteria. The core purpose of this filtering was to exclude low-depth sites—when the total number of reads is less than 3, random fluctuations in sequencing data may lead to biases in the calculation of AF values ​​(for example, with only one read, AF can only be 0 or 1, failing to accurately reflect the genotype), thus interfering with the accuracy of subsequent window homozygosity (HR).

[0040] This step retains only SNP sites with a total read count of ≥3 and their corresponding AF values. For example, a site may have 1 wild-type read and 2 mutant reads (total read count = 3), or 3 wild-type reads and 1 mutant read (total read count = 4). This generates a dataset of SNP sites and their AF values ​​that meet the criteria, ensuring that the sites included in subsequent analyses have sufficient sequencing depth to support the results and reducing the impact of random errors.

[0041] Integrating the above processing results, the allele mutation frequency calculation and screening results are generated, including the screened SNP loci and their allele mutation frequency values. The results include a high-quality SNP locus list screened by vcftools (ensuring locus reliability), genotype information and read count for each SNP locus (from the vcf dataset after GATK re-alignment), allele mutation frequencies (AF values, reflecting the proportion of mutant alleles) calculated using R software, and the total read count for each locus (the sum of wild-type and mutant reads, with loci having been screened for a total read count ≥3).

[0042] By integrating this information, the resulting complete result not only includes the basic characteristics of each SNP locus (such as chromosomal location and genotype), but also covers quantified mutation frequency data and quality indicators (total number of reads), ensuring data traceability and reliability. This result directly provides the core input for calculating the homozygosity rate (HR) after genome windowing in subsequent steps—the calculation of the HR value depends on the number (n) of SNP loci with AF=0 or AF=1 in each window and the total number of SNP loci (N). The high-quality SNP loci and their AF values ​​selected in this step are the basis for accurately calculating n and N.

[0043] S104 divides the entire genome into 0.1Mb ​​windows, counts the number of valid SNP sites with a total read count ≥3 in each window and excludes low-quality windows, calculates the HR value of the retained windows, and generates the window homozygosity calculation results.

[0044] In one implementation, the R package GenomicRanges is used to systematically divide the whole genome into windows, generating non-overlapping windows of size 0.1Mb, and all windows are globally numbered across chromosomes (e.g., bin1, bin2...binN). Each window is identified by a unique number to ensure accurate location in subsequent analyses.

[0045] Statistical processing is performed on SNP loci within each window to filter out valid SNP loci with a total read count ≥ 3, generating information on the number of valid SNP loci for each window. Specifically, each SNP locus within each window is statistically analyzed, and valid SNP loci with a total read count (wild-type reads + mutant reads) ≥ 3 are filtered out, excluding low-depth loci to reduce errors. For example, if a window contains 12 SNP loci, and 8 of them have a total read count ≥ 3, then the number of valid SNP loci for that window is 8, and an information table containing the number of valid SNP loci for each window is generated.

[0046] Based on the number of valid SNP loci, low-quality windows with less than 5 valid SNP loci are excluded (these windows may lead to HR value calculation errors due to insufficient loci), and only windows with 5 or more valid SNP loci are retained. For example, the window with 8 valid SNP loci meets the criteria and is retained, while the window with 3 valid SNP loci is excluded, ultimately generating a high-quality window dataset.

[0047] For high-quality windows, use the formula Calculate the HR value and generate HR value data for each window, where N≥5, n is the number of SNP sites in the window with AF=0 or AF=1 (AF=0 indicates that all reads are wild-type, i.e., wild-type homozygous; AF=1 indicates that all reads are mutant, i.e., mutant homozygous, and both types of sites are considered homozygous sites), and N is the total number of SNP sites in the window (N≥5 is required, i.e., the window retains ≥5 valid SNP sites after screening to ensure statistical reliability).

[0048] During the calculation, the homozygosity of a genomic region is quantified by the ratio of homozygous sites (n) to the total effective SNP sites (N) within a statistical window. For example, if there are 12 total effective SNP sites (N) within a window, with 5 sites having AF=0 and 3 sites having AF=1 (n=5+3=8), then the HR value for this window is approximately 8 / 12 ≈ 0.67. This value reflects the proportion of homozygous sites within the 0.1 Mb window; the closer the value is to 1, the higher the homozygosity of the region, and the more likely it is to contain loss of heterozygosity (AOH).

[0049] By calculating the HR value for each high-quality window, HR value data for each window across the entire genome is generated, providing key quantitative indicators for subsequent construction of nonlinear weighted references based on normal samples and identification of AOH regions using the improved CBS algorithm.

[0050] Integrate high-quality windows and their HR values ​​to generate window homozygosity calculation results, where the key is the window identifier and the value is the HR value of the corresponding window. When generating window homozygosity calculation results, it is necessary to integrate the key information of all selected and retained high-quality windows (with ≥5 effective SNP loci), specifically including: the unique identifier of each window (i.e., a global cross-chromosome number, such as bin1, bin2…binN, generated by the R package GenomicRanges when dividing windows), and information obtained through the formula… The calculated corresponding window HR value (e.g., 0.7, 0.85, etc.).

[0051] The integrated results are presented in key-value pair format, where the "key" is the unique identifier of the window, ensuring that each window can be accurately located and traced; the "value" is the HR value of that window, quantifying the homozygosity of the 0.1Mb ​​genomic region. This structured result not only achieves a one-to-one correspondence between windows and HR values, but also provides direct data support for subsequent steps: on the one hand, when constructing a nonlinear weighted reference based on normal samples, the mean (μ_control) of the HR values ​​of each window in the normal samples needs to be calculated, and then weights are generated through the tanh function; on the other hand, when the improved CBS algorithm identifies AOH regions, the HR value serves as the core input data, combined with the weighted reference list to calculate the significant changes in HR values ​​between windows, thereby achieving the division of significant AOH segments.

[0052] S105. Select normal sample CNV-seq data verified by CMA for processing, obtain the mean HR of each window, construct a non-linear weight reference and fill the missing value window with the median, integrate to obtain the weight reference list, and generate the non-linear weight reference result.

[0053] In one implementation, CNV-seq data of selected CMA-validated normal samples are processed to generate HR value data for each window, with each window having a unique identifier. CNV-seq data from 15 CMA-validated normal samples are selected, and the same processing procedure (including SNP screening, window division, HR value calculation, etc.) is followed as for the samples to be tested to generate HR value data for each sample. Each window is uniquely identified by a globally cross-chromosome number (e.g., bin1, bin2…binN) to ensure a one-to-one correspondence with the windows of the samples to be tested.

[0054] The mean HR value for each window is calculated, and missing value windows are imputed using the median mean HR value of normal samples, generating a complete window HR mean dataset. For each window (matched by unique identifier), the mean HR value (μ_control) of 15 normal samples is calculated to form the HR mean data for each window. For missing windows where no HR value is found in all normal samples (e.g., certain genomic regions where no valid SNPs are detected in normal samples), imputed using the median mean HR value of all windows in normal samples, ultimately generating a complete window HR mean dataset.

[0055] Through formula Calculate weight values ​​for each window to generate window weight data, where μ_controlj is the mean HR of that window, and the output range of the tanh function is [0, tanh(1)≈0.76], therefore the weight value range is [0.24, 1]. For example, if the mean HR of a window is μ_control=0.5, then the weight value = 1-tanh(0.5)≈1-0.462=0.538. This weight data provides the basis for the subsequent construction of the weighted likelihood ratio statistic. This weight calculation principle can be achieved through... Figure 2 (A diagram illustrating the weight values ​​of linear and nonlinear mappings) This visually demonstrates that the smooth decay characteristics of the tanh function in the high HR region can reduce interference in the normal high HR region.

[0056] By integrating window identifiers with corresponding weight values, a weight reference list is generated, which improves the weighted likelihood ratio statistic in the CBS algorithm. The calculation provides key inputs and generates nonlinear weighted reference results, where the key is the window identifier and the value is the weight value of the corresponding window and the basic parameters required for the weighted likelihood ratio calculation. This is a specific value for a particular window k in the weight reference list. Ensure that the unique identifier of each window (global cross-chromosome number, such as bin1, bin2...binN) corresponds one-to-one with the previously divided 0.1Mb ​​genome windows, guaranteeing a precise association between the weight value and the window. The list is presented in a structured key-value pair format, where the "key" is the unique number of the window, used for quickly locating specific genomic regions in subsequent analysis; the "value" contains two core pieces of information: one is the value obtained through the formula... The calculated window weight value (e.g., a window weight value of 0.538) and the basic parameters required for calculating the weighted likelihood ratio statistic include the variance of the whole-genome HR value. — This variance is obtained through statistical analysis of the HR value distribution of normal samples and is used to quantify the overall homozygosity fluctuation level of the genome.

[0057] Weight values ​​for a specific window k It serves as the weighted likelihood ratio statistic in the improved CBS algorithm. The core input directly affects the size of the statistic: the higher the weight value, the lower the HR value (low homozygosity) of the window in the normal sample, and the greater the contribution of the change in its HR value to AOH identification; conversely, the influence of windows with lower weight values ​​(such as regions with naturally high HR in normal samples) on the statistic is weakened, thereby reducing false positives.

[0058] The nonlinear weight reference results generated through this integration method achieve both precise mapping between windows and weights and provide standardized basic parameters for the CBS algorithm, ensuring that the algorithm can dynamically adjust the contribution of different genomic regions. The construction logic of this weight reference list can be achieved through... Figure 2 (Illustrative diagram of weight values ​​for linear and nonlinear mapping) The diagram clearly shows how the tanh function uses nonlinear mapping to distribute weights to different HR mean windows, avoiding information loss caused by extreme weights in linear mapping.

[0059] S106. Input the HR data into the target CBS algorithm and introduce the weight reference list. By calculating the significant changes in HR values ​​between windows and recursively segmenting significant AOH segments, generate a table containing the chromosome position, start and end positions, mean HR value, and number of effective bins in the region of the AOH segment to complete the identification of loss of heterozygosity.

[0060] In one implementation, the window HR value data of the samples (e.g., HR values ​​of 0.8, 0.92, etc. for each window) are input into the improved Cyclic Binary Segmentation (CBS) algorithm, along with the previously constructed weight reference list Weight_ref. The HR value of each window is associated with a corresponding weight value (e.g., 0.538) through a unique identifier (e.g., bin123), ensuring that the algorithm can dynamically adjust the contribution of different windows during segmentation—the HR value changes of high-weight windows (regions with low HR values ​​in normal samples) have a greater impact on the segmentation result, while the impact of low-weight windows (regions with high HR values ​​in normal samples) is moderately weakened.

[0061] The significance of the change in HR values ​​between windows was calculated using the CBS algorithm, based on the weighted likelihood ratio statistic. Recursive segmentation is performed to divide significant AOH contiguous regions. Specifically, the algorithm first takes the HR values ​​of each window across the entire genome as input, and combines them with the corresponding window weight values ​​from the weight reference list Weight_ref. Construct a weighted likelihood ratio statistic .in, This represents the mean HR from the first window to the kth window. This represents the mean HR from the (k+1)th window to the last window, and the difference between the two reflects the overall difference in HR values ​​on both sides of the cutoff point k. The variance of the whole genome HR value is obtained by statistical analysis of the HR value distribution of normal samples and is used to standardize the difference to eliminate the interference of overall genome fluctuations on local differences.

[0062] As the weight value of window k (from the formula) The algorithm calculates and weights the statistic. Windows with lower HR values ​​(higher weight) in the normal sample contribute more to the statistic from changes in HR values, while windows with higher HR values ​​(lower weight) have a weaker contribution, thus reducing interference from the normal high HR region. Secondly, the algorithm determines the significance of the statistic through a test: when... When the threshold value exceeds the set threshold (α=0.01, adjusted by Bonferroni to control for multiple test error), it indicates that there is a statistically significant difference in the HR values ​​on both sides of the split point k. At this point, the genome is split into two sub-regions.

[0063] Finally, repeat the above process (calculating potential segmentation points) for each of the segmented sub-regions. The algorithm iterates through several steps, including significance testing and sub-region segmentation, until no significant difference in HR values ​​is found within any of the sub-regions, ultimately identifying continuous significant AOH regions. This recursive segmentation method can accurately identify the boundaries of AOH regions. For example, after processing, sample C1105 was divided into multiple significant AOH regions across the entire genome (corresponding to...). Figure 3 The distribution of whole-genome HR values ​​was consistent with the whole-genome AOH results verified by CMA.

[0064] Information is extracted from the segmented AOH regions to generate tabular data containing chromosome location, start and end positions, mean HR value, and the number of effective bins within the region. For significant AOH regions, information is extracted to generate tabular data containing the following core information: chromosome location (e.g., chr2), start and end positions of the region (e.g., 100001-89200001), mean HR value of the region (e.g., 0.96), and the number of effective windows within the region (num.mark, e.g., 25 0.1Mb ​​windows). This information accurately reflects the genomic location and homozygosity of the AOH. For example, two AOH regions were detected on chromosome 2 of sample C814, with mean HR values ​​of 0.96 and 0.97, respectively.

[0065] The AOH region tabular data is visualized to present the distribution characteristics of AOH regions in the genome. This visualization is achieved through genome-wide HR value distribution maps or specific chromosome HR value maps. For example, the genome-wide HR value map of sample C1105 (…). Figure 3 The results showed multiple regions on all chromosomes with HR values ​​close to 1, suggesting genome-wide AOH. Figure 4 This is a visualization of the whole-genome HR value of the negative control sample C20230189, showing the distribution characteristics of AOH regions through the whole-genome HR value distribution. The figure shows that the overall HR value of sample C20230189 is at a low level throughout the whole genome, with no regions showing multiple HR values ​​close to 1, and no significant AOH segment distribution, consistent with the characteristics of a negative control sample, i.e., no loss of heterozygosity was detected. This visualization result is consistent with the expectation that the negative control sample has no AOH, further validating the effectiveness of this calculation method in distinguishing between normal samples and samples with AOH. (HR value map of chromosome 6 of sample C749 is also shown.) Figure 5 The results clearly show the distribution of high HR values ​​in specific segments of the chromosome, consistent with the CMA validation results.

[0066] By integrating tabular data and visualization results, the identification of loss of heterozygosity is completed. The key represents the AOH segment identifier, and the value represents detailed information about the corresponding segment and related data from the visualization results. The final result for identifying loss of heterozygosity is formed by integrating the tabular data and visualization results of the AOH segments. The result is presented in key-value pairs—the key is the unique identifier of the AOH segment (e.g., chr2_1), and the value is detailed information about the segment (chromosomal location, start and end positions, mean HR value, number of effective bins) and related data from the visualization chart. For example, two AOH segments were identified on chromosome 16 of sample C822, with overlap rates of 98.46% and 99.13% with the CMA results, respectively, confirming the accuracy of the identification and thus completing the precise identification of loss of heterozygosity.

[0067] This application enables accurate and efficient identification of loss of heterozygosity (AOH) without additional CMA testing. Genomic DNA is extracted from prenatal samples, fragmented, and end-repaired before sequencing to complete whole-genome library construction. The sequencing data is processed using trim_galore to obtain a clean FASTQ file, followed by genome alignment to obtain a BAM file. GATK is then used for variant identification to generate sequencing data analysis results. The sample BAM files are re-aligned to obtain genotypes and generate a VCF dataset. High-quality SNP loci are screened, and effective loci are calculated and selected to obtain allele mutation frequencies. The whole genome is divided into 0.1Mb ​​windows, and the homozygosity rate (HR) of the retained windows is calculated after excluding low-quality windows. A non-linear weighted reference is constructed using CMA-validated normal sample data, and missing values ​​are filled to create a weighted reference list. The HR data is input into an improved CBS algorithm, combined with the weighted reference list, to segment significant AOH regions and generate a table containing information such as chromosome location, completing the identification. Through multi-step data processing and algorithm optimization, the accuracy and efficiency of AOH identification are improved.

[0068] In one implementation, such as Figure 6 As shown, this application also provides a computational device for accurately identifying loss of heterozygosity in prenatal diagnosis, comprising:

[0069] The acquisition module 601 is used to extract genomic DNA from prenatal samples and perform fragmentation, end repair, A-tailing, adapter ligation, PCR amplification and magnetic bead purification before sequencing to complete the whole genome library construction.

[0070] Processing module 602 is used to remove impurities from the offline data using trim_galore to obtain a pure FASTQ file. After genome alignment, deduplication, Indel re-alignment, and base correction, a BAM file is obtained. Depth and coverage are calculated, and GATK is used to perform variant identification, filtering, and annotation to generate sequencing data analysis results. The module then re-aligns the sample BAM file, inserts missing regions to obtain genotypes, and generates a VCF dataset. After screening high-quality SNP sites using VCFtools, the proportion of mutant alleles at the SNPs is calculated, and sites with a total read count ≥3 are selected, generating allele mutation frequency calculation and screening results. The entire genome is divided into 0.1Mb ​​windows, and the total read count within each window is calculated. The number of valid SNP loci with ≥3 reads was counted, and low-quality windows were excluded. The HR values ​​of the retained windows were calculated, and the homozygosity calculation results were generated. CNV-seq data of normal samples verified by CMA were selected for processing to obtain the mean HR value of each window. A nonlinear weighted reference was constructed and the median was used to fill the missing value windows. The weighted reference list was integrated to generate the nonlinear weighted reference results. The HR data was input into the target CBS algorithm and introduced into the weighted reference list. By calculating the significant change of HR values ​​between windows and recursively segmenting significant AOH segments, a table of information containing the chromosome position, start and end positions, mean HR value, and number of valid bins in the region of AOH segments was generated to complete the identification of heterozygosity loss.

[0071] The computer-readable storage medium provided in the above embodiments of this application and the calculation method for accurately identifying heterozygosity loss in prenatal diagnosis provided in the embodiments of this application are based on the same inventive concept and have the same beneficial effects as the methods adopted, run or implemented by the application programs stored therein.

[0072] The various embodiments in this application are described in a related manner. Similar or identical parts between embodiments can be referred to mutually. Each embodiment focuses on describing the differences from other embodiments. In particular, the embodiments of a calculation method, electronic device, electronic device, and readable storage medium for evaluating the accurate identification of lost heterozygosity in prenatal diagnosis are basically similar to the embodiments of the calculation method for accurately identifying lost heterozygosity in prenatal diagnosis described above, and are therefore described relatively simply. Relevant parts can be referred to in the description of the embodiments of the calculation method for accurately identifying lost heterozygosity in prenatal diagnosis described above.

Claims

1. A computational method for accurately identifying loss of heterozygosity in prenatal diagnosis, characterized in that, include: Genomic DNA was extracted from prenatal samples and subjected to fragmentation, end repair, A-tailing, adapter ligation, PCR amplification, and magnetic bead purification before sequencing to complete the whole genome library construction. The sequencing data was cleaned by removing impurities using the trim_galore method to obtain a pure FASTQ file. After genome alignment, deduplication, indel re-alignment, and base correction, a bam file was obtained. Sequencing depth and coverage were calculated. Then, GATK was used to identify, filter, and annotate variants, generating sequencing data analysis results. This process included processing the sequencing data using trim_galore to remove end sequencing adapters, amplification primers, and low-quality bases, generating a pure FASTQ file; aligning the pure FASTQ file with a reference genome, removing duplicate reads, performing indel re-alignment, and base quality correction to generate a bam file; calculating the sequencing depth and genome coverage of the bam file to obtain depth and coverage information; using GATK software to identify, filter, and annotate variants in the bam file, generating a dataset containing variant site information; and integrating the above processing results to generate complete sequencing data analysis results. The inserted and missing regions in the sample BAM file were re-aligned to obtain genotypes and generate a VCF dataset. After screening high-quality SNP sites using VCFtools, the proportion of mutated alleles at the SNPs was calculated. Sites with a total read count ≥3 were screened, and allele mutation frequency calculation and screening results were generated. This included using GATK software to re-align the inserted and missing regions in the sample BAM file, obtain sample genotypes, and generate a VCF format dataset. Each variant site has corresponding genotype information. The VCF dataset was filtered using vcftools to retain high-quality SNP loci, generating a filtered SNP locus dataset. The filtered SNP loci were then calculated using R software based on the formula... The allele mutation frequency of each SNP locus was obtained, generating a dataset of SNP locus allele mutation frequency values. This represents the proportion of mutated alleles at this SNP. This is a mutation-type read. This is a wild-type reading passage; The dataset of SNP loci allele mutation frequency values ​​was filtered, and loci with ≥3 wild-type reads + mutant reads were retained to generate a dataset of SNP loci and their allele mutation frequency values ​​that met the criteria. The above processing results were integrated to generate allele mutation frequency calculation and filtering results containing the filtered SNP loci and their allele mutation frequency values. The whole genome was divided into 0.1Mb ​​windows. The number of effective SNP sites with a total read count ≥3 in each window was counted and low-quality windows were excluded. The HR value of the retained windows was calculated, and the window homozygosity calculation results were generated. Normal sample CNV-seq data verified by CMA were selected for processing to obtain the mean HR of each window. Non-linear weight references were constructed and missing value windows were filled with the median. The weight reference list was then integrated to generate non-linear weight reference results. HR data is input into the target CBS algorithm and a weighted reference list is introduced. By calculating the significant changes in HR values ​​between windows and recursively segmenting significant AOH segments, a table of information containing the chromosome position, start and end positions, mean HR value, and number of effective bins in the region of AOH segments is generated to complete the identification of loss of heterozygosity. The method described is a non-disease diagnosis method.

2. The method as described in claim 1, characterized in that, The entire genome was divided into 0.1Mb ​​windows. The number of valid SNP sites with a total read count ≥3 within each window was counted, and low-quality windows were excluded. The HR value of the retained windows was calculated, and the window homozygosity calculation results were generated, including: The whole genome was divided into windows using the R package GenomicRanges, generating windows of 0.1Mb ​​in size with global cross-chromosome numbering, each window having a unique identifier; Statistical processing was performed on the SNP sites in each window to filter out the effective SNP sites with a total number of reads ≥3, and the number of effective SNP sites in each window was generated. Low-quality windows are excluded based on the number of valid SNP sites, and windows with ≥5 valid SNP sites are retained to generate a high-quality window dataset. For high-quality windows, use the formula Calculate the HR value and generate HR value data for each window, where N≥5, n is the number of SNP sites with AF=0 or AF=1 in the window, and N is the total number of SNP sites in the window. Integrate high-quality window data and their HR values ​​to generate window homozygosity calculation results, where the key is the window identifier and the value is the HR value of the corresponding window.

3. The method as described in claim 1, characterized in that, Normal sample CNV-seq data validated by CMA were selected and processed to obtain the mean HR for each window. A non-linear weighted reference was constructed, and missing value windows were filled with the median. The weighted reference list was then integrated to generate the non-linear weighted reference results, including: The CNV-seq data of selected normal samples verified by CMA are processed to generate HR value data for each window, and each window has a unique identifier. Calculate the HR mean for each window, and fill in the missing value windows with the median of the HR mean of normal samples to generate a complete window HR mean dataset. Through public Calculate the weight value for each window to generate window weight data, which provides the basic weight data for the subsequent construction of the weighted likelihood ratio statistic. By integrating window identifiers with corresponding weight values, a weight reference list is generated, which improves the weighted likelihood ratio statistic in the CBS algorithm. The calculation provides key inputs and generates nonlinear weighted reference results, where the key is the window identifier and the value is the weight value of the corresponding window and the basic parameters required for the weighted likelihood ratio calculation. This is a specific value for a particular window k in the weight reference list. This represents the mean HR from the first window to the kth window. This represents the mean HR from the (k+1)th window to the last window, and the difference between the two reflects the overall difference in HR values ​​on both sides of the cutoff point k. The variance of whole-genome HR values ​​is represented by .

4. The method as described in claim 3, characterized in that, The HR data is input into the improved CBS algorithm. After introducing a weighted reference list, the significant changes in HR values ​​between windows are calculated and recursively segmented to divide significant AOH regions. Finally, a table containing information such as chromosome location is output and visualized to complete the identification of loss of heterozygosity, including: The window HR value data of the sample is input into the improved CBS algorithm, and a weight reference list is introduced, in which the HR value of each window is associated with the corresponding weight value. The significance of the change in HR values ​​between windows was calculated using the CBS algorithm, based on the weighted likelihood ratio statistic. Perform recursive segmentation to divide significant AOH continuous segments; Information is extracted from the segmented AOH regions to generate tabular data containing chromosome location, start and end positions, mean HR value, and number of effective bins within the region; Visualize the AOH region tabular data to present the distribution characteristics of AOH regions in the genome; By integrating tabular data and visualization results, the identification of missing heterozygosity is completed, where the key is the AOH segment identifier and the value is the detailed information of the corresponding segment and the associated data of the visualization results.

5. A computational device for accurately identifying loss of heterozygosity in prenatal diagnosis, characterized in that, The apparatus for implementing the method of claim 1 includes: The acquisition module is used to extract genomic DNA from prenatal samples and perform fragmentation, end repair, A-tailing, adapter ligation, PCR amplification, and magnetic bead purification before sequencing to complete the whole genome library construction. The processing module removes impurities from the sampled data using `trim_galore` to obtain a clean FASTQ file. After genome alignment, deduplication, Indel realignment, and base correction, a BAM file is generated. Depth and coverage are calculated, and GATK is used to identify, filter, and annotate variants, generating sequencing data analysis results. The module then re-aligns the sample BAM file, inserts missing regions to obtain genotypes, and generates a VCF dataset. After screening for high-quality SNP sites using VCFtools, the proportion of mutant alleles at the SNPs is calculated, and sites with a total read count ≥3 are selected, generating allele mutation frequency calculation and selection results. The entire genome is divided into 0.1Mb ​​windows, and effective reads with a total read count ≥3 are counted within each window. The number of SNP loci is counted and low-quality windows are excluded. The HR value of the retained windows is calculated, and the homozygosity of the windows is calculated. CNV-seq data of normal samples verified by CMA are selected for processing to obtain the mean HR value of each window. A nonlinear weighted reference is constructed and the median is used to fill the missing value windows. The weighted reference list is integrated to generate the nonlinear weighted reference result. The HR data is input into the target CBS algorithm and introduced into the weighted reference list. By calculating the significant change of HR value between windows and recursively segmenting significant AOH segments, a table of information containing the chromosome position, start and end positions, mean HR value, and number of effective bins in the region of AOH segments is generated to complete the identification of heterozygosity loss. The method described is a non-disease diagnosis method.

6. An electronic device, characterized in that, include: First processor; and memory for storing executable instructions of the first processor; The first processor is configured to execute the computational method for accurately identifying loss of heterozygosity in prenatal diagnosis as described in any one of claims 1 to 4 by executing the executable instructions.

7. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the second processor, it implements the calculation method for accurately identifying loss of heterozygosity in prenatal diagnosis as described in any one of claims 1 to 4.

Images