CN108319817B - Method and device for processing circulating tumor DNA repetitive sequence - Google Patents

Abstract

The invention discloses a method and device for processing the repeating sequence of circulating tumor DNA. The method includes: acquiring sequencing data of circulating tumor DNA to be detected and a reference genome sequence, wherein the sequencing data is data obtained by high-throughput sequencing of circulating tumor DNA to be detected; comparing the sequencing data with the reference genome sequence, Obtaining a first alignment result, wherein the first alignment result at least includes: the genomic positions, base sequences and corresponding base quality value sequences of multiple pairs of double-ended sequences; based on the first alignment result, at least one sequence set is obtained ; Perform network clustering on at least a pair of double-ended sequences in each sequence set to obtain at least one first network; Obtain the largest number of double-ended sequences in each first network to obtain independent sequences corresponding to each first network . The invention solves the technical problem of low accuracy in the prior art sequencing data processing method to delete or mark repetitive sequences for sample sequencing.

Description

Method and device for processing circulating tumor DNA repeats

technical field

The present invention relates to the field of genetic engineering, in particular, to a method and device for processing repeated sequences of circulating tumor DNA.

Background technique

During the process of division and proliferation, tumor cells will undergo apoptosis, death, and necrosis, and will also actively release DNA fragments carrying tumor mutations into body fluids, that is, circulating tumor DNA (ctDNA), which is mostly found in the body fluids. In body fluids such as blood, synovial fluid and cerebrospinal fluid, especially in plasma cell-free DNA (cfDNA). Through the sequencing of ctDNA, the detection of the genomic region of the base sequence change (mutation) on the DNA molecule of tumor cells can effectively reflect the patient's response to treatment; after the detection of drug response, the tumor may become resistant to drug treatment, ctDNA detection can also track the generation of drug-resistant mutations, qualitatively and quantitatively; detect whether there is residual tissue after surgery, judge the prognosis and early tumor screening. Different from conventional genomic DNA, ctDNA fragments are short, usually only 100-400 bp, and the content in blood is small, so the amount of ctDNA that can be extracted in practice is very low.

Due to the short and low content of ctDNA fragments, when the amount of extraction is small, multiple rounds of polymerase chain reaction (PCR) need to be carried out in the library construction stage to expand the content of the original extracted DNA. High-throughput sequencing (NGS sequencing for short) and subsequent bioinformatics analysis are performed to generate enough DNA molecules. Since PCR amplification leads to multiple mirror copies of a molecule, resulting in duplicated reads, these invalid duplicate data are very easy to introduce artificial errors for the detection of variants. For the most ideal NGS data analysis process, it is necessary to remove all the sequencing data obtained by PCR as much as possible and restore it to a state without PCR.

In the prior art, two methods for deduplication of repetitive sequences are provided, samtools rmdup and Picard's MarkDuplicates. Among them, the working principle of samtools rmdup is: the sequence (read) obtained by NGS sequencing is compared with the human reference genome (mapping) to obtain the alignment position of this read. If different reads are aligned to the same genomic position, then This part of the reads is considered to be multiple repetitive sequences generated by PCR, only the reads with the highest mapping quality are retained, and the rest of the repetitive sequences are deleted. For PE reads, if the read ratio at both ends is too large on different chromosomes in the genome or the distance between the two is too long (that is, not Proper Paired), it will not be considered. The basic idea of Picard'sMarkDuplicates is the same as that of samtools rmdup. By comparing the mapping positions at the 5' end of the reads, for sequences with the same 5' position, the reads with the highest sequencing quality are selected as the only reads retained after deduplication, and PE reads are not The case of Proper Paired will also be deduplicated.

However, at the same position in the genome, there may often be multiple original molecules. These original molecules are not PCR repeats in the usual sense, and there may be meaningful mutations (for example, tumor-related mutations in ctDNA), but in the above-mentioned In the deduplication method, samtools rmdup and Picard's MarkDuplicates will mistakenly think that they are the same original molecule, and only one pair of reads will be retained, resulting in excessive deduplication and wasting part of the meaningful data. And because the PCR process may be accompanied by random errors, these errors are likely to result in false positive mutation calls, and samtools rmdup and Picard's MarkDuplicates have no way to correct them. Both Samtools rmdup and Picard's MarkDuplicates only consider a certain quality value of the read, but do not consider the specific sequence differences on the read, resulting in the occurrence of unique reads that have been duplicated or incorrectly selected in the past.

Aiming at the problem of low accuracy of performing repetitive sequence deletion or labeling on sample sequencing in the prior art sequencing data processing methods, no effective solution has been proposed so far.

SUMMARY OF THE INVENTION

Embodiments of the present invention provide a method and device for processing repeating sequences of circulating tumor DNA, so as to at least solve the technical problem of low accuracy of repeating sequence deletion or labeling for sample sequencing in prior art sequencing data processing methods.

According to an aspect of the embodiments of the present invention, there is provided a method for processing circulating tumor DNA repeat sequences, including: acquiring sequencing data and reference genome sequence of circulating tumor DNA to be detected, wherein the sequencing data is the high-resolution analysis of circulating tumor DNA to be detected. The data obtained by flux sequencing, the sequencing data includes: multiple pairs of double-ended sequences; the sequencing data and the reference genome sequence are compared to obtain a first alignment result, wherein the first alignment result at least includes: multiple pairs of double-ended sequences The genome position, base sequence and corresponding base quality value sequence of ; based on the first alignment result, at least one sequence set is obtained, wherein each sequence set includes: at least one pair of double-ended sequences, the The genomic positions of the double-ended sequences are the same; perform network clustering on at least one pair of double-ended sequences in each sequence set to obtain at least one first network; obtain the largest number of double-ended sequences in each first network, and obtain each The independent sequence corresponding to the first network.

Further, performing network clustering on at least a pair of double-ended sequences in each sequence set to obtain at least one first network includes: comparing at least a pair of double-ended sequences in each sequence set to obtain at least one node and the number of members of each node, where the base sequences of the double-ended sequences corresponding to the same node are the same, and the number of members is used to represent the number of double-ended sequences corresponding to the same node; obtain any sequence in each sequence set. Edit distance between two nodes; determine whether the edit distance between any two nodes is less than the preset distance; if the edit distance between any two nodes is less than the preset distance, determine any two nodes belong to the same network, obtain at least one first network.

Further, obtaining the double-ended sequence with the largest number in each first network, and obtaining the independent sequence corresponding to each first network includes: Step A, obtaining the first node with the largest number of members in each first network, and calculating the first node. The sum of the base quality values of all bases contained in each pair of paired-end sequences corresponding to a node is obtained, the sum of the base qualities of each pair of paired-end sequences corresponding to the first node is obtained, and the maximum base quality is added to the corresponding The double-ended sequence is used as an independent sequence corresponding to each first network; in step B, the first node and the nodes adjacent to the first node are deleted from each first network to obtain at least one second network; Step C, judging whether any second network includes a node; Step D, if any second network includes a node, use any second network as the first network, and execute Step A and Step B cyclically.

Further, before it is determined that any two nodes belong to the same network and at least one first network is obtained, the above method further includes: judging whether the number of members of any two nodes satisfies a preset condition; If the number of members satisfies the preset condition, it is determined that any two nodes belong to the same network, and at least one first network is obtained.

Further, judging whether the number of members of any two nodes satisfies the preset condition includes: obtaining a preset multiple of the number of members of the first node in any two nodes, and obtaining the first value; obtaining the first value and the preset The difference between the values; determine whether the number of members of the second node in any two nodes is greater than the difference; if the number of members of the second node is greater than the difference, determine that the number of members of any two nodes satisfies the preset condition .

Further, after obtaining the largest number of double-ended sequences in each first network, and obtaining the independent sequence corresponding to each first network, the above method further includes: dividing the first comparison result by the number corresponding to each first network. The alignment results of other paired-end sequences other than independent sequences are deleted to obtain a second alignment result.

Further, after deleting the alignment results of other double-ended sequences except the independent sequence corresponding to each first network in the first alignment result, after obtaining the second alignment result, the above method also includes: according to each The genomic position of the independent sequence, the second alignment result is sorted, and the third alignment result is obtained.

Further, after sorting the second alignment result according to the genomic position of each independent sequence to obtain the third alignment result, the above method further includes: acquiring a capture sequencing interval; Perform single nucleotide variation detection and indel detection to obtain the detection results.

