CN111326212A - A method for detecting structural variation - Google Patents

Abstract

The invention discloses a method for detecting structural variation, which simultaneously supports the identification of the structural variation of RNA and DNA, and has the advantages of high sensitivity, high specificity, high speed and low resource consumption. The invention also provides a complete system or device, a computer readable storage medium and equipment established based on the method.

Description

A method for detecting structural variation

technical field

The invention belongs to the technical field of gene detection, and in particular relates to a detection method for structural variation and a related system, device, computer-readable storage medium and equipment.

Background technique

Originating from structural variation within the genome, including deletions, inversions, duplications within the same chromosome, and abnormal connections between different chromosomes. No matter what kind of event is caused, the result often shows that the different parts of the two genes are physically linked together. After transcription, a new transcript consisting of parts of the transcripts of two different genes can be obtained at the transcriptional level. Book. These structurally variable genes have important scientific significance in the occurrence and development of cancer, and have extremely important medical value for studying the mechanism of tumor occurrence and development, and for tumor treatment and monitoring. For example, BCR->ABL is widely found in hematological tumors and bladder cancer, lung cancer, malignant glioma and other tumors, FGFR3->TACC3 mainly occurs in bladder cancer, cervical squamous cell carcinoma and cervical adenocarcinoma, EML4->ALK mainly occurs in in lung cancer.

Structural variation detection technology based on next-generation sequencing technology has emerged for a long time, mainly by sequencing the target region or the whole genome, and analyzing the sequence obtained to determine the occurrence of structural variation.

The detection of the DNA level is mainly by aligning the data obtained by sequencing with the genome, and according to whether the aligned reads have a break alignment, that is, the two parts of the read are aligned to different positions in the genome, to collect possible Evidence that supports the occurrence of structural variation. If there is a break alignment, the two parts of the read with the break alignment are further analyzed. According to the position of the alignment of the two parts, the chain has always reversed the cause of the structural variation, and calculated the structure. result of mutation. In the case of paired-end sequencing, two reads will be generated from the same template after sequencing, and abnormal read pairs that support structural variation can be collected according to whether the alignment of the two paired reads is abnormal. Normally, two reads One segment should be aligned to the positive strand of the genome and the other to the negative strand of the genome, and the transcription direction of the nucleic acid should be consistent, and the length of the inserted fragment should be within a reasonable distribution range. If the two reads come from the structure For the two parts of the mutated gene, the direction is abnormal, or the length of the implied insert fragment is abnormal. However, many methods currently released have shortcomings such as long calculation time, low sensitivity, high false positive, and no annotation module.

The detection of RNA level often requires the alignment of the two reference sequences, the genome and the transcriptome, through the transcriptome alignment, and then the coordinates are mapped to the genome, and then through the alignment features of the read segments, referring to the judgment method of DNA calculation, deduce The mechanism of the occurrence of structural variation events, and to calculate the types and genes of structural variation. This technical solution not only consumes too much resources, but also suffers from intron interference, which often leads to inaccurate calculations, and there are many false positives and false negatives. Calculation results It is also often lack of annotations, which is extremely inconvenient to use.

SUMMARY OF THE INVENTION

One of the objectives of the present invention is to provide a structural variation detection method that supports both RNA and DNA structural variation identification, with high sensitivity, high specificity, high speed and low resource consumption, aiming at the defects of the prior art.

In order to achieve the above object, the present invention provides a detection method for structural variation, the detection method comprises the following steps:

1) aligning the sequencing data to a reference genome sequence or a reference master transcript sequence;

2) Look for the normal aligned reads (reads), the reads with fragmented alignments, and the inconsistent aligned reads;

3) Categorize the read segment in which the fragmentation alignment occurred and the inconsistently aligned read pair;

4) Grouping different types of fragmented alignment reads and inconsistent alignment read pairs respectively, and grouping reads supporting the same structural variation event into the same set;

5) For the structural variation event determined by the fragmented alignment reads, a conserved sequence is formed by assembling the reads supporting the structural variation event;

6) Determine precise breakpoint positions based on conserved sequences;

7) Merge the structural variants supported by the break alignment reads and the structural variants supported by the inconsistent alignment reads respectively, and the merging refers to merging structural variant events with similar breakpoints and the same type into the same structural variant event ;

8) Merging the structural variants supported by the break alignment reads with similar breakpoints and the same type with the structural variants supported by the discordant alignment read pairs;

9) Deletion of conserved sequences that can be completely and continuously matched to a sequence on the genome or structural variation events that can generate multiple consistent alignments;

10) Calculate the frequency of structural variation events.

In a specific embodiment, in step 1), if the sequencing data is DNA data, it is aligned to the reference genome sequence; if the sequencing data is RNA data, it is aligned to the reference master transcript sequence.

