WO2019132010A1 - Method, apparatus and program for estimating base type in base sequence - Google Patents

Abstract

Provided are a method, apparatus and program for estimating base type in a base sequence. This method involves the following steps: acquiring a lead sequence derived from each nucleic acid molecule (S11); mapping the acquired lead sequence on a reference sequence (S12); classifying the lead sequence on the basis of mapping results (S13); determining a reliability score corresponding to the consensus lead and each base in the consensus lead, on the basis of the classification results (S14); arbitrarily selecting a base type at the specific position from the consensus lead, acquiring the number of consensus leads containing each base at the selected position, with respect to the total number of consensus leads, constructing an error model according to a parametric error model using a reliability score corresponding to the base type at the selected specific position (S15); estimating the incidence according to analysis error for the base type at the specific positon with the set reliability, and estimating, on the basis of incidence, whether the specific base type at the specific position is present in the clonal nucleic acid mixture at the set reliability (S16).

Description

Method, apparatus and program for estimating base type in base sequence

The present invention relates to a technique for estimating a base type in a base sequence.

Cancer cells are said to be generated by mutation of base sequences in oncogenes and tumor suppressor genes (cancer related genes) in normal cells. In recent years, development of a method to distinguish normal cells from cancer cells has been promoted by applying genetic abnormalities of these cancer-related genes. For example, Non-Patent Document 1 and Non-Patent Document 2 disclose that a trace amount of abnormal DNA derived from cancer cells present in blood or stool could be detected. Furthermore, with the advent of molecular targeting drugs that target abnormal proteins translated by mutations in cancer-related genes, it is required to determine the base sequence with high accuracy.

As a result of comprehensive profiling of mutations in tumor specimens by next-generation sequencers, it has been reported that in one tumor, there are multiple genomic clones derived from cancer cells, which have genomes of sequences different from normal cells. This phenomenon is called intratumoral heterogeneity. Since this phenomenon causes the tumor itself to be in a polyclonal mixed state, there is a need for exhaustive detection including low frequency mutations.

In current next generation sequencing technology, errors occur at a certain frequency during nucleic acid amplification and sequencing. Therefore, detection of low frequency mutations of 1 to 5% or less is difficult because the base sequences of genomic DNA subjected to sequence analysis include false positive mutations caused by technical errors. However, in recent years, it has been reported that the application of molecular barcode technology suppresses errors in low frequency mutation detection.

For example, although non-patent document 3 discloses smCounter as a technique for enhancing detection accuracy by statistical processing of base sequence information obtained from a sequence analysis device, higher detection sensitivity and accuracy such as liquid biopsy can be obtained. In the case of lack of precision, it is insufficient.

Armengol G et al., J Mol Diagn. 2016 Jul; 18 (4): 471-9. Pmid: 27155048 Diehl F et al., Nat Med. 2008; 14 (9): 985-90. Pmid: 18670422 Chang Xu et al., BMC Genomics, Vol. 18, page e5, 2017 (doi: 10.1186 / s12864-016-3425-4))

An object of the present invention is to provide a method, an apparatus, and a computer program for carrying out the method for estimating a base type in a base sequence.

The present inventor determines the reliability score for each molecular barcode from the large amount of base sequence data (read information) of the next-generation sequencer obtained from the DNA fragment prepared by adding the molecular barcode to the template DNA, We found that parametric analysis can detect mutations with higher sensitivity than conventional methods.

Also, in order to estimate the base type at a specific position, which exists infrequently from the sequence data, the present inventor calculates an index for efficiently removing an error generated by nucleic acid amplification or sequencing, and the true base type We have developed a statistical method that determines with high accuracy.

The technique for enhancing detection accuracy of the present invention can be expected, for example, to detect variant nucleic acids at an earlier disease stage.

In a first aspect of the present disclosure, a computer is used to provide a first estimation method for estimating the presence of a specific base type present at a specific position in the base sequence of a nucleic acid molecule constituting a clonal nucleic acid mixture.
The first estimation method is
a) obtaining a lead sequence derived from each nucleic acid molecule,
b) mapping the read sequence obtained in step a onto a reference sequence to obtain a mapping result;
c) obtaining a clustering result of the lead sequence based on the mapping result obtained in step b,
d) obtaining consensus lead information consisting of a consensus lead and a reliability score corresponding to each base in the consensus lead based on the clustering result,
e) From the sequence information of each consensus lead contained in the set of consensus reads obtained from each nucleic acid molecule constituting the clonal nucleic acid mixture through steps a to d, the base type present at the position corresponding to the specific position on the reference sequence And the process of selecting their reliability score,
f) obtaining a ratio of the number of consensus reads containing each specific base species at the particular position selected in step e to the total number of consensus leads including the particular position used in step e,
g) Estimate the occurrence rate of appearance of a specific base due to an analysis error at a specified position selected in step e at a set threshold or significance level by a parametric error model using the reliability score obtained in step d Process,
h) When the ratio obtained in step f is significantly higher than the appearance rate obtained in step g, the threshold value set in step g is a nucleic acid molecule having a specific base type at a specific position selected in step e. Or estimating at a significant level as present in the clonal nucleic acid mixture,
including.

In a second embodiment of the present disclosure, a second estimation is performed using a computer to estimate the proportion of a specific base species present at a specific position in a base sequence of a nucleic acid molecule constituting a clonal nucleic acid mixture in the clonal nucleic acid mixture. A method is provided.
The second estimation method is
a) obtaining a lead sequence derived from each nucleic acid molecule,
b) mapping the read sequence obtained in step a onto a reference sequence to obtain a mapping result;
c) obtaining a clustering result of the lead sequence based on the mapping result obtained in step b,
d) obtaining consensus lead information comprising a consensus lead and a reliability score corresponding to each base in the consensus lead based on the clustering result,
e) From the sequence information of each consensus lead contained in the set of consensus reads obtained from each nucleic acid molecule constituting the clonal nucleic acid mixture through steps a to d, the base type present at the position corresponding to the specific position on the reference sequence And the process of selecting their reliability score,
f) obtaining a ratio of the number of consensus reads containing each base species specific to the particular position selected in step e to the total number of consensus leads including the particular position used in step e,
g) Estimate the occurrence rate of a distinctive base due to an analysis error at a specific position selected in step e at a set threshold or significance level, using a parametric error model using the reliability score obtained in step d Process,
h) When the ratio obtained in step f is significantly higher than the appearance rate obtained in step g, the threshold value set in step g is a nucleic acid molecule having a specific base type at a specific position selected in step e. Or a step of presuming presence in clonal nucleic acid mixture at a significant level, and i) containing the base species at the specific position obtained in step f, when it is judged that the base species selected in step e is present in step h Calculating the ratio of the specific base species at the specific position selected in step e at the threshold or significance level set in step g from the number of consensus leads, and adopting it as the ratio present in the clonal nucleic acid mixture,
including.

In a third aspect of the present disclosure, there is provided a third estimation method for estimating a base sequence of a specific region of each nucleic acid molecule constituting a clonal nucleic acid mixture using a computer.
The third estimation method is
a) obtaining a lead sequence derived from each nucleic acid molecule,
b) mapping the read sequence obtained in step a onto a reference sequence to obtain a mapping result;
c) obtaining a clustering result of the lead sequence based on the mapping result obtained in step b,
d) obtaining consensus lead information consisting of a consensus lead and a reliability score corresponding to each base in the consensus lead based on the clustering result,
e) From the sequence information of each consensus lead contained in the set of consensus reads obtained from each nucleic acid molecule constituting the clonal nucleic acid mixture through steps a to d, the base type present at the position corresponding to the specific position on the reference sequence And the process of selecting their reliability score,
f) obtaining a ratio of the number of consensus reads containing each specific base species at the particular position selected in step e to the total number of consensus leads including the particular position used in step e,
g) The base at the specific position confirmed in step e at a set threshold or significance level by the parametric error model using the reliability score obtained in step d corresponding to the base type at the specific position selected in step e Estimating the rate of occurrence due to analysis errors of the species,
h) When the ratio obtained in step f is significantly higher than the appearance rate obtained in step g, the threshold value set in step g is a nucleic acid molecule having a specific base type at a specific position selected in step e. Or judging that the base species at the specific position is unknown if it is present in the clonal nucleic acid mixture at a significant level or lower than the occurrence rate,
obtained by applying i) the result obtained in step h to each consensus lead of the base to create secondary sequence information, and j) changing the base position in step e and repeating steps e to i. Integrating all secondary sequence information of the specific region and adopting it as the base sequence of the specific region of each nucleic acid molecule constituting the clonal nucleic acid mixture,
including.

In a fourth aspect of the present disclosure, there is provided a computer readable program for causing a computer to execute any of the above first to third estimation methods.

In a fifth aspect of the present disclosure, there is provided an estimation device for estimating the presence of a specific base species at a specific position in the base sequence of a nucleic acid molecule constituting a clonal nucleic acid mixture.
The estimation device is
A lead sequence acquisition unit for acquiring a lead sequence derived from each nucleic acid molecule;
A mapping information acquisition unit that supplies the acquired lead array to a mapping device and acquires the mapping result on the reference array by the mapping device;
A clustering result acquisition unit for acquiring a clustering result of the mapped lead sequence;
A consensus lead information acquisition unit comprising a consensus lead and a reliability score corresponding to each base in the consensus lead based on the clustering result;
Equipped with

In a sixth aspect of the present disclosure, there is provided an estimation device for estimating the presence of a specific base species at a specific position in the base sequence of a nucleic acid molecule constituting a clonal nucleic acid mixture.
The estimation device is
A lead sequence acquisition unit for acquiring a lead sequence derived from each nucleic acid molecule;
A mapping information acquisition unit that supplies the acquired lead array to a mapping device and acquires the mapping result on the reference array by the mapping device;
A first clustering result acquisition unit for acquiring a clustering result of the mapped lead sequence;
A second clustering result acquisition unit that executes a local assembly using a lead sequence included in each clustering and obtains a clustering result of the assembly sequence;
A consensus lead information acquisition unit comprising a consensus lead and a reliability score corresponding to each base in the consensus lead based on the clustering result of the second clustering result acquisition unit;
Equipped with

According to the present invention, it is possible to detect infrequent mutations occurring on a part of nucleic acids derived from a body sample of a subject, and to estimate the base type in the base sequence.

The figure which shows the structure of the system which detects the variation | mutation of the DNA sequence which concerns on embodiment of this indication. Figure showing the procedure of a novel mutation detection method in DNA sequences Diagram illustrating FASTQ format, which is a file format of base sequence data Figure showing the specific procedure of mutation analysis Diagram illustrating SAM / BAM format, which is a file format of mapping data Diagram showing the outline of UMI classification work Diagram showing creation of de Bruijn Graph for use in local assembly Diagram to explain statistical model for detecting mutations Figure for explaining the statistical model to detect Large Indel mutation The figure which shows the comparison of the error rate calculated from Phred quality score in Example 1, and the error rate calculated from UQ. FIG. 6 shows the evaluation results of sensitivity for detecting eight types of mutations in a nucleic acid mixture in Example 1. FIG. 8 shows the abundance ratio comparison of mutant molecules estimated from the read ratio in the nucleic acid mixture in Example 2. Block diagram showing internal configuration of mutation detection device Flow chart showing processing executed by control unit of mutation detection apparatus Flow chart showing processing executed by control unit of mutation detection apparatus