Further, before obtaining at least one sequence set based on the first alignment result, the method further includes: sorting the first alignment result according to the genomic positions of the multiple pairs of double-ended sequences to obtain a fourth alignment result, and establishing an index for the fourth alignment result; filtering the fourth alignment result to obtain a fifth alignment result; and obtaining at least one sequence set based on the fifth alignment result.

Further, the sequencing data is compared with the reference genome sequence, and obtaining the first comparison result includes: obtaining the matching degree of each sequence in the multiple pairs of double-ended sequences and each sequence in the reference genome sequence; obtaining the highest matching degree corresponding to At least one sequence of each sequence is obtained, and the matching sequence of each sequence is obtained; according to the matching sequence of each sequence, the genomic position of each sequence is determined.

According to another aspect of the embodiments of the present invention, there is also provided an apparatus for processing circulating tumor DNA repeat sequences, including: a first acquisition module for acquiring sequencing data and reference genome sequences of circulating tumor DNA to be detected, wherein the sequencing The data is the data obtained by high-throughput sequencing of the circulating tumor DNA to be detected. The sequencing data includes: multiple pairs of double-ended sequences; an alignment module for aligning the sequencing data with the reference genome sequence to obtain the first alignment result, Wherein, the first alignment result at least includes: genomic positions, base sequences and corresponding base quality value sequences of multiple pairs of double-ended sequences; a processing module, configured to obtain at least one sequence set based on the first alignment result, wherein , each sequence set includes: at least a pair of double-ended sequences, and the genomic positions of the double-ended sequences in the same sequence set are the same; a clustering module is used to perform network clustering on at least a pair of double-ended sequences in each sequence set. class, to obtain at least one first network; a second acquisition module, used to obtain the largest number of double-ended sequences in each first network, and obtain an independent sequence corresponding to each first network.

According to another aspect of the embodiments of the present invention, a storage medium is also provided, the storage medium includes a stored program, wherein when the program is run, the device where the storage medium is located is controlled to execute the above-mentioned method for processing circulating tumor DNA repeats.

According to another aspect of the embodiments of the present invention, a processor is also provided, and the processor is configured to run a program, wherein the above-mentioned method for processing circulating tumor DNA repeats is executed when the program runs.

In the embodiment of the present invention, the sequencing data of the circulating tumor DNA to be detected and the reference genome sequence are obtained, the sequencing data and the reference genome sequence are compared to obtain a first comparison result, and based on the first comparison result, at least one sequence is obtained set, and further perform network clustering on at least a pair of double-ended sequences in each sequence set to obtain at least one first network, obtain the largest number of double-ended sequences in each first network, and obtain the corresponding Independent sequences, thus realizing the deduplication of repetitive sequences. It is easy to notice that, after aligning the sequencing data with the reference genome sequence, it is necessary to perform network clustering on the double-ended sequences with the same genomic position, and select the largest number of double-ended sequences in each network to obtain independent sequences. Therefore, while considering the sequence quality value, the differences in specific sequences are considered, so as to retain more original molecules, reduce manual errors, and improve the technical effect of processing accuracy, thereby solving the processing method of sequencing data in the prior art. Repeated sequence deletion or labeling for sample sequencing is a technical problem with low accuracy.

Description of drawings

The accompanying drawings described herein are used to provide further understanding of the present invention and constitute a part of the present application. The exemplary embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute an improper limitation of the present invention. In the attached image:

1 is a flowchart of a method for processing circulating tumor DNA repeats according to an embodiment of the present invention;

2 is a schematic diagram of a first alternative method for processing circulating tumor DNA repeats according to an embodiment of the present invention;

3 is a schematic diagram of a second alternative method for processing circulating tumor DNA repeats according to an embodiment of the present invention;

4 is a schematic diagram of a third alternative method for processing circulating tumor DNA repeats according to an embodiment of the present invention;

FIG. 5 is a flowchart of an alternative method for processing circulating tumor DNA repeats according to an embodiment of the present invention; and

Fig. 6 is a schematic diagram of a processing device for circulating tumor DNA repeats according to an embodiment of the present invention.

Detailed ways

In order to make those skilled in the art better understand the solutions of the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only Embodiments are part of the present invention, but not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by persons of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

It should be noted that the terms "first", "second" and the like in the description and claims of the present invention and the above drawings are used to distinguish similar objects, and are not necessarily used to describe a specific sequence or sequence. It is to be understood that the data so used may be interchanged under appropriate circumstances such that the embodiments of the invention described herein can be practiced in sequences other than those illustrated or described herein. Furthermore, the terms "comprising" and "having" and any variations thereof, are intended to cover non-exclusive inclusion, for example, a process, method, system, product or device comprising a series of steps or units is not necessarily limited to those expressly listed Rather, those steps or units may include other steps or units not expressly listed or inherent to these processes, methods, products or devices.

Below, the technical terms that appear in the embodiments of the present invention are first explained as follows:

NGS sequencing: High-throughput sequencing (High-throughput sequencing), also known as "Next-generation" sequencing ("Next-generation" sequencing, abbreviated as NGS), is relative to traditional Sanger Sequencing (Sanger Sequencing). It is marked by the ability to sequence hundreds of thousands to millions of DNA molecules in parallel and the general short read length. NGS technology has been widely used in many fields such as diagnosing genetic diseases, measuring gene expression levels, constructing evolutionary trees, distinguishing morphologically similar species, de novo sequencing of non-model organisms, and constructing their reference genome sequences. application. NGS sequencing can be divided into single-end (Single-end, referred to as SE) and pair-end (Paird-end, PE) sequencing. At present, most of them use PE sequencing technology. PE refers to the two ends of the same DNA molecule. Test one end first, then reverse and test the other end. The two sequences tested are called Paired ends Reads (PE Reads), and in many cases, they are also called "Mate Paires".

Capture and sequencing: the genomic region of interest is customized into specific probes and hybridized with genomic DNA on a sequence capture chip (or solution), and the DNA fragments of the target genomic region are enriched and then sequenced using NGS sequencing technology. Research strategy . Due to its short capture area, it greatly reduces the cost of sequencing, and is widely used in various fields. The strategy of capture and sequencing is also used in the field of tumor detection.

ctDNA: Circulating tumor DNA, a DNA fragment that has undergone gene mutation actively secreted into body fluids by tumor cells during the process of division and proliferation.

PCR: The polymerase chain reaction, a technique used to amplify specific DNA fragments.

Repeated sequences: The consequence of multiple mirror copies of a molecule due to PCR amplification.

Edit distance: Refers to the minimum number of edit operations required to convert two strings from one to the other. Permitted editing operations include replacing one character with another, inserting a character, deleting a character, and the smaller the edit distance, the more similar the two strings are.

reads: Sequencing read length, the genome or transcriptome sequence fragments measured by the sequencer.

fasta: A text-based format for representing nucleotide or amino acid sequences.

fastq: A common high-throughput sequencing file type, usually raw sequencing data are stored in this file type.

bwa: an alignment method software used to find the position of the sequenced sequence in the human gene reference sequence, and can output the result file in bam format.

sam: a sequence alignment format used to store the results of sequencing sequences backed into reference genomes

bam: The binary compression format of the sam file, which is used to store the results of the sequencing sequence back to the reference genome.

samtools: A tool for working with bam/sam files.

picard: A tool for processing high-throughput sequencing data, which can be used to process alignment result files such as sam/bam.

Alignment quality: used to quantify the likelihood of an alignment to the wrong location, with higher values indicating lower likelihood.

CIGAR: Brief Alignment Information Expression, which is based on a reference sequence and uses data plus letters to indicate alignment results.

Example 1

According to an embodiment of the present invention, an embodiment of a method for processing circulating tumor DNA repeats is provided. It should be noted that the steps shown in the flowchart of the accompanying drawings may be implemented in a computer system such as a set of computer-executable instructions. and, although a logical order is shown in the flowcharts, in some cases the steps shown or described may be performed in an order different from that herein.

FIG. 1 is a flowchart of a method for processing circulating tumor DNA repeats according to an embodiment of the present invention. As shown in FIG. 1 , the method includes the following steps:

Step S100, obtaining sequencing data and reference genome sequences of the circulating tumor DNA to be detected, wherein the sequencing data is data obtained by high-throughput sequencing of the circulating tumor DNA to be detected, and the sequencing data includes: multiple pairs of double-ended sequences.

Specifically, the above-mentioned circulating tumor DNA to be detected can be the gene sequence extracted from the patient's blood, lymph fluid, interstitial fluid, cerebrospinal fluid and other body fluids. In the embodiment of the present invention, the ctDNA extracted from blood is used as an example. Note; the above-mentioned sequencing data can be the fastq data of ctDNA sample capture sequencing obtained by NGS sequencing of the ctDNA to be detected; the above-mentioned reference genome sequence can be the human reference genome fasta data downloaded from the public database NCBI and other websites.

Step S102, compare the sequencing data with the reference genome sequence to obtain a first comparison result, wherein the first comparison result at least includes: the genome positions, base sequences and corresponding base quality values of multiple pairs of double-ended sequences sequence.

Specifically, the above-mentioned genome position can be the position in which each pair of PE reads is aligned with the reference genome sequence, and different PE reads can be aligned to the same position; the above-mentioned base sequence can be each pair of double-ended sequences in each pair. The base type at the base position. The DNA sequence contains four types of bases, namely G, C, T, and A. During the NGS sequencing process, the base type at each base position can be determined, and Obtain the base quality value of the base type; the above-mentioned base quality value can be the sequencing quality obtained by NGS sequencing, which is used to measure the accuracy of the base type measurement at each base position. The larger the base quality value is , indicating that the accuracy of the base type measurement is higher.

In an optional solution, the human reference genome fasta data and the ctDNA sample capture sequencing fastq data can be obtained, and the genome alignment tool bwa mem can be used to perform sequence alignment to obtain an alignment result file (.bam), that is, to obtain The above-mentioned first comparison result, the comparison result file is in bam format, including the name of each pair of PE reads, position information, SAM marker, comparison quality information, CIGAR string, mate pair information, fragment sequence, sequencing quality, etc.) .

It should be noted that the base sequences and base quality values of the multiple pairs of double-ended sequences in the first alignment result are data directly inherited from the sequencing data of NGS sequencing, and the first alignment result includes the genome alignment position. In addition to the information on the comparison situation, the base sequences and base quality values of multiple pairs of double-ended sequences are also stored, which is convenient for subsequent analysis, and the fastq file is no longer used.

Step S104, obtaining a first sequence set based on the first alignment result, wherein each sequence set includes: at least a pair of double-ended sequences, and the double-ended sequences in the same sequence set have the same genomic position.

Step S106: Perform network clustering on at least a pair of double-ended sequences in each sequence set to obtain at least one first network.

Step S108: Acquire the largest number of double-ended sequences in each first network, and obtain an independent sequence corresponding to each first network.

In an optional solution, since the same DNA molecule fragment (fragment) generates multiple identical fragments through PCR without PCR and sequencing errors, all the fragments can be aligned to the same position in the reference genome, And the sequence is the same; PCR and sequencing errors are small probability events that occur randomly, so the base sequence generated by them accounts for less than the correct base sequence; and DNA molecules from different fragments, although they may be aligned to the same genome. However, because they may belong to different DNA molecules (for example, DNA extracted from blood contains two types, one is ctDNA molecules containing tumor information; the other is self-free DNA in blood, mostly Released from the rupture of cells or white blood cells in the body, it is generally considered to be harmless, and will be cleaned up by the body itself in a short time. The information carried by the two types of DNA is different; for example, different ctDNA molecules), so its sequence may not be same. Based on the above basic facts, the base sequences of all PE reads aligned to the same position in the genome can be clustered by network. It is more likely that the data from two other different fragments are divided into two different classes, so as to obtain at least one first network, and further, the most selected pair of sequences in each first network can be selected as the final representative unique reads, and the rest are considered repetitive sequences. The above scheme also reduces the possibility of the final selection of sequences with PCR or sequencing errors while retaining more original molecules, and improves the accuracy of the results.

It should be noted that, before the deduplication is carried out by the deduplication method provided by the present invention, or after deduplication is carried out by the deduplication method provided by the present invention, all bam files contain the following information: the name of each pair of PE reads, Position information, SAM marker, alignment quality information, CIGAR string, mate pair information, fragment sequence, sequencing quality, etc.

According to the above embodiment of the present invention, the sequencing data of the circulating tumor DNA to be detected and the reference genome sequence are obtained, the sequencing data and the reference genome sequence are compared to obtain a first comparison result, and based on the first comparison result, at least one sequence is obtained set, and further perform network clustering on at least a pair of double-ended sequences in each sequence set to obtain at least one first network, obtain the largest number of double-ended sequences in each first network, and obtain the corresponding Independent sequences, thus realizing the deduplication of repetitive sequences. It is easy to notice that, after aligning the sequencing data with the reference genome sequence, it is necessary to perform network clustering on the double-ended sequences with the same genomic position, and select the largest number of double-ended sequences in each network to obtain independent sequences. Therefore, while considering the sequence quality value, the differences in specific sequences are considered, so as to retain more original molecules, reduce manual errors, and improve the technical effect of processing accuracy, thereby solving the processing method of sequencing data in the prior art. Repeated sequence deletion or labeling for sample sequencing is a technical problem with low accuracy.

Optionally, in the above embodiment of the present invention, in step S106, network clustering is performed on at least a pair of double-ended sequences in each sequence set to obtain at least one first network including:

Step S1061, compare the base sequences of at least a pair of double-ended sequences in each sequence set to obtain at least one node and the number of members of each node, wherein the bases of the double-ended sequences corresponding to the same node are obtained. The base sequences are the same, and the number of members is used to characterize the number of double-ended sequences corresponding to the same node.

Specifically, PE reads can be used as nodes, PE reads with the same sequence can be determined as a node, and the number of members of the node can be obtained according to the number of PE reads with the same sequence, so that the number of members in each sequence set can be further analyzed. Nodes are clustered in the network.

Step S1062: Obtain the edit distance between any two nodes in each sequence set.

Specifically, the edit distance between any two nodes can represent the similarity between any two nodes. The smaller the edit distance, the higher the similarity, and the more likely any two nodes are from the same fragment; the higher the edit distance, the higher the similarity. The larger the similarity, the more likely it is from two different fragments.

It should be noted that the edit distance between any two nodes can be obtained by using the existing edit distance processing method, which is not specifically limited in the present invention.

Step S1063, it is judged whether the edit distance between any two nodes is smaller than the preset distance.

Specifically, the above-mentioned preset distance may be a threshold value set according to actual needs, for example, the preset distance may be 1.

Step S1064, if the edit distance between any two nodes is less than the preset distance, it is determined that any two node columns belong to the same network, and at least one first network is obtained.

In an optional solution, the embodiment of the present invention provides the first network clustering method, that is, the Cluster algorithm, which can align all the sequences in the same position in the genome and the editing distance between sequences is less than a threshold (1 by default). The sequence is connected into a network, with the sequence as the node, the number of members of each node represents the number of times the sequence appears at this position, the entire network can be considered to come from the same original fragment, and within the same group, select A sequence with the largest number of members is determined, and the sequence with the highest base quality and the highest in this sequence is determined as the final representative unique reads.

For example, as shown in Figure 2, for the 6 kinds of PE reads shown in Figure 2, namely AGTC, CGTC, TGTC, TGAC, TAAC, TGAG, it is assumed that the genome positions of the 6 kinds of PE reads are the same, and the editing distance is less than 1, The six PE reads are connected into a network, and further according to the number of members of each node, that is, the number of members of AGTC is 2, the number of members of CGTC is 2, the number of members of TGTC is 365, and the number of members of TGAC is 62, The number of members of TAAC is 80, and the number of members of TGAG is 1. It can be determined that the final representative unique reads are the base quality and the highest TGTC.

Optionally, in the above-mentioned embodiment of the present invention, in step S108, obtaining the double-ended sequence with the largest number in each first network, and obtaining the independent sequence corresponding to each first network includes:

Step A: Obtain the first node with the largest number of members in each first network, calculate the sum of the base quality values of all bases contained in each pair of double-ended sequences corresponding to the first node, and obtain the corresponding first node. The base mass sum of each pair of paired-end sequences, and the maximum base mass and the corresponding paired-end sequence are taken as independent sequences corresponding to each first network.