In a specific embodiment, step 2) further comprises the step of calculating the length of the insert by using the normal alignment reads, and calculating the main parameter of the length distribution of the insert; the main parameter is preferably a maximum value, a minimum value and/or a mean.

In a specific embodiment, in step 3), the classification of the reads with fragmentation alignment can be performed based on the following indicators:

Whether the two parts of the split alignment are aligned to the same chromosome, aligned to different orientations of the genome, and/or whether the splicing position is upstream of the aligned position.

In a specific embodiment, in step 3), the reads with fragmentation alignment can be classified according to the following criteria:

In a specific embodiment, in step 3), the classification of discordantly aligned read pairs can be performed based on the following indicators:

Whether to align to the same chromosome, to different orientations of the genome and/or insert size.

In a specific embodiment, in step 3), the discordantly aligned read pairs can be classified according to the following criteria:

In a specific embodiment, SA tags can be used to find reads that have split alignments and/or discordantly aligned read pairs.

In a specific embodiment, when identifying reads that are split aligned, parts that align elsewhere are not counted.

In a specific embodiment, when identifying reads of a break alignment, a break alignment record from the same aligned read can be considered as one entity.

In a specific embodiment, in step 4), the grouping can be performed by cluster analysis.

Preferably, the criteria for clustering the break-aligned reads include the type of structural variation, the reference sequence name for the two-part alignment of the break-aligned reads, the first breakpoint position and/or the second breakpoint position.

Further, reads with the same type of structural variation, the same name of the alignment reference sequence, and the first and second breakpoint positions are within m bases apart as a class, where m is a natural number within 30, preferably 10.

Further, if the read support number in a class is higher than a preset threshold (the threshold is selected from a natural number above 1), then the average breakpoint position is obtained by using the average value of all the read breakpoint position information in the class. .

Preferably, the criteria for clustering discordantly aligned read pairs include the type of structural variation, the reference sequence name for the read alignment, the reference sequence name for the paired read alignment, the read alignment position and/or the paired read comparison location.

Further, the reads with the same name of the alignment reference sequence and the alignment position difference within the range of the maximum insert size are regarded as one type.

Further, if the read support number in a class is higher than a preset threshold (the threshold is selected from a natural number above 1, preferably 2), then the breakpoint range determined by the read pairs in the class is used to determine the breakpoint range. Estimating the breakpoint position; preferably, according to the alignment start position and end position of the read pair, a progressive method is used to continuously narrow the interval where the breakpoint position is located to estimate the breakpoint position.

In a specific embodiment, in step 5), the assembly can be performed by multiple sequence alignment.

In a specific embodiment, step 5) further includes the step of parsing and reconstructing the sequence involved in the assembly, the parsing and reconstructing comprising: extracting a short insertion sequence near the breakpoint and/or adjusting the sequence direction The 5' end of the read is in the same direction as the reference sequence.

In a specific embodiment, in step 6), the conserved sequence is aligned with the reference breakpoint sequence, and the exact breakpoint position is determined according to the position where the breakpoint alignment occurs; the reference breakpoint sequence comprises two parts, one part The other portion is derived from the reference sequence spanning one breakpoint of the structural variation event, and the other portion is derived from the reference sequence spanning the other breakpoint of the structural variation event.

In a specific embodiment, in step 7), structural variation events with higher support numbers are retained when merging.

In a specific embodiment, in step 7), for the structural variation events supported by the break alignment reads, the conditions for merging include that the breakpoint distance is within the length of any one of the two homologous sequences or less than n bases , and is the same type of structural variation, where n is a natural number within 30, preferably 10.

In a specific embodiment, in step 7), for the structural variation events supported by the discordantly aligned read pairs, the conditions for merging include that the breakpoint distance is within the range of the maximum insertion sequence length, and is the same type of structural variation.

In a specific embodiment, in step 8), the conditions for merging include a breakpoint of a structural variation event supported by a break alignment read and a breakpoint of a structural variation event supported by another discordant alignment read of the same type. The point distance is within the maximum insertion sequence length.

In a specific embodiment, the alignment can be performed using the BWT algorithm in step 9).

In a specific embodiment, in step 9), if the conserved sequence can be completely aligned to the reference genome sequence or a continuous region on the reference master transcript, the set of records is recorded and/or filtered out; if the conserved sequence The break alignment can occur on the reference genome sequence or the reference main transcript sequence, and the breakpoint position determined by the break alignment is within 10bp of the previously calculated breakpoint position, and there are more than 2 groups of such break alignments. This group of records is recorded and/or filtered; if such a break alignment is 1 group, the breakpoint is further refined according to the breakpoint position; if such a break alignment is less than 1 group, no marking and/or modification is performed.

In a specific embodiment, in step 10), the structural variation event frequency is calculated by the reference type count and variant type count of two breakpoints of the structural variation event, and the reference type is the same as the reference genome sequence or the reference master transcript. The type of sequence consistency, the variant type is the type consistent with the sequence of the structural variation event; preferably, the reference type with higher support numbers on both sides of the breakpoint is used as the reference type count of the structural variation event.