Hereinafter, embodiments of the present disclosure will be described with reference to the accompanying drawings.

<Definition of terms>
The "clonal nucleic acid mixture" in the present disclosure refers to a mixture of multiple nucleic acids generated by repeating biological replication from one nucleic acid molecule. The sample from which such clonal nucleic acid is obtained is not particularly limited by species, but includes, for example, human biological samples, tissue, blood, plasma, serum, urine, saliva, mucus drainage, sputum, stool and Tears are illustrated. In addition to the FFPE (formalin fixed paraffin embedded) sample of the tissue, it may be a cell population by culture in vitro. The FFPE sample-derived DNA is often damaged such as short fragments, nicks and gaps, and deamination (urasylation) of cytosine. Since it directly contributes to poor analysis, it is desirable to repair the DNA using a commercially available reagent such as NEBNext FFPE DNA Repair Mix (manufactured by NEB).

The "read sequence" in the present disclosure is base sequence data output from a sequencer.

The “infrequently occurring base species at a specific position” in the present disclosure is not particularly limited, and includes, for example, “infrequent mutations”. Here, the term "low frequency mutation" is intended to mean (a sudden) mutation of a base in a nucleic acid generated in vivo and satisfying the following two conditions a and b.
a) In a nucleic acid molecule, it appears at a frequency of 1 × 10 ^-3 / base or less (ie, a frequency of 1 or less per 1,000 bases).
b) In a sample containing nucleic acid molecules, the proportion of nucleic acid molecules having different sequences of bases at specific positions is 10% or less of the total number of nucleic acid molecules in the sample.

Hereinafter, the present disclosure will be described by taking “low frequency mutation” as an example as “base type at a specific position which exists at low frequency” as an example, but the present disclosure is not limited to detection of low frequency mutations. Absent.

The (sudden) mutation in the present disclosure may be any of substitution, insertion and deletion. Also, the (sudden) mutation may be known or novel.

Embodiment
1. Configuration for Mutation Detection FIG. 1 is a diagram showing a configuration of a mutation detection system (an example of a estimation system) that implements a method of estimating a base sequence according to an embodiment of the present invention. The mutation detection system 100 detects mutations in DNA sequences. The mutation detection system 100 analyzes the fragmented DNA with the sequencer 50, and inputs the lead sequence data obtained as a result of the analysis into the mutation detection apparatus 10 (an example of an estimation apparatus). The mutation detection device 10 acquires lead sequence data from the sequencer 50 and analyzes it to detect mutations in the DNA sequence.

2. Mutation Detection Method FIG. 2 is a diagram showing the procedure of a novel mutation detection method in the DNA sequence of the present embodiment. The procedure of the method for detecting mutations in DNA sequences will be described with reference to FIG.

First, genomic DNA and RNA (clonal nucleic acid) are extracted from a biological sample, fragmented so as to contain DNA (for example, 1000 base or less) of a length suitable for a sequencer, and fragmented clonal nucleic acid is generated (S1). If the biological sample is already fragmented DNA (cell-free DNA, FFPE sample, etc.) or short RNA (non-coding RNA, etc.), the fragmentation step may not be performed.

The present invention is particularly suitably used for analysis of nucleic acids which may contain low frequency mutations. The clonal nucleic acid to be analyzed is DNA or RNA, preferably DNA, more preferably genomic DNA.

Depending on the model of the next-generation sequencer, the sequence length that can be analyzed may be as short as several hundred bp. Therefore, for example, when the clonal nucleic acid is genomic DNA, it is necessary to fragment it into a length that can be analyzed by a sequencer. There is no limitation on the method of fragmentation, and fragmentation may be performed by a known method. Although the present invention is not limited in any way, for example, DNA Shearing System M220 / ME220 (manufactured by Covaris), which physically cleaves using ultrasound, QIAseq FX Library Kit (Qiagen), which cleaves using an enzyme. Commercially available instruments and reagents such as those manufactured by Preferably, physical cutting using ultrasound is used.

Next, a nucleic acid called "molecular barcode" is added to each fragmented clonal nucleic acid (S2). The template nucleic acid thus obtained is copied later, but the molecular barcode can be added to identify which clonal nucleic acid the copy is derived from. In general, when a copy of each nucleic acid molecule constituting a clonal nucleic acid is produced and a lead sequence of the copy is obtained, it is impossible to specify the origin of the original nucleic acid molecule unit only from the base information of the lead sequence. To address this, the molecular barcode is used. In the present disclosure, the molecular barcode given to the clonal nucleic acid may be one or more unique. Also, the sequence length and sequence pattern are not particularly limited, but at least 5 bases long (4 ⁵ = 1024) are required to distinguish 1000 types of clonal nucleic acids, and the sequence is usually 6 to 10 bases long. used. The molecular barcode has various names, and may be called "unique molecular identifier (UMI)", "unique molecular tag (UMT)", or simply "molecular tag". There is no particular limitation on the method of adding the molecular barcode to the clonal nucleic acid, but it can be prepared by binding with an enzyme such as ligase.

Hereinafter, a clonal nucleic acid to which a molecular barcode is added may be referred to as a template nucleic acid or a nucleic acid to be analyzed.

It is easy to perform the below-mentioned nucleic acid amplification procedure by adding a sequence that can be aligned with a primer for PCR (polymerase chain reaction) amplification in addition to the molecular barcode to the template nucleic acid. For preparation of such a template nucleic acid, for example, a commercially available product, ThruPLEX (registered trademark) Tag-Seq Kit (manufactured by Takara Bio Inc.) can be used without limitation. In addition, when the nucleic acid region to be analyzed is determined in advance, addition of the molecular barcode to the clonal nucleic acid and amplification can be performed simultaneously by using a nucleic acid in which the PCR primer capable of amplifying the region and the molecular barcode are linked. it can. Alternatively, template nucleic acid including the target region may be concentrated by the so-called probe capture method or the like. The concentration operation may be performed after obtaining a copy.

The molecular bar coat may further contain a sequence called "stem" or the like. The stem sequence is a sequence that marks the start of the molecular barcode. Although the sequence length and sequence pattern are not particularly limited, for example, a sequence of 8 to 10 bases in length is used. The stem sequence is placed at a position between the clonal nucleic acid and the molecular barcode.

After addition of the molecular barcode to the clonal nucleic acid, the resulting template nucleic acid is amplified to create a library (S3). In the present embodiment, when the sequencer 50 is a next-generation sequencer, a nucleic acid molecule containing a low frequency mutation also increases the amount of nucleic acid, that is, copies are obtained because it requires a certain amount of nucleic acid to perform sequence analysis. There is a need.

Although the method for producing the copy of the template nucleic acid is not particularly limited, a nucleic acid amplification method may be mentioned, and it is preferable to carry out by a method based on PCR from the viewpoint of ease of operation and the like. The sample to be decoded by the sequencer obtained in this manner is called a library.

Mutations in the DNA sequence are analyzed using the nucleotide sequence data obtained by sequencing this library. In the present embodiment, amplification of the template nucleic acid is performed by a PCR-based method, but is not limited thereto. The method for producing the library is not particularly limited. In the library preparation process, the region to be analyzed may be concentrated by PCR method or probe capture method, but there is no limitation on the concentration method in practice.

The library is input to the sequencer 50, and the base sequence is analyzed (S4). That is, the lead sequence can be obtained by nucleotide sequence sequencing by a known sequencing method. The platform for the sequencer is not particularly limited, but next generation sequencers are desirable when the amount of information to be analyzed is enormous, such as genomic DNA. As next-generation sequencing, for example, HiSeq (manufactured by illumina Inc.) (R) and MiSeq ®, Ion Proton ^TM and Ion ^PGM TM (manufactured by Thermo Fisher Scientific Inc.), PacBio ® RS II and PacBio (registered TM) ^Sequel TM (manufactured by Pacific Biosciences, Inc.), MinION, GridION X5, PromethION , and SmidgION (Oxford Nanopore Technologies Inc.) and the like, the base treated from tens of millions of DNA fragments hundreds of millions concurrently Devices which determine the sequence are preferred.

When the analysis by the sequencer 50 is completed, a mutation in the DNA sequence is detected using the analysis result from the sequencer 50 (S5). This step (S5) is performed by the mutation detection device 10.

The file format of the base sequence data (read sequence) obtained from the sequencer 50 is FASTQ format, but may be other format. The FASTQ format is a format known in the art and is as shown in FIG. 3 (A). The meaning of each column of data is as shown in FIG. In the FASTQ format sequence data, a Phred quality score is output as an index indicating the accuracy for base call (base designation) of each base. In the present embodiment, the error rate by the sequence per base is indicated by the error frequency (P _error ). The relationship between Phred quality score (a) and error frequency (P _error ) is expressed by the following equation.
P _error = 10 ^{-a / 10} (/ base)
a = ASCII-33
ASCII code values are output in FASTQ format data.

Phred quality score is calculated automatically by the next generation sequencer. Although the conversion formula to error frequency differs depending on each platform, there is no limited factor in this embodiment.

In the step of incorporating a lead sequence derived from each nucleic acid molecule, so-called lead filtering (also referred to as lead cleaning) processing is performed to remove low quality lead sequences, that is, sequence data itself containing bases with low basecall accuracy (Phred quality score). You may A known technique can be used for this process, and examples thereof include known software such as Cutadapt (Dortmund University of Technology) and Trimmomatic (Aachen University of Technology).

In addition, when the base of a part of the lead sequence is of low quality, so-called trimming processing may be performed to remove only the base sequence of the corresponding part. The trimming process is sometimes called a masking process. A known technique can be used for the treatment, and examples thereof include known software such as Trimmomatic and sickle (https://github.com/najoshi/sickle).

FIGS. 4 (a) and 4 (b) are diagrams showing a more specific procedure of the mutation analysis process (S5) of FIG. Among these, FIG. 4 (a) shows the procedure by data analysis in consideration of PCR errors and sequence errors as described in detail later, and FIG. 4 (b) considers lead alignment. It shows the procedure by the data analysis. As is clear from both figures, FIG. 4 (b) differs from FIG. 4 (a) only in that it includes the steps of "local assembly" (S17). Hereinafter, each step of the mutation analysis process (S5) will be specifically described with reference to both FIG. 4 (a) and FIG. 4 (b).

1) Acquisition of read sequence (S11)
First, read sequence data is acquired from the sequencer 50 (S11).

2) Mapping to reference sequence (S12)
Since the lead sequence is a collection of sequence information of copy molecules derived from various fragments derived from each nucleic acid molecule constituting the clonal nucleic acid, first, mapping processing is performed to classify the lead sequence by comparison with a reference sequence. Here, the reference sequence is a reference sequence, and any sequence can be used. Although the present invention is not limited in any way, for example, sequences registered in NCBI GenBank, DDBJ, UCSC Genome Browser, etc. can be used. For example, when the analysis target is human genomic DNA, Genome Reference Consortium human build 38 (GRCh38) or UCSC human genome 19 (hg19) can be used as reference sequences (these versions are updated as needed, and as necessary) Choose the appropriate version). Further, the region to be the reference sequence can be appropriately selected according to the purpose. It may be the entire length of the sequence included in the record, or only the desired region. When the reference sequence is compared with a consensus read described later, and a base different from the reference sequence is in the consensus read, the base is a candidate for mutation.