Step B, delete the first node and the nodes adjacent to the first node from each first network to obtain at least one second network.

Step C, judging whether any second network includes a node.

Specifically, when the second network does not contain nodes, it may be determined that all nodes in the first network are completely traversed.

In step D, if any second network includes a node, any second network is used as the first network, and steps A and B are executed cyclically.

In an optional solution, the embodiment of the present invention provides a second network clustering method, that is, the Adjacency algorithm, which can align all the sequences that are aligned to the same position in the genome and whose editing distance between sequences is less than a threshold (1 by default). Sequences are connected into a network, with sequences as nodes. The number of members of each node represents the number of times the sequence appears at this position. The sequences in the entire network may come from multiple original fragments. Further, uniqe reads The selection rule includes the following steps: first select the node with the largest number of members from the network, and select the base quality and the highest paired-end sequence from the node as the first uniqe reads; The node directly connected to the node; further select the node with the largest number of members among the remaining nodes, select the base quality and the highest paired-end sequence from the node as the second representative reads, and then delete all directly connected to it. nodes; repeat the above method until all nodes are traversed, so as to obtain all independent sequences.

For example, as shown in Figure 3, for the 6 PE reads shown in Figure 3, which are AGTC, CGTC, TGTC, TGAC, TAAC, and TGAG, assuming that the genomic positions of the 6 PE reads are the same, and the editing distance is less than 1, The six PE reads are connected into a network, and further according to the number of members of each node, that is, the number of members of AGTC is 2, the number of members of CGTC is 2, the number of members of TGTC is 365, and the number of members of TGAC is 62, The number of members of TAAC is 80, the number of members of TGAG is 1, and the sequence with the maximum number of members and the highest base quality can be obtained, that is, the TGTC with the highest base quality and the highest quality can be obtained as an independent sequence, and the sequence directly connected to TGTC can be further deleted. Sequence, that is, delete AGTC, CGTC and TGAC, the remaining sequences are TAAC and TGAG, obtain the sequence with the largest number of members and the highest base quality, that is, obtain the base quality and the highest TAAC as an independent sequence, because TAAC and TGAG is not directly connected, therefore, the base quality and the highest TAAC and the base quality and the highest TGAG can be regarded as independent sequences. Therefore, three independent sequences can be obtained, namely the base quality and the highest TGTC, the base quality and the highest TGAG. Base quality and highest TAAC and base quality and highest TGAG.

Optionally, in the above embodiment of the present invention, in step S1064, before it is determined that any two node columns belong to the same network and at least one first network is obtained, the method further includes:

Step S1065, it is judged whether the number of members of any two nodes satisfies a preset condition.

Specifically, the above-mentioned preset conditions are used to represent the conditions that any two nodes can be connected to form a network.

Step S1066, if the number of members of any two nodes satisfies a preset condition, it is determined that any two nodes belong to the same network, and at least one first network is obtained.

In an optional solution, the embodiment of the present invention provides a third network clustering method, that is, the Direct algorithm, which can be aligned to the same position in the genome and the editing distance between sequences is less than a threshold (default is 1) and the same All sequences of adjacent node a and node b that meet certain conditions are connected into a network, with the sequence as the node, the number of members of each node represents the number of times this sequence appears at this position, so that at this position One or more networks can be formed at the location, the sequences in each network are derived from an original fragment, the sequence with the largest number of members is selected, and the sequence with the highest base quality and the highest base quality in this sequence is determined as the final representative.

For example, as shown in Figure 4, for the six PE reads shown in Figure 4, they are AGTC, CGTC, TGTC, TGAC, TAAC, TGAG, where AGTC and CGTC are adjacent, CGTC and TGTC are adjacent, and TGTC and TGTC are adjacent to each other. AGTC is adjacent, TGTC and TGAC are adjacent, TGAC and TAAC are adjacent, and TGAC and TGAG are adjacent. It is assumed that the genomic positions of the six PE reads are the same, and the editing distance is less than 1, and the AGTC, CGTC, TGTC, TGAC and TGAG are similar. If two adjacent nodes meet certain conditions, AGTC, CGTC, TGTC, TGAC and TGAG can be connected into a network, and TAAC is a separate network. Further according to the number of members of each node, that is, the number of members of AGTC is 2, the number of members of CGTC is 2, the number of members of TGTC is 365, the number of members of TGAC is 62, the number of members of TAAC is 80, and the number of members of TGAG is 80. When the number is 1, the sequence with the largest number of members and the highest base quality in each network can be obtained, that is, the base quality and the highest TGTC and the base quality and the highest TAAC can be obtained as independent sequences.

Optionally, in the above-mentioned embodiment of the present invention, step S1065, judging whether the number of members of any two nodes satisfies a preset condition includes:

Step S10652: Obtain a preset multiple of the number of members of the first node in any two nodes to obtain a first value.

Step S10654: Obtain the difference between the first value and the preset value.

Specifically, the above-mentioned preset value can be set according to actual needs, for example, it can be 1.

Step S10656, determine whether the number of members of the second node in any two nodes is greater than the difference.

Step S10658, if the number of members of the second node is greater than the difference, it is determined that the number of members of any two nodes meets the preset condition.

In an optional solution, the number of members of adjacent nodes a and b can be aligned to the same position in the genome and the editing distance between sequences is less than a certain threshold (default is 1) and the number of members of adjacent nodes a and b satisfies: node The number of members of a ≥ 2 * the number of members of node b - 1 all nodes are connected to form a network.

Optionally, in the above-mentioned embodiment of the present invention, in step S108, after obtaining the double-ended sequence with the largest number in each first network, and obtaining the independent sequence corresponding to each first network, the method further includes:

In step S110, the alignment results of other double-ended sequences except for the independent sequences corresponding to each first network in the first alignment result are deleted to obtain a second alignment result.

In an optional solution, all repetitive sequences except independent sequences in all PE reads can be deleted, and only independent sequences are retained to avoid artificial errors caused by repetitive sequences.

Optionally, in the above-mentioned embodiment of the present invention, in step S110, the comparison results of other double-ended sequences other than the independent sequences corresponding to each first network in the first comparison result are deleted to obtain a second comparison. After evaluating the results, the method also includes:

Step S112: Rank the second alignment result according to the genomic position of each independent sequence to obtain a third alignment result.

In an optional solution, due to network clustering based on edit distance, the sequence position in the alignment result file (.bam) changes. In order to make PE reads at the same position adjacent to each other, it is convenient for subsequent PE reads to be removed. For reprocessing, the Picard's SortSam module can be called to sort the comparison result file (.bam) (that is, the above-mentioned second comparison result) according to the comparison position to obtain the third comparison result.

Optionally, in the above embodiment of the present invention, in step S112, the second alignment result is sorted according to the genomic position of each independent sequence, and after obtaining the third alignment result, the method further includes:

Step S114, acquiring the capture sequencing interval.

Step S116, according to the capture sequencing interval, perform single nucleotide variation detection and indel detection on each independent sequence to obtain a detection result.

In an optional solution, a Bed file of the capture sequencing interval can be obtained, and the varscan2mpileup2snp module can be called to detect single nucleotide variation (SNV), and the mpileup2indel module can be used to detect insertion deletion (INDEL), where the single nucleotide variation refers to the reference A change in base type occurs at a certain position in the genome. Indels refer to the insertion of a small new sequence or the deletion of a certain sequence in a certain sequence of the reference genome.

Optionally, in the above-mentioned embodiment of the present invention, in step S104, before obtaining the first sequence set based on the first alignment result, the method further includes:

Step S118: Rank the first alignment results according to the genomic positions of the pairs of double-ended sequences to obtain a fourth alignment result, and establish an index for the fourth alignment result.

In an optional solution, Picard's SortSam module can be called to sort the comparison result file (.bam) (that is, the first comparison result above) according to the comparison position, and at the same time create an index file (.bai) of the bam file. ). The comparison result files are sorted according to the comparison positions, so that the PE reads in the same position are adjacent to each other, which is convenient for subsequent deduplication processing of PE reads.

Step S120, filtering the fourth comparison result to obtain the fifth comparison result.

In an optional solution, since the same PE reads may be compared to multiple genomic locations, before deduplication processing, the result file (.bam) needs to be filtered first. Specifically, the samtools view module can be called to The sorted bam files are filtered to obtain the fifth alignment result.