Further, the frequency of the structural variation event=the number of molecules supporting the structural variation/(the number of molecules supporting the structural variation+the number of molecules supporting the reference type).

Preferably, if two reads of the same template support a certain structural variation event, it is counted as one molecule.

The detection method of the present invention may further comprise the step of annotating the fusion gene and/or outputting the detection result; preferably, the annotated information items include but are not limited to the fusion gene name, fusion gene breakpoint position, exon, intron position, Whether the connection is correct, the sequence of structural variants and/or the number of structural variants supported; preferably, the result output is in tab-separated mode or binary compression mode.

Structural variation detection software in the prior art usually does not have a fusion gene annotation module, and often requires third-party annotation software. The third-party annotation software often cannot fully utilize the data of structural variation calculation software, and the annotation is often insufficient or even fails. In order to solve this problem, the software developed based on the method of the present invention can have a built-in annotation module, after calculating the structural variation events, the structural variation genes can be directly annotated.

In order to achieve the purpose of personalized output, various standards can be set for the scope of the structural variation report, such as the limit of the number of supports, the limit of the reported range, the limit of the direction of structural variation, etc. The report format can be divided into two categories: one The class is a tab-separated pattern, which is easy to read by humans, and the class is a binary compression format, which is easy for computer processing. You can also set up a companion tool module to further mine the results, and perform operations such as filtering based on existing results, regenerating reports, and merging support background pools.

The output information can be selected as needed, including but not limited to: fusion genes related to structural variation, number of molecules supporting structural variation, number of molecules supporting reference type, frequency of structural variation, whether structural variation occurs in COSMIC or ONCOKB, fusion genes Status, distance between two genes involved in fusion, chromosomal information, breakpoint location, type of structural variation, size of structural variation, insertion sequence near breakpoint, fusion mask and/or conserved sequence near breakpoint, etc.

The detection method of the present invention may further comprise the step of sequencing and/or obtaining sequencing data.

The technical solutions of the present invention can be used in various diagnostic and non-diagnostic application scenarios of cancer. The technical solution of the present invention can be applied to tumors of any stage, such as very early stage tumors, early stage tumors, mid-stage tumors, and late stage tumors; it is preferably used for early stage tumors or very early stage tumors.

Another object of the present invention is to provide a system or device for detecting structural variation, the system or device comprising:

1) a sequencing module and/or a sequencing data acquisition module; and

2) An identification module for performing the structural variation detection method of the present invention with respect to the sequencing data.

In a specific embodiment, the system or device further comprises:

3) An annotation module for annotating fusion genes; and/or

4) The output module is used to output the detection result.

The present invention also provides a computer-readable storage medium, the computer-readable storage medium comprising a stored computer program, the computer program comprising a program for executing the structural variation detection method of the present invention.

The present invention also provides an apparatus comprising a processor, a memory, and a computer program stored in the memory, the computer program comprising a program for performing the structural variation detection method of the present invention.

The beneficial effects of the present invention include at least the following aspects:

(1) The detection method of the present invention supports the identification of structural variation of RNA and DNA at the same time, with extremely high sensitivity, high specificity, high speed and low resource consumption. As long as there is one read/read pair that supports variation, the detection can be realized, while the existing technology usually requires at least 3 or 3 pairs. Such sensitivity can achieve a qualitative leap, because in practical sequencing applications, especially in the sequencing of tumor samples with low frequency and high heterogeneity, it is often difficult to obtain enough (at least 3) long enough samples that just span the breakpoint. Reads meet the requirements of prior art methods.

(2) For the calculation of structural variation based on RNA sequencing, the present invention adopts a novel primary transcript alignment calculation method, which completely avoids the influence of intron alignment and calculation, and achieves excellent results.

(3) When using the method of the present invention, when classifying the read segment of the fracture alignment and the read segment with inconsistent alignment, the classification can be completed by only visiting the same read segment once, and does not depend on other comparisons of the read segment. Records, and the requirements for comparison records are more stringent. If the fractured comparison part is compared to other places, it will be regarded as a repeating area and will not be counted.

(4) When classifying the read segments, the prior art method often sorts according to the breaking position, which is easy to cause false negatives and false positives, and the method of the present invention can effectively avoid this problem.

(5) When adopting the method of the present invention, when grouping different types of break alignment reads and inconsistent alignment read pairs, the use of complex calculation models is avoided, and the calculation is simplified, and the method of the prior art determines breakpoints complex nodes, trees and clustering will be constructed, which is not only complex in operation, but also prone to false positives.

(6) Using the method of the present invention to assemble the sequence can effectively improve the accuracy of breakpoint judgment.

