A full-site bioinformatics analysis method for detecting ultra-low-frequency mutations in blood ctDNA

Abstract

The present application discloses a bioinformatics analysis method for detecting ultra-low frequency mutations in ctDNA in the blood of tumor patients, which covers data splitting, quality control of second-generation sequencing data, comparison with the human reference genome, UMI deduplication correction of the comparison results, left-alignment of the comparison results and correction of the base quality value, statistical base composition of all sites in the detection area, determination of background value based on negative samples, clinical sample background polishing and first round of filtering of variation results, annotation and second round of filtering of variation results. Unlike the common variation detection software on the market, the method provided in this application can accurately detect ultra-low frequency ctDNA mutations in the blood of tumor patients by making detailed statistics on the base composition of all sites in the tumor panel detection area, and formulating background polishing algorithms and variation screening methods. The lowest detection frequency can reach 5 / 10,000, which has important clinical application prospects in early cancer screening, disease monitoring, treatment effect evaluation and prognosis monitoring.

Description

Technical Field

[0001] The present invention belongs to the field of bioinformatics, and in particular relates to a full-site bioinformatics method for detecting ultra-low frequency mutations in ctDNA in the blood of tumor patients. Background Art

[0002] ctDNA (circulating tumor DNA) refers to tumor DNA that exists in body fluids such as blood and cerebrospinal fluid and is free outside the cells. It is about 170bp in length. When tumor cells grow, divide or die, the DNA fragments they release can enter the blood circulation and form circulating tumor DNA. As a potential biomarker, circulating tumor DNA (ctDNA) has attracted much attention in the field of tumor diagnosis and monitoring in recent years. With the rise of the concept of personalized treatment, liquid biopsy technology represented by ctDNA has gradually become a research hotspot as a non-invasive tumor diagnosis and monitoring method. ctDNA carries the genetic information of tumor cells, including tumor-specific mutations, gene rearrangements and other variation information. By detecting ctDNA in blood samples, non-invasive and real-time tumor genetic information can be obtained, providing a new means for early diagnosis of tumors, monitoring of treatment effects, and prognosis evaluation.

[0003] Traditional methods for detecting gene mutations primarily rely on PCR and its derivatives, such as digital PCR and digital titration PCR (ddPCR). However, these methods are typically used to detect known mutation sites and are inefficient for multiple sites or large-scale screening. In contrast, capture-based next-generation sequencing (NGS) enriches circulating tumor DNA (ctDNA) and utilizes high-throughput sequencing technologies (such as Illumina sequencing) to sequence entire exomes or specific regions, enabling comprehensive detection of mutations in multiple genes. This approach, with its high throughput and high sensitivity, plays an increasingly important role in personalized cancer treatment. However, analyzing ctDNA high-throughput sequencing data presents challenges. Tumor mutations are typically low in frequency, requiring high-throughput data. However, as data volume increases, background noise, such as sequencing errors and PCR-induced false positives, inevitably increases. In this context, it is crucial to ensure that the influence of background noise is mitigated to accurately identify low-frequency tumor mutation signals. Summary of the Invention

[0004] In order to solve the above problems, the present invention provides a full-site bioinformatics analysis method for detecting ultra-low-frequency mutations in blood ctDNA. By detailed statistics of the base composition of all sites in the tumor panel sequencing area and formulating a background polishing algorithm and variation screening method, it can accurately detect ultra-low-frequency ctDNA mutations in the blood of tumor patients, with the lowest detection frequency reaching 0.5%, which has important clinical application prospects in early cancer screening, disease monitoring, treatment effect evaluation and prognosis monitoring.

[0005] The present application provides a full-site bioinformatics analysis method for detecting ultra-low-frequency mutations in blood ctDNA, comprising the following steps:

[0006] 1. After the sequencing data is downloaded, the data is split according to the label sequence (barcode) of each sample to obtain the second-generation sequencing data of each sample in Fastq format.

[0007] 2. Quality control of second-generation sequencing data, including removal of adapter sequences, removal of sequences with base quality values ​​below 30, removal of read sequences with 1 or more N bases, removal of one base at the end of the read, removal of low-complexity sequences, and removal of read sequences with a length of less than 50bp after control.

[0008] 3. Compare the quality-controlled data to the human reference genome, store the comparison results in BAM format files, and sort and index the BAM files.

[0009] 4. Use UMI deduplication to correct the comparison results, including the following steps:

[0010] (1) Reads aligned to the reference genome are grouped according to chromosome, start, and end positions, and reads with the same alignment position are grouped together.