Step S122, based on the fifth alignment result, obtain at least one sequence set.

In an optional solution, the same position of the genome in the fifth alignment result can be divided into a sequence set, so as to facilitate subsequent network aggregation of the base sequences of all PE reads in each sequence set kind.

Optionally, in the above-mentioned embodiment of the present invention, in step S102, the sequencing data is compared with the reference genome sequence, and obtaining the first comparison result includes:

Step S1022, obtaining the matching degree of each sequence in the multiple pairs of double-ended sequences and each sequence in the reference genome sequence.

Step S1024: Obtain at least a sequence corresponding to the highest matching degree, and obtain a matching sequence of each sequence.

Specifically, the above-mentioned preset similarity may be set according to actual detection requirements, which is not specifically limited in the present invention.

Step S1026: Determine the genomic location of each sequence according to the matching sequence of each sequence.

In an optional solution, the matching degree between each read in each pair of PE reads and the human reference genome sequence can be calculated, and whether each read is from a certain segment of the human reference genome sequence can be judged by the matching degree. The higher the probability that each read comes from the sequence in the human reference genome sequence, the more likely each read is to be aligned to the sequence with the highest matching degree, so that the genomic position of the read can be obtained according to the position of the sequence.

It should be noted that, in this embodiment of the present invention, a comparison algorithm provided in the prior art may be used for comparison, which is not specifically limited in the present invention.

FIG. 5 is a flowchart of an optional method for processing circulating tumor DNA repeats according to an embodiment of the present invention. A preferred embodiment of the present invention will be described in detail below with reference to FIG. 5 . As shown in Figure 5, the method may include the following steps: input cfDNA samples to capture and sequence fastq files and human reference genome fasta files, and use bwa mem software to perform genome alignment; call Picard software to sort reads; call samtools software to filter reads; According to data characteristics or accuracy requirements, a suitable algorithm can be selected from the Cluster, Adjacency and Direct algorithms provided in the above-mentioned embodiments of the present invention for deduplication; Picard software is called to sort reads to obtain the bam file of the cfDNA sample after deduplication; input Capture the Bed file of the sequencing interval, and call samtools mpileup to display the alignment and quality values of all reads by position in the labeled bam file; call the varscan2mpileup2snp module to identify SNVs, and the mpileup2indel module to identify INDELs.

It should be noted that the above-mentioned cfDNA samples may also be other body fluid samples containing ctDNA.

The input file of the present invention includes: a sequencing data file (bam format, including the name, SAM marker, position information, alignment quality information, CIGAR string of each sequencing fragment) generated after the samples to be tested are compared, sorted, filtered, etc. , mate pair information, fragment sequence, sequencing quality, etc.), human reference genome sequence (fastq format);

The output file of the present invention includes: a comparison result file (bam format) after the marking of the sample to be tested is repeated, and a vcf format file of the detected SNV and INDEL.

Through the above scheme, the network clustering method can not only retain more original molecules, but also eliminate the influence of most random errors, thereby reducing the influence of random errors in the final mutation detection in a certain procedure, making the mutation detection more accurate. ,reliable. For samples or sequencing solutions with severe fragmentation of DNA molecules, small coverage of the genome, and multiple rounds of PCR, especially the captured sequencing data of plasma ctDNA samples, more original molecules can be retained, the base sequence can be effectively used, and the accuracy of the original data can be improved. utilization, and ultimately the accuracy of variant detection.

The above embodiment is verified by the single base variation (SNV) gradient dilution cell line test experimental verification below.

1. Cell line culture

Cell lines HCT116, KYSE450, NCI-H1573, NCI-H1975, NCI-H441, PC-9, SK-HEP-1, SW48, THP-1, BEAS-2B were purchased from Nanjing Kebai Biotechnology Co., Ltd. The cell culture was carried out according to the instructions, that is, 10% fetal bovine serum was added to the RPMI-1640 medium, and the culture was carried out under the condition of 37 degrees.

2. Cell DNA extraction

After collecting the cell suspension, centrifuge at 300g at room temperature for 5 minutes, discard the supernatant, resuspend the cells with 200uL PBS, and then use QIAamp DNA Mini Kit (Cat. No. 51304; Qiagen, Germany) for genomic DNA extraction. After lysis, it was purified by column, and finally the DNA was eluted with low-TE buffer.

3. Determine the theoretical VAF of the mutation sites in the above cell lines by ddPCR

Using the genomic DNA extracted from the cells as a template, ddPCR experiments were performed, and the theoretical VAFs of the mutation sites in the above cell lines are shown in Table 1. ddPCR uses Biohler's instruments, commercial probes and reaction systems. The reaction system consists of: 10ulddPCR supermix for probes (no dUTP), 1ul mutant probe, 1ul wild-type probe, and 20ng DNA to be tested. After the reaction system was prepared, chyle was generated according to the method of using the instrument, and the chyle was sucked into a 96-well PCR plate, and then heat-sealed with Pierceable Foil Heat Seal. The conditions of PCR reaction were: enzyme activation at 95°C for 8 min; melting at 94°C for 30s, annealing and extension at 55°C for 1 min, a total of 39 cycles; enzyme inactivation at 98°C for 10 min; incubation at 4°C. After PCR amplification, Biole's droplet reader reads the number of fluorescent droplets in each reaction well. Ultrapure water was used as a negative control for each batch of reactions. Three replicate wells were made for each DNA to be tested as technical replicates.

Table 1

4. Preparation of samples containing 11 mutation sites

The 10 cell lines in the above table were mixed according to the mass percentages in Table 2 below to prepare 1 sample, and the expected VAF value was calculated.

Table 2

5. ddPCR results of samples

The VAF value of each site in the above list was detected by ddPCR experiment. As shown in Table 3, 20ng of sample DNA was added to each reaction system, and three replicates were made for each sample as technical replicates.

table 3

Gene mutation DDPCR VAF KRAS G13D 0.53 PIK3CA H1047R 1.06 EGFR G719S 0.88 NRAS Q61K 1.80 EGFR L858R 1.26 EGFR T790M 1.52 KRAS G12V 1.43 EGFR E746_A750del 4.76 BRAF V600E 0.92 EGFR S768I 2.42 NRAS G12D 4.48

6. Library construction, capture and sequencing of samples

The DNA of the mixed cell line samples was firstly broken into DNA fragments of about 200bp with covaris ultrasound. After qubit fluorescence quantification, as shown in Table 4, different starting amounts of DNA were used. The hyperpreparation kit (Roche) was used to construct the library. After end repair, polyA addition to the 3' end, ligation of sequencing adapters, unbiased amplification, and purification, the library was obtained.

Table 4

sample Input DNA (ng) PCR cycle number Sample 1 20 6 Sample 2 5 8 Sample 3 5 8

Details are as follows:

1) The ends are flush and A is added at the 3' end: the reaction system is shown in Table 5 below:

table 5

reagent volume Fragmented, double-stranded DNA 50μL End Repair&A-Tailing Buffer 7μL End Repair&A-Tailing Enzyme Mix 3μL total capacity 60μL

Buffer and enzyme should be mixed in the EP tube in advance, and the DNA should be vortexed and mixed according to the following reaction. The reaction steps are shown in Table 6 below:

Table 6

This step sets the PCR tube lid temperature to 85°C instead of 105°C. If the next experiment is performed immediately after this operation, the termination temperature should be set to 20°C instead of 4°C.

2) Ligation adapters: According to the instructions of the library construction instructions, 20ng DNA should use 7.5uM adapters. Prepare the reaction system as shown in Table 7 below:

Table 7

reagent volume reaction product 60μL Joint volume 5μL Ultra-pure water 5μL connection buffer 30μL DNA ligase 10μL total capacity 110μL

Buffer and enzyme should be mixed in an EP tube in advance, vortexed, centrifuged, and incubated at 20°C for 15 minutes.

3) Purification after ligation: Add 88 ul of Agencourt AMPure XP purification magnetic beads to the reaction system (110 ul) of the previous step.