In the mapping technology, for biological species including human whose genome sequence has already been decoded, the genome sequence can be generally obtained in FASTA format as a reference sequence. In the present disclosure, it is possible to determine the presence or absence of a mutation by mapping a lead sequence analyzed from a template nucleic acid to a reference sequence, and then comparing it with a next generation sequence. The presence or absence can be detected.

In this embodiment, the read sequence is mapped to the reference sequence acquired above. For example, but not limited to, BWA (Wellcome Trust Sanger Institute), Bowtie (Johns Hopkins University), DNASTAR (registered trademark) SeqMan NGen (made by DNASTAR), Hisat 2 (Johns Hopkins University), Novoalign (NOVACRAFT) Company software etc. can be used for mapping work.

In the present embodiment, the mapping result is output in a known file format called a SAM / BAM format file as shown in FIG. BAM format files are binary files of SAM format files. The file format for outputting the mapping result is not limited to the BAM format. FIG. 5A is a SAM / BAM format file to which the result of mapping to the reference sequence (S12) is output, and FIG. 5B is added with the result of molecular bar code clustering (S13) It is a SAM / BAM format file. The meaning of each column of the SAM / BAM format data is as shown in FIG.

Next, a clustering process is performed on the thus mapped lead sequences to classify the lead sequences derived from the individual template nucleic acids. When the above-mentioned molecular barcode is added to the template nucleic acid, it is possible to classify based on the sequence information of the molecular barcode region present in the lead sequence. Hereinafter, the molecular barcode may be referred to as "UMI".

It is important how to use the obtained lead sequence for analysis, as the Bayesian estimation becomes higher in analysis accuracy as the number of lead sequences used for analysis increases.
In the existing analysis methods, a read sequence having a low Phred quality score is not used for the subsequent analysis because it contributes to poor analysis. Therefore, there has been a problem that while obtaining a large amount of sequence data, the number of read sequences available for analysis decreases.

On the other hand, in the method of the present disclosure, “N” (a base where the base type can not be determined) is regarded as “N” (a base having a low Phred quality score) and “completely matched base” when compared with the reference sequence. It could be used for the subsequent analysis by assuming it as

Even if the mapping results in one group, different molecular barcode read sequences may be mixed. In that case, by grouping based on molecular barcode sequences, it is possible to reorganize into correct groups.

In the existing analysis method, the unaligned read sequence contributes to analysis failure and is not used for the subsequent analysis. Therefore, there has been a problem that while obtaining a large amount of sequence data, the number of read sequences available for analysis decreases.

However, in the present invention, the unmapped (unmapped) lead sequences are left for judgment and used for analysis in comparison with the generated family in the grouping step based on molecular barcode sequences. Determine if you can. In this way, the number of lead sequences that can be used for analysis is increased (the number of non-analyzable lead sequences is reduced).

3) Molecular barcode (UMI) clustering (S13)
In the present embodiment, after obtaining a series of information including mapping information (eg, lead sequence, UMI sequence, mapping position relative to reference sequence) of lead sequences via a SAM / BAM format file, each lead sequence is obtained. UMI sequencing step (clustering of molecular barcodes) is performed, in which the UMI sequences assigned to are determined and the lead sequence groups derived from the same clonal nucleic acid are put together. FIG. 6A shows an example of a state after the lead arrangement is grouped by the mapping start point and end point on the reference arrangement. FIG. 6A exemplifies three groups (Position family) classified based on the mapping position on the reference sequence. FIG. 6 (B) shows an example of the state where it is further divided by the sequence of UMI (molecular barcode). FIG. 6 (B) illustrates two groups (UMI family) classified according to the sequence of UMI (molecular barcode).

In the step of classification of UMI (grouping by arrangement of molecular barcodes), the following steps ad are performed.

Step a: Classify the lead sequence group independently of the UMI sequence based on the position of each lead sequence consisting of the mapping position to the reference sequence. The classified groups are managed as Position family. That is, a read sequence group having the same mapping start point and end point on the reference sequence described in the SAM / BAM format file is managed as a Position family.
Step b: Count the number of lead sequences in each Position family. A flag is set for Position family whose aggregate value is less than the threshold.
Step c: Classify the lead sequences based on the similarity of UMI sequences within each Position family. A group classified based on the similarity of UMI sequences is called UMI family. In addition, when two or more types of UMI exist on the same library, all UMI sequences are combined. Specifically, in the example of FIG. 6, L-UMI and R-UMI are combined. The read sequences in the Position family are further classified based on the similarity of UMI sequences. By setting a similarity threshold and allowing a certain degree of mismatch, mismatched UMIs within the threshold have the same sequence. Depending on the number of occurrences of UMI sequences in the Position family, a read sequence having a UMI sequence having a smaller number of occurrences than the total number of read sequences is flagged on the SAM / BAM file.
Step d: Based on the number of lead sequences in the Position family, reclassification of the final UMI family is performed with limitation to the Position family located in the vicinity. That is, it is re-judged based on the degree of similarity of the UMI arrangement whether the Position family which has set the flag and which has a small number of reads can be integrated into the nearby Position family. If it is determined that integration is possible, the flagged position family with a small number of leads is integrated into the neighboring position family.

In this embodiment, determination of sequence similarity as to whether UMI sequences are identical uses editing distance such as Hamming distance. However, determination of sequence similarity is not limited to this, and, for example, similarity of Zscore-based vectors assuming (similarity) network analysis based on motif appearance frequency, phylogenetic distance and evolution distance used in molecular phylogenetic analysis An indicator such as is also selectable. Furthermore, the criteria for determining the degree of similarity may be changed depending on the length of the UMI sequence and the copy number of the template nucleic acid.

In the work of next-generation sequencing, UMI sequences themselves may have errors due to PCR amplification and sequencing. Therefore, in the present embodiment, the following steps can be performed with respect to the possibility of the occurrence of the error.
Step a: Bases with low Phred quality scores in the sequence data are masked with N (a symbol indicating an arbitrary base) and excluded from comparison targets of sequence similarity.
Step b: Adjust the edit distance threshold to allow for a mismatch between UMI sequences. That is, the read sequence group in the Position family is further classified based on the degree of similarity (such as Hamming distance) of the UMI sequence. By setting a similarity threshold and allowing a certain degree of mismatch, mismatched UMIs satisfying the threshold condition have the same sequence.

Furthermore, WO 2016/176091 suggests a UMI substitution phenomenon that replaces a sequence different from the UMI sequence initially given in the PCR amplification and sequencing process (eg UMI Jumping, PCR Jumping [UMI -Tools: modeling sequencing errors in Unique Molecular Identifiers to improve quantification accuracy, Tom Smith et al. Http://dx.doi.org/ 10.1101 / gr.209601.116.], Index swapping [Characterization and remediation of ample index swaps by non-redundant dual indexing on massively parallel sequencing platforms, Maura Costello et al. (http://dx.doi.org/10.1101/200790)], Characterization and remediation of sample index swaps by non-redundant dual indexing on massively parallel sequencing platforms, Maura Costello et al. (http://dx.doi.org/10.1101/2007 Report / suggestion already at 0)). In order to cope with this possibility, in the case of a library configuration in which there are a plurality of UMI sequences, if at least one UMI is completely matched within the range within the Position family, the step is classified as a lead sequence having the same UMI. Is included. That is, in the case of a library configuration in which a plurality of UMIs exist, if at least one UMI completely matches in the Position family, they are put together as a read sequence having the same UMI.

Although the present disclosure is not limited in any way, for example, when the molecular barcode consists of 6 bases, a mismatch of 2 bases or less in one or both molecular barcodes, or a lead sequence of the molecular barcode region at either end If the six bases match, it is regarded as a lead sequence to be subjected to clustering. This enables clustering that takes into account the possibility of errors present in the UMI sequences themselves.

Here, the error present in the UMI sequence itself refers to a mismatched base, a base call bad base, a mapping error and the like. In the existing analysis methods, only mismatched bases were considered. On the other hand, in the present invention, the analysis accuracy can be improved by considering the base call defect base, the mapping error and the like.

Although the present invention is not limited at all, the presence rate of errors present in the molecular barcode sequence itself is 50% or less, preferably 40% or less, more preferably 35% or less, particularly preferably the molecular barcode base length. Is less than or equal to 30%, clustering is possible in consideration of the possibility of the error.

When molecular barcodes are added to both ends of a clonal nucleic acid, the error present in the UMI sequence itself may be present in only one molecular barcode or in both molecular barcodes according to the invention. In the method, clustering is possible in consideration of the possibility of the error.

　本実施の形態では、リードのアライメントパターンのステータス（ＸＰ：ｓｔａｔｕｓ）情報は、上述のように、０、１、２、３の値を採り、このうちＸＰ＝０の場合が「正常」のステータスを表すものとされている。また、ＸＰ＝１、２、３のときのデータを取り扱うことにより、別染色体にマッピングされるリードの場合や、片方のリードが参照配列上にマッピングされない場合なども、ＵＭＩの分類結果として出力することができる。 In the present embodiment, the classification result of UMI is added to the data of the SAM / BAM format file as described above. That is, as shown in the SAM / BAM format (output format) of FIG. 5 (B), the name of the lead sequence in the first column is [chromosome number]-[mapping start point]-[mapping end point]-[alignment pattern] -[UMI sequence after classification]-[UMI sequence before classification] Furthermore, the information obtained by adding necessary information to BC: Z: sequence: sequence, XD: i: count, XP: i: status is output as option information of the SAM / BAM format file (see Table 1 below). The form of appending is not limited to this form.

In the present embodiment, as described above, the status (XP: status) information of the alignment pattern of the lead takes the values of 0, 1, 2, 3 and, among these, the status of “normal” when XP = 0 Is supposed to represent. Also, by handling data when XP = 1, 2, 3, UMI's classification result is output even in the case of reads mapped to another chromosome, or when one of the reads is not mapped onto the reference sequence, etc. be able to.

In this embodiment, mutation detection is performed based on a SAM / BAM format file to which UMI classification results are added.

4) Local assembly (S17)
The step of this local assembly is a step included in the mutation analysis process shown in FIG. 4 (b). In the step of local assembly, processing is performed to determine a sequence region not smaller than the reference sequence but longer than the sample-specific read sequence length. By comparing this sequence with the reference sequence, it is possible to confirm base substitution, insertion or deletion having a sequence size equal to or greater than the length of the lead sequence. The local assembly of this embodiment is described in Rimmer et al. It is coded with reference to Platypus's algorithm published in (Nature Genetics, 2014), and a two-colored de Bruijn Graph is constructed from a reference sequence and a lead sequence to construct an assembly sequence. The options available for local assembly and their examples are shown in Table 2 below.

In de Bruijn graph, the base sequence is first cut into base sequences of length k (k-mer) to make it a Node (apex). Furthermore, each Node (vertex) is connected by Edge (edge). A k-mer Node-Edge is generated with a plurality of read sequences or reference sequences, and a graph in which nodes having the same k-mer and the same Node are combined becomes a de Bruijn graph (see FIG. 7). In this embodiment, the origin of the base sequence (derived from the reference sequence or the lead sequence, presence or absence of the molecular barcode and its sequence information) and the reliability score (the Phred quality score of the sequence) for Node and Edge of de Bruijn graph Can be given a weight calculated from

(4-1) Creation of de Bruijn graph In this embodiment, de Bruijn graph targets the partial length of the analyzed gene and can change the sequence length (which can be changed by the window size option in Table 2). Create using the read array mapped to the same region as the reference sequence corresponding to the specified target region. Further, the read sequences to be used include the read sequences themselves mapped to the same region and paired reads (that is, unmapped reads not mapped to the reference sequence or reads mapped to the outside of the target region). Also, in this embodiment, the color of the Edge of the graph is set to 0 for the reference array side, 1 for the lead array side, and 0.5 for the portion where both exist is derived from the Edge (edge). It is possible to distinguish.

(4-2) Restructuring of Graph In this embodiment (in the case of FIG. 7, for example), the length of the k-mer is searched in the 5 'end of the constructed graph and circularized. Can be stretched to create a graph again.

(4-3) Detection of Shortest Path In the present embodiment, Dijkstra's algorithm, which is an algorithm based on the best priority search for solving the single-source shortest path problem when the weight of an edge is a non-negative number, in path search on a graph. Adopt (priority queue) (https://en.wikipedia.org/wiki/Dijkstra method or https://en.wikipedia.org/wiki/Dijkstra%27s_algorithm). That is, individual search paths become assemble sequences, and further, Edges specific to the reference sequence or edge specific to the read sequence can be discriminated by the difference in color of Edges on the paths, so that mutation candidates can be identified. However, the route search method and algorithm are not limited to the above.

5) Consensus lead (S14)
As described above, especially in "3) Molecular Barcode (UMI) clustering (S13)", based on UMI sequences, taking into account PCR errors and sequence errors, the same UMI sequences (same UMI family) A common read sequence, ie, a consensus read, is output from the read sequence having. The consensus lead may be simply referred to as a consensus, a consensus sequence, or a consensus lead sequence. FIG. 8A is a view for explaining the processing steps of the consensus read process for obtaining a consensus lead. (See “(a) Data Analysis Taking Account of PCR Error and Sequence Error” below).
On the other hand, in the analysis using the local assembly, the same UMI sequence (the same UMI family as described in the above-mentioned "4) local assembly (S17)" is taken into consideration, taking into account the sequence error and the mapping distance of the lead. An assemble sequence is constructed from the read sequence having the above-mentioned), and based on the above-mentioned "(4-3) detection of the shortest path", among these, the highly reliable assemble sequence is regarded as a consensus read. FIG. 9A is a view for explaining the processing steps for obtaining a consensus lead. (See “(b) Data analysis considering lead alignment (Large Indel detection)” below).
Furthermore, there is another method in “(a) Data analysis in consideration of PCR error and sequence error” (see “(a ′) Another method of data analysis in consideration of PCR error and sequence error”).