(7) The present invention utilizes the SA tag to search for the read segment and/or the inconsistently aligned read segment pair, which can complete the classification by only accessing the same read segment once, while in the prior art, it is necessary to obtain a read segment at the same time The calculation can be carried out only if the two comparison records are obtained, and the method of the present invention not only saves the calculation time, but also saves the storage space of the comparison file.

(8) The method of the present invention has stricter requirements on the recording, and when identifying the read segment of the fracture alignment, if part of the alignment is in other places, it is regarded as a repeating area and will not be counted.

(9) When identifying the reads of the break alignment, the break alignment records from the same aligned read are considered as one entity. However, in the prior art, the sequence is often based on the position of breakage, and does not distinguish whether the source of the same read is not. The methods of the present invention help to avoid false negatives and false positives.

Description of drawings

FIG. 1 is a flow chart of the structural variation detection method of the present invention.

Detailed ways

Unless otherwise specified, the terms used in the present invention have the usual meanings in the art, and the reagents used are all conventional commercial reagents in the art.

The term "reference sequence" in the present invention refers to a sequence derived from a DNA sequence file of a reference genome or an RNA sequence file of a transcriptome.

In the present invention, the term "normally aligned read (read)" refers to a continuous region on the genome that can be continuously aligned in the process of aligning a read with the genome. If the read length is m, then the read Able to align to a continuous closed region [a,b] on the genome, allowing mismatches and short indels, but no splicing at the beginning and end.

In the present invention, the term "a read segment with a split alignment" refers to the process of aligning a read segment with the genome, and two parts of the read segment are aligned to a discontinuous region of the genome. For example, the read length is m, 10<n<m ,Read [0,n] is normally aligned to genome [a,b], read [n+1,m] is normally aligned to genome [c,d],[a,b],[c,d] The two intervals either belong to different chromosomes or the distance is greater than 100.

The term "inconsistently aligned read pair" in the present invention refers to a pair of reads obtained by paired-end sequencing of a template. When alignment occurs, the alignment results do not support the same normal template, and the normal template refers to the upper part of the genome. A continuous region [a,b] of .

Regarding the term "insertion sequence" in the present invention, the insertion sequence has different meanings in different contexts. When calculating the template length distribution, the insertion sequence of a read refers to the length of the template, that is, the minimum length of a pair of read alignments. The absolute value of the difference between the starting position and the maximum ending position; when calculating the structural variation, it means that there is an extra sequence at a certain position in the sample relative to the reference genome.

In the present invention, the terms "largest insert fragment" and "largest insert sequence" refer to the insert sequence (template) with the longest length.

In the present invention, the term "support" means a causal relationship, A supports B, and a phenomenon indicates that B has occurred.

The term "structural variation event" in the present invention refers to a large change in the genome of a sample relative to a normal reference genome sequence, specifically including inversion, insertion, deletion, duplication, and transposition.

The term "conserved sequence" in the present invention refers to a collection of a series of sequences, through multiple sequence alignment, in the alignment result, the base with the support number greater than three sequences and the highest support number is extracted as the conserved base at a certain position , and connect all conserved bases according to their relative positions in the original sequence to obtain a conserved sequence.

In the present invention, the term "main transcript" refers to the artificially selected, among several transcripts of the same gene, the longest, the most thoroughly studied and the most clinically valuable transcript.

The term "breakpoint position" in the present invention refers to the start and end position of a structural variation event, and the position is a position on the genome or transcriptome.

The term "reference type" in the present invention refers to the same genotype as the reference genome or reference transcript.

The term "variant type" in the present invention refers to the same variant type as a structural variant event.

The term "molecular counting" in the present invention refers to counting according to the number of template strips.

The term "template" in the present invention refers to the DNA molecule of the original sample used for sequencing.

In the present invention, the terms "smaller alignment position" and "larger alignment position" are a pair of relative concepts, and are a measure of the alignment position of a read on the reference genome sequence or reference master transcript sequence, close to the 5' The position is smaller, and the position close to 3' is larger. The terms "smaller coordinates" and "larger coordinates" in the present invention mean that the coordinates on the same chromosome or transcript are smaller at the position close to 5', and larger at the position close to 3'.

In the present invention, the terms "upstream splicing" and "downstream splicing" refer to when a sequence is aligned with the reference genome sequence or the reference master transcript sequence, a part of which is completely aligned, and the other part is Shear out, upstream shearing if shearing occurs upstream of the fully aligned section, and downstream shearing if shearing occurs downstream of the fully aligned section.

In the present invention, the terms "large chromosome" and "small chromosome" refer to the fact that chromosomes are mapped into natural numbers in the BAM file, and the larger chromosomes are mapped to the larger natural numbers, and the smaller chromosomes are small chromosomes.

The technical content of the present invention will be described in detail below with reference to the accompanying drawings and specific embodiments. Those skilled in the art will understand that the following examples are only used to illustrate the present invention and should not be construed as limiting the scope of the present invention.