Vortex well and centrifuge gently. Adsorb at room temperature for 5-15 minutes to fully bind DNA and magnetic beads to the EP tube and place it on a magnetic stand to adsorb until the liquid is clear. Slowly aspirate the supernatant from the EP tube and discard. Add 200 μL of 80% ethanol to the EP tube and incubate for 30 seconds. Slowly aspirate the ethanol in the EP tube and discard. Repeat the ethanol wash of the magnetic beads. Dry the EP tube at room temperature for 3-5 minutes until the ethanol is completely evaporated. Remove the EP tube from the magnetic stand, add 22 μL of ultrapure water, vortex, and incubate for 2 minutes at room temperature to elute DNA. Place the EP tube on the magnetic stand to absorb the liquid until the liquid is clear, transfer the supernatant to a new EP tube, and take 1 μL The DNA concentration in the supernatant was measured, and the remainder was amplified.

4) PCR amplification: The PCR system was prepared as shown in Table 8 below.

Table 8

reagent volume KAPA HiFi HotStart ReadyMix(2X) 25μL KAPA Library Amplification Primer Mix(10X)* 5μL Linker Libraries 20μL total capacity 50μL

After sufficient shaking, the samples were centrifuged quickly, and the PCR reaction was carried out according to the conditions shown in Table 9 below.

Table 9

5) Purification after amplification: add Agencourt AMPure XP purification magnetic beads (50 μl) in the same volume as the PCR reaction system.

Fully vortexed, centrifuged gently, and adsorbed at room temperature for 5-15 minutes to fully bind the DNA to the magnetic beads. Place the EP tube on a magnetic stand to absorb the liquid until the liquid is clear, slowly aspirate the supernatant from the EP tube and discard. Add 200 μL of 80% ethanol to the EP tube and incubate for 30 seconds, slowly aspirate the ethanol in the EP tube and discard. Repeat the ethanol wash of the magnetic beads. Dry the EP tube at room temperature for 3-5 minutes until the ethanol is completely evaporated. Remove the EP tube from the magnetic stand, add 52 μL of ultrapure water, vortex, and centrifuge gently. Incubate at room temperature for 2 minutes to elute the DNA, place the EP tube on a magnetic stand to adsorb until the liquid is clear, transfer the supernatant to a new EP tube, and take 1 μL of the supernatant to measure the DNA concentration.

6) Using the probe capture method before sequencing, the target region containing 11 mutation sites was enriched and further amplified with the Roche NimbleGen probe to obtain a library of the target region. On-board sequencing was performed after q-PCR quantification.

7. Process the fastq data from the computer as an input file that can be used by each software.

After the data gets off the computer, first process the off-camera data from the fastq file into a bam file. The specific software and steps are as follows:

7.1 Comparison

Call bwa-0.7.12mem to align each pair of fastq files as PE reads to the hg19 human reference genome sequence. Except for the -M parameter and the ID of the specified Reads Group, the other parameter options are not used to generate the initial bam file.

7.2 Sorting

Call the SortSam module of picard-2.1.0, sort the initial bam file according to the chromosome position, and set the parameter as "SORT_ORDER=coordinate".

7.3 Screening

Call samtools-1.3view to filter the sorted bam files, using "-F 0x900" as the parameter.

7.4 Indexing

Call the index module of samtools-1.3 to index the final generated bam file, and generate a bai file paired with the filtered bam file.

8. Mark duplicates

8.1 Use Picard’s MarkDuplicates module to mark duplicates. During subsequent variant detection, these duplicate sequences will be automatically filtered and then analyzed.

8.2 According to the method provided by the above-mentioned embodiment of the present invention, the "Cluster", "Adjacency" and "Direct" methods are respectively invoked to remove duplicate sequences from the filtered bam file, and a duplicate-removed bam file is generated.

8.3 Statistical comparison:

Call the flagstat module of samtools-1.3 to count the final generated bam file, and generate the alignment file of the bam file after deduplication, including the number of total reads, the number of duplicate reads, the number of reads aligned to the reference genome, The number of paired reads data, the number of read1, the number of read2, the number of reads that are perfectly matched to the reference sequence (properly paired), the number of pairs of reads that are aligned to the reference sequence, and only one of a pair of reads matches the reference sequence. The number of matches, the number of pairs of reads that align to different chromosomes, the number of pairs of reads that align to different chromosomes and have an alignment quality value greater than 5, etc.

8.4 Results comparison:

The statistical results of the data volume of the algorithms provided by the above embodiments of the present invention and the Picard method are shown in Table 10 below. It can be seen from Table 10 that the algorithms provided by the present invention all retain more data volume than the Picard method, which improves the data volume. The effective utilization of , and satisfies Adjacency>Direct>Cluster>Picard.

Table 10

sample Picard Cluster Adjacency Direct Sample 1 24872747 51683842 51930648 51698546 Sample 2 13687626 46207524 46730986 46247784 Sample 3 14290322 44097044 44541914 44129370

9. Variation detection

9.1 Stacking

Call samtools-1.3mpileup to display the alignment and quality values of all reads by position in the labeled bam file. The parameter is set to "q=1". The result file of mpileup (mpileup file) contains chromosomes, genome positions, and reference genomes. Base type, sequencing depth at this site, alignments and quality values that all cover the reads at this site.

Due to the limited number of positive sites verified by ddPCR, mpileup processing was performed only on the following intervals, and the parameter "-lpositive.bed" was used. The positive.bed file is shown in Table 11.

Table 11

chromosome starting point end position Gene chr1 115256527 115256530 NRAS chr1 115258745 115258748 NRAS chr3 178952083 178952086 PIK3CA chr12 25398279 25398282 KRAS chr12 25398282 25398285 KRAS chr7 140453134 140453137 BRAF chr7 55241706 55241709 EGFR chr7 55242414 55242513 EGFR chr7 55249003 55249006 EGFR chr7 55249069 55249072 EGFR chr7 55259513 55259516 EGFR

9.2 Statistical average sequencing depth of positive.bed interval

Use a simple script or bash command to count the average sequencing depth in the positive.bed interval for different deduplication methods based on the mpileup file. The results are shown in Table 12.

Table 12

sample Picard Cluster Adjacency Direct Sample 1 1625.370 4111.270 4130.270 4112.240 Sample 2 533.496 3251.800 3302.600 3258.130 Sample 3 627.380 3084.870 3121.860 3087.800

The methods provided by the present invention all have higher average depth in the positive.bed interval than the Picard method, and satisfy Adjacency>Direct>Cluster>Picard.

9.3 Variation detection

Call the varscan2mpileup2snp module to detect single nucleotide variants (SNV), the mpileup2indel module to detect insertion and deletion markers (INDEL), parameter settings "--min-coverage 100--min-reads2 2--min-var-freq0.001--p -value 0.05--min-avg-qual 20".

The statistical variation results of the ddPCR-positive sites of the above three samples after different deduplication methods are shown in the following Tables 13 to 15 (the values in the table are mutation frequencies), among which, Table 13 shows the variation results of sample 1 , Table 14 shows the mutation results of sample 2, and Table 15 shows the mutation results of sample 3.

Table 13

Gene Aachange Picard Cluster Adjacency Direct NRAS p.Q61K 0 1 1.05 1 PIK3CA p.H1047R 0.96 1.24 1.23 1.24 BRAF p.V600E 0.83 0.81 0.8 0.81 NRAS p.G12D 3.87 4.22 4.24 4.21 EGFR p.G719S 0.88 0.77 0.77 0.77 EGFR p.L858R 1.64 1.98 1.97 1.98 EGFR p.S768I 2.15 2.36 2.34 2.36 KRAS p.G13D 0.6 0.52 0.51 0.52 EGFR p.745_750del 3.05 2.79 2.78 2.79 KRAS p.G12V 1.02 1.16 1.15 1.15 EGFR p.T790M 1.39 1.1 1.09 1.1

Table 14

Gene Aachange Picard Cluster Adjacency Direct NRAS p.Q61K 4.22 4.3 4.28 4.3 PIK3CA p.H1047R 0 1.68 1.67 1.68 BRAF p.V600E 0 1.03 1.01 1.02 NRAS p.G12D 0 1.55 1.52 1.55 EGFR p.G719S 0.93 1.14 1.15 1.14 EGFR p.L858R 2.3 2.23 2.23 2.26 EGFR p.S768I 1.04 0.93 0.91 0.93 KRAS p.G13D 1.07 0.88 0.87 0.88 EGFR p.745_750del 2.92 2.05 2.05 2.05 KRAS p.G12V 1.34 1.85 1.88 1.85 EGFR p.T790M 0.96 0.95 0.93 0.95