[0011] (2) For each of the above groups, cluster them according to the UMI in the group. Each UMI cluster includes at least 2 read sequences. If the data does not have UMI, skip this step.

[0012] (3) If the reads are paired, find the overlapping regions of each pair of reads and score the bases in the overlapping regions based on the consistency and quality of the bases.

[0013] (4) For each base position in each UMI cluster, calculate the score of each base A / T / C / G / N.

[0014] (5) Correct the site based on the frequency of occurrence. If there is a dominant base at a certain position, use this base as the corrected base for that position. If there are one or more reads that are consistent with a high-quality reference genome, or if all reads at that position are of low quality, use the corresponding base in the reference genome as the corrected base for that position.

[0015] 5. Before variant detection, it is necessary to left-align the insertions and deletions (indels) in the resulting BAM file and correct the base quality in the BAM file to eliminate base quality deviations introduced during the sequencing process.

[0016] 6. Mutation detection is performed by statistically analyzing the base composition of all sites within the detection region. The principles are as follows:

[0017] (1) Count the number of A / T / C / G / N and indels with a base quality value greater than 30 at each site in the tumor panel detection area;

[0018] (2) For sites without indels, the base with the largest number, excluding the bases identical to the reference genome, is the variant base at that position. The frequency of the variant base is equal to the number of variant bases divided by the total depth of the site.

[0019] (3) For sites with indels, count the number of reads with indels at the site, and divide this number by the total depth of the site to obtain the frequency of the indel;

[0020] (4) If the remaining three bases and the number of insertions and deletions except the bases identical to the reference genome are all 0, it indicates that no base mutation has occurred at this site;

[0021] (5) If a site has a point mutation and the number of insertions and deletions is not unique, the site will not be counted.

[0022] 7. The present invention statistically analyzed the base composition of all sites in 100 negative samples and further determined the background value of sequencing data, which specifically includes the following steps:

[0023] (1) Count the base composition of all sites in the panel detection area of ​​100 negative samples, and calculate the mutation frequency of the variant base if there is a mutation;

[0024] (2) Count the number of reads at the variant sites of 100 negative samples and calculate the 99% confidence interval of the number of reads at each variant site. The upper limit of the interval is taken as the read background value.

[0025] (3) Count the mutation frequencies of the variant sites in 100 negative samples and calculate the 99% confidence interval of the mutation frequency of each variant site. The upper limit of the interval is taken as the background value of the mutation frequency;

[0026] (4) Variants with mutation frequencies above 30% are not included in the calculation of confidence intervals in (2) and (3).

[0027] 8. Background polishing of clinical samples refers to the actual mutation read count of a site minus the background read count of that site, i.e., the mutation read count of a clinical sample variant site after background polishing. Similarly, the mutation frequency of a clinical sample variant site after background polishing is equal to the actual mutation frequency of the site minus the background mutation frequency of that site.

[0028] 9. The first round of filtering of variant sites, the filtering principles are as follows:

[0029] (1) Remove mutations with chain preference, which means that the mutation site accounts for more than 90% of a single chain (positive or negative chain);

[0030] (2) The mutation sites were distinguished based on whether they fell into the COSMIC database. As shown in Table 1, different thresholds were set based on the total depth of the site, the number of reads, and the mutation frequency to screen the mutation sites after background polishing.

[0031] Table 1. Variant site filtering conditions

[0032] Total depth Number of reads Mutation frequency non-cosmic >＝500 >＝5 >＝0.10％ cosmic >＝500 >＝2 >＝0.05％

[0033] (3) Read number and mutation frequency backtracking: after filtering, the actual read number and mutation frequency before background polishing are output.

[0034] 10. Variant site annotation and second round of filtering follow the following principles:

[0035] (1) SnpEff, SnpSift, and ANNOVAR were used to annotate variant sites;

[0036] (2) Remove mutations that fall in UTR, intron, and are annotated as synonymous_variant;

[0037] (3) Remove mutations with a population frequency of more than 5%. If a population frequency of more than 5% is found in one of the three databases, 1000Genomes Project, ExAC, and genomAD, then the mutations will be removed.

[0038] (4) Remove mutation sites annotated as benign by ClinVar or InterVar and with a mutation frequency of 45% < 55% or a mutation frequency of > 95%;

[0039] (5) Use the hotspot database built by this application to annotate the mutation sites after filtering.

[0040] Compared with the prior art, this application has the following beneficial effects:

[0041] 1. The UMI correction and background polishing algorithms provided by this invention remove the impact of background noise, such as sequencing errors and false positives introduced by PCR, on variant detection, thereby reducing false positives.