(A) Data Analysis Considering PCR Error and Sequence Error In the present embodiment, first, based on the UMI sequence, having the same UMI sequence in consideration of PCR error and sequence error (that is, within the same UMI family) The common lead sequence (consensus lead) is determined from the lead sequence.

In addition, in order to determine whether a specific base type is due to an error in sequence analysis or a base present in a clonal nucleic acid, a consensus read and a reliability score (UQ) for each base in the sequence are used. Do. The reliability score (UQ) is calculated for each base in each UMI family, and the reliability score (UQ) is expressed by the following equation.

When calculating the reliability score (UQ), statistical processing is performed instead of simple majority voting. Consensus leads can be obtained using commercial products such as Connor (https://github.com/umich-brcf-bioinf/Connor).

Before starting the consensus lead process, the following steps a and b are performed.
Step a) Extract UMI family (No = 1, 2, 3 ... j) mapped to each position (No = 1, 2, 3 ... k),
Step b) Assign numbers (No = 1, 2, 3 ... i) to all the read sequences in each UMI family, select the type of base on the corresponding position, the number of read sequences of each base, and the error frequency of the sequence ( Get the data of P _error ).

These processes are performed on data of SAM / BAM file after UMI classification. In the case of the example shown in FIG. 8 (B), in the specific, ie, the j-th UMI family (No = j) mapped to the 38,930th base pair (bp) of chromosome 3, the total number of read sequences is 7 There is a book (the number of the lead arrangement No = 1 ... 7). There are four lead sequences in which the base at position 38, 930 in chromosome 3 is A, and three lead sequences in which C is present. Each error frequency (P _error ) can be obtained from the ASCII information of the 11th column of the SAM / BAM format file (see FIGS. 5B and 5C).

Although the reliability score (UQ) for each base in each UMI family is not particularly limited, for example, in the Bayesian approach, it can be estimated as a posterior distribution. Specifically, in each UMI family, the likelihood for each target base (each base) at the corresponding place is estimated. The target bases are a subset of the whole set θ, and all insertions or deletions detected in all four types of base pairs (A, T, G, C) and all UMI families mapped onto the same reference sequence It is composed of (other obs. Alleles). In the example of FIG. 8B, two kinds of bases A and C are detected at the corresponding place (chr3: 38, 930 bp) of UMI family (No = j), but the configuration of the set θ is A, There are 4 types of T, G, C (insertion and defective bases are not included because they are not detected).

The statistical model of the well-known software SmCounter is a specification with low detection sensitivity to insertion or deletion bases. On the other hand, in the statistical model of the present invention, by taking into consideration the insertion or deletion base at θ, the detection sensitivity also increases for base insertion and base deletion in addition to base substitution.

In the reliability score (UQ) formula, P (UMI _j | each base) is a likelihood function, and P (each base) is a prior distribution showing the appearance rate of the target base. The prior distribution can adopt a uniform distribution such as 1 / θ. However, if there is any information about θ or P (each base) from prior data etc., it is possible to change.

The likelihood function P (UMI _j | each base) is expressed by the following equation.

The likelihood function P (UMI _j | each base) expresses the probability of the generation of the read sequence data of each UMI family from the prior distribution of the target bases, in consideration of the sequence error and the PCR amplification error. For example, considering the example shown in FIG. 8B, when the target base is A (each base = A), the target base is A when the original A is correct in the UMI family. A likelihood P (UMI _j | A) is calculated, in which seven lead sequences having four lead sequences and three C lead sequences are generated.

The likelihood function expresses a sequence error (S _error ) and a PCR amplification error. In the likelihood function equation, the first term indicates the frequency of sequence errors (S _error ) and the second term indicates the frequency of occurrence of PCR amplification errors. In the example of FIG. 8B, when the target base is A (each base = A), the sequence error (error from A to C) is S _error frequency in the case of three read sequences where the target base is C. It occurs in On the other hand, in the case of four read sequences in which the target base is A, the sequence is correctly performed with the frequency of (1-S _error ), and the likelihood under this condition is determined.

In the first paragraph, it is assumed that no PCR amplification error occurs, so corC _p (probability of no PCR amplification error) is multiplied. Furthermore, n _y indicates the number (redundancy) of the read sequences of each constituent base (Y) of the set θ. In the example of FIG. 8B, in the case of Y = A, n _y is four. On the other hand, for T and G which could not actually be detected in the set θ, n _y will be zero, and the calculation of the first term will not be performed. That is, only the likelihood (second term) due to PCR amplification error is calculated.

In the method of the present disclosure, the sequence error (S _error ) is calculated based on P _error calculated from the above-mentioned Phred quality score. However, there is also a report that the possibility that insertion and deletion bases will occur as a sequencing error is lower than single nucleotide substitution (SNV) (Marinier E et al., BMC Bioinformatics. 2015 Jan 16; 16: 10. doi: 10 Schirmer M et al., Nucleic Acids Res. 2015 Mar 31; 43 (6): e37. Doi: 10.1093 / nar / gku1341). Therefore, in the present embodiment, when the target base (each base) is an insertion / deletion base, detection is performed by multiplying P _error by correction value m (1/100) to reduce the error rate, and detection is performed. I am improving the performance.

The second term of the likelihood function is a code based on the frequency of PCR amplification errors (misC _p ). misC _p is calculated from PCR _error (x).

　前者の変異出現の確率はポアソン分布に従う。ポアソン分布のパラメーターはＰＣＲ酵素のエラー出現率（μ）と、ＰＣＲ酵素の増幅効率（λ）とにより期待値として求めることができる。後者の確率は二項分布の期待値（最大尤度）として算出する。よってＰＣＲ_{ｅｒｒｏｒ}（ｘ）は下記式で求められる。

The PCR _error (x) is determined by the product of the probability that one mutation occurs in the PCR cycle of x cycles and the probability that the PCR product of x cycles is included in the PCR product after the final cycle (n).

The probability of the appearance of mutations in the former follows the Poisson distribution. The parameters of Poisson distribution can be obtained as expected values from the error appearance rate (μ) of the PCR enzyme and the amplification efficiency (λ) of the PCR enzyme. The latter probability is calculated as the expected value (maximum likelihood) of the binomial distribution. Therefore, PCR _error (x) can be obtained by the following equation.

In calculating PCR _error (x), it is important to estimate the number of PCR cycles in which a mutation appears, but in this embodiment, this estimation is performed based on the number of lead sequences of each base in the UMI family. . In the example of FIG. 8 (B), since the case where A is a target base corresponds to the case where C appears due to a PCR amplification error, the ratio is 3/7 (43.8%). The proportion of the number of lead sequences is approximately shown by the following formula of the proportion of PCR products derived from PCR amplification errors contained in all PCR products.
Sm (1 + λ _m ) ^{n x} / So (1 + λ _w ) ⁿ
Here, Sm is an initial nucleic acid copy number of PCR amplification error and is 1. So is the first initial nucleic acid copy number, which is 1. λ m is the PCR amplification efficiency for the nucleic acid product of PCR error. λ w is the PCR amplification efficiency for nucleic acid products without PCR errors.

Assuming that the PCR amplification efficiencies of the PCR error product and the normal product are the same (λ _m = λ _w ), since the PCR cycle number n is known, the PCR cycle number x at which the error occurred can be estimated. .

In the second term of the likelihood function P (UMI _j | each base), the error appearance rate (μ) of the PCR enzyme and the amplification efficiency (λ) of the PCR enzyme become important setting parameters.

In this embodiment, in the likelihood function, only two types of PCR amplification error and sequence error are considered and treated as mutually exclusive events, but the effect of interaction (both PCR amplification error and sequence error Cases that arise may apply) and other factors may also be added to the equation. Other factors include mapping errors (in the case where the read sequence is mapped to different places and different places) and base substitution by oxidation reaction.

Through the above steps, the reliability score (UQ) is calculated for each base in each UMI family. The numbers in FIG. Taking the UMI family of j as an example, the reliability score for each of A, T, G and C is calculated for the 38,930th base, and the base with the highest reliability score among them is calculated. Sequence of consensus reads. In the example of FIG. 8 (B), since the reliability score for A is the highest, A is taken as the consensus lead sequence. This method is different from the conventional method in that the reliability score is calculated for the consensus lead.

(B) Data analysis considering lead alignment (Large Indel detection)
In this embodiment, in analysis using a local assembly, a common read sequence, that is, a consensus read is determined from a read sequence having the same UMI sequence (the same UMI family) in consideration of the mapping errors of the sequence and the reads. . FIG. 9A is a diagram showing processing steps for obtaining a consensus lead.