Table 15

Gene Aachange Picard Cluster Adjacency Direct NRAS p.Q61K 0 1.1 1.09 1.1 PIK3CA p.H1047R 0 0.46 0.46 0.46 BRAF p.V600E 0.99 0.98 0.97 0.98 NRAS p.G12D 5.45 5.86 5.8 5.85 EGFR p.G719S 0 1.26 1.27 1.25 EGFR p.L858R 0.76 0.94 0.97 0.94 EGFR p.S768I 1.66 1.55 1.56 1.55 KRAS p.G13D 0 0.3 0.31 0.3 EGFR p.745_750del 2.56 2.52 2.49 2.52 KRAS p.G12V 1.54 2.06 2.04 2.06 EGFR p.T790M 0 0.95 0.93 0.95

The mutation frequency detected by Picard at multiple positive loci was 0 (frequency>0 is positive, frequency=0 is negative), while the four methods of Adjacency, Direct and Cluster detected positive at all 11 loci. To sum up, it can be seen that more positive sites can be detected by using the present invention than Picard deduplication, and while ensuring positive detection, the Direct and Cluster methods remove more repetitive sequences.

Example 2

According to an embodiment of the present invention, an embodiment of an apparatus for processing circulating tumor DNA repeats is provided.

FIG. 6 is a schematic diagram of an apparatus for processing circulating tumor DNA repeats according to an embodiment of the present invention. As shown in FIG. 6 , the apparatus includes:

The first acquisition module 60 is configured to acquire sequencing data and reference genome sequences of the circulating tumor DNA to be detected, wherein the sequencing data is data obtained by high-throughput sequencing of the circulating tumor DNA to be detected, and the sequencing data includes: multiple pairs of double-ended sequences .

Specifically, the above-mentioned circulating tumor DNA to be detected can be the gene sequence extracted from the patient's blood, lymph fluid, interstitial fluid, cerebrospinal fluid and other body fluids. In the embodiment of the present invention, the ctDNA extracted from blood is used as an example. Explain; the above-mentioned sequencing data can be the fastq data of ctDNA sample capture and sequencing obtained by NGS sequencing of the ctDNA to be detected; the above-mentioned reference genome sequence can be the human reference genome fasta data downloaded from the public database NCBI and other websites, in the embodiment of the present invention Take the reference genome sequence extracted from blood as an example.

The alignment module 62 is configured to compare the sequencing data with the reference genome sequence to obtain a first alignment result, wherein the first alignment result at least includes: the genomic positions, base sequences and corresponding pairs of double-ended sequences. A sequence of base quality values.

Specifically, the above-mentioned genome position can be the position in which each pair of PE reads is aligned with the reference genome sequence, and different PE reads can be aligned to the same position; the above-mentioned base sequence can be each pair of double-ended sequences in each pair. The base type at the base position. The DNA sequence contains four types of bases, namely G, C, T, and A. During the NGS sequencing process, the base type at each base position can be determined, and Obtain the base quality value of the base type; the above-mentioned base quality value can be the sequencing quality obtained by NGS sequencing, which is used to measure the accuracy of the base type measurement at each base position. The larger the base quality value is , indicating that the accuracy of the base type measurement is higher.

In an optional solution, the human reference genome fasta data and the ctDNA sample capture sequencing fastq data can be obtained, and the genome alignment tool bwa mem can be used to perform sequence alignment to obtain an alignment result file (.bam), that is, to obtain The above-mentioned first comparison result, the comparison result file is in bam format, including the name of each pair of PE reads, position information, SAM marker, comparison quality information, CIGAR string, mate pair information, fragment sequence, sequencing quality, etc.) .

It should be noted that the base sequences and base quality values of the multiple pairs of double-ended sequences in the first alignment result are data directly inherited from the sequencing data of NGS sequencing, and the first alignment result includes the genome alignment position. In addition to the information on the comparison situation, the base sequences and base quality values of multiple pairs of double-ended sequences are also stored, which is convenient for subsequent analysis, and the fastq file is no longer used.

The processing module 64 is configured to obtain at least one sequence set based on the first alignment result, wherein each sequence set includes: at least one pair of double-ended sequences, and the double-ended sequences in the same sequence set have the same genomic position.

The clustering module 66 is configured to perform network clustering on at least a pair of double-ended sequences in each sequence set to obtain at least one first network.

The second obtaining module 68 is configured to obtain the largest number of double-ended sequences in each first network, and obtain an independent sequence corresponding to each first network.

In an optional solution, since the same DNA molecule fragment (fragment) generates multiple identical fragments through PCR without PCR and sequencing errors, all the fragments can be aligned to the same position in the reference genome, And the sequence is the same; PCR and sequencing errors are small probability events that occur randomly, so the base sequence generated by them accounts for less than the correct base sequence; and DNA molecules from different fragments, although they may be aligned to the same genome. However, because they may belong to different DNA molecules (for example, DNA extracted from blood contains two types, one is ctDNA molecules containing tumor information; the other is self-free DNA in blood, mostly Released from the rupture of cells or white blood cells in the body, it is generally considered to be harmless, and will be cleaned up by the body itself in a short time. The information carried by the two types of DNA is different; for example, different ctDNA molecules), so its sequence may not be same. Based on the above basic facts, the base sequences of all PE reads aligned to the same position in the genome can be clustered by network. It is more likely that the data from two other different fragments are divided into two different classes, so as to obtain at least one first network, and further, the most selected pair of sequences in each first network can be selected as the final representative unique reads, and the rest are considered repetitive sequences. The above scheme also reduces the possibility of the final selection of sequences with PCR or sequencing errors while retaining more original molecules, and improves the accuracy of the results.

It should be noted that, before the deduplication is carried out by the deduplication method provided by the present invention, or after deduplication is carried out by the deduplication method provided by the present invention, all bam files contain the following information: the name of each pair of PE reads, Position information, SAM marker, alignment quality information, CIGAR string, mate pair information, fragment sequence, sequencing quality, etc.

According to the above embodiment of the present invention, the sequencing data of the circulating tumor DNA to be detected and the reference genome sequence are obtained, the sequencing data and the reference genome sequence are compared to obtain a first comparison result, and based on the first comparison result, at least one sequence is obtained set, and further perform network clustering on at least a pair of double-ended sequences in each sequence set to obtain at least one first network, obtain the largest number of double-ended sequences in each first network, and obtain the corresponding Independent sequences, thus realizing the deduplication of repetitive sequences. It is easy to notice that, after aligning the sequencing data with the reference genome sequence, it is necessary to perform network clustering on the double-ended sequences with the same genomic position, and select the largest number of double-ended sequences in each network to obtain independent sequences. Therefore, while considering the sequence quality value, the differences in specific sequences are considered, so as to retain more original molecules, reduce manual errors, and improve the technical effect of processing accuracy, thereby solving the processing method of sequencing data in the prior art. Repeated sequence deletion or labeling for sample sequencing is a technical problem with low accuracy.

Example 3

According to an embodiment of the present invention, an embodiment of a storage medium is provided, the storage medium includes a stored program, wherein when the program is run, the device where the storage medium is located is controlled to execute the above-mentioned method for processing circulating tumor DNA repeats.

Example 4

According to an embodiment of the present invention, an embodiment of a processor is provided, and the processor is configured to run a program, wherein the above-mentioned method for processing circulating tumor DNA repeats is executed when the program runs.

The above-mentioned serial numbers of the embodiments of the present invention are only for description, and do not represent the advantages or disadvantages of the embodiments.

In the above-mentioned embodiments of the present invention, the description of each embodiment has its own emphasis. For parts that are not described in detail in a certain embodiment, reference may be made to related descriptions of other embodiments.

In the several embodiments provided in this application, it should be understood that the disclosed technical content can be implemented in other ways. The device embodiments described above are only illustrative. For example, the division of the units may be a logical function division. In actual implementation, there may be other division methods. For example, multiple units or components may be combined or may be Integration into another system, or some features can be ignored, or not implemented. On the other hand, the shown or discussed mutual coupling or direct coupling or communication connection may be through some interfaces, indirect coupling or communication connection of units or modules, and may be in electrical or other forms.

The units described as separate components may or may not be physically separated, and components shown as units may or may not be physical units, that is, may be located in one place, or may be distributed to multiple units. Some or all of the units may be selected according to actual needs to achieve the purpose of the solution in this embodiment.