Example 1

15,000 samples from patients with solid tumors and hematological tumors were collected, DNA was extracted for paired-end second-generation capture and library construction, and RNA was extracted from some samples for capture and sequencing. For each sample library obtained by sequencing, for DNA library, use BWA alignment software to align with the human reference genome HG19; for RNA library, use BWA alignment software to align with the personalized constructed transcriptome. The transcriptome is constructed by:

Select one transcript for each gene as the master transcript. The master transcript must be preferentially included in the COSMIC and ONCOKIB databases. If it is not included, it must be the one recorded in the master transcript database in UCSC. If the main transcript is determined, the longest one is used as the main transcript.

The results obtained by the comparison are removed from PCR duplicates and corrected for errors, and a BAM file for structural variation detection is generated. The BAM file is used as the input file of the method of the present invention. The minimum number of seeds for inconsistent comparison reads is 2, the minimum number of seeds for broken comparison reads is 1, the running memory is set to 8G, the thread is set to 8, and the cluster is delivered for analysis.

Example 2

Traverse the alignment result file, take 1 million normal aligned reads (the extraction standard is that no fragmentation alignment occurs, and the read pair conforms to normal reads, and does not support any structural variation) to count the length of the inserted fragment and calculate the length of the inserted fragment The main parameters of the distribution, maximum, minimum, mean.

Traverse the alignment result file again, access all the reads, find out all the reads with broken alignment and the read pairs with abnormal alignment, and analyze the reads with broken alignment and the abnormal alignment according to the following characteristics. Classification of read pairs: The reads of the break alignment are mainly classified according to the alignment of the two parts of the break alignment. The classification criteria include whether they are aligned to the same chromosome, whether they are aligned to different directions of the genome, and whether they are cut or not. The cut positions are all upstream of the alignment record, and the detailed classification table is shown in Table 1; the classification criteria of the reads with abnormal alignment mainly include whether the alignment is on the same chromosome, whether the alignment is in different directions of the genome, and the length of the inserted fragment. Size, detailed classification table see Table 2.

Table 1 Breakdown alignment read classification table

Table 2 Discordant alignment read pair classification table

Clustering analysis is performed on different types of break alignment reads and inconsistent alignment reads. The purpose of clustering is to put reads that support the same structural variation event into a set to obtain more accurate breakpoints. Minimize random errors in breakpoints caused by alignment and sequencing errors.

The method of clustering the reads of the break alignment is as follows: Rank the reads of the break alignment, and the ranking criteria include "structural variant type, the reference sequence name of the two-part alignment of the break alignment reads, and the first break point. position, the second breakpoint position"; the reads with the same type of structural variation, the same name of the alignment reference sequence, and both break positions within 10 bases are regarded as a class. If the read support number in a class is higher than a preset threshold (specifically set to 1 in this example, considering that probe design and sequencing interruption may not be able to perfectly detect many reads spanning the breakpoint, the The minimum threshold can be selected as 1, which can achieve extremely high sensitivity), then the average breakpoint position is obtained by taking the average value of the breakpoint position information of all reads in this class, and a structural variation event supported by the broken read is generated, which is the same as the average value. The breakpoint maximum error and minimum error are used as the breakpoint error range.

The clustering method for inconsistently aligned read pairs is as follows: Rank the inconsistently aligned reads, and the ranking criteria include "structural variation type, read alignment reference sequence name, paired read alignment reference sequence name, read Alignment position, the position of the paired read alignment"; the read pairs with the same alignment reference sequence name and the alignment position difference within the range of the maximum insert size are regarded as a class; if the read pairs in a class support The number is higher than the preset threshold (specifically set to 2 in this example, the value can be adjusted according to the actual sequencing quality, the minimum is 1, which can achieve a great sensitivity, it is recommended to set it to 2, because the reads are inconsistently compared For the total number, it is often impossible to obtain a higher value than the initial evidence collection stage, and if the setting is too low, it will result in a large amount of computation), then use the breakpoint range determined by the read pairs in this class to estimate the breakpoint position, and generate a Structural variation events supported by discordantly aligned read pairs.

For the structural variation event determined by the fragmented alignment reads, a conserved sequence is obtained by assembling the reads supporting the structural variation event, and the assembly is performed by multiple sequence alignment technology. The conserved sequence can reduce the small size caused by sequencing errors. The effect of indels or single nucleotide mutations on breakpoint judgment. In the process of sequence assembly, the sequences involved in the assembly were analyzed and reconstructed, mainly including: extracting short insertion series near the breakpoint, the insertion sequence will affect the alignment with the reference sequence and affect the accurate judgment of the breakpoint; Adjust the direction to ensure that the 5' end of the read is in the same direction as the reference sequence as much as possible, which is convenient for unified processing and breakpoint calculation in the later stage.