That is, according to the following procedure, the maximum posterior probability P is applied to a route having 10 base or more of Insert and / or deletion from the shortest path (assembly sequence, haplotype candidate, H) of each local assembly and mutation list derived from each. Determine {each haplotype | UMI _j } and link it to mutation detection (described later). The maximum posterior probability P in the present embodiment corresponds to the reliability score (UQ) described in (a) above.

［数７］の（２）において、前半部分はリードのマッピング（アライメント）距離の考慮部分を、後半部分はシーケンスのエラーの考慮部分を示す。ここで、前半部分の「Ｐ{dist_k|μ，σ}：Ｉｎｓｅｒｔサイズ(正規分布μ，σ)に対するリード間の距離の確率」は、同じＵＭＩ配列（同じＵＭＩ Ｆａｍｉｌｙ、ＵＭＩ_ｊ）を有する全リード配列（ｋ∈ＲｅａｄＰａｉｒ）をハプロタイプの配列に再アラインメントして、リード間の距離（ｄｉｓｔ_ｋ）に対する確率を表すものである。なお、本実施の形態における母集団は、全リードの参照配列上のリード間の距離の平均値（μ）と分散（σ^２）によって定義される正規分布を、Ｉｎｓｅｒｔサイズの分布(正規分布)とする。 A description of each of the above mentioned variables is given below. Please refer to FIG. 9 (B).

In (2) of [Equation 7], the first half shows the consideration portion of the mapping (alignment) distance of the lead, and the second half shows the consideration portion of the error of the sequence. Here, “P {dist _k | μ, σ}: Probability of the distance between the leads to the insert size (normal distribution μ, σ)” in the first half is all with the same UMI sequence (same UMI Family, UMI _j ) The read sequence (k∈ReadPair) is realigned to the sequence of the haplotype to represent the probability for the distance between the reads (dist _k ). The population in this embodiment is a normal distribution defined by the mean value (μ) and the variance (σ ² ) of the distances between the leads on the reference sequence of all the leads, and the distribution of Insert size (normal distribution) I assume.

(B-1) Extraction of Lead Sequence The lead sequence used in the present embodiment can be selected by the mapping information to the reference sequence. For example, in this embodiment, a read having I / D / S / H (unmapped read not mapped to the reference sequence) in Cigar (see FIG. 5C) of the read sequence which can be obtained from the mapping result of the SAM / BAM format Extract). Furthermore, the leads having the same UMI are also extracted (even when the above conditions are not met).

(B-2) and re-align the leads of the realignment particular UMI (UMI _j) haplotype (h _i), to obtain the information of the distance matched / mismatched bases between the leads. In the realignment in this embodiment, the Smith-Waterman method is employed. The realignment algorithm may be different from the Smith-Waterman method.

(B-3) Determination of maximum a posteriori probability P {each haplotype | UMI _j } (haplotype estimation by MAP estimation)
Seeking posterior probability of each haplotype relative to the particular UMI (UMI _j), it is assumed that the specific UMI (UMI _j) is generated from (among all haplotypes) haplotypes showing a maximum a posteriori _{(h i).} That is, P {each haplotype | UMI _j } indicates the probability that a particular UMI (UMI _j ) is derived from a particular haplotype (hi). In the subsequent processing, mutation detection is performed by Poisson distribution based on the maximum posterior probability.

(A ′) Another Method of Data Analysis Taking Account of PCR Error and Sequence Error Further, another method of the above-mentioned “(a) Data analysis taking into account PCR error and sequence error” will be described.

この数９は、上述の「（ａ）ＰＣＲエラーやシーケンスエラーを考慮したデータ解析」における［数２］と対応するが、ＰＣＲエラーのＣ_ｅの部分が異なる。 (A'-1) Detection of SNV, MNP, and Small Indels The same as the above-described Large Indel (refer to the above-mentioned "(b) data analysis in consideration of lead alignment"), Bayesian estimation is performed. P {A _i | UMI _k } is to obtain a posteriori probability that the specific UMI is derived from the allele Ai (Equation 9). The likelihood at that time is composed of two conditions of 1) when there is no PCR error and only a base call error occurs, and 2) when a PCR error occurs and there is no bale call error (Equation 10).

Although this number 9 corresponds to [Equation 2] in the above-mentioned “(a) Data analysis in consideration of PCR errors and sequence errors”, the part of C _e of PCR errors is different.

例えば、「Ｓ_{ｅｒｒｏｒ（ｙ）ｊ}」は、数２と同じ定義である。 The explanation of each variable mentioned above is as follows.

For example, "S _{error (y) j} " has the same definition as the number 2.

(A'-2) Extraction of Leads In the lead sequences mapped to the corresponding Positions, leads having mutations at the corresponding positions are selected based on Cigar and NM tags. Furthermore, the same UMI reads are also extracted regardless of mutations. All mutations in the corresponding part are represented by a set θ, and an element Ai composed of A, T, G, C and all mutations (obs Allele) in the corresponding part belongs.

(A'-3) Estimation of PCR error (C _e ) It is assumed that 7 “A” and 3 “T” and “T” are caused by PCR error for all 10 leads having a specific UMI. In this case, "C _e " is calculated as follows.

これはポアソン分布を使った割合である。個々の変数の意味は以下の通りである。
n:ＰＣＲサイクル数
λ:ＰＣＲ増幅効率∈[0,1]
μ:ＰＣＲエラー率（ＰＣＲの１サイクルで1塩基に対して生じるエラー頻度）
G: テンプレートの長さ（本実施の形態では１ｂｐである） Moreover, the definition and calculation method of each value were as follows.
Percentage of PCR products without PCR errors (noFq):

This is a ratio using Poisson distribution. The meaning of each variable is as follows.
n: PCR cycle number λ: PCR amplification efficiency ∈ [0, 1]
μ: PCR error rate (error frequency occurring for one base in one cycle of PCR)
G: Template length (1 bp in this embodiment)

Number of PCR products without PCR errors (no PCR):

Number of PCR products with PCR errors (ePCR):

θ:集合θの要素数 (A,T,C,G,obsAllele) Percentage of PCR products with PCR errors (eFq):

θ: Number of elements of set θ (A, T, C, G, obsAllele)

(A'-4) Determination of maximum a posteriori probability P {each base | UMI _j } An allyl showing the largest posterior probability (among all alleles) by calculating the posterior probability of each allyl for a specified UMI (UMI _j ) Suppose that a specific UMI (UMI _j ) is generated from (each base). That is, P {each base | UMI _j } indicates the probability that a particular UMI (UMI _j ) is derived from a particular allyl (each base). In the subsequent processing, mutation detection is performed by Poisson distribution based on the maximum posterior probability.

6) Error model construction (S15)
Subsequently, in the present embodiment, a base type at a specific position is estimated using consensus reads obtained from all UMI families and a reliability score (UQ) for each base (FIG. 8 (A), FIG. 9 (A )reference).

In order to conduct the analysis, first, a specific base type present at a specific position is arbitrarily selected from a set of consensus leads. Here, "arbitrary" means "all cases" and "all cases". The "specific position" means that the base sequence of the reference sequence and the consensus lead is different. The “specific base type” means a base of a consensus lead different from the base of the reference sequence.

The position and base type to be confirmed in the clonal nucleic acid may be set (selected) arbitrarily. When it is desired to confirm the presence of a known mutant base, the site and base type for the mutation may be selected in advance. In addition, when it is desired to confirm an unknown mutation site, a position where a base type different from the reference sequence is present may be extracted by comparing the reference sequence with the consensus lead. In addition to the base substitution type (SNP as mutation), the base type may also be a base type corresponding to the insertion site or a base type equivalent to deletion (in this case, treatment with no such equivalent base type). In addition, for the position where there is no consensus lead different from the reference sequence, the subsequent estimation step may not be performed.

Next, the ratio of the number of consensus reads containing each specific base species at a specific position to the total number of consensus reads used for selection of base types is obtained.

On the other hand, a parametric error model using a reliability score estimates the occurrence rate of appearance of a specific base by analysis error at a specific position at a set threshold (or significance level).

Then, when the ratio based on the above consensus read number is significantly higher than the appearance rate due to analysis error, it is estimated that the specific base type at the specific base position is present in the clonal nucleic acid mixture at the above threshold (or significant level). .

The basic model of base type estimation builds a parametric error model by Poisson distribution, and performs base type estimation at arbitrary places. Depending on the embodiment, it is possible to change the error model. Besides Poisson distribution, for example, tests with multi-valued variables (prior distribution is Dirichlet distribution), tests with beta binomial distribution with binary variables, supervised learning models, etc. can be used as parametric error models of the present invention.

Base type estimation includes the following steps a to d.
Step a) Select a target location (chr3: 38, 930 in FIG. 8 (B)).
Step b) Extract all UMI families mapped to the corresponding part (in FIG. 8 (B), extract j UMI families).
Process c) Acquire the score of the consensus and UQ in the corresponding part of each UMI family.
Step d) Select only the consensus lead and UQ score different from the reference sequence (A) having the base type to be confirmed.
Step e) An error model (Poisson distribution) is constructed from the information obtained through steps a to d (pattern of consensus lead having selected base species and UQ score).

The parameter λ of the Poisson distribution can preferably be expressed as follows.

ここで、ａｒｇｍｉｎＰ｛Ｘ│ＵＭＩ｝は、塩基Ｘとなったｎ個のＵＭＩ　ｆａｍｉｌｙのＵＱ（Ｘ）の中の最小値を示す。 If m UMI families are obtained at a certain position, and a consensus read of bases X (= A, T, G or C) is obtained with n UMI families among them, Poisson distribution parameter λ of Expected value E [λ] can be calculated.

Here, argminP {X | UMI} indicates the minimum value among UQ (X) of n UMI families that have become base X.

For example, when 1000 UMI families are obtained at chr3: 38, 930, and a consensus read of C is obtained for 5 UMI families among them, the parameter λ of Poisson distribution can be calculated according to the following equation.
λ = 1000 * (1−argminP {C | UMI})
Here, argminP {C | UMI} indicates the minimum value among five UQ scores UQ (c) (P {C | UMI}) for five consensus reads where the base is C.

By obtaining λ as described above, an error model using a Poisson distribution can be constructed.

7) Mutation detection (S16)
Mutations are detected using the error model (Poisson distribution) constructed as described above. That is, when m UMIs are obtained from the error model, it is possible to estimate the probability that n consensus leads can be obtained by the errors.

For example, when 1,000 UMIs are obtained, it is possible to estimate the probability that an error will yield 5 consensus leads. Here, the above-mentioned threshold is also called a significance level, and is generally set to 0.05 (5%), but if necessary, 0.01 (1%) or 0.001 (0.1%) Set to). Then, the probability estimated from the error model is compared with a predetermined threshold (such as 1% level or 5% level). If the estimated probability is lower than a predetermined threshold value, the hypothesis due to an error will be rejected, so it is determined to be the selected base type. This makes it possible to estimate the presence of a base type in which only one of 1000 molecules is present. Thus, the present disclosure can be suitably used for detection of low frequency mutations.

Depending on the embodiment, the error model can be changed as appropriate.