[0042] 2. The mutation detection method provided by the present invention is different from the common mutation detection software on the market. This method can accurately detect ultra-low-frequency ctDNA mutations in the blood of tumor patients by comprehensively counting the base composition of all sites in the tumor panel sequencing area and screening out reliable mutation results. The minimum detection frequency can reach 0.5%, thereby improving the sensitivity of detection.

[0043] 3. This invention provides detailed annotation of variant sites and filtering principles, which helps to gain a deeper understanding of the pathogenesis of tumors, guide individualized treatment, evaluate prognosis, discover new drug targets, and determine drug sensitivity, which is of great significance for improving the treatment effect and quality of life of cancer patients. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] Figure 1 The present invention provides a flow chart of a bioinformatics analysis method for detecting ultra-low frequency mutations in blood ctDNA.

[0045] Figure 2 Flowchart for determining background read and mutation frequencies.

[0046] Figure 3 Flowchart of clinical sample background polishing and first round of filtration.

[0047] Figure 4 Variant site annotation and second-round filtering flowchart. DETAILED DESCRIPTION

[0048] In order to more clearly illustrate the contents disclosed by the present invention, including the purpose, technical solutions and advantages, we will use the drawings and detailed description to illustrate. The purpose of providing these embodiments is to make the understanding of the present invention more thorough and comprehensive.

[0049] Unless otherwise defined, all technical and scientific terms used herein are those commonly understood by those skilled in the art in the art of the present invention. The terms used in the specification of the present invention herein are only for the purpose of describing specific embodiments and are not intended to limit the present invention.

[0050] Example 1

[0051] like Figure 1 As shown, a specific embodiment of the present invention provides a bioinformatics analysis method for detecting ultra-low frequency mutations in blood ctDNA, comprising the following steps:

[0052] 1. Data splitting. In this paper, bcl2fastq v1.2.6 software was used to split the next-generation sequencing data to obtain the sequencing data in fastq format for each sample. The --barcode-mismatches parameter was set to 1, which means that the barcode sequence allowed a one-base mismatch.

[0053] 2. Second-generation sequencing data quality control. In this invention, fastp v0.23.2 is used to perform quality control on the data after splitting, including removing adapter sequences, removing sequences with base quality values ​​less than 30, removing read sequences with 1 or more N bases, removing 1 base at the end of the read, removing low-complexity sequences, and removing read sequences with a length less than 50bp after control. For the entire sample, a qualified sample must meet the following three conditions: the Q30 base data volume is greater than or equal to 90%, the adapter contamination ratio is less than 1%, and the effective data volume after quality control must be greater than 80%. In addition, in this step, the UMI length is specified by the --umi_len parameter, and the UMI position is specified by the --umi_loc parameter to facilitate subsequent analysis.

[0054] 3. In the present invention, bwa v0.7.17-r1188 was used to align the quality-controlled fastq sequences with the reference genome sequences, and samtools v1.6 was used to sort and index the alignment results in BAM format.

[0055] 4. Use UMI deduplication to correct the alignment results. In the present invention, UMI deduplication correction is performed on the alignment results using gencore v0.17.2 to more accurately identify and correct deviations in the PCR amplification and sequencing processes, thereby improving the quality and reliability of cfDNA sequencing data. The parameters are set to -s 2 --umi_prefix=UMI --ratio_threshold=0.9.

[0056] 5. In the present invention, the BAM format alignment results were preprocessed before variant detection using GATK 4.0.5.1 software. First, the insertions and deletions (indels) in the BAM file were left-aligned using the LeftAlignIndels command. Then, the BaseRecalibrator command was used to correct the base quality deviation in the BAM file based on the known variant site information. By analyzing the known variant sites and the base quality scores in the sequencing data, BaseRecalibrator can calculate the correction value of each base quality. Finally, the base quality correction value generated by BaseRecalibrator was applied to the original BAM file using the ApplyBQSR command to correct the base quality in the BAM file and eliminate the base quality deviation introduced during the sequencing process.

[0057] 6. The present invention uses bam-readcount v1.0.1 in combination with a self-written Python script to count the number of A / T / C / G / N and indels with a base quality value greater than 30 at each site in the tumor panel detection area of ​​the BAM file; the bam-readcount parameters are set to -b 30 -w 0; for sites without indels, the base with the largest number, excluding the bases identical to the reference genome, is the variant base at that position; a variant base of 0 indicates that the site has not mutated, and the mutation frequency of the variant base is equal to the number of variant bases divided by the total depth of the site; for sites with indels, the number of reads with indels at the site is counted, and this number divided by the total depth of the site is the frequency of the indel; if the number of the other three bases, excluding the bases identical to the reference genome, and the number of indels are all 0, it indicates that the site has not undergone base mutation; if the number of point mutations and indels at a site is not unique, the site is not counted.