After the conserved sequence is obtained, it is compared with the reference sequence, and the exact breakpoint position is determined according to the position where the break alignment occurs (if no suitable break alignment is generated, the structural variation event is not a reliable event and is directly discarded. Subsequent analyses were not performed, which were often false positives due to sequencing errors or excessive alignment errors). After obtaining the breakpoint, take a short sequence spanning 10 bases upstream and downstream of the breakpoint as a probe (when calculating the frequency of structural variation events later, the original BAM file will be traversed again, and more spans will be extracted from the original BAM file. Reads with breakpoints. These reads may be seed reads that have undergone break alignment and have been used to calculate structural variation events, or may be short breaks and failed to generate BAM in the previous alignment. By performing a fragmentation alignment, by aligning with the probe, it is possible to know whether these fragmented reads support the corresponding structural variation events), and then have a more accurate calculation of the frequency. For cases involving inserts, two copies of the probe are kept, one with the insert and one without the insert.

The reference sequence consists of two parts, one from the reference sequence spanning one breakpoint of the structural variation event, and the other from the reference sequence spanning the other breakpoint of the structural variation event, and the manner in which the two are connected (including which one is in the 5' end, which is at the 3' end, whether the forward sequence participates in the ligation or the reverse complement participates in the ligation) is determined by the alignment of the corresponding fragmented reads that support the structural variation event. The details are as follows, the reference sequence of the type I inversion class structural variation event, the upstream part comes from the reference sequence spanning the breakpoint with smaller coordinates, the downstream part comes from the reference sequence spanning the breakpoint with larger coordinates, the upstream is positively involved in the link, the downstream part is Reverse complement is involved in ligation; the reference sequence for type II inversion class structural variation events, the upstream part comes from the reference sequence spanning the smaller breakpoint, the downstream part comes from the reference sequence spanning the larger breakpoint, and the upstream is reverse complementary. The reference sequence of the deletion-type structural variation event, the upstream part comes from the reference sequence spanning the smaller coordinate breakpoint, and the downstream part comes from the reference sequence spanning the larger coordinate breakpoint, both upstream and downstream are both. The forward sequence is involved in the connection; the reference sequence of the replicative structural variation event, the upstream part comes from the reference sequence that spans the larger coordinate, and the downstream part comes from the reference sequence that spans the smaller coordinate breakpoint, and both upstream and downstream sides participate in the connection with the forward sequence ; the reference sequence of the insertion structural variation event is directly taken from the positive strand sequence of a reference sequence spanning the insertion site; the structural variation event that connects the 5' and 5' ends of the two chromosomes, the upstream part comes from the larger spanning chromosome number The reference sequence of the breakpoint, the downstream part comes from the reference sequence of the breakpoint spanning the smaller chromosome number, the upstream forward is involved in the link, and the downstream reverse complement is involved in the link; the structural variation of the 3' and 3' ends of the two chromosomes is connected Events, the upstream part comes from the reference sequence spanning the breakpoint with the larger chromosome number, the downstream part comes from the reference sequence spanning the breakpoint with the smaller chromosome number, the upstream is involved in the connection in a reverse complementary manner, and the downstream is involved in the forward connection; Structural variation events connecting the 5' end of the large chromosome to the 3' end of the small chromosome, the upstream part comes from the reference sequence spanning the breakpoint with the larger chromosome number, and the downstream part comes from the reference sequence spanning the breakpoint with the smaller chromosome number , the upstream and downstream are connected in a forward direction; the structural variation event connecting the 3' end of the large chromosome and the 5' end of the small chromosome, the upstream part comes from the reference sequence spanning the larger breakpoint of the chromosome number, and the downstream part comes from the spanning chromosome The reference sequences of the smaller breakpoints are connected upstream and downstream in a forward direction.

After the refinement of the breakpoint locations, there are some structural variation events with similar breakpoints and the same type, which themselves are events from the same structural variation, but the sequencing error or alignment error caused two events to occur. At this time, only You need to merge and retain the higher number of supports. Structural variants supported by the split reads were merged if the breakpoint distance was within the length of either homologous sequence or less than 10 bases. For the structural variants supported by the inconsistent reads in the alignment, if the breakpoint distance is within the range of the maximum insertion sequence length, and the structural variants of the same type are merged, the higher supporters are retained.

Structural variation events determined by fragmentation alignment reads are often supported by inconsistent alignment reads. In the previous calculation of this example, these two events were separated. To strengthen the evidence, the two events were merged at this time. The criterion for merging is that if the breakpoint of a structural variation event supported by a fragmented read is just within the maximum insertion sequence length from the breakpoint of a structural variation event supported by another discordant read of the same type, the merge will be performed. Structural variant events supported by discordantly aligned read pairs that could not be merged were retained directly.