The ratio that the specific base species at the specific base position estimated as described above is present in the clonal nucleic acid mixture can be calculated by the following formula. The ratio is also called mutation frequency.
Ratio (mutation frequency) = (number of UMIs of consensus reads with mutations mapped to mutation sites) / (total UMI number of consensus reads mapped to mutations sites)

The calculated ratio can be regarded as the percentage (low abundance) of low frequency mutations present in the clonal nucleic acid mixture.

Within each UMI family, it can be presumed that the specific base type at the specific base position can be present in the clonal nucleic acid mixture for all bases in the specific region to be analyzed, or the base type at the specific base position is unknown By verifying whether or not it can be estimated and aggregating information in which the specific base type is present, the aggregated sequence information can be the sequence of low frequency variants present in the clonal nucleic acid mixture. That is, it can be estimated that the specific base type at the specific base position can be estimated to be present in the clonal nucleic acid mixture, or the base type at the specific base position can be estimated unknown, the result obtained by verifying It applies to the read to create secondary sequence information, and further integrates all secondary sequence information of a specific region obtained by repeating the above-mentioned process while changing the base position, and each diffusion constituting the clonal nucleic acid mixture It can be employed as the base sequence of a specific region of the molecule. In the position where there is no consensus read different from the reference sequence (the reference sequence and all the consensus reads are the same), it is judged that the specific base type does not exist, and the base at that position of the clonal nucleic acid is identical to the reference sequence. I assume.

If necessary, mutation filtering can be further performed. For example, P.I. Flaherty et al., Nucleic Acids Res. Vol. 40, No. 1, page e2, 2012 (DOI: 10.1093 / nar / gkr861), using the method described above, the mutations detected in normal samples are excluded from the mutations detected in tumor samples You may Thereby, mutations specific to the tumor sample can be extracted (filtered).

3. EXAMPLE Hereinafter, an example to which the method of the present embodiment is applied will be described. The application of the present disclosure is not limited to the following examples.

(1) Example 1: Detection of Infrequent Mutations in Fragmented Nucleic Acid Mixture In order to investigate the detection sensitivity of rare mutations present in the nucleic acid mixture, 4 types of EGFR, KRAS, NRAS and PIK3CA are used. Multiplex I cfDNA Reference Standard Set (horizon) as a human-derived clonal nucleic acid having, at different mixing ratios, a copy of a wild-type gene and a copy of eight mutant genes in which rare mutations occur in the same gene in a human gene group. Company company) was used. In this example, two types of clonal nucleic acids, in which the copy number of the mutant gene was 1% (Part No .: HD778) and 0.1% (Part No .: HD779), were used as the mixing ratio.

The Multiplex I cfDNA Reference Standard Set is a standard product of clonal nucleic acid artificially fragmented into about 160 bp by resembling cell-free nucleic acid (cell-free DNA) in plasma.

The following procedure was performed using 50 ng of each clonal nucleic acid. Since human genomic DNA is composed of about 3 billion base pairs and the average molecular weight of one base is 330, about 3000 copies of genomic DNA are present in 10 ng of clonal nucleic acid. When the mixing ratio of mutant genes is 1%, 30 copies of mutant genes are present, and when the mixing ratio is 0.1%, 3 copies of mutant genes are present.

Using a protocol of ThruPLEX Tag-seq 6S (12) Kit (manufactured by Takara Bio Inc.), a Stem-Loop adapter was added to both ends of 50 ng of clonal nucleic acid by a ligation method to obtain a template nucleic acid. Then, using the obtained template nucleic acid as a template, PCR amplification was performed with a universal primer having an adapter sequence to prepare an amplified nucleic acid product (adapter + template nucleic acid) for Illumina's sequencer. In the ThruPLEX technology, UMI is added to each template nucleic acid before amplification because UMI of random sequence consisting of 6 bases is contained in the Stem-Loop adapter (total 2 UMI is included) ).

In order to sequence only the template nucleic acid derived from the eight gene regions that become the target region, using a custom SureSelect capture probe (manufactured by Agilent), a specific target region including eight gene regions from the amplified nucleic acid Was concentrated. The specific operation was performed according to the manual recommended by ThruPLEX Tag-seq Kit.

Using a massively parallel sequencing approach using Illumina's next-generation sequencer HiSeq, sequences at both ends 100 bp for each sample were performed to obtain base sequence data (read sequence data).

Low quality read sequence data was removed from the obtained read sequence data using a Trimmomatic (version: 0.32).

The read sequence data was mapped to the human reference sequence (hg19) using BWA-MEM (version: 0.7.15). From the resulting SAM / BAM format file, mapping information for the reference sequence and UMI sequence information were obtained for each read sequence.

Molecular barcode clustering was performed on a specific target area based on the mapping position information and UMI sequence information according to the method described above. The conditions for classification are as follows.
Condition 1: Two or less mismatched bases are identical UMI.
Condition 2: Even if the condition 1 is not satisfied, if at least one UMI sequence is identical, the same UMI is used.

Consensus lead generation was performed according to the method described above. Thereafter, the bases of 178, 935, 997-178, 936, 121 bp corresponding to the exon 10 of the PIK3CA gene located on the 3rd chromosome are calculated from the Phred quality score given to the read sequence of the Illumina sequencer of 1) The error rate Perror and the error rate 1-UQ (base) calculated from 2) UQ were calculated, and two types of error rates were compared. The results are shown in FIG. In FIG. 10, using the SAM / BAM file after molecular bar code clustering, the above two types of error rates are calculated for the base with the largest number of UMI in each location, and the Phred scale (−10 log ₁₀ (error It plotted with the rate).

As shown in FIG. 10, the error rate (× in FIG. 10) calculated from the Phred quality score is between 30 and 40 on the Phred scale, but the error rate based on UQ (FIG. 10) ◆) indicates a value over 50, and the error rate was reduced to about 1/10 times over the entire area.