[0058] 7. Such as Figure 2 As shown, the present invention counted the base composition of all sites of 100 negative samples through a self-written Python script, and further determined the background value. First, the base composition of all sites in the panel detection area of ​​100 negative samples was counted. If there was a mutation, the mutation frequency of the variant base was calculated, and mutations with a frequency of more than 30% were removed; then the number of reads of the variant sites of the 100 negative samples was counted, and the 99% confidence interval of the number of reads of each variant site was calculated, and the upper limit of the interval was taken as the read background value; finally, the mutation frequency of the variant sites of the 100 negative samples was counted, and the 99% confidence interval of the mutation frequency of each variant site was calculated, and the upper limit of the interval was taken as the mutation frequency background value.

[0059] 8. Such as Figure 3As shown, in the present invention, by writing a self-written Python script, the actual mutation read number of all sites is subtracted from the read background value of the corresponding site, and the resulting value is the mutation read number of the clinical sample variant site after background polishing; similarly, the mutation frequency of the clinical sample variant site after background polishing is equal to the actual mutation frequency of the site minus the mutation frequency background value of the site.

[0060] 9. In the present invention, two rounds of filtering of variants were performed by self-written Python scripts, such as Figure 3 As shown, in the first round of filtering, mutations with a single-stranded fraction exceeding 90% are removed. Then, the mutation sites are classified based on whether they are included in the COSMIC database. For mutations included in the COSMIC database, the thresholds are set to total depth >= 500, mutation read count >= 2, and mutation frequency >= 0.05%. For mutations outside the COSMIC database, the thresholds are set to total depth >= 500, mutation read count >= 5, and mutation frequency >= 0.1%. Finally, the mutation read count and frequency are back-tested to output the actual read count and mutation frequency before background polishing.

[0061] 10. The present invention first uses SnpEff, SnpSift and ANNOVAR combined with self-written Python scripts to annotate the variant sites; then Figure 4 As shown, a second round of filtering is performed on the sites after annotation. The filtering principles are as follows:

[0062] (1) Remove mutations that fall in UTR, intron, and are annotated as synonymous_variant;

[0063] (2) Remove mutations with a population frequency of more than 5%. If a population frequency of more than 5% is found in one of the three databases, 1000Genomes Project, ExAC, and genomAD, then the mutations will be removed.

[0064] (3) Remove mutation sites annotated as benign by ClinVar or InterVar and with a mutation frequency of 45% < 55% or a mutation frequency of > 95%;

[0065] (4) Finally, the hotspot database built by the present invention is used to annotate the mutation sites after filtering.

[0066] Example 2

[0067] To evaluate the effectiveness of the bioinformatics analysis method for detecting ultra-low-frequency mutations in blood ctDNA provided in Example 1, the present invention targeted six genes, namely EGFR, KRAS, BRAF, PIK3CA, TP53, and BRCA1, and mixed DNA fragments containing mutant and wild-type sequences in a certain proportion to form standards with different mutation frequencies. To test the detection sensitivity, the standards were designed with five frequency gradients of 0.03%, 0.05%, 0.1%, 0.15%, and 0.2%, with a total of 10 replicates. Table 2 below lists the 13 known mutation sites used in the standards.

[0068] Table 2. 6 genes and 13 mutation sites in the standard

[0069]

[0070]

[0071] In the present invention, the variation detection method provided by the present invention was tested and compared in parallel with the commonly used variation detection software GATK, Vardict, and VarScan2 on the market. The four methods performed variation detection on the same alignment results. Tables 3, 4, 5, and 6 show the detection of standards at different frequency gradients. It can be clearly seen that the variation detection method provided by the present invention has significantly higher sensitivity than the other three software, especially in the detection of indels, with a detection rate significantly higher than the other three software.

[0072] Table 3. Detection of 0.03% standard products

[0073]

[0074] Table 4. Detection of 0.05% standard products

[0075]

[0076]

[0077] Table 5. Detection of 0.10% standard products

[0078]

[0079] Table 5. Detection of 0.15% standard products

[0080]

[0081]

[0082] Table 6. Detection of 0.2% standard products

[0083]

Claims