In addition, each functional unit in each embodiment of the present invention may be integrated into one processing unit, or each unit may exist physically alone, or two or more units may be integrated into one unit. The above-mentioned integrated units may be implemented in the form of hardware, or may be implemented in the form of software functional units.

The integrated unit, if implemented in the form of a software functional unit and sold or used as an independent product, may be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention is essentially or the part that contributes to the prior art, or all or part of the technical solution can be embodied in the form of a software product, and the computer software product is stored in a storage medium , including several instructions for causing a computer device (which may be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes: U disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), mobile hard disk, magnetic disk or optical disk and other media that can store program codes .

The above are only the preferred embodiments of the present invention. It should be pointed out that for those skilled in the art, without departing from the principles of the present invention, several improvements and modifications can be made. It should be regarded as the protection scope of the present invention.

Claims (11)

1. A method for processing circulating tumor DNA repeat sequences, characterized in that, comprising: Obtaining sequencing data and a reference genome sequence of the circulating tumor DNA to be detected, wherein the sequencing data is data obtained by high-throughput sequencing of the circulating tumor DNA to be detected, and the sequencing data includes: multiple pairs of double-ended sequences; Comparing the sequencing data with the reference genome sequence to obtain a first alignment result, wherein the first alignment result at least includes: the genomic positions, base sequences and corresponding pairs of pairs of double-ended sequences The base quality value sequence of ; Based on the first alignment result, at least one sequence set is obtained, wherein each sequence set includes: at least one pair of double-ended sequences, and the genomic positions of the double-ended sequences in the same sequence set are the same; performing network clustering on at least a pair of double-ended sequences in each sequence set to obtain at least one first network; Obtaining the largest number of double-ended sequences in each first network, and obtaining an independent sequence corresponding to each first network; Wherein, performing network clustering on at least a pair of double-ended sequences in each sequence set to obtain at least one first network includes: performing a network clustering on the base sequences of at least a pair of double-ended sequences in each sequence set By comparison, at least one node and the number of members of each node are obtained, wherein the base sequences of the double-ended sequences corresponding to the same node are the same, and the number of members is used to characterize the double-ended sequences corresponding to the same node. number; obtain the edit distance between any two nodes in each sequence set; determine whether the edit distance between any two nodes is less than the preset distance; if the distance between any two nodes is If the edit distance is less than the preset distance, then it is determined that the any two nodes belong to the same network, and the at least one first network is obtained; The obtaining the double-ended sequence with the largest number in each first network, and obtaining the independent sequence corresponding to each first network includes: Step A, obtaining the first node with the largest number of members in each first network , calculate the sum of the base quality values of all bases contained in each pair of double-ended sequences corresponding to the first node, and obtain the sum of the base qualities of each pair of double-ended sequences corresponding to the first node, Taking the maximum base quality and the corresponding double-ended sequence as the independent sequence corresponding to each first network; Step B, from each first network, the first node sum and the first The nodes adjacent to the nodes are deleted to obtain at least one second network; Step C, determine whether any second network contains a node; Step D, if any second network contains a node, then the Any second network is used as the first network, and the step A and the step B are executed cyclically. 2 . The method according to claim 1 , wherein before it is determined that the any two nodes belong to the same network and the at least one first network is obtained, the method further comprises: 3 . Judging whether the number of members of any two nodes satisfies a preset condition; If the number of members of the any two nodes satisfies the preset condition, it is determined that the any two nodes belong to the same network, and the at least one first network is obtained. 3. The method according to claim 2, wherein judging whether the number of members of the any two nodes satisfies a preset condition comprises: Obtain the preset multiple of the number of members of the first node in the any two nodes, and obtain the first value; obtaining the difference between the first numerical value and the preset value; Judging whether the number of members of the second node in any two nodes is greater than the difference; If the number of members of the second node is greater than the difference, it is determined that the number of members of any two nodes satisfies the preset condition. 4. The method according to claim 1, wherein after obtaining the double-ended sequence with the largest number in each first network and obtaining the independent sequence corresponding to each first network, the method further comprises: A second alignment result is obtained by deleting alignment results of other double-ended sequences in the first alignment result except the independent sequences corresponding to each of the first networks. 5. The method according to claim 4, wherein, in the first alignment result, the alignment results of other double-ended sequences except the independent sequences corresponding to each first network are deleted, After obtaining the second comparison result, the method further includes: The second alignment result is sorted according to the genomic position of each independent sequence to obtain a third alignment result. 6. The method according to claim 5, wherein, after the second alignment result is sorted according to the genomic position of each independent sequence to obtain the third alignment result, the method further comprises: Obtain the capture sequencing interval; According to the capture sequencing interval, single nucleotide variation detection and insertion deletion detection are performed on each independent sequence to obtain the detection result. 7. The method according to claim 1, wherein before obtaining at least one sequence set based on the first alignment result, the method further comprises: According to the genomic positions of the multiple pairs of double-ended sequences, the first alignment result is sorted to obtain a fourth alignment result, and an index is established for the fourth alignment result; The fourth comparison result is filtered to obtain the fifth comparison result; Based on the fifth alignment result, the at least one sequence set is obtained. 8. The method according to claim 1, wherein the sequencing data is compared with the reference genome sequence, and obtaining the first comparison result comprises: Obtain the matching degree of each sequence in the multiple pairs of double-ended sequences and each sequence in the reference genome sequence; Obtain at least one sequence corresponding to the highest matching degree, and obtain the matching sequence of each sequence; The genomic location of each sequence is determined based on the matching sequence of each sequence. 9. A device for processing circulating tumor DNA repeats, comprising: The first acquisition module is used to acquire the sequencing data of the circulating tumor DNA to be detected and the reference genome sequence, wherein the sequencing data is data obtained by high-throughput sequencing of the circulating tumor DNA to be detected, and the sequencing data includes: a plurality of pairs of double-ended sequence; an alignment module for aligning the sequencing data with the reference genome sequence to obtain a first alignment result, wherein the first alignment result at least includes: the genomic positions of the multiple pairs of double-ended sequences , the base sequence and the corresponding base quality value sequence; a processing module, configured to obtain at least one sequence set based on the first alignment result, wherein each sequence set includes: at least one pair of double-ended sequences, and the genomic positions of the double-ended sequences in the same sequence set are the same; a clustering module, configured to perform network clustering on at least a pair of double-ended sequences in each sequence set to obtain at least one first network; a second acquisition module, configured to acquire the largest number of double-ended sequences in each first network, and obtain an independent sequence corresponding to each of the first networks; Wherein, performing network clustering on at least a pair of double-ended sequences in each sequence set to obtain at least one first network includes: performing a network clustering on the base sequences of at least a pair of double-ended sequences in each sequence set By comparison, at least one node and the number of members of each node are obtained, wherein the base sequences of the double-ended sequences corresponding to the same node are the same, and the number of members is used to characterize the double-ended sequences corresponding to the same node. number; obtain the edit distance between any two nodes in each sequence set; determine whether the edit distance between any two nodes is less than the preset distance; if the distance between any two nodes is If the edit distance is less than the preset distance, then it is determined that the any two nodes belong to the same network, and the at least one first network is obtained; The obtaining the double-ended sequence with the largest number in each first network, and obtaining the independent sequence corresponding to each first network includes: Step A, obtaining the first node with the largest number of members in each first network , calculate the sum of the base quality values of all bases contained in each pair of double-ended sequences corresponding to the first node, and obtain the sum of the base qualities of each pair of double-ended sequences corresponding to the first node, Taking the maximum base quality and the corresponding double-ended sequence as the independent sequence corresponding to each first network; Step B, from each first network, the first node sum and the first The nodes adjacent to the nodes are deleted to obtain at least one second network; Step C, determine whether any second network contains a node; Step D, if any second network contains a node, then the Any second network is used as the first network, and the step A and the step B are executed cyclically. 10 . A storage medium, characterized in that the storage medium comprises a stored program, wherein when the program runs, a device where the storage medium is located is controlled to execute the circulating tumor according to any one of claims 1 to 8 Methods of processing DNA repeats. 11 . A processor, wherein the processor is configured to run a program, wherein when the program is run, the method for processing circulating tumor DNA repeats according to any one of claims 1 to 8 is executed.

Images