It was verified by the method of the present disclosure whether eight mutant genes could be detected. Furthermore, as comparison with the conventional method, 1) LoFreq (version: 2.1.2 http://csb5.github.io/lofreq/) and 2) SmCounter (version: https://github.com/xuchang116/ Mutation detection was performed using smCounter 170823 acquisition). The results are shown in FIG. In FIG. 11, the gray hatched mutations indicate the detected mutations, REF indicates the number of lead sequences of the reference sequence, ALT indicates the number of read sequences having mutant sequences, and AF is calculated from the ratio of read sequences. Allele frequency is shown.

Here, LoFreq is a mutation detection software reported in a paper (Wilm A et al., Nucleic Acids Res. 2012 40 (22): 11189-11201), and Bernoulli trial is performed based on the sequence error rate calculated from Phred quality score. It is software that is implemented to detect mutations. In FIG. 11, after acquiring consensus leads using Connor (version: 0.5.1 https://github.com/umich-brcf-bioinf/Connor) for SAM / BAM files mapped by BWA Mutation detection was performed.

Here, SmCounter is software that performs error correction based on molecular barcodes by the Bayesian approach to detect mutations (Non-Patent Document 3). FIG. 11 shows the result of mutation detection using the SmCounter, with the SAM / BAM format file subjected to molecular barcode clustering (the disclosed method) as an input file.

As shown in FIG. 11, among the 8 mutations with a 1.0% mixing ratio in LoFreq, the deletion mutation (Deletion) of Δ746-750 of the EGFR gene, the single nucleotide substitution (SNV) of T790M, and the NRAS gene are indicated by hatching. 4 types of SNVs of Q61K and A59T were detected, and the detection sensitivity was 50% (4 types / 8 types). On the other hand, in SmCounter, among 8 types of mutations at a mixing ratio of 1.0%, 6 types of mutations other than insertion mutations of Δ746-750 and V769-D770, which are indicated by hatching, are detected by hatching, detection sensitivity Was 75% (6 types / 8 types). However, when examined with a 0.1% mixing ratio mutation, neither of the eight types of mutant genes could be detected in either software.

As described in FIG. 11, the analysis code of the present invention (in the figure: “unique pipeline”) detects all eight mutant genes with a mixing ratio of 1.0% indicated by hatching, and the detection sensitivity is It was 100% (8 types / 8 types). Furthermore, while the mixing ratio is 0.1%, five kinds of mutations are shown: hatched SNV of EGFR gene L858R and Deletion of Δ746-750, SN12 of NRAS gene G12D and A59T, and SNK of PIK3CA E545K. The detection sensitivity was 62.5% (5 types / 8 types). The p-values of the detected mutations were all on the order of E-06 (10 ^-6 ), which were lower than the threshold (1% level).

As described above, according to the results of Example 1, the detection method of the present embodiment showed a detection sensitivity superior to that of the conventional method, and furthermore, the mutation could be detected even in a very low mutation which could not be detected by the conventional method.

(2) Example 2: Quantification of mutation frequency in a nucleic acid mixture In order to investigate the mixing ratio (mutation frequency) of rare mutations present in a nucleic acid mixture, in 19 gene groups, a copy of a wild-type gene, further Tru-Q7 (1.3% Tier) Reference Standard (Part No. :) with copies of 35 mutant genes in which rare mutations occur in the same gene at various mixing ratios of 1.0% -30%. HD734) (manufactured by Horizon) was used as a template nucleic acid.

As in Example 1, a library for Illumina sequencer was prepared using ThruPLEX Tag-seq 6S (12) Kit, and sequence reads were obtained by HiSeq.

As a result of performing mutation detection in the same manner as in Example 1, 34 types of mutations could be detected out of 35 types of mutations (detection sensitivity 97.1%). The results are shown in FIG. The frequency of these detected mutations was calculated from SAM / BAM file after clustering of molecular barcodes. The calculation method is as follows.

The mutation frequency is calculated by the following equation.
Mutation frequency = (number of UMIs of consensus reads with mutations mapped to mutation sites) / (total UMI number of consensus reads mapped to mutations sites)

Calculated by comparing the mutation frequency calculated by the above equation (the Observed Allele Frequency) and Droplet Digital ^PCR TM mutation frequency obtained by (manufactured by Bio-Rad) (Expected Allele Frequency, the stated value in the product information) did. The results are shown in FIG. As a result, the frequency estimated by the method of the present disclosure was very close to the value determined by Digital PCR, and the correlation coefficient R2 was 0.9945.

As described above, according to the result of Example 2, UMI classification can be realized with high accuracy by clustering of molecular barcodes according to the method of the present disclosure. In addition, low-frequency mutations could be detected with high accuracy even with degraded (fragmented) low-quality genomic DNA.

4. Mutation Detection Device FIG. 13 is a block diagram showing a specific internal configuration of the mutation detection device 10 for performing the mutation detection method described above. The mutation detection device 10 is configured by an information processing device such as a personal computer. The mutation detection device 10 includes a control unit 11 that controls the overall operation, a display unit 13 that performs screen display, an operation unit 15 that a user performs operations, a data storage unit 17 that stores data and programs, and an external device. And an interface unit 18 which communicates with the network, and a RAM (Random Access Memory) 16 which is a main memory for the control unit 11 to perform processing.

The display unit 13 is configured of, for example, a liquid crystal display or an organic EL display. The operation unit 15 includes a keyboard, a mouse, a touch panel, and the like.

The interface unit 18 can connect various devices (printer, communication device, input device, etc.) conforming to an interface such as USB, HDMI (registered trademark), SCSI, etc., and data between the connected device and the mutation detection device 10 And enable communication of control commands. In particular, the mutation detection apparatus 10 acquires the read sequence data from the next-generation sequencer 50 via the interface unit 18.

The control unit 11 controls the entire operation of the mutation detection device 10, and is configured by a CPU or an MPU that realizes a predetermined function by executing a program. The program executed by the control unit 11 can be provided via a communication line, or a recording medium such as a CD, a DVD, a memory card or the like. Alternatively, the control unit 11 may be configured by a dedicated hardware circuit (FPGA, ASIC, etc.) designed to realize a predetermined function.

The program uses, for example, a computer language such as JAVA (registered trademark), C, C ++, C #, OBJECTIVE-C, SWIFT, or a script language such as PERL, PYTHON, BIO-RUBY, or any suitable language. For example, it may be implemented as software code using conventional or object oriented techniques.

The software code may be stored as a sequence of instructions on a computer readable medium for storage and / or transmission, such as RAM 16. Alternatively, it can be transmitted using a carrier signal that is encoded and adapted for transmission over wired, optical, and / or wireless communication lines that conform to various protocols, including the Internet.

The data storage unit 17 is a device that stores data and programs, and can be configured by, for example, a hard disk (HDD), an SSD (SOLID STATE DRIVE), a semiconductor memory device, or an optical disk. The data storage unit 21 stores a control program executed by the control unit 11, read sequence data, and the like.

FIGS. 14A and 14B are flowcharts showing processing executed by the control unit 11 of the mutation detection device 10 of FIG. Among them, the flowchart of FIG. 14A shows a process related to data analysis in consideration of PCR errors and sequence errors. The flowchart of FIG. 14B shows a process related to data analysis in consideration of lead alignment.

First, processing relating to data analysis in consideration of PCR errors and sequence errors shown in FIG. 14A will be described.

The mutation detection device 10 acquires the read sequence data from the sequencer 50 via the interface unit 18 (S51). The control unit 11 of the mutation detection device 10 performs mapping processing of each read sequence data to a reference sequence (S52). The control unit 11 classifies the read sequence data based on the position on the reference sequence and the UMI according to the method described above (S53). Thereafter, the control unit 11 performs a consensus read process to create a consensus read for each UMI family (S54).

The control unit 11 constructs an error model (S55). Specifically, the control unit 11 calculates a reliability score UQ, obtains a parameter λ of an error model (Poisson distribution) from the pattern of consensus leads and the reliability score UQ, and constructs an error model.

Then, the control unit 11 detects a mutation with reference to the error model (S56). Specifically, the control unit 11 selects the position and base type for which the mutation is to be confirmed, and estimates the occurrence probability of the mutation with reference to the error model. Then, when the estimated value is smaller than the threshold value, the control unit 11 determines that no mutation has occurred, and determines that the selected base type is correct. In the flow chart shown in FIG. 14 (a), small sized mutations can be detected.

As an optional step, if necessary, the control unit 11 performs mutation filtering (S57) and excludes a mutation that satisfies a predetermined condition from the mutation detection result.

Next, processing related to data analysis in consideration of lead alignment shown in FIG. 14B will be described. The mutation detection device 10 acquires the read sequence data from the sequencer 50 via the interface unit 18 (S51). The control unit 11 of the mutation detection device 10 performs mapping processing of each read sequence data to a reference sequence (S52). The control unit 11 classifies the read sequence data based on the position on the reference sequence and the UMI according to the method described above (S53).

In the flowchart (process) shown in FIG. 14B, the control unit 11 sets the reference sequence to “clustering of molecular barcodes (S53)” and (following) “consensus read processing (S54)”. A local assembly process (S58) is performed to determine a sequence region longer than the sample-specific read sequence length.

Thereafter, the control unit 11 performs a consensus read process to create a consensus read for each UMI family (S54). Next, the control unit 11 constructs an error model (S55). Specifically, the control unit 11 calculates a reliability score UQ, obtains a parameter λ of an error model (Poisson distribution) from the pattern of consensus leads and the reliability score UQ, and constructs an error model. Furthermore, the control unit 11 detects a mutation with reference to an error model (S56). Specifically, the control unit 11 selects the position and base type for which the mutation is to be confirmed, and estimates the occurrence probability of the mutation with reference to the error model. Then, when the estimated value is smaller than the threshold value, the control unit 11 determines that no mutation has occurred, and determines that the selected base type is correct. Furthermore, if necessary, the control unit 11 performs mutation filtering (S57) and excludes a mutation that satisfies a predetermined condition from the mutation detection result.

In the flowchart (process) shown in FIG. 14 (b), the “local assembly process (S58)” is between the “molecular bar code clustering (S53)” and the (following) “consensus read process (S54)”. The processing is the same as the processing of FIG. However, in the flow chart shown in FIG. 14 (b), mutations of large size can be detected.

As shown in the flow chart of FIG.
Read arrangement processing (S51),
Mapping of reference sequence (S52), and
· Molecular bar code clustering (classification of UMI) (S53)
After the processing of
Consensus lead processing (S54a),
Error model construction (S55a), and
· Mutation detection (small mutation) (S56a)
Process flow for the detection of small mutations,
・ Local assembly processing (S58)
Consensus lead processing (S 54 b),
Error model construction (S55 b), and
-Mutation detection (large mutation) (S56b)
And the processing flow for the detection of large mutations in parallel, and then output data of the union of the detection of “small mutation” and the detection of “large mutation” is created (“combination processing (S 56 c "," May be a process flow.

As described above, the mutation detection apparatus 10 can detect a mutation in the DNA sequence based on the data of the read sequence from the sequencer 50.

(This disclosure)
The above embodiments disclose the following technical ideas.

(1) A method of using a computer to estimate the presence of a specific base type present at a specific position in the base sequence of a nucleic acid molecule constituting a clonal nucleic acid mixture,
a) obtaining a lead sequence derived from each nucleic acid molecule,
b) mapping the read sequence obtained in step a onto a reference sequence to obtain a mapping result;
c) obtaining a clustering result of the lead sequence based on the mapping result obtained in step b,
d) obtaining consensus lead information consisting of a consensus lead and a reliability score corresponding to each base in the consensus lead based on the clustering result,
e) From the sequence information of each consensus lead contained in the set of consensus reads obtained from each nucleic acid molecule constituting the clonal nucleic acid mixture through steps a to d, the base type present at the position corresponding to the specific position on the reference sequence And the process of selecting their reliability score,
f) obtaining a ratio of the number of consensus reads containing each specific base species at the particular position selected in step e to the total number of consensus leads including the particular position used in step e,
g) Estimate the occurrence rate of appearance of a specific base due to an analysis error at a specified position selected in step e at a set threshold or significance level by a parametric error model using the reliability score obtained in step d Process,
h) When the ratio obtained in step f is significantly higher than the appearance rate obtained in step g, the threshold value set in step g is a nucleic acid molecule having a specific base type at a specific position selected in step e. Or estimating at a significant level as present in the clonal nucleic acid mixture,
Method including.

(2) A method of using a computer to estimate the proportion of a specific base species present at a specific position in a base sequence of a nucleic acid molecule constituting a clonal nucleic acid mixture in the clonal nucleic acid mixture,
a) obtaining a lead sequence derived from each nucleic acid molecule,
b) mapping the read sequence obtained in step a onto a reference sequence to obtain a mapping result;
c) obtaining a clustering result of the lead sequence based on the mapping result obtained in step b,
d) obtaining consensus lead information consisting of a consensus lead and a reliability score corresponding to each base in the consensus lead based on the clustering result,
e) From the sequence information of each consensus lead contained in the set of consensus reads obtained from each nucleic acid molecule constituting the clonal nucleic acid mixture through steps a to d, the base type present at the position corresponding to the specific position on the reference sequence And the process of selecting their reliability score,
f) obtaining a ratio of the number of consensus reads containing each base species specific to the particular position selected in step e to the total number of consensus leads including the particular position used in step e,
g) Estimate the occurrence rate of a distinctive base due to an analysis error at a specific position selected in step e at a set threshold or significance level, using a parametric error model using the reliability score obtained in step d Process,
h) When the ratio obtained in step f is significantly higher than the appearance rate obtained in step g, the threshold value set in step g is a nucleic acid molecule having a specific base type at a specific position selected in step e. Or a step of presuming presence in clonal nucleic acid mixture at a significant level, and i) containing the base species at the specific position obtained in step f, when it is judged that the base species selected in step e is present in step h Calculating the ratio of the specific base species at the specific position selected in step e at the threshold or significance level set in step g from the number of consensus leads, and adopting it as the ratio present in the clonal nucleic acid mixture,
Method including.

(3) A method of estimating the base sequence of a specific region of each nucleic acid molecule constituting a clonal nucleic acid mixture, using a computer,
a) obtaining a lead sequence derived from each nucleic acid molecule,
b) mapping the read sequence obtained in step a onto a reference sequence to obtain a mapping result;
c) obtaining a clustering result of the lead sequence based on the mapping result obtained in step b,
d) obtaining consensus lead information comprising a consensus lead and a reliability score corresponding to each base in the consensus lead based on the clustering result,
e) From the sequence information of each consensus lead contained in the set of consensus reads obtained from each nucleic acid molecule constituting the clonal nucleic acid mixture through steps a to d, the base type present at the position corresponding to the specific position on the reference sequence And the process of selecting their reliability score,
f) obtaining a ratio of the number of consensus reads containing each specific base species at the particular position selected in step e to the total number of consensus leads including the particular position used in step e,
g) The base at the specific position confirmed in step e at a set threshold or significance level by the parametric error model using the reliability score obtained in step d corresponding to the base type at the specific position selected in step e Estimating the rate of occurrence due to analysis errors of the species,
h) When the ratio obtained in step f is significantly higher than the appearance rate obtained in step g, the threshold value set in step g is a nucleic acid molecule having a specific base type at a specific position selected in step e. Or judging that the base species at the specific position is unknown if it is present in the clonal nucleic acid mixture at a significant level or lower than the occurrence rate,
obtained by applying i) the result obtained in step h to each consensus lead of the base to create secondary sequence information, and j) changing the base position in step e and repeating steps e to i. Integrating all secondary sequence information of the specific region and adopting it as the base sequence of the specific region of each nucleic acid molecule constituting the clonal nucleic acid mixture,
Method including.

(4) In the method of (3), when the reference sequence adopted in step b and the base sequence on the consensus lead are the same, the base of the secondary information prepared in step i as it is without passing through steps e to h You may adopt as a seed.

(5) In the method according to any one of (1) to (3), the specific base position selected in step e may be a position different from the reference sequence adopted in step b and the base type on the consensus read .

(6) In the method of any of (1) to (3), the nucleic acid molecule may have a molecular barcode at at least one end.

(7) In any of the methods (1) to (6), each nucleic acid molecule may have a different molecular barcode.

In the method of (8) (7), each nucleic acid molecule may have one or more different molecular barcodes.

(9) In any of the methods (1) to (8), the clustering in step c may be performed using the arrangement of the molecular barcode region in the lead sequence as an index.

(10) In the method of (9), the molecular barcode may be a random sequence of base sequences.

(11) In the method of (10), the mismatch in the lead sequence of the molecular barcode region may be within the allowable range of the sequence similarity between the lead sequences of the molecular barcode region.