1. A bioinformatics analysis method for detecting ultra-low frequency mutations in ctDNA in the blood of tumor patients, characterized in that: The following steps are involved: S-1. Data splitting to obtain Fastq format second-generation sequencing data; S-2. Quality control of next-generation sequencing data; S-3. Align the genome to the human reference genome and store the results in a BAM format file. S-4. Use UMI (Unique Molecular Identifier) ​​to de-duplicate and correct the alignment results; The UMI deduplication and correction comparison results include the following steps: (1) Group the reads aligned to the reference genome according to chromosome, start, and end positions, and group reads with the same alignment position into one group; (2) For each group in (1), cluster according to the UMI in the group. Each UMI cluster includes at least 2 read sequences. If the data does not have UMI, skip this step. (3) If the reads are paired, find the overlapping regions of each pair of reads and score the bases in the overlapping regions based on the consistency and quality of the bases; (4) For each base position in each UMI cluster, calculate the score of each base A / T / C / G / N; (5) Correct the site based on the frequency of occurrence. If there is a dominant base at a certain position, use this base as the base after correction at that position. If there are one or more reads that are consistent with the high-quality reference genome, or if all reads at that position are of low quality, use the corresponding base in the reference genome as the base after correction at that position. S-5. Left-align the alignment results and correct the base quality values; S-6. Detect mutations by statistically analyzing the base composition of all sites within the detection region; S-7. Determine the background value based on the negative sample; The method of determining the background value based on the negative sample comprises the following steps: (1) Count the base composition of all sites in the panel detection area of ​​100 negative samples, and calculate the mutation frequency of the variant base if there is a mutation; (2) Count the number of reads at the variant sites of 100 negative samples and calculate the 99% confidence interval of the number of reads at each variant site. The upper limit of the interval is taken as the read background value. (3) Count the mutation frequencies of the variant sites in 100 negative samples and calculate the 99% confidence interval of the mutation frequency of each variant site. The upper limit of the interval is taken as the background value of the mutation frequency; (4) Variants with mutation frequencies above 30% are not included in the calculation of confidence intervals in (2) and (3); S-8. Clinical sample background polishing; The clinical sample background polishing refers to the actual mutation read number of the site minus the background read value of the site, that is, the mutation read number of the clinical sample variant site after background polishing; similarly, the mutation frequency of the clinical sample variant site after background polishing is equal to the actual mutation frequency of the site minus the background mutation frequency of the site; S-9. First round of filtering of mutation results; S-10. Mutation result annotation and second round of filtering.

2. The bioinformatics analysis method for detecting ultra-low frequency mutations in ctDNA in the blood of tumor patients according to claim 1, characterized in that: The principle of performing variation detection by statistically analyzing the base composition of all sites in the detection region in step S-6 is as follows: (1) Count the number of A / T / C / G / N and indels with a base quality value greater than 30 at each site in the tumor panel detection area; (2) For sites without indels, the base with the largest number, excluding the bases identical to the reference genome, is the variant base of the site. The mutation frequency of the variant base is equal to the number of variant bases divided by the total depth of the site; (3) For sites with indels, count the number of reads with indels at the site, and divide this number by the total depth of the site to obtain the frequency of the indel; (4) If the remaining three bases and the number of insertions and deletions except the bases identical to the reference genome are all 0, it indicates that no base mutation has occurred at this site; (5) If a site has a point mutation and the number of insertions and deletions is not unique, the site will not be counted.

3. The bioinformatics analysis method for detecting ultra-low frequency mutations in ctDNA in the blood of tumor patients according to claim 1, characterized in that: In the first round of filtering of variant sites described in step S-9, chain-biased mutations are first removed. Chain bias refers to mutation sites that account for more than 90% of a single chain. Then, the mutation sites are distinguished based on whether they fall into the COSMIC database. As shown in Table 1, different thresholds are set based on the total depth of the site, the number of reads, and the frequency to screen the mutation sites after background polishing. Finally, the number of reads and the frequency of mutations are back-tested to output the actual number of reads and the frequency of mutations before background polishing. Table 1. Variant site filtering conditions 。 4. A bioinformatics analysis method for detecting ultra-low frequency mutations in ctDNA in the blood of tumor patients according to claim 1, characterized in that: The variant site annotation and second round of filtering described in step S-10 are specifically performed as follows: (1) SnpEff, SnpSift, and ANNOVAR were used to annotate variant sites; (2) Remove mutations that fall in UTR, intron, and are annotated as synonymous_variant; (3) Remove mutations with a population frequency of more than 5%. If a population frequency of more than 5% is found in one of the three databases, 1000Genomes Project, ExAC, and genomAD, then the mutations will be removed. (4) Remove mutation sites annotated as benign by ClinVar or InterVar and with a mutation frequency of 45% < 55% or a mutation frequency of > 95%; (5) Use the hotspot database built by this application to annotate the mutation sites after filtering.

Images