For the structural variation events supported by the fragmented alignment reads, some false positives come from the reads in the repetitive region. These reads cannot be identified from a certain position in the genome, and their conserved sequences can often be aligned. To multiple locations in the genome, and the early determination of breakpoints was performed in a small range of alignments, the accuracy was limited to a certain extent. In order to further refine the location of the breakpoints, a global alignment with the reference genome was performed to mark the repetitive regions. structural variation events.

The alignment is performed using the BWT algorithm. If the conserved sequence can be completely and continuously matched to a sequence on the genome, or multiple consistent alignments can occur, the structural variation event is not a true structural variation event. The operation scheme is as follows: if the conserved sequence can be completely and contiguously matched on the reference genome sequence or the reference master transcript sequence, the shuffling record is recorded and filtered; if the conserved sequence can occur on the reference genome sequence or the reference master transcript sequence Break alignment, and the breakpoint position determined by break alignment is within 10 bp from the previously calculated breakpoint position, and such break alignment exceeds 2 groups, the record of this group is recorded and filtered; if such break If the alignment is 1 group, the breakpoint is further refined according to the position of the breakpoint; if such a break alignment is less than 1 group, no marking and modification are performed.

After obtaining the structural variation events with further refined breakpoints, the frequency of structural variation events is calculated. To calculate the frequency of a structural variation event, two related indicators need to be calculated. One is the reference count of the two breakpoints of the structural variation event. , and the other is the variant count of structural variant events. In order to avoid the high count caused by using the number of reads as the variant count in the traditional practice, in this embodiment, the number of molecules is counted. Molecular counts are also used for the counts of the reference type, and the reference type with higher support numbers on both sides of the breakpoint is used as the reference type count for the structural variation event.

Whether the reads and read pairs spanning the breakpoint support the reference type is determined according to the classification of reads of different structural variants (see Table 1 and Table 2), and any variants that do not support are counted as the reference type.

The structural variant frequency is equal to the number of molecules supporting the structural variant divided by the sum of the number of molecules supporting the structural variant and the number of molecules supporting the reference type.

Example 3

The structural variation detection method of the present invention can support the identification of RNA and DNA structural variation at the same time, while the prior art methods seldom support both identifications at the same time. For example, both DELLY and GeneFuse only support the identification of DNA structural variation, and TophatFusion only supports Supports RNA structural variant identification. To test the performance of the method of the present invention, it was compared with FusionMap, a commercial software capable of supporting both DNA and RNA structural variant identification.

269 positive samples were extracted from the sequencing data of Example 1, and each sample had 1-2 fusions that were experimentally verified to exist. The method of DNA library building and sequencing analysis and FusionMap software were used to detect structural variation. The results showed that the method of Example 2 could achieve a 100% detection rate (see Table 3), and FusionMap software could also be used to detect 100%. , indicating that the software can detect positive fusions well, and the number of fusions detected in each sample is very small, and the specificity is highly guaranteed. Although FusionMap can also detect 100% of the fusions of these samples, it will report For other unreliable fusions (more than 3 fusions are detected on average per sample in FusionMap), it can be seen that the method of the present invention has higher specificity than FusionMap.

Table 3 Positive sample analysis results of the detection method of the present invention

Number of test samples 269 The detection rate 100% Report the number of fusions/sample Average of less than 2 fusions/sample

A total of 18 samples with poor quality raw data and samples with more complex libraries were used for further testing. All 18 difficult samples were identified by the method of Example 2, but were not identified in the FusionMap software (see Table 4). . It can be seen that compared with the commercial paid software FusionMap, the method of the present invention can better detect samples with poor original data quality or complex libraries.

Table 4 Analysis results of complex or low frequency positive samples

Number of test samples 18 pcs The number of samples detected by the method of embodiment 2 18 (100% detected) Number of samples detected by FusionMap 0

The test scale was further expanded, 8000 samples were extracted from the sequencing data of Example 1 for testing, and FusionMap was used for detection first, and a total of 98 credible fusion samples were missed. The samples missed by FusionMap were detected by the method of Example 2, and all 98 missed samples were detected (Table 5).

Table 5 Analysis results of clinical samples

Number of test samples 8000 FusionMap missed detection and trusted fusion 98 Detection of missed samples by the method of embodiment 2 98

In addition, the reason why the method of the present invention can obtain high-sensitivity detection results can be at least partially based on the principle of the method of the present invention, when searching for the initial evidence of possible structural variation events, only one fragmentation alignment read or only one fragment can be achieved. Calculations in the case of the support of inconsistently aligned read pairs, while the FusionMap software as a control requires at least 2 break alignment reads to identify structural variants, and another commonly used open source software, DELLY, requires at least 3 break ratios. Structural variation calculation can only be started from the reads or three pairs of discordantly aligned read pairs, so the detection method of the present invention has higher sensitivity than traditional methods.