(12) In any of the methods (1) to (11), the reliability score may be obtained by the Bayesian approach from the information on the clustered lead sequences.

(13) In the method of (12), the lead sequence may be a lead sequence of the amplification product of the nucleic acid molecule.

In the method of (14) (13), an amplification product may be obtained by polymerase chain reaction.

(15) In the method of (13) or (14), the likelihood function of the Bayesian approach may be a function taking into account sequencer errors and amplification errors.

(16) In the method of (15), if it is determined that different bases in the lead sequence to the reference sequence correspond to insertions or deletions in the mapping step in step b, a function with further correction for sequencer error is adopted. May be

(17) In any of the methods (1) to (3), the parametric error model of step g may be based on a Poisson distribution.

(18) In the method of (1) or (2), the specific base type at the specific position may be a target base of low frequency mutation.

(19) In the method of any one of (1) to (3), step c may proceed directly to step d.

(20) In the method of any of (1) to (3), after the step c,
k) performing local assembly using the read sequences included in each clustering obtained in step c, and obtaining clustering results of the assemble sequences;
You may

(21) A method comprising the step of determining a specific base species present at a specific position in the base sequence of the diffusion molecule constituting the clonal nucleic acid mixture based on the results obtained by the methods of (19) and (20). May be

(22) A computer readable program for causing a computer to execute the method according to any one of (1) to (21).

(23) An apparatus for estimating the presence of a specific base type at a specific position in the base sequence of a nucleic acid molecule constituting a clonal nucleic acid mixture,
A lead sequence acquisition unit for acquiring a lead sequence derived from each nucleic acid molecule;
A mapping information acquisition unit that supplies the acquired lead array to a mapping device and acquires the mapping result on the reference array by the mapping device;
A clustering result acquisition unit for acquiring a clustering result of the mapped lead sequence;
A consensus lead information acquisition unit comprising a consensus lead and a reliability score corresponding to each base in the consensus lead based on the clustering result;
An estimation device equipped with

(24) An apparatus for estimating the presence of a specific base type at a specific position in the base sequence of a nucleic acid molecule constituting a clonal nucleic acid mixture,
A lead sequence acquisition unit for acquiring a lead sequence derived from each nucleic acid molecule;
A mapping information acquisition unit that supplies the acquired lead array to a mapping device and acquires the mapping result on the reference array by the mapping device;
A first clustering result acquisition unit for acquiring a clustering result of the mapped lead sequence;
A second clustering result acquisition unit that executes a local assembly using a lead sequence included in each clustering and obtains a clustering result of the assembly sequence;
A consensus lead information acquisition unit comprising a consensus lead and a reliability score corresponding to each base in the consensus lead based on the clustering result of the second clustering result acquisition unit;
An estimation device equipped with

DESCRIPTION OF SYMBOLS 10 Mutation detection apparatus 11 Control part 13 Display part 15 Operation part 16 RAM
17 data storage unit 18 interface unit 50 sequencer 100 mutation detection system

Claims (24)

A method of using a computer to estimate the presence of a specific base species present at a specific position in the base sequence of a nucleic acid molecule constituting a clonal nucleic acid mixture,
a) obtaining a lead sequence derived from each nucleic acid molecule,
b) mapping the read sequence obtained in step a onto a reference sequence to obtain a mapping result;
c) obtaining a clustering result of the lead sequence based on the mapping result obtained in step b,
d) obtaining consensus lead information consisting of a consensus lead and a reliability score corresponding to each base in the consensus lead based on the clustering result,
e) From the sequence information of each consensus lead contained in the set of consensus reads obtained from each nucleic acid molecule constituting the clonal nucleic acid mixture through steps a to d, the base type present at the position corresponding to the specific position on the reference sequence And the process of selecting their reliability score,
f) obtaining a ratio of the number of consensus reads containing each specific base species at the particular position selected in step e to the total number of consensus leads including the particular position used in step e,
g) Estimate the occurrence rate of appearance of a specific base due to an analysis error at a specified position selected in step e at a set threshold or significance level by a parametric error model using the reliability score obtained in step d Process,
h) When the ratio obtained in step f is significantly higher than the appearance rate obtained in step g, the threshold value set in step g is a nucleic acid molecule having a specific base type at a specific position selected in step e. Or estimating at a significant level as present in the clonal nucleic acid mixture,
Method including. A method of using a computer to estimate the proportion of a specific base species present at a specific position in a nucleotide sequence of a nucleic acid molecule constituting a clonal nucleic acid mixture in the clonal nucleic acid mixture,
a) obtaining a lead sequence derived from each nucleic acid molecule,
b) mapping the read sequence obtained in step a onto a reference sequence to obtain a mapping result;
c) obtaining a clustering result of the lead sequence based on the mapping result obtained in step b,
d) obtaining consensus lead information consisting of a consensus lead and a reliability score corresponding to each base in the consensus lead based on the clustering result,
e) From the sequence information of each consensus lead contained in the set of consensus reads obtained from each nucleic acid molecule constituting the clonal nucleic acid mixture through steps a to d, the base type present at the position corresponding to the specific position on the reference sequence And the process of selecting their reliability score,
f) obtaining a ratio of the number of consensus reads containing each base species specific to the particular position selected in step e to the total number of consensus leads including the particular position used in step e,
g) Estimate the occurrence rate of a distinctive base due to an analysis error at a specific position selected in step e at a set threshold or significance level, using a parametric error model using the reliability score obtained in step d Process,
h) When the ratio obtained in step f is significantly higher than the appearance rate obtained in step g, the threshold value set in step g is a nucleic acid molecule having a specific base type at a specific position selected in step e. Or a step of presuming presence in clonal nucleic acid mixture at a significant level, and i) containing the base species at the specific position obtained in step f, when it is judged that the base species selected in step e is present in step h Calculating the ratio of the specific base species at the specific position selected in step e at the threshold or significance level set in step g from the number of consensus leads, and adopting it as the ratio present in the clonal nucleic acid mixture,
Method including. A method of estimating the base sequence of a specific region of each nucleic acid molecule constituting a clonal nucleic acid mixture, using a computer, comprising:
a) obtaining a lead sequence derived from each nucleic acid molecule,
b) mapping the read sequence obtained in step a onto a reference sequence to obtain a mapping result;
c) obtaining a clustering result of the lead sequence based on the mapping result obtained in step b,
d) obtaining consensus lead information comprising a consensus lead and a reliability score corresponding to each base in the consensus lead based on the clustering result,
e) From the sequence information of each consensus lead contained in the set of consensus reads obtained from each nucleic acid molecule constituting the clonal nucleic acid mixture through steps a to d, the base type present at the position corresponding to the specific position on the reference sequence And the process of selecting their reliability score,
f) obtaining a ratio of the number of consensus reads containing each specific base species at the particular position selected in step e to the total number of consensus leads including the particular position used in step e,
g) The base at the specific position confirmed in step e at a set threshold or significance level by the parametric error model using the reliability score obtained in step d corresponding to the base type at the specific position selected in step e Estimating the rate of occurrence due to analysis errors of the species,
h) When the ratio obtained in step f is significantly higher than the appearance rate obtained in step g, the threshold value set in step g is a nucleic acid molecule having a specific base type at a specific position selected in step e. Or judging that the base species at the specific position is unknown if it is present in the clonal nucleic acid mixture at a significant level or lower than the occurrence rate,
obtained by applying i) the result obtained in step h to each consensus lead of the base to create secondary sequence information, and j) changing the base position in step e and repeating steps e to i. Integrating all secondary sequence information of the specific region and adopting it as the base sequence of the specific region of each nucleic acid molecule constituting the clonal nucleic acid mixture,
Method including. When the reference sequence adopted in step b and the base sequence on the consensus read are the same, it is characterized in that it is adopted as the base species of the secondary sequence information prepared in step i without passing through steps e to h. A method for determining a base sequence according to claim 3. The method according to any one of claims 1 to 3, wherein the specific base position selected in step e is a position different from the reference sequence adopted in step b and the base type on the consensus lead. A method according to any of claims 1 to 3, characterized in that the nucleic acid molecule has a molecular barcode at at least one end. 7. The method of claim 6, wherein each nucleic acid molecule has a different molecular barcode. 8. The method of claim 7, wherein each nucleic acid molecule has one or more different molecular barcodes. The method according to any one of claims 1 to 8, wherein the clustering in step c is performed using the arrangement of molecular barcode regions in the lead sequence as an indicator. The method according to claim 9, wherein the molecular barcode is a random sequence of base sequences. The method according to claim 10, wherein the mismatch in the lead sequence of the molecular barcode region is within the allowable range of the sequence similarity between the lead sequences of the molecular barcode region. The method according to any one of claims 1 to 11, wherein the confidence score is obtained by the Bayesian approach from the information of the clustered lead sequences. The method according to claim 12, wherein the lead sequence is a lead sequence of an amplification product of a nucleic acid molecule. The method according to claim 13, wherein the amplification product is obtained by polymerase chain reaction. The method according to claim 13 or 14, wherein the likelihood function of the Bayesian approach is a function taking into account sequencer errors and amplification errors. The method according to claim 15, characterized in that if it is determined in the mapping step in step b that different bases in the read sequence correspond to insertions or deletions with respect to the reference sequence, a function further corrected for sequencer error is employed. Method. The method according to any of the preceding claims, wherein the parametric error model of step g is based on a Poisson distribution. The method according to claim 1 or 2, wherein a specific base type at a specific position is a target base of low frequency mutation. A method according to any of the preceding claims, which proceeds directly from step c to step d. Next to step c
k) performing local assembly using the read sequences included in each clustering obtained in step c, and obtaining clustering results of the assemble sequences;
The method according to any one of claims 1 to 3, characterized in that A method comprising the step of determining a specific base species present at a specific position in the base sequence of the diffusion molecule constituting the clonal nucleic acid mixture, based on the results obtained by the method according to claim 19 and claim 20. A computer readable program for causing a computer to perform the method according to any of the preceding claims. An apparatus for estimating the presence of a specific base type at a specific position in the base sequence of a nucleic acid molecule constituting a clonal nucleic acid mixture, comprising:
A lead sequence acquisition unit for acquiring a lead sequence derived from each nucleic acid molecule;
A mapping information acquisition unit that supplies the acquired lead array to a mapping device and acquires the mapping result on the reference array by the mapping device;
A clustering result acquisition unit for acquiring a clustering result of the mapped lead sequence;
A consensus lead information acquisition unit comprising a consensus lead and a reliability score corresponding to each base in the consensus lead based on the clustering result;
Estimation device equipped with An apparatus for estimating the presence of a specific base type at a specific position in the base sequence of a nucleic acid molecule constituting a clonal nucleic acid mixture, comprising:
A lead sequence acquisition unit for acquiring a lead sequence derived from each nucleic acid molecule;
A mapping information acquisition unit that supplies the acquired lead array to a mapping device and acquires the mapping result on the reference array by the mapping device;
A first clustering result acquisition unit for acquiring a clustering result of the mapped lead sequence;
A second clustering result acquisition unit that executes a local assembly using a lead sequence included in each clustering and obtains a clustering result of the assembly sequence;
A consensus lead information acquisition unit comprising a consensus lead and a reliability score corresponding to each base in the consensus lead based on the clustering result of the second clustering result acquisition unit;
Estimation device equipped with

Images