Method and device for detecting copy number variation types in thalassemia patients
By constructing a baseline database of copy number of thalassemia-related genes and performing sliding window analysis, combined with softclip technology and comparison with a self-built database, the problem of detecting complex copy number variations in thalassemia patients was solved, achieving efficient and accurate detection results.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- NANODIGMBIO (NANJING) BIOTECHNOLOGY CO LTD
- Filing Date
- 2023-11-24
- Publication Date
- 2026-05-22
Smart Images
Figure CN117649873B_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of bioinformatics, and specifically relates to a method and apparatus for detecting copy number variation types in thalassemia patients. Background Technology
[0002] Thalassemia is a common single-gene disorder caused by pathogenic mutations in the α or β globin genes, leading to an imbalance between the α and β globin chains. Traditional detection methods for thalassemia (and other hemoglobinopathies) begin with hematological tests, such as measuring mean corpuscular hemoglobin (Hb) and mean corpuscular volume (MCV), followed by Hb electrophoresis, and finally genetic testing to identify the specific pathogenic mutation. While this process is effective in most cases, it requires significant laboratory work and involves incremental decision-making, making large-scale testing quite cumbersome. In rare cases, some disease mutations may not present with positive hematological characteristics, leading to potential false-negative hematological screening results.
[0003] Novel mutations or complex structural changes in globin genes pose a major challenge to genetic counseling and prenatal diagnosis. Currently, carrier screening and prenatal testing are considered the preferred protective measures to prevent thalassemia and related diseases. Traditional mutation detection methods require the use of multiple technologies to analyze various types of gene mutations. However, NGS technology based on targeted enrichment can significantly improve the efficiency of mutation detection. Compared to multiplex amplicon assays, probe hybridization capture is more efficient and has the following advantages: ① it can detect both common and uncommon mutations; ② it is more cost-effective in large-scale population screening; ③ it can reduce the risk of misdiagnosis or missed diagnosis of patient mutations.
[0004] The α-globin gene cluster is located in the 16p13.3 region of chromosome 16, a region possessing a series of unique characteristics, including homologous sequences, scattered repetitive sequences, and various types of genetic variation. Due to this sequence homology, fragments from NGS sequencing often map to multiple different locations rather than a single location. This not only leads to missed detections of point mutations and insertion / deletion mutations (InDels) but also results in inaccuracies in sequencing depth, which is crucial for detecting structural variants (SVs). Furthermore, complex structural variants are quite common in thalassemia-endemic regions. For example, --SEA / -α3.7 is a complex heterozygous genotype that causes Hb H disease and is widespread in Hong Kong and some other Southeast Asian regions. However, there are few commercially available methods for detecting complex heterozygous genes, especially structural variants. Therefore, accurate detection of these complex genotypes using bioinformatics analysis tools is essential for the comprehensive implementation of molecular diagnosis and carrier screening for thalassemia. Summary of the Invention
[0005] To address the above problems, this invention discloses a method for detecting copy number variation types in thalassemia patients, comprising the following steps:
[0006] S1. Collect variation information related to thalassemia copy number from the HbVar, LOVD, and Itha databases, and build a thalassemia copy number variation database after removing redundancy.
[0007] S2. Capture the nucleic acid sequences of genes related to thalassemia in healthy individuals using hybridization probes and perform high-throughput sequencing to construct a baseline database of copy numbers of thalassemia-related genes in healthy individuals;
[0008] S3. Capture the nucleic acid sequence in the sample to be identified by hybridization probe, perform high-throughput sequencing, compare it with the human reference genome, and use softclip to determine the breakpoints in the sample genome where copy number variations may occur;
[0009] S4. Using the sequencing data of the sample to be identified in step S3, compare it with the baseline database in step S2, and select the appropriate copy number reference baseline for the sample;
[0010] S5. By comparing the sequencing data of the sample to be identified in step S3 with the copy number reference baseline determined in step S4, the suspected breakpoint region of copy number change is obtained;
[0011] S6. Combining the regional breakpoint information in step S3 and the suspected breakpoint region of copy number change in step S5, determine the range of the region where the copy number of the sample to be tested changes and the type of copy number change in that region.
[0012] S7. Using the range of regions where the copy number of the sample to be tested changes and the copy number in that region from step S6, and combining this with the self-built thalassemia copy number variation database from step S1, determine the type of thalassemia copy number variation to which the sample belongs.
[0013] Further, in step S2, the number of healthy individuals used to construct the copy number baseline database is X. Capture sequencing is performed on each of the X samples, and the depth of each point within the capture interval is calculated. Using W bp as the standard sliding window size and S bp as the step length, the average depth within each standard sliding window of each sample's capture interval is calculated, resulting in X reference depth sets. X' samples are randomly selected from the reference depth sets, requiring X' to be greater than 1 / 2X and less than X-1, repeated M times to form M candidate reference depth libraries. The average depth of all samples in each candidate library within each of the aforementioned standard sliding windows is calculated to form the average depth set within the standard sliding window. The overall average depth of the capture interval for all samples is also calculated, resulting in the M average depth sets within the standard sliding windows and their corresponding overall average depths of the capture intervals, thus forming the copy number baseline database.
[0014] Further, in step S3, the sequencing results of the sample to be identified are captured and compared with the reference genome. Reads containing softclips are extracted from the alignment result file. Reads with softclip lengths greater than a certain value are selected, and the location where the softclip occurs is recorded as breakpoint A. Simultaneously, the number of reliable reads at point A is AR = AR + 1. The softclip sequence of this read is extracted, and alignment is performed again on the chromosome aligned to this read. The aligned location is recorded as breakpoint B, and the number of reliable reads at point B is BR = BR + 1. When both the number of reliable reads AR and BR are greater than a certain value, points A and B are considered to be the two ends of the region where copy number variation occurs. Points A and B are paired to form an AB paired point set. When only one read in AR or BR has a read count greater than a certain value, the location of point A or point B is recorded separately to form an AB single point set. The copy number variation type of each paired point within the AB paired point set is preset, with the preset value in the form FR-R'-F', where F and R are integers ranging from 0 to 4, R = R', F = F', and |RF| ≤ 2. The copy number mutation type of each point in the AB single-point set is preset, and the preset value is in the form of FR, where the values of F and R are both integers from 0 to 4 and |RF|≤2.
[0015] Further, in step S4, using the sequencing results of the sample to be identified, the depth of each point within the capture interval is calculated. Using W bp as the standard sliding window size and S bp as the step length, the average depth within each standard sliding window of each sample and the overall average depth of the capture interval are calculated, yielding the average depth of each standard sliding window of the sample to be tested and the average depth of the sample to be tested. Here, W and S are consistent with those used when constructing the copy number baseline database. The difference between the average depth of the sample to be tested and the overall average depth of each capture interval in the copy number baseline database is calculated. The overall average depth of the capture interval with the smallest difference and the corresponding average depth within the standard sliding window are used as the copy number reference baseline for the sample to be tested.
[0016] Further, in step S5, the copy number of each standard sliding window region of the test sample is calculated using the average depth of the test sample and the copy number reference baseline of the test sample. Using the GC content of the standard sliding window within a certain range as a standard, the positions of the standard sliding windows that meet this standard are recorded. The copy number of the standard sliding window regions corresponding to these standard sliding window positions is used to calculate the mean, and the difference is calculated with the copy number at GC correction anchor point 2. The result is the correction offset, which is used to correct the copy number of the standard sliding window regions of the test sample. The corrected copy number replaces the original copy number. Furthermore, the correction offset is calculated for each chromosome separately, and each chromosome is corrected separately.
[0017] Further, in step S5, the standard sliding windows are sorted by position. When the difference in the corrected copy number of adjacent standard sliding windows is within a certain range, the adjacent standard sliding windows are merged into bin windows, and the average copy number within the bin window is recalculated. When the difference in the average copy number of adjacent bin windows exceeds a certain range, the integer part of the average copy number of the first and last standard sliding windows within the bin window is recorded. When the integer part of the average copy number of the first and last standard sliding windows is different, if the absolute value of the difference between the average copy number of the 5 bin windows upstream of the first standard window and the average copy number of the first window is less than 0.5, then the copy number of the first standard window is true. At the same time, if the absolute value of the difference between the average copy number of the 5 bin windows downstream of the last standard window and the average copy number of the last window is less than 0.5, then the copy number of the last standard window is true. When the copy numbers of both the first and last standard windows are true, the adjacent bin windows are merged into usable bin windows. The region C' of the standard windows connected when the bin windows are merged is recorded, with its midpoint being C. Furthermore, when there are multiple C, a set of C points is formed, corresponding to multiple sets of C'.
[0018] Further, in step S5, the copy number variation type within the available bin sliding window is calculated. If the average copy number of the first standard sliding window is close to F, and the average copy number of the last standard sliding window is close to R, then the copy number variation type corresponding to C is FR. Specifically, the values of F and R are integers from 0 to 4, and |RF| ≤ 2. If the average copy number of the standard sliding window is greater than 4, it is considered 4; if it is less than 0, it is considered 0. The average copy number of the standard sliding window being close to F or R means that the absolute value of the difference between this copy number and F or R is less than 0.5. FR is the copy number variation type of region C'.
[0019] Further, in step S6, after pairing all points in the obtained AB single-point set with all points in the C point set, and combining this with the AB paired-point set, a paired endpoint set for the variant region to be tested is obtained. More specifically, the pairing rule is: single points with copy number change type FR are paired with all single points with copy number change type R'-F', where numerically F = F' and R = R'. After pairing, the copy number change type of the paired endpoint is updated to FR-R'-F'.
[0020] Further, in step S6, the standard sliding window number between each paired endpoint in the set of paired endpoints of the variant region to be tested is calculated as W_all, and the sequence length between paired endpoints is calculated as L. When F is greater than R in the copy number change type FR-R'-F' of the paired endpoints of the variant region to be tested, the standard sliding window number with the average copy number between each paired endpoint in the range of 0 to (R+r) is calculated as W_support; when F is less than R, the standard sliding window number with the average copy number between each paired endpoint in the range of (Rr) to 4 is calculated as W_support. Here, r is a decimal in the range of 0 to 1.
[0021] Further, in step S6, when W_support > 20 and L ≥ 100000, W_support / W_all ≥ 95%, and the type of change in the paired endpoint and the paired copy number is true; when W_support > 20 and 20000 ≤ L < 100000, W_support / W_all ≥ 85%, and the type of change in the paired endpoint and the paired copy number is true; when W_support > 20 and 5000 ≤ L < 20000, W_support / W_all ≥ 80%, and the type of change in the paired endpoint and the paired copy number is true. The endpoint and the type of copy number change for the paired endpoint are true. When W_support > 20 and L < 5000, W_support / W_all ≥ 70%, and the type of copy number change for the paired endpoint and the paired endpoint is true. If the copy number change types of two paired endpoints that are true are consistent, and the regions covered by the endpoints overlap, the paired endpoints can be merged. The minimum value of the new paired endpoint is the minimum value of the merged front point, and the maximum value of the new paired endpoint is the maximum value of the merged front point. The type of copy number change is consistent with the type of copy number change before merging. This yields the range of regions where the copy number of the tested sample changes and the type of copy number change in that region.
[0022] Further, in step S7, the region of copy number variation in each test sample must achieve a 90% similarity to the variation region indicated by the entry in the self-built thalassemia copy number variation database. That is, the overlap length between the two regions must be more than 90%, and the type of copy number variation in the test sample must be consistent with the entry in the database. If so, the sample is considered to contain this thalassemia copy number variation type. Regions of copy number variation that do not reach 90% similarity with the database, along with the type of copy number variation in those regions, are considered newly discovered copy number variation regions and new copy number variation types.
[0023] The present invention also discloses an apparatus for detecting copy number variation types in patients with thalassemia, using the method described in any of the above-mentioned embodiments;
[0024] The device includes:
[0025] Thalassemia Copy Number Variance Database Unit: This unit is configured to contain a database of various thalassemia copy number variant types and their descriptions.
[0026] Baseline construction unit for thalassemia-related genes in healthy individuals: This unit is designed to construct a baseline for the sequencing depth of thalassemia-related genes using sequencing data from healthy individuals.
[0027] Data cut-off points and verification units: set to use reads containing softclips in sequencing reads to determine the location and reliability of large genomic variants;
[0028] Baseline selection unit for test sample: It is set to select the most suitable baseline for the sample for subsequent analysis by jointly analyzing the sequencing depth information of the test sample with the baseline database;
[0029] Suspected breakpoint scanning unit: It is set to use the depth information of the sample to be tested and the selected baseline depth information to calculate the copy number of each interval of the sample to be tested, and extract the intervals where there may be breakpoints by using the copy number of each interval.
[0030] Copy number variation interval verification module: set to verify the locations of large-scale variations and intervals where there may be breakpoints, to determine the intervals and types of copy number variations that have occurred;
[0031] The copy number variation annotation module is configured to compare the defined copy number variation ranges and types with the database to complete the annotation of copy number variation ranges and report newly discovered copy number variation types.
[0032] The present invention has the following beneficial effects:
[0033] Traditional methods for detecting copy number variants in thalassemia patients target only individual genes or variant types. This invention provides a method for detecting copy number variants in thalassemia patients that extends beyond individual genes or variant types. It can clearly detect known structural variants associated with thalassemia in samples, and is particularly effective in accurately identifying atypical thalassemia-related genomic structural variants. Validated and applied according to standards, this method can be used to detect known thalassemia copy number variants and discover new ones. Attached Figure Description
[0034] Figure 1 This diagram illustrates the calculation of the standard sliding window depth after using the DepthOfCoverage command in the GATK software to calculate the depth of each point in the capture interval.
[0035] Figure 2 A schematic diagram for constructing a copy number baseline database. Detailed Implementation
[0036] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Unless otherwise specified, the following embodiments and features can be combined with each other.
[0037] The detailed description of the embodiments of the present invention provided below is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.
[0038] Example 1
[0039] Blood samples were extracted from 20 healthy individuals and blood samples from patients who may have thalassemia. After high-throughput sequencing using probes targeting thalassemia-related genes, the resulting sequencing data were used for subsequent analysis.
[0040] The thalassemia structural variation detection for this sample includes the following steps:
[0041] Information on thalassemia copy number variations was collected from the HbVar, LOVD, and Itha databases. After redundancy removal, a thalassemia copy number variation database was constructed. This database contains type information, location information, and related gene information for common Mediterranean copy number variations such as SEA, alpha3.7, and alpha4.2, as shown in Table 1.
[0042] Table 1. Examples of information on copy number variation in the Mediterranean Sea.
[0043]
[0044] Twenty healthy human samples were captured and sequenced. The depth of each point within the capture interval was calculated using the DepthOfCoverage command in the GATK software, with a standard sliding window size of 75 bp and a step length of 10 bp. The average depth within each standard sliding window of each sample's capture interval was calculated, resulting in 20 reference depth sets. Fifteen samples were randomly selected from the reference depth sets, repeated 10 times, to form 10 candidate reference depth libraries. The average depth of all samples in each candidate library within each of the aforementioned standard sliding windows was calculated to form the average depth set within the standard sliding window. The overall average depth of the capture interval for all samples was also calculated, resulting in the 10 average depth sets within the standard sliding windows and their corresponding overall average depths of the capture intervals, forming the copy number baseline database.
[0045] The overall average depths of the capture intervals of the 10 reference depth candidate libraries are 2670, 2684, 2666, 2798, 2803, 2850, 2855, 2903, 2950, and 3076, respectively.
[0046] Furthermore, the sequencing results of the sample to be identified were captured and compared with the human reference genome hg38. Reads containing softclips were extracted from the alignment results file. Reads with softclip lengths greater than 20 bp were selected, and the location of the softclip was recorded as breakpoint A, in the form of chr16:165397. Simultaneously, the confidence read count AR = AR+1 for point A was set. The softclip sequence of this read was extracted, and the alignment was performed again on the chromosome to which this read was aligned. The aligned location was set as breakpoint B, in the form of chr16:184783. Simultaneously, the confidence read count BR for point B was set as BR+1. When both the confidence read counts AR and BR are greater than 100, points A and B are considered to be the two ends of the region where copy number variation occurred. The positions of points A and B are paired to form the AB pairing point set, in the form of: pairing chr:165397-184783. When only one read count in AR or BR is greater than a certain value, the position of point A or point B is recorded separately to form an AB single-point set, in the form of: single-point chr: 177277. The copy number mutation type of each paired point in the AB paired-point set is preset, with the preset value in the form FR-R'-F', where F and R are integers from 0 to 4, R = R', F = F', and |RF| ≤ 2. For example, if the copy number mutation type of the above paired point is 2-1-1-2, it is recorded as (2-1-1-2, chr16, 165339, 184783). The copy number mutation type of each point in the AB single-point set is preset, with the preset value in the form FR, where F and R are integers from 0 to 4 and |RF| ≤ 2. For example, if a single point in this sample has a copy number mutation type of 2-1.
[0047] Furthermore, using the sequencing results of the samples to be identified, the depth of each point within the capture interval was calculated using the DepthOfCoverage command in the GATK software. With a standard sliding window size of 75 bp and a step length of 10 bp, the average depth within each standard sliding window and the overall average depth of the capture interval for each sample were calculated, yielding the average depth of each standard sliding window and the average depth of the sample as 2864. The difference between the average depth of the sample and the overall average depth of each capture interval in the copy number baseline database was calculated. The overall average depth of the capture interval with the smallest difference was determined to be 2855, and the set of average depths within the corresponding standard sliding window was used as the copy number reference baseline for the sample.
[0048] Furthermore, using the average depth of the sample under test and the copy number reference baseline, the copy number of each standard sliding window region of the sample under test was calculated. Using a GC content of 0.3-0.43 for the standard sliding window as the standard, the positions of standard sliding windows conforming to this standard were recorded. The copy numbers of the standard sliding window regions corresponding to these standard sliding window positions were used to calculate the mean of each chromosome, and the difference was calculated with the GC correction anchor point 2 copy number. The correction offsets for each chromosome were: chr2->0.29293889395426; chr6->0.473280205151707; chr11->0.518167043398662; chr16->0.267824595285057; chr19->0.230622071720058. These were used to correct the copy number of the standard sliding window region of each chromosome in the sample under test, and the corrected copy number replaced the original copy number.
[0049] Furthermore, based on the standard sliding window positions, when the corrected copy number difference between adjacent standard sliding windows is between 0 and 0.6, the adjacent standard sliding windows are merged into bin windows, and the average copy number within the bin window is recalculated. When the average copy number difference between adjacent bin windows is greater than 0.6, the integer part of the average copy number of the first and last standard sliding windows within the bin window is recorded. When the integer part of the average copy number of the first and last standard sliding windows is different, if the absolute value of the difference between the average copy number of the 5 bin windows upstream of the first standard window and the average copy number of the first window is less than 0.5, then the copy number of the first standard window is true. At the same time, if the absolute value of the difference between the average copy number of the 5 bin windows downstream of the last standard window and the average copy number of the last window is less than 0.5, then the copy number of the last standard window is true. When the copy numbers of both the first and last standard windows are true, the adjacent bin windows are merged into usable bin windows, and the region C' of the standard windows connected when the bin windows are merged is recorded, with its midpoint being C. Furthermore, when there are multiple C, a set of C points is formed, corresponding to multiple sets of C'.
[0050] Furthermore, the copy number variation type within the bin sliding window can be calculated. If the average copy number of the first standard sliding window is close to F, and the average copy number of the last standard sliding window is close to R, then the copy number variation type corresponding to C is FR. Specifically, the values of F and R are integers from 0 to 4, and |RF| ≤ 2. If the average copy number of the standard sliding window is greater than 4, it is considered 4; if it is less than 0, it is considered 0. The average copy number of the standard sliding window being close to F or R means that the absolute value of the difference between this copy number and F or R is less than 0.5. FR represents the copy number variation type of region C'.
[0051] Furthermore, after pairing all points in the obtained AB single-point set with all points in the C point set, and combining this with the AB paired-point set, we obtain the paired endpoint set of the variant region to be tested. Even further, the pairing rule is: single points with copy number change type FR are paired with all single points with copy number change type R'-F', where numerically F = F' and R = R'. After pairing, the copy number change type of the paired endpoint is updated to FR-R'-F'. For example, pairing 2-1 with 1-2 yields 2-1-1-2.
[0052] Further, the standard sliding window number between each paired endpoint in the set of paired endpoints of the variant region to be tested is calculated as W_all, and the sequence length between paired endpoints is calculated as L. When F is greater than R in the copy number variation type FR-R'-F' of the paired endpoints of the variant region to be tested, the standard sliding window number with the average copy number between each paired endpoint in the range of 0 to (R+r) is calculated as W_support; when F is less than R, the standard sliding window number with the average copy number between each paired endpoint in the range of (Rr) to 4 is calculated as W_support. Here, r is a decimal in the range of 0 to 1, preferably 0.5.
[0053] Furthermore, when W_support > 20 and L ≥ 100000, W_support / W_all ≥ 95%, and the type of change in the paired endpoint and paired copy number is true; when W_support > 20 and 20000 ≤ L < 100000, W_support / W_all ≥ 85%, and the type of change in the paired endpoint and paired copy number is true; when W_support > 20 and 5000 ≤ L < 20000, W_support / W_all ≥ 80%, and the type of change in the paired endpoint and paired copy number is true. The type of copy number change for the paired endpoints is true. When W_support > 20 and L < 5000, W_support / W_all ≥ 70%, and the type of copy number change for the paired endpoints is true. If the copy number change types of two paired endpoints with true copy number change types are consistent, and the regions covered by the endpoints overlap, the paired endpoints can be merged. The minimum value of the new paired endpoint is the minimum value of the merged front point, and the maximum value of the new paired endpoint is the maximum value of the merged front point. The type of copy number change is consistent with the type of copy number change before merging. This yields the range of regions where the copy number of the tested sample changes and the type of copy number change in that region. For example, in (2-1-1-2, chr16, 165339, 184783), the standard sliding window number between paired endpoints is 260, the length L is 19444, and W_sopport = 237. It satisfies the condition that W_support / W_all ≥ 80% when W_support > 20 and 5000 ≤ L < 20000, therefore the pairing endpoint and pairing copy number change type of this entry are true.
[0054] Furthermore, the region of copy number variation in each test sample must achieve a 90% similarity to the variation region indicated by the entry in the self-built thalassemia copy number variation database. This means the overlap length between the two regions must be more than 90%, and the type of copy number variation in the test sample must match the entry in the database. Only then is the sample considered to contain this type of thalassemia copy number variation. Regions of copy number variation that do not reach 90% similarity to the database, along with their corresponding copy number variation types, are considered newly discovered copy number variation regions and new copy number variation types.
[0055] Example 2
[0056] An apparatus for detecting copy number variation types in patients with thalassemia, the apparatus comprising:
[0057] Thalassemia Copy Number Variance Database Unit: This unit is configured to contain a database of various thalassemia copy number variant types and their descriptions.
[0058] Baseline construction unit for thalassemia-related genes in healthy individuals: This unit is designed to construct a baseline for the sequencing depth of thalassemia-related genes using sequencing data from healthy individuals.
[0059] Data cut-off points and verification units: set to use reads containing soffclip in sequencing reads to determine the location and reliability of large genomic variants;
[0060] Baseline selection unit for test sample: It is set to select the most suitable baseline for the sample for subsequent analysis by jointly analyzing the sequencing depth information of the test sample with the baseline database;
[0061] Suspected breakpoint scanning unit: It is set to use the depth information of the sample to be tested and the selected baseline depth information to calculate the copy number of each interval of the sample to be tested, and extract the intervals where there may be breakpoints by using the copy number of each interval.
[0062] Copy number variation interval verification module: set to verify the locations of large-scale variations and intervals where there may be breakpoints, to determine the intervals and types of copy number variations that have occurred;
[0063] The copy number variation annotation module is configured to compare the defined copy number variation ranges and types with the database to complete the annotation of copy number variation ranges and report newly discovered copy number variation types.
[0064] The method of using the above-mentioned device includes the following steps:
[0065] Input the thalassemia copy number variation database into the thalassemia copy number variation database unit so that subsequent analyses can call the thalassemia copy number variation database information in this unit.
[0066] Blood samples were extracted from 20 healthy individuals. After high-throughput sequencing using probes targeting thalassemia-related genes, the sequencing data were input into a baseline construction unit for thalassemia-related genes in healthy individuals. This unit will use the sequencing data from healthy individuals to construct a sequencing depth baseline for thalassemia-related genes.
[0067] Furthermore, this unit uses the DepthOfCoverage command in GATK software to calculate the depth of each point within the capture interval, using a standard sliding window size of 75 bp and a step length of 10 bp. The step length is used to calculate the average depth within each standard sliding window of each sample capture interval (e.g., ...). Figure 1 As shown in the figure, 20 reference depth sets were obtained. Fifteen samples were randomly selected from each reference depth set, repeated 10 times, to form 10 reference depth candidate libraries. The average depth of all samples in each candidate library within each of the above standard sliding windows was calculated to form the average depth set within the standard sliding window. The overall average depth of the capture interval for all samples was also calculated. The overall average depths of the capture intervals for the 10 reference depth candidate libraries were obtained as follows: 2670, 2684, 2666, 2798, 2803, 2850, 2855, 2903, 2950, and 3076. Simultaneously, 10 average depth sets within the standard sliding windows were generated, corresponding to the reference depth candidate libraries, thus forming a copy number baseline database (e.g., ...). Figure 2 (As shown).
[0068] Blood samples were extracted from patients with thalassemia. After high-throughput sequencing using probes targeting thalassemia-related genes, the sequencing data was input into a data cut-off and verification unit. This unit used the captured sequencing data to align with the human reference genome hg38, extracting reads containing softclips from the alignment results. Reads with softclip lengths greater than 20 bp were selected, and their softclip locations were recorded as cut-off point A, in the form of chr16:165397. The confidence read count AR at point A was set to AR+1. The softclip sequence of this read was then extracted and re-aligned on the chromosome to which it was aligned. The aligned location was set to cut-off point B, in the form of chr16:184783. The confidence read count BR at point B was set to BR+1. When both the confidence read counts AR and BR were greater than 100, points A and B were considered the two ends of the region where copy number variation occurred. Points A and B were paired to form an AB pairing set, in the form of: pairing chr:165397-184783. When only one read count in AR or BR is greater than a certain value, the position of point A or point B is recorded separately to form an AB single-point set, in the form of: single-point chr: 177277. The copy number mutation type of each paired point in the AB paired-point set is preset, with the preset value in the form FR-R'-F', where F and R are integers from 0 to 4, R = R', F = F', and |RF| ≤ 2. For example, if a pair of paired points in this sample has a copy number mutation type of 2-1-1-2, it is recorded as (2-1-1-2, chr16, 165339, 184783). The copy number mutation type of each point in the AB single-point set is preset, with the preset value in the form FR, where F and R are integers from 0 to 4 and |RF| ≤ 2. For example, if a single point in this sample has a copy number mutation type of 2-1. The resulting file is shown in Table 2.
[0069] Table 2 Examples of Single-Point Copy Number Variation
[0070] Pairing / Single Point type chromosome Breakpoint A Breakpoint B single point 2-1 chr16 177277 ... ... ... ... single point 1-2 chr16 177277 pair 2-1-1-2 chr16 165397 184783 ... ... ... ... ...
[0071] The capture sequencing data of the sample to be tested is input into the baseline selection unit. This unit uses the DepthOfCoverage command in GATK software to calculate the depth of each point within the capture interval, with a standard sliding window size of 75 bp and a step length of 10 bp. The step length is used to calculate the average depth within each standard sliding window of each sample and the overall average depth of the capture interval, resulting in the average depth of each standard sliding window and the average depth of the sample to be tested, which is 2864. The average depth of the sample to be tested is then subtracted from the overall average depth of each capture interval in the copy number baseline database. The overall average depth of the capture interval with the smallest difference is determined to be 2855, and the average depth set within the corresponding standard sliding window is used as the copy number reference baseline for the sample to be tested.
[0072] The suspected breakpoint scanning unit uses the average depth of the sample under test and the copy number reference baseline to calculate the copy number of each standard sliding window region of the sample under test. Using a GC content of 0.3-0.43 for the standard sliding window as a standard, the positions of standard sliding windows conforming to this standard are recorded. The copy numbers of the standard sliding window regions corresponding to these standard sliding window positions are used to calculate the mean, and the difference is taken with the copy number at the GC correction anchor point 2. The correction offsets for each chromosome are: chr2->0.29293889395426; chr6->0.473280205151707; chr11->0.518167043398662; chr16->0.267824595285057; chr19->0.230622071720058. This is used to correct the copy number of the standard sliding window region of each chromosome in the sample under test, and the corrected copy number replaces the original copy number.
[0073] Furthermore, based on the standard sliding window positions, when the difference in corrected copy number between adjacent standard sliding windows is within 0 to 0.6, the adjacent standard sliding windows are merged into bin windows, and the average copy number within the bin window is recalculated. When the difference in the average copy number between adjacent bin windows is greater than 0.6, the integer part of the average copy number of the first and last standard sliding windows within the bin window is recorded. When the integer part of the average copy number of the first and last standard sliding windows is different, if the absolute value of the difference between the average copy number of the 5 bin windows upstream of the first standard window and the average copy number of the first window is less than 0.5, then the copy number of the first standard window is true. At the same time, if the absolute value of the difference between the average copy number of the 5 bin windows downstream of the last standard window and the average copy number of the last window is less than 0.5, then the copy number of the last standard window is true. When the copy numbers of both the first and last standard windows are true, the adjacent bin windows are merged into usable bin windows, and the region C' of the standard windows connected when the bin windows are merged is recorded, with its midpoint being C. Furthermore, when there are multiple C, a set of C points is formed, corresponding to multiple sets of C'.
[0074] Furthermore, the copy number variation type within the available bin sliding window is calculated. If the average copy number of the first standard sliding window is close to F, and the average copy number of the last standard sliding window is close to R, then the copy number variation type corresponding to C is FR. Specifically, the values of F and R are integers from 0 to 4, and |RF| ≤ 2. If the average copy number of the standard sliding window is greater than 4, it is considered 4; if it is less than 0, it is considered 0. The average copy number of the standard sliding window being close to F or R means that the absolute value of the difference between this copy number and F or R is less than 0.5. FR is the copy number variation type for region C'. The resulting file is shown in Table 3.
[0075] Table 3 Data Example
[0076] type chromosome Point C 2-1 chr16 159725 ... ... 1-2 chr16 184746 ... ...
[0077] The copy number variation interval verification module pairs all points in the AB single-point set obtained by the above module with all points in the C point set, and then combines this with the AB paired point set to obtain the paired endpoint set of the variation region to be tested. Furthermore, the pairing rule is: single points with copy number variation type FR are paired with all single points with copy number variation type R'-F', where numerically F = F' and R = R'. After pairing, the copy number variation type of the paired endpoint is updated to FR-R'-F'. For example, 2-1 is paired with 1-2, resulting in 2-1-1-2. The resulting file is shown in Table 4:
[0078] Table 4 Data Example
[0079] type chromosome Pair 1 Pair 2 2-1-1-2 chr16 177277 159725 ... ... ... 2-1-1-2 chr16 177277 184746 2-1-1-2 chr16 165397 184783 ... ... ... ...
[0080] Further, the standard sliding window number between each paired endpoint in the set of paired endpoints of the variant region to be tested is calculated as W_all, and the sequence length between paired endpoints is calculated as L. When F is greater than R in the copy number variation type FR-R'-F' of the paired endpoints of the variant region to be tested, the standard sliding window number with the average copy number between each paired endpoint in the range of 0 to (R+r) is calculated as W_support; when F is less than R, the standard sliding window number with the average copy number between each paired endpoint in the range of (Rr) to 4 is calculated as W_support. Here, r is a decimal in the range of 0 to 1, preferably 0.5.
[0081] Furthermore, when W_support > 20 and L ≥ 100000, W_support / W_all ≥ 95%, and the type of change in the paired endpoint and paired copy number is true; when W_support > 20 and 20000 ≤ L < 100000, W_support / W_all ≥ 85%, and the type of change in the paired endpoint and paired copy number is true; when W_support > 20 and 5000 ≤ L < 20000, W_support / W_all ≥ 80%, and the type of change in the paired endpoint and paired copy number is true. The type of copy number change for the paired endpoints is true. When W_support > 20 and L < 5000, W_support / W_all ≥ 70%, and the type of copy number change for the paired endpoints is true. If the copy number change types of two paired endpoints with true copy number change types are consistent, and the regions covered by the endpoints overlap, the paired endpoints can be merged. The minimum value of the new paired endpoint is the minimum value of the merged front point, and the maximum value of the new paired endpoint is the maximum value of the merged front point. The type of copy number change is consistent with the type of copy number change before merging. This yields the range of regions where the copy number of the tested sample changes and the type of copy number change in that region. For example, in (2-1-1-2, chr16, 165339, 184783), the standard sliding window number between paired endpoints is 260, the length L is 19444, and W_sopport = 237. It satisfies the condition that W_support / W_all ≥ 80% when W_support > 20 and 5000 ≤ L < 20000. Therefore, the pairing endpoint and pairing copy number change type of this entry are true and will be output to the subsequent module.
[0082] The variation interval annotation module compares the paired endpoints and points where the paired copy number change type is true with the variation range indicated by the entries in the self-built thalassemia copy number variation database. If it finds that the similarity with the --SEA entry in the database reaches 90%, that is, the overlap length of the two regions accounts for more than 90% of the two regions and the copy number change type of the sample to be tested is consistent with this entry in the database, then the sample is considered to contain a --SEA type variation.
[0083] As can be seen from the above embodiments, the present invention provides a method for detecting copy number variants in thalassemia patients. This method can clearly detect known structural variants associated with thalassemia carried by the sample and accurately identify atypical thalassemia-related genomic structural variants. Verification has shown that, when used according to standards, this method can be used to detect known thalassemia copy number variant types and discover new thalassemia copy number variant types.
[0084] Finally, it should be noted that the scope of protection of the present invention is not limited to the above embodiments. Those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein; and all such modifications or substitutions will fall within the scope of protection of the present invention.
Claims
1. A method for detecting copy number variation types in patients with thalassemia, characterized in that, Includes the following steps: S1. Collect variation information related to thalassemia copy number from the HbVar, LOVD, and Itha databases, and build a thalassemia copy number variation database after removing redundancy. S2. Capture the nucleic acid sequences of genes related to thalassemia in healthy individuals using hybridization probes and perform high-throughput sequencing to construct a baseline database of copy numbers of thalassemia-related genes in healthy individuals; S3. Capture the nucleic acid sequence in the sample to be identified by hybridization probe, perform high-throughput sequencing, compare it with the human reference genome, and use softclip to determine the breakpoints in the sample genome where copy number variations may occur; S4. Using the sequencing data of the sample to be identified in step S3, compare it with the baseline database in step S2, and select the appropriate copy number reference baseline for the sample; S5. By comparing the sequencing data of the sample to be identified in step S3 with the copy number reference baseline determined in step S4, the suspected breakpoint region of copy number change is obtained; S6. Combining the regional breakpoint information in step S3 and the suspected breakpoint region of copy number change in step S5, determine the range of the region where the copy number of the sample to be tested changes and the type of copy number change in that region. S7. Using the range of regions where the copy number of the sample to be tested changes and the copy number in that region from step S6, and combining this with the self-built thalassemia copy number variation database from step S1, determine the type of thalassemia copy number variation to which the sample belongs.
2. The method according to claim 1, characterized in that, Step S2 includes the following steps: constructing a copy number baseline database: The number of healthy people used to construct the copy number baseline database is X. X samples are captured and sequenced, and the depth of each point in the capture interval is calculated. W bp is the standard sliding window size and S bp is the step length. The step length is used to calculate the average depth in each standard sliding window within the capture interval of each sample, and X reference depth sets are obtained. X' samples are randomly selected from the reference depth set, where X' is greater than 1 / 2X and less than X-1. This process is repeated M times to form M reference depth candidate libraries. The average depth of all samples in each candidate library within each of the above standard sliding windows is calculated to form the average depth set within the standard sliding window. The overall average depth of the capture interval of all samples is also calculated to obtain the average depth sets within the M standard sliding windows and their corresponding overall average depths within the capture interval, thus forming the copy number baseline database.
3. The method according to claim 1, characterized in that, Step S3 includes the following steps: using softclip to determine breakpoints in regions where copy number variations may occur: using the sample to be identified to capture sequencing results, comparing them with the reference genome, extracting reads containing softclips from the alignment result file, selecting reads with softclip lengths greater than a certain value, recording the location of the softclip as breakpoint A, and simultaneously setting the number of reliable reads at point A AR = AR + 1, extracting the softclip sequence of this read, and performing another alignment on the chromosome aligned to this read, the aligned location is set as breakpoint B, and simultaneously setting the number of reliable reads at point B BR = BR + 1; when both the number of reliable reads AR and BR are greater than a certain value, points A and B are considered to be the two ends of the region where copy number variations have occurred, and the positions of points A and B are paired to form the AB paired point set; when only one read in AR or BR has a read count greater than a certain value, the position of point A or the position of point B is recorded separately to form the AB single point set; The copy number mutation type of each paired point in the AB paired point set is preset, with the preset value in the form FR-R'-F', where F and R are both integers in the range of 0 to 4, R=R', F=F', and |RF|≤2; the copy number mutation type of each point in the AB single point set is preset, with the preset value in the form FR, where F and R are both integers in the range of 0 to 4 and |RF|≤2.
4. The method according to claim 1, characterized in that, Step S4 includes the following steps: Selecting a copy number reference baseline applicable to the sample to be identified; using the sequencing results captured from the sample to be identified, calculating the depth of each point within the capture interval; using W bp as the standard sliding window size and S bp as the step length; calculating the average depth within each standard sliding window of each sample and the overall average depth of the capture interval; obtaining the average depth of each standard sliding window of the sample to be tested and the average depth of the sample to be tested, wherein W and S are consistent with W and S in claim 2; The average depth of the sample to be tested is subtracted from the overall average depth of each capture interval in the copy number baseline database. The overall average depth of the capture interval with the smallest difference and the set of average depths within the corresponding standard sliding window are used as the copy number reference baseline for the sample to be tested.
5. The method according to claim 1, characterized in that, Step S5 includes the following steps: determining the suspected breakpoint region of copy number change, and calculating the copy number of each standard sliding window region of the sample under test using the average depth of the sample under test and the copy number reference baseline of the sample under test; Using the GC content of a standard sliding window within a certain range as a standard, the positions of standard sliding windows conforming to this standard are recorded. The copy number of the standard sliding window region corresponding to these standard sliding window positions is used to calculate the mean, and the difference is calculated with the copy number at the GC correction anchor point 2. The result is the correction offset, which is used to correct the copy number of the standard sliding window region of the test sample. The corrected copy number replaces the original copy number. The standard sliding windows are sorted by position. When the difference in the corrected copy number of adjacent standard sliding windows is within a certain range, adjacent standard sliding windows are merged into bin sliding windows, and the average copy number within the bin sliding window is recalculated. When the difference in the average copy number of adjacent bin sliding windows exceeds a certain range, the integer part of the average copy number of the first and last standard sliding windows within the bin sliding window is recorded. When the integer part of the average copy number of the first and last standard sliding windows is different, if the first standard sliding window has 5 bins upstream... If the absolute value of the difference between the average copy number of the sliding window and the average copy number of the first sliding window is less than 0.5, then the copy number of the first standard sliding window is true. Simultaneously, if the absolute value of the difference between the average copy number of the five downstream bins of the last standard sliding window and the average copy number of the last sliding window is less than 0.5, then the copy number of the last standard sliding window is true. When both the copy numbers of the first and last standard sliding windows are true, the adjacent bins are merged to form a usable bin. The region C' of the standard sliding window connected when the preceding and following bins are merged is recorded, with its midpoint as C. The copy number change type within the usable bin is calculated. If the average copy number of the first standard sliding window is close to F, and the average copy number of the last standard sliding window is close to R, then the copy number change type corresponding to C is FR. The average copy number of the standard sliding window being close to F or R means that the absolute value of the difference between this copy number and F or R is less than 0.5; FR is the copy number change type of region C'.
6. The method according to claim 1, characterized in that, Step S7 includes the following steps: determining the thalassemia copy number variant type of the sample by combining the self-built thalassemia copy number variant database: For each sample to be tested, the range of regions where the copy number changes must have a 90% similarity to the range of variation indicated by the entry in the self-built thalassemia copy number variation database. That is, the overlap length between the two regions must be more than 90% and the copy number variation type of the sample to be tested must be consistent with the entry in the database. In this case, the sample is considered to contain this thalassemia copy number variation type. The remaining regions of copy number variations that do not reach 90% similarity with the database and the copy number variation type of those regions are considered newly discovered copy number variations and new copy number variation types.
7. A device for detecting copy number variation types in patients with thalassemia, characterized in that, Apply the method as described in any one of claims 1-6; The device includes: Thalassemia Copy Number Variance Database Unit: This unit is configured to contain a database of various thalassemia copy number variant types and their descriptions. Baseline construction unit for thalassemia-related genes in healthy individuals: This unit is designed to construct a baseline for the sequencing depth of thalassemia-related genes using sequencing data from healthy individuals. Data cut-off points and verification units: set to use reads containing softclips in sequencing reads to determine the location and reliability of large genomic variants; Baseline selection unit for test sample: It is set to select the most suitable baseline for the sample for subsequent analysis by jointly analyzing the sequencing depth information of the test sample with the baseline database; Suspected breakpoint scanning unit: It is set to use the depth information of the sample to be tested and the selected baseline depth information to calculate the copy number of each interval of the sample to be tested, and extract the intervals where there may be breakpoints by using the copy number of each interval. Copy number variation interval verification module: set to verify the locations of large-scale variations and intervals where there may be breakpoints, to determine the intervals and types of copy number variations that have occurred; The copy number variation annotation module is configured to compare the defined copy number variation ranges and types with the database to complete the annotation of copy number variation ranges and report newly discovered copy number variation types.