The aforementioned test items have partially demonstrated the high specificity of the method of the present invention. Next, the scale of the test is further expanded to analyze all 15,000 samples in Example 1. Using the method of Example 2, 1-2 samples were obtained for each sample. However, there are a large number of more than 5 structural variants in the test results of the commercial software FusionMap, and DELLY has several test results with dozens or hundreds of structural variants. Generally speaking, 1-2 structural variants are credible, and the results of additional structural variants are often unreliable or even clinically meaningful. This result indicates that the method of the present invention has higher specificity.

In order to test the speed and resource consumption of the method of the present invention, the method of Example 2 was used to test 15,000 samples in Example 1. The average analysis time is 3 minutes/sample (see Table 6), the average memory usage is 6G/sample, and the thread is 8/sample. General open source software or commercial software takes at least half an hour or longer, and the resource consumption is also more. It can be seen that compared with the prior art, the method of the present invention has higher speed and less resource consumption.

Table 6 Resource consumption analysis

testing sample 15000 samples Average memory consumption 6G/sample Threads 8/sample Average analysis time 3 minutes/sample

Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, but not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand: It is still possible to modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements to some or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the technical solutions of the embodiments of the present invention. range.

Claims (10)

1. A method for detecting structural variation, the method comprising:

1) aligning the sequencing data to a reference genomic sequence or a reference master transcript sequence;

2) searching a normal comparison read (read), a read subjected to fracture comparison and an inconsistent comparison read pair;

3) classifying the read with fracture comparison and the read with inconsistency comparison;

4) grouping the fracture comparison reads of different types and the inconsistent comparison reads respectively, and classifying the reads supporting the same structural variation event into the same set;

5) for structural variation events determined by the fragmentation alignment reads, conserved sequences are formed by assembling the reads that support the structural variation events;

6) determining an exact breakpoint location based on the conserved sequence;

7) merging the structural variation supported by the fracture comparison read and the structural variation supported by the inconsistent comparison read, wherein merging refers to merging the structural variation events with similar breakpoints and the same type into the same structural variation event;

8) merging the structural variation supported by the fracture comparison read with similar breakpoints and the same type with the structural variation supported by the inconsistent comparison read;

9) the deletion of the conserved sequence can be completely and continuously matched with a segment of sequence on the genome or can generate a plurality of structural variation events which are aligned in a consistent way;

10) calculating the frequency of structural variation events.

2. The detection method according to claim 1, wherein the step 2) further comprises the steps of counting the lengths of the inserts by the normal alignment reads and calculating the main parameters of the distribution of the lengths of the inserts; the main parameter is preferably a maximum, a minimum and/or a mean.

3. The detection method according to any one of claims 1-2, wherein the SA tags are used to find the reads that are aligned for fragmentation and/or the read pairs that are aligned for inconsistency.

4. The method according to any one of claims 1 to 3, wherein in step 8), the merging condition comprises that the breakpoint of one split aligned read supports a structural variant event is within the maximum insertion sequence length of the breakpoint of another non-aligned read of the same type to the supported structural variant event.

5. The detection method according to any one of claims 1 to 4, wherein in step 9), if the conserved sequence can be completely and continuously matched with a segment of sequence on the genome, the record of the group is recorded and/or filtered; if the conserved sequence can be subjected to fragmentation alignment on a reference genome sequence or a reference major transcript sequence, the breakpoint position determined by the fragmentation alignment is within 10bp of the breakpoint position calculated before, and the fragmentation alignment exceeds 2 groups, recording and/or filtering the group of records; if the fracture ratio is 1 group, further accurate breakpoint is carried out according to the breakpoint position; if such fragmentation ratio is less than 1 group, no labeling and/or modification is performed.

6. The method according to any one of claims 1 to 5, wherein the structural variant event frequency is calculated in step 10) by counting the number of reference types and the number of variant types at two breakpoints of the structural variant event, wherein the reference type is of the type corresponding to the reference genomic sequence or the reference master transcript sequence, and the variant type is of the type corresponding to the structural variant event sequence; preferably, the reference pattern with higher support number on both sides of the breakpoint is used as the reference pattern count of the structural variation event.

7. The detection method according to claim 6, wherein the structural variation event frequency is the number of molecules supporting structural variation/(number of molecules supporting structural variation + number of molecules supporting reference type).

8. A system or apparatus for detecting structural variations, the system or apparatus comprising:

1) a sequencing module and/or a sequencing data acquisition module; and

2) an identification module for performing a method of detecting a structural variation according to any one of claims 1-7 on sequencing data.

9. A computer-readable storage medium, characterized in that the computer-readable storage medium comprises a stored computer program containing a program for executing the method for detecting a structural variation according to any one of claims 1 to 7.

10. An apparatus comprising a processor, a memory, and a computer program stored in the memory, the computer program comprising a program for performing the method of detecting a structural variation according to any one of claims 1-7.

Images