Sequence data-based sample analysis method and device and storage medium

Abstract

The invention discloses a method and equipment for analyzing samples based on sequence data and a storage medium. According to the invention, density-based spatial clustering analysis is carried out on a data set obtained based on sequence data from a sample, the real condition of the sample is judged according to a clustering cluster obtained through analysis, and / or the non-real condition or abnormal condition of the sample is judged according to discrete data relative to the clustering cluster obtained through analysis. Wherein the minimum number k of the data points contained in the clustering cluster is obtained through calculation instead of being directly set or preset according to experience. According to the method, factors of subjective design clustering distance and clustering quantity are eliminated, and the sensitivity and specificity of a detection result are obviously improved.

Description

Technical Field

[0001] The present invention relates to the field of sequence data analysis, and in particular to a method, device and storage medium for analyzing samples based on sequence data. Background Art

[0002] Based on the second-generation sequencing (NGS), it is possible to detect fusion genes at the DNA or RNA level. Due to its short length, it is usually necessary to splice the sequencing data and then compare it, looking for reads that span the fusion breakpoints to identify the fusion event. In addition, the third-generation / fourth-generation long-read sequencing technology does not need to break the RNA, and can obtain longer RNA sequence information. For reads with fusion breakpoints, more intuitive and clear characterization results can be obtained. The sequencing data can be directly compared to the reference genome to obtain the aligned BAM file to find reads containing SA tags to directly obtain the full length of the fusion gene. The second-generation sequencing data is spliced ​​and then compared, and the reads that span the fusion breakpoints are found to identify the fusion event. If the fusion position is in the genomic repetitive sequence, it will cause multiple comparisons and make the detection results inaccurate.

[0003] Software for fusion gene detection on sequencing data, such as LongGF, Genion, JAFFAL, FusionSeeker, etc., all detect a large number of false positive results after testing, and the detection sensitivity of Genion and FusionSeeker is also poor. At present, it is urgent to develop a new fusion detection method that takes into account both detection sensitivity and accuracy to better support the needs of fusion gene detection.

[0004] The information in the background technology is only for illustrating the general background of the present invention and should not be regarded as admitting or suggesting in any form that such information constitutes the prior art known to a person skilled in the art. Summary of the invention

[0005] In order to solve at least some of the technical problems in the prior art, the present invention provides a method for judging the true or false situation of a sample based on density-based spatial clustering analysis. In an exemplary embodiment, the method of the present invention can greatly reduce the false positive of fusion gene detection during sequencing data analysis, while maintaining a high detection rate. Specifically, the present invention includes the following contents.

[0006] In a first aspect, the present invention provides a method for analyzing a sample based on sequence data, comprising:

[0007] A density-based spatial clustering analysis is performed on a data set obtained based on sequence data from the sample, and the true situation of the sample is judged based on the clusters obtained from the analysis, and / or the non-true situation or abnormal situation of the sample is judged based on discrete data relative to the clusters obtained from the analysis, wherein the minimum number k of data points contained in the clusters is obtained by calculation.

[0008] In certain embodiments, according to the method for analyzing samples based on sequence data of the present invention, the data set is obtained by weighted reconstruction based on depth;

[0009] And / or, the cluster-dependent distance parameter e is obtained by calculation,

[0010] Preferably, the cluster-dependent distance parameter e is the median, mean or percentile of the neighbor distance between any two data points in the data set, more preferably, the median, and more preferably, the neighbor distance is the Euclidean distance;

[0011] And / or, the minimum number k is the median, mean or percentile of the number of supports for each data point in the data set, preferably the median.

[0012] In certain embodiments, the support number may be the number of aligned fragment pairs supporting each candidate fusion gene.

[0013] And / or, the neighbor distance between any two data points in the data set is Euclidean distance or Manhattan distance, preferably Euclidean distance.

[0014] In certain embodiments, according to the method for analyzing a sample based on sequence data of the present invention, wherein the sample is a biological sample,

[0015] And / or, the sequence data is selected from the sequencing data of any one or more of polynucleotides, polypeptides, polysaccharides, and lipids, preferably polynucleotides,

[0016] And / or, the density-based spatial clustering is selected from DBSCAN, OPTICS, or DENCLUE, preferably DBSCAN,

[0017] and / or, the data set obtained based on the sequence data from the sample is a data set that has been subjected to normalization and / or denoising processing by non-density-based clustering analysis,

[0018] And / or, each data point in the data set obtained based on the sequence data from the sample includes at least two dimensions, preferably, the at least two dimensions include a time dimension and / or a space dimension.

[0019] In certain embodiments, according to the method for analyzing a sample based on sequence data of the present invention, the actual situation of the sample includes the actual situation of gene fusion in the biological sample.

[0020] And / or, each data point in the data set obtained based on the sequence data from the sample is the breakpoint position coordinates x of two fragments of the same sequencing read of the sample in two different genes corresponding to the same reference sequence or reference genome. i and Y j .

[0021] In certain embodiments, according to the method for analyzing a sample based on sequence data of the present invention, the method for obtaining a data set based on sequence data from a sample comprises:

[0022] Aligning the reads in the sequencing data from the sample to a reference sequence or a reference genome, and obtaining aligned fragments, wherein the aligned fragments are portions of the continuous sequence regions in the reads that can be aligned to the reference sequence or the reference genome;

[0023] Screening according to the alignment lengths of the alignment fragments to obtain a set of candidate alignment fragments;

[0024] Determine the gene information corresponding to each alignment fragment in the candidate alignment fragment set according to the reference sequence or the reference genome annotation file, and combine two alignment fragments from the same read segment and corresponding to different genes into an alignment fragment pair, thereby obtaining initial fusion information, which includes the gene information and breakpoint coordinate information corresponding to each of the first alignment fragment and the second alignment fragment;

[0025] Perform consistency analysis on the fusion breakpoints in the initial fusion information according to the breakpoint coordinate information, exclude inconsistent alignment fragment pairs, and obtain an alignment fragment set after filtering the fusion information, that is, the data set;

[0026] The data set is weightedly reconstructed based on the total sequencing depth at the breakpoint position, which includes: calculating the weight coefficient fn based on the larger sequencing depth Dn near the two breakpoints of each alignment fragment pair, where n is the number of the candidate alignment fragment set, and scaling the coordinate matrix of the fusion information of each alignment fragment pair based on fn, wherein:

[0027]

[0028] The stretching is to multiply fn by the coordinate matrix of the fusion information of the corresponding aligned fragment pairs;

[0029] Preferably, the coordinate matrix of the fusion information of each aligned fragment pair is simplified, and more preferably, at least the last three digits of the genomic coordinate values ​​of the coordinates are retained after the simplification.

[0030] In certain embodiments, according to the method for analyzing samples based on sequence data of the present invention, the step of obtaining the aligned fragments further comprises the step of filtering and removing the aligned fragments with multiple alignments.

[0031] And / or, the screening of the candidate alignment fragment set according to the alignment length of the alignment fragment includes, after preliminarily confirming the alignment fragment corresponding to the potential fusion gene according to the SA tag, further screening the candidate alignment fragment set according to the alignment length of the alignment fragment,

[0032] and / or, the calculation of the alignment length in the candidate alignment segment set obtained by screening according to the alignment length of the alignment segment includes analyzing the sequence in the alignment segment according to the CIGAR information to obtain a string containing at least one of M, I, D, N, H and S, counting the number of M in the string in order until H or S appears in the string, and terminating the counting, and calculating the sum of the lengths corresponding to each M as the alignment length,

[0033] and / or, the screening according to the comparison length of the comparison fragments to obtain the candidate comparison fragment set includes removing the comparison fragments whose comparison length is less than a threshold, and simultaneously removing the comparison fragments whose coordinates of the sequences of the M region in the first comparison fragment and the second comparison fragment overlap by a length greater than a threshold,

[0034] and / or, when |the minimum value between the ending coordinate position of gene X and the ending coordinate position of the M region of the aligned fragment|-|the maximum value between the starting coordinate position of gene X and the starting coordinate position of the M region of the aligned fragment| is greater than a threshold value, the aligned fragment is determined to correspond to gene X;

[0035] And / or, performing consistency analysis on the fusion breakpoints in the initial fusion information based on the breakpoint coordinate information includes counting the distribution of the breakpoint coordinates of the first aligned fragment in each aligned fragment pair and the breakpoint coordinates of the second aligned fragment in each aligned fragment pair among multiple aligned fragment pairs supporting the same fusion, and removing the aligned fragments corresponding to the breakpoint coordinates that deviate from the threshold.

[0036] In certain embodiments, the method for analyzing samples based on sequence data according to the present invention further includes a step of determining whether the fusion gene has a complete gene structure. If so, the result of the actual situation of the gene fusion in the sample is determined to be credible; if not, the result of the actual situation of the gene fusion in the sample is determined to be unreliable.

[0037] In a second aspect, the present invention provides a method for analyzing fusion genes based on sequencing data, comprising: performing cluster analysis on a data set obtained based on sequencing data from a sample, and judging the true situation and / or false situation or abnormal situation of gene fusion in the sample based on the analysis, wherein:

[0038] Each data point in the data set is the breakpoint position coordinates X of two fragments of the same sequencing read of the sample in two different genes corresponding to the same reference sequence or reference genome. i and Y j ; The data set is obtained by weighted reconstruction based on the total sequencing depth at the breakpoint position.

[0039] In most cases, if the original breakpoint information is directly used for clustering analysis, there will be more false positive information. In order to reduce false positives as much as possible, the original breakpoint information is deeply weighted to obtain a new weighted value that balances the depth and the number of fusion genes. The new weighted value can be used to normalize the background noise as much as possible. The new data set obtained is more conducive to estimating the true fusion gene information in the subsequent clustering algorithm.

[0040] In certain embodiments, the method for obtaining the dataset comprises:

[0041] Aligning the reads in the sequencing data from the sample to a reference sequence or a reference genome, and obtaining aligned fragments, wherein the aligned fragments are portions of the continuous sequence regions in the reads that can be aligned to the reference sequence or the reference genome;

[0042] Screening according to the alignment lengths of the alignment fragments to obtain a set of candidate alignment fragments;

[0043] Determine the gene information corresponding to each alignment fragment in the candidate alignment fragment set according to the reference sequence or the reference genome annotation file, and combine two alignment fragments from the same read segment and corresponding to different genes into an alignment fragment pair, thereby obtaining initial fusion information, which includes the gene information and breakpoint coordinate information corresponding to each of the first alignment fragment and the second alignment fragment;

[0044] The fusion breakpoints in the initial fusion information are analyzed for consistency according to the breakpoint coordinate information, and inconsistent alignment fragment pairs are excluded to obtain an alignment fragment set after filtering the fusion information, thereby obtaining the data set.

[0045] In certain embodiments, the data set is weighted and reconstructed based on the total sequencing depth at the breakpoint position, comprising:

[0046] Based on the larger sequencing depth Dn near the two breakpoints of each alignment fragment pair (such as 20 bp upstream and downstream), n is the number of the candidate alignment fragment set, the weight coefficient fn is calculated, and the coordinate matrix of the fusion information of each alignment fragment pair is scaled based on fn, where:

[0047]

[0048] The stretching is to multiply fn by the coordinate matrix of the fusion information of the corresponding aligned fragment pairs;

[0049] Preferably, the coordinate matrix of the fusion information of each aligned fragment pair is simplified, and more preferably, at least the last three digits of the genomic coordinate values ​​of the coordinates are retained after the simplification.

[0050] In certain embodiments, the cluster analysis uses density-based spatial cluster analysis;

[0051] Preferably, the true situation of the sample is judged based on the clusters obtained by the analysis, and / or the untrue situation or abnormal situation of the sample is judged based on the discrete data relative to the clusters obtained by the analysis.

[0052] Preferably, the density-based spatial clustering is selected from DBSCAN, OPTICS, or DENCLUE, more preferably, DBSCAN.

[0053] In certain embodiments, the cluster-dependent distance parameter e is obtained by calculating,

[0054] Preferably, the cluster-dependent distance parameter e is the median, mean or percentile of the neighbor distance between any two data points in the data set. More preferably, the neighbor distance is the Euclidean distance.

[0055] And / or, the minimum number k is the median, mean or percentile of the number of supports for each data point in the data set.

[0056] In a third aspect, the present invention provides a device for analyzing a sample based on sequence data, comprising:

[0057] a data acquisition unit configured to acquire sequence data;

[0058] A data processing unit is configured to be able to perform the method described in the first aspect or the second aspect.

[0059] A fourth aspect of the present invention provides a device for analyzing samples based on sequence data, which includes a memory, a processor, and a program stored in the memory and capable of running on the processor, wherein the program is configured to implement the method described in the first aspect or the second aspect.

[0060] According to a fifth aspect of the present invention, a computer-readable storage medium is provided, wherein the computer-readable storage medium stores a program for analyzing samples based on sequence data, and when the program is executed by a processor, the method according to the first aspect or the second aspect is implemented.

[0061] The density-based spatial clustering analysis of the method of the present invention does not pre-limit the minimum number k of data points contained in the cluster, or the distance parameter e on which the cluster depends and the minimum number k of data points contained in the cluster. It can improve the accuracy of the results, reduce unnecessary debugging processes or subsequent verifications due to pre-limitation, reduce the amount of calculation, and improve the efficiency of result judgment, thereby improving the accuracy of judgment on sample information.

[0062] In an exemplary embodiment, the method of the present invention is used to analyze gene sequencing data, which is applicable to third-generation sequencing platforms including QNome of Qitan Technology, Oxford Nanopore and Pacific Biosciences, and is also applicable to second-generation sequencing platforms such as BGI and Illumina. After full-length transcriptome sequencing or targeted capture and full-length transcript sequencing based on DNA or cDNA libraries, there is no need to assemble the sequencing fragments, and the full-length transcript can be directly obtained for fusion gene analysis results, which significantly improves the sensitivity and specificity of the fusion gene detection results. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 The figure is an exemplary flow chart for illustrating the method of the present invention by taking fusion gene analysis as an example.

[0064] Figure 2 This is a flowchart for comparing and screening the sequence fragments obtained by exemplary sequencing with the reference genome.

[0065] Figure 3 The figure shows the results of clustering analysis of exemplary fusion gene breakpoints.

[0066] Figure 4-7 To exemplify different fusion types to determine whether the fusion gene has a complete gene structure. DETAILED DESCRIPTION

[0067] Various exemplary embodiments of the present invention will now be described in detail. This detailed description should not be considered as limiting the present invention, but should be understood as a more detailed description of certain aspects, features, and embodiments of the present invention.

[0068] It should be understood that the terms described in the present invention are only for describing special embodiments and are not intended to limit the present invention. In addition, for the numerical range in the present invention, it should be understood that the upper and lower limits of the scope and each intermediate value therebetween are specifically disclosed. Each smaller range between the intermediate value in any stated value or stated range and any other stated value or intermediate value in the described range is also included in the present invention. The upper and lower limits of these smaller ranges can be independently included or excluded in the scope.

[0069] Unless otherwise indicated, all technical and scientific terms used herein have the same meanings as commonly understood by those skilled in the art to which the invention belongs. Although the present invention describes only preferred methods and materials, any methods and materials similar or equivalent to those described herein may also be used in the implementation or testing of the present invention. All documents mentioned in this specification are incorporated by reference to disclose and describe the methods and / or materials related to the documents. In the event of a conflict with any incorporated document, the content of this specification shall prevail.

[0070] Herein, the term "sequence data" refers to a data set arranged in a certain order, and each element in the sequence data is called an item of the sequence. The sequence data can be finite or infinite. Sequence data usually has a temporal or spatial association, and their order is important because the position information of each item is of great significance for understanding and processing data. Examples of sequence data of the present invention include but are not limited to text data, time series data, audio signals, videos, DNA sequences, and rows or columns of images. Among them, text data includes sentences, paragraphs, articles, etc., which are all text sequence data. In natural language processing, text data is usually divided into words or characters to form a sequence. Time series data refers to data recorded in chronological order, such as stock prices, temperatures, electrocardiograms, etc. The data at each time point can be regarded as an item of the sequence. An audio signal can be regarded as a continuous sequence, and the amplitude value at each time point constitutes a part of the sequence. A video can be regarded as a series of image frames, and each frame is an item of the sequence. The rows or columns of an image, when processing an image, each row or column of the image can be regarded as a sequence. A nucleic acid sequence refers to a sequence composed of a series of bases, including a DNA sequence and an RNA sequence. Amino acid sequence refers to a sequence composed of a series of amino acids. Nucleic acid sequence and amino acid sequence are naturally occurring sequence data in biology.

[0071] Herein, the term "sample" refers to an object to be analyzed containing sequence data, including natural samples and artificially generated samples, examples of which include but are not limited to biological samples such as any tissue containing nucleic acids or proteins, and file samples such as audio, video or text files. Examples of biological samples include but are not limited to polynucleotides, polypeptides, polysaccharides and lipids.

[0072] Herein, the term "read" refers to a sequence fragment of a nucleotide (such as DNA or RNA) sample or reference, and is sometimes written as "read(s)" herein. Herein, a "read" may be a sample read and / or a reference read. Typically, a read represents a short sequence of adjacent base pairs in a sample or reference. A read may be represented by a base pair sequence (i.e., ATCG format) of a sample or reference fragment. The read may be obtained directly from a sequencing instrument or device or indirectly from stored sequence information about a sample. The read may be stored in a memory device and may be processed as needed to determine whether the read matches a reference sequence or meets other criteria. The length of the read is generally more than 10 bp, such as more than 50 bp, more than 100 bp, more than 200 bp, more than 500 bp, more than 800 bp, more than 1000 bp, more than 2000 bp, more than 5000 bp, or even more than 1 Mbp, more than 2 Mbp, more than 5 Mbp, etc. Exemplarily, the read may be directly obtained from a second-generation sequencing instrument or a third-generation sequencing instrument. Exemplary sequencing technologies include, but are not limited to, sequencing by synthesis, pyrosequencing, ion semiconductor technology, single molecule real-time sequencing, and sequencing by ligation.

[0073] As used herein, the term "alignment" refers to the process of comparing a read segment with a reference sequence to obtain detailed alignment information of the read segment and determine whether the read segment is within the reference sequence or corresponds to the position of the reference sequence.

[0074] Herein, the term "reference sequence" or "reference genome" refers to an identified sequence that can be used as a reference and is from a subject, which can be any specific known genome sequence of any organism. Exemplarily, the reference genome can be a reference genome of humans or many other organisms collected by the National Center for Biotechnology Information. "Genome" refers to the complete genetic information of an organism or virus expressed in a nucleic acid sequence. The genome includes both genes and non-coding sequences of DNA. The reference sequence may be larger than the read segment aligned with it. In some examples, the reference genome sequence is the sequence of the full-length human genome. In some other examples, the reference genome sequence is limited to a specific human chromosome. In some specific implementations, the reference chromosome is a chromosome sequence from the human genome version hg19. In a specific implementation, the reference genome can be a consensus sequence or other combination derived from multiple individuals. However, in some cases, the reference sequence can also be taken from a specific individual.

[0075] As used herein, the term "sequence" refers to a string or chain of nucleotides that include or represent nucleotides coupled to each other. Nucleotides may be DNA or RNA. A sequence may include or may be divided into a plurality of subsequences. A read or reference genome herein may consist of one or more sequences. In order to distinguish between various different sequences, different labels (e.g., reference sequences, etc.) may be given to the sequence.

[0076] Herein, the term "alignment fragment" refers to the portion of a read that can be aligned to a continuous sequence region of a reference genome. Wherein "aligned to" means that a portion of the sequence in the read is substantially consistent with the sequence in the reference genome. Substantial consistency is generally understood to mean that a portion of the sequence in the read has more than 50%, preferably more than 60%, more than 70%, more than 80%, more than 90%, more than 95%, or even 100% sequence identity with the sequence in the reference genome. Herein, a read may contain one or more alignment fragments. In the case where the read contains a fusion site, the read typically contains at least two alignment fragments.

[0077] Herein, the term "sequence identity" refers to the sequence identity at the nucleotide level between genes, which is a measure at the nucleotide level. Sequence identity can be determined by comparing the nucleotide sequence at a given position in each sequence when comparing sequences. The method for sequence alignment is well known in the art, and such methods include GAP, BESTFIT, BLAST, FASTA and TFASTA. The BLAST algorithm calculates sequence identity percentage and performs statistical analysis on the similarity between two sequences. The software for performing BLAST analysis can be publicly obtained by, for example, the National Center for Biotechnology Information (NCBI) website.

[0078] Herein, the term "aligned fragment pair" refers to two aligned fragments consisting of two aligned fragments from the same read segment and corresponding to different genes. In order to distinguish the two aligned fragments in the aligned fragment pair, one of the aligned fragment pair can be called the first aligned fragment and the other can be called the second aligned fragment.

[0079] Herein, the terms "first comparison fragment" and "second comparison fragment" represent two different comparison fragments, wherein "first" and "second" are used only for the purpose of distinction, and do not have the meaning of order or direction such as before or after or before and after. Exemplarily, the first comparison fragment and the second comparison fragment in this article are comparison fragments of the same read segment and different regions of the reference genome, respectively. For example, the first comparison fragment is aligned to a partial sequence of gene A, and the second comparison fragment is aligned to a partial sequence of gene B. Since the first and second have no meaning of before or after or before and after order, the fusion gene can be gene A-gene B or gene B-gene A at this time. Wherein, "-" here represents a covalent bond. Therefore, the first comparison fragment can correspond to the 5' end of the fusion gene, or it can correspond to the 3' end of the fusion gene. Similarly, the second comparison fragment can correspond to the 3' end of the fusion gene, or it can correspond to the 5' end of the fusion gene.

[0080] Herein, the term "coordinate" refers to the position of a specific nucleotide in a reference genome, in particular a specific nucleotide in a sequence or read or a fragment thereof corresponding to the position of the reference genome. The coordinate is usually expressed as a distance away from a fixed position (such as the 5' start point of the reference genome).

[0081] Herein, the term "alignment length" refers to the length of the sequence in the alignment fragment that is completely consistent with the nucleotide at the corresponding position of the reference genome. The alignment length is usually expressed as the number of nucleotides, usually in bp.

[0082] Herein, the term "candidate alignment fragment set" refers to a set of multiple alignment fragments obtained by screening and removing potential noise fragments according to the alignment length and coordinates of the alignment fragments. The alignment fragments in the candidate alignment fragment set come from multiple different reads, which may contain at least one alignment fragment pair. The candidate alignment fragment set can be further used to obtain different alignment fragment pairs.

[0083] Herein, the term "gene fusion" refers to the fusion of part or all of the sequences of two different genes due to a certain mechanism (such as genomic variation) to form a new gene. Exemplarily, the fusion gene is produced by chromosome rearrangement, including chromosome translocation, insertion, inversion, deletion, etc. The generation of fusion genes changes the protein coding sequence or regulatory sequence of the gene, causing changes in gene function, which has a greater impact on the body. In general, gene fusion refers to fusion at the genome level, but fusion may also occur at the transcriptome level, mainly because the RNA produced by transcription of two different genes is fused together for some reason to form a new fusion RNA. The RNA may encode protein or may be non-coding. The fusion gene generated at the genome level may be expressed or not expressed (such as destruction of the promoter region or other reasons) depending on the fusion situation. Gene fusion is closely related to the occurrence and development of various diseases, especially cancer, and is even a direct cause of some cancers, so gene fusion has also become an important research content in the current omics big data analysis. At present, it has been reported that the occurrence of many cancers is closely related to gene fusion. These fusion genes may also be potential drug targets. It is very necessary to conduct in-depth research on them, which is of great significance for the prevention, treatment and comprehensive understanding of such tumor diseases. In actual analysis, the information of predicted fusion gene sites can be used to obtain evidence involving transcription sequence mutations and structural rearrangements, and then find changes in the cancer genome related to its function. The identification of gene fusion can be based on whole genome sequencing data (WGS), or transcriptome sequencing data (RNA-seq), or a combination of the two technologies is better. The gene fusion identified by whole genome sequencing can basically determine that it is caused by a certain mutation at the genome level, but if there is no transcriptome sequencing data, it is impossible to accurately determine whether the new gene generated after the fusion can be expressed, or the expression level. The gene fusion identified by transcriptome sequencing data can clearly be a gene fusion that can be expressed, but it is impossible to completely determine whether it is caused by genomic mutations or RNA fusion after transcription of two different genes. Therefore, if conditions permit, combining whole genome sequencing and transcriptome sequencing to identify gene fusion can obtain more accurate identification results.

[0084] Herein, the term "initial fusion information" includes any information indicating that a gene is a potential fusion gene, including but not limited to the name of the fusion gene, the chromosome number of the alignment fragment corresponding to the reference genome, the breakpoint coordinates of the alignment fragment, etc. Exemplarily, the initial fusion information is provided in the form of a unit corresponding to the alignment fragment pair. For example, for an alignment fragment pair, the corresponding initial fusion information can be expressed as (gene A--gene B, the chromosome number corresponding to the first alignment fragment, the breakpoint coordinates of the first alignment fragment, the chromosome number corresponding to the second alignment fragment, the breakpoint coordinates of the second alignment fragment).

[0085] Methods for analyzing samples based on sequence data

[0086] In a first aspect, the present invention provides a method for analyzing samples based on sequence data, which includes but is not limited to: performing density-based spatial clustering analysis on a data set obtained based on sequence data from the sample, judging the true situation of the sample based on the clusters obtained by the analysis, and / or judging the untrue situation or abnormal situation of the sample based on discrete data relative to the clusters obtained by the analysis.

[0087] In the present invention, outlier data is filtered out by density-based spatial cluster analysis to obtain a set of clustered alignment fragment pairs. This step includes calculating the density around each data point to determine the cluster structure without pre-specifying the number of clusters, thereby automatically discovering at least one cluster and discrete data (i.e., noise points) relative to the cluster.

[0088] In the present invention, the data set is deeply weighted reconstructed, and the minimum number k of data points contained in the cluster and the distance parameter e on which the cluster depends are obtained by calculation rather than being directly set or preset based on experience.

[0089] In certain embodiments, the data set is obtained by recording the left and right breakpoints of all candidate fusion genes in the sample and storing them as a matrix A:

[0090] A=[[F1:[L1,L2,L3,L4...],[R1,R2,R3,R4...]],

[0091] [F2:[L1,L2,L3,L4...],[R1,R2,R3,R4...]],

[0092] [F3:[L1,L2,L3,L4...],[R1,R2,R3,R4...]],

[0093] ...,

[0094] [Fn:[L1,L2,L3,L4...],[R1,R2,R3,R4...]]]

[0095] Among them, F represents the fusion gene number, L represents the coordinate value of the left breakpoint after processing, and R represents the coordinate value of the right breakpoint after processing. The L and R coordinate processing method is to take the last 3 digits of the genome coordinate value, in order to limit the value to hundreds to reduce the amount of calculation, such as chr22:23290411, the processed coordinate value is 411, chr9:130854065, the processed coordinate value is 65.

[0096] In certain embodiments, the dataset A is weighted and reconstructed based on the total sequencing depth at the breakpoint coordinates. The specific method is to obtain the average depth of 20 bp near the breakpoint of each candidate fusion gene pair F1, F2...Fn, and define the larger value of the two as D1, D2..Dn. The weight coefficient of each fusion gene Fn is calculated based on D1..Dn, and according to the weight coefficient fn, the coordinate matrix of each pair of candidate fusion genes is expanded or contracted according to the fn coefficient by copying or deleting, where the calculation method of fn is:

[0097]

[0098] The method for scaling the fusion coordinate matrix based on fn is:

[0099] A'=[[F1:[L1,L2,L3,L4...]*f1,[R1,R2,R3,R4...]*f1],

[0100] [F2:[L1,L2,L3,L4...]*f2,[R1,R2,R3,R4...]*f2],

[0101] [F3:[L1,L2,L3,L4...]*f3,[R1,R2,R3,R4...]*f3],

[0102] ...,

[0103] [Fn:[L1,L2,L3,L4...]*fn,[R1,R2,R3,R4...]*fn]]

[0104] The calculation rule for multiplying a list by fn is:

[0105] Then, according to the number of lists of F'n, randomly select values ​​from Fn until the value that satisfies the number of lists of F'n is obtained. If the number of F'n is less than 5, only the first value is taken as noise. For example, Fn = [[2,3,2,2], [10,8,8,9]], the number of elements is 4, indicating that the fusion gene has obtained 4 detected breakpoints on the upstream and downstream genes. Assuming fn = 1.5, the number of F'n after scaling is 4 times 1.5, that is, 6. F'n randomly selects values ​​from the elements of Fn, and needs to take 6 values, then: F'n = [[2,3,2,2,2,3], [10,8,8,9,10,8]]. Assuming fn = 0.6, 4 times 0.6 is rounded up to 3, F'n = [[2,3,2], [10,8,8]]. The reconstructed matrix A is:

[0106] A'=[[F'1:[L1,L2,L3,L4...],[R1,R2,R3,R4...]],

[0107] [F'2:[L1,L2,L3,L4...],[R1,R2,R3,R4...]],

[0108] [F'3:[L1,L2,L3,L4...],[R1,R2,R3,R4...]],

[0109] ...,

[0110] [F'n:[L1,L2,L3,L4...],[R1,R2,R3,R4...]]]

[0111] In certain embodiments, the neighbor distance e, a parameter required for clustering, is obtained by calculating the distance d between two points based on the reconstructed matrix A' using the coordinates of the upstream and downstream breakpoints of each fusion gene reconstructed F'n, and calculating the median of all breakpoints d of all fusion genes.

[0112] in x, y correspond to the coordinate values ​​L and R of the upstream and downstream breakpoints of F'n respectively.

[0113] In certain embodiments, the minimum number k of data points contained in the cluster is the number of fusion gene breakpoints from F'1 to F'n in the reconstructed matrix A', sorted from small to large, to form a list B = [b1, b2...bn]. When the number in the B list is greater than or equal to 3, k is the value in the middle of the list B; when B is less than 3, k is the value in the middle of the list B multiplied by 0.6. The specific position is calculated by int(len(B) / 2)+1, where len represents the length of B and int represents rounding down the value.

[0114] The density-based spatial clustering of the present invention is not limited, and can be, for example, DBSCAN, OPTICS, or DENCLUE. In certain preferred embodiments, it is based on DBSCAN and further optimized, wherein the definition based on the algorithm includes: object domain Neps(q): the neighborhood within a given object radius eps (Euclidean distance) is called the neighborhood of the object, and Neps(q) is defined as: Neps(q):{p belongs to D||distance(p,q)<=eps}, p, q is a point in D, and D is all points in the data set. distance(p,q) refers to the Euclidean distance from p to q. Neps(q) contains two elements: the domain radius eps and the core object of the point in the domain: if the neighborhood of an object contains at least the minimum number of MinPts (MinPts is the minimum number of domain points for a given point to become a core object in the domain) of objects, then the object is called a core object. |Neps(p)|>=MinPts. Boundary object: the neighborhood of an object is less than MinPts objects, but is in the neighborhood of a core object. Noise point: The neighborhood of the object is less than MinPts objects and it is not a boundary point. Direct density reachable: Given a set of objects D, if p is in the neighborhood of q and q is a core object, then object p is directly density reachable from object q. Density reachable: If there exists a chain of objects p1, p2, ..., pn such that for pi belongs to D (1≤i≤n-1), pi+1 is directly density reachable from pi with respect to eps and MinPts, then object pn is density reachable from p1 with respect to eps and MinPts.

[0115] In the optimization process of this step, the input parameters are: the dataset of the breakpoint coordinates of the first and second alignment segments: dataset and the domain radius (epslion) set based on empirical values. The clustering process includes: calculating the minimum domain point number k of the core object and clustering based on the obtained k value and epslion as input.

[0116] In the present invention, the specific procedure flow of clustering is as follows:

[0117] 1. If the number of objects contained in Neps(p) is less than k: mark object p as a boundary point or a noise point, indicating that the density around the point is low, and skip processing object p.

[0118] 2. If the number of objects contained in Neps(p) is greater than or equal to k: mark object p as a core point, indicating that the density around this point is high. Create a new cluster Cluster and add object p to the cluster Cluster.

[0119] 3. For all unprocessed objects q in Neps(p): check the domain Neps(q) of object q. If the number of objects contained in Neps(q) is greater than or equal to k: add the objects in Neps(q) that are not included in any cluster to cluster Cluster, and expand the members of cluster Cluster. Recursively process the domain of object q and continue to explore more objects that are density-reachable from the core point object p.

[0120] 4. When all objects in the dataset have been processed, the clustering process ends. At the same time, cluster labels are added to each point in the dataset. If the value of the cluster label is -1, it means it is unclassified, that is, it is a noise point and is filtered. The clusters whose labels are not -1 are the collections of credible fusion gene breakpoints. The clustering-based method can further filter the fusion information of the alignment fragments.

[0121] In certain embodiments, the dataset obtained based on sequence data from a sample is a dataset that has been subjected to normalization and / or denoising by non-density-based clustering analysis.

[0122] In certain embodiments, the method based on sequence data analysis sample of the present invention is a method for analyzing gene fusion, which is included in the comparison fragment set, and further filtering outlier data based on density cluster analysis to obtain clustered comparison fragment set. This step can obtain highly credible comparison fragments supporting fusion, thereby further reducing false positives and false negatives. Due to nanopore sequencing, if the length of the fused transcript is long, insertion or deletion caused by sequencing will be introduced, resulting in a large fluctuation in the position of the fusion breakpoint. If algorithm processing is not performed, a large amount of false positives and false negatives will be produced. Therefore, by designing a highly sensitive clustering algorithm, the breakpoint range can be more accurately reduced.

[0123] In a preferred embodiment, the aligned fragment set of the present invention is obtained by a method comprising the following steps:

[0124] (1) aligning the reads in the sequencing data to a reference genome, and screening to obtain aligned fragments, wherein the aligned fragments are portions of the reads that can be aligned to a continuous sequence region of the reference genome;

[0125] (2) screening to obtain a set of candidate alignment fragments according to the alignment lengths of the alignment fragments;

[0126] (3) determining the gene information corresponding to each alignment fragment in the candidate alignment fragment set according to the reference genome annotation file, and combining two alignment fragments from the same read segment and corresponding to different genes into an alignment fragment pair, thereby obtaining initial fusion information, which includes the gene information and breakpoint coordinate information corresponding to each of the first alignment fragment and the second alignment fragment;

[0127] (4) Perform consistency analysis on the fusion breakpoints in the initial fusion information according to the breakpoint coordinate information, exclude inconsistent alignment fragment pairs, and obtain an alignment fragment set after filtering the fusion information.

[0128] Those skilled in the art will appreciate that the above numbers (1)-(4) etc. are only for the purpose of distinguishing different steps and do not represent the order of the steps. As long as the purpose of the present invention can be achieved, the order of the above steps is not particularly limited. In addition, two or more of the above steps may be combined and performed simultaneously. In addition, those skilled in the art will also appreciate that other steps or operations may be included before or after the above steps (1)-(4), or between these arbitrary steps, for example, to further optimize and / or improve the method described in the present invention.

[0129] The above step (1) of the present invention is a step of obtaining a comparison fragment, which includes comparing the reads in the sequencing data to the reference genome, thereby obtaining a comparison fragment. The comparison can be performed in a manner known in the art, which is not particularly limited. For example, existing software such as minimap2 can be used for comparison. In certain embodiments, step (1) further includes a step of screening the obtained comparison fragments. Exemplarily, the screening here includes filtering and removing the comparison fragments of multiple comparisons. Here, "multiple comparisons" means that a certain sequence fragment can be compared to multiple positions on the reference genome.

[0130] Step (2) of the present invention is a step of screening the alignment fragments, which includes screening to obtain a set of candidate alignment fragments according to the alignment lengths of the alignment fragments.

[0131] In the present invention, the exemplary calculation method of the alignment length is based on CIGAR information. In CIGAR, M indicates that the alignment fragment sequence is optimally aligned with the reference genome, I indicates that the alignment fragment sequence has an insertion in the reference genome, D indicates that the alignment fragment sequence has a deletion in the reference genome, N indicates that the alignment fragment skips a certain region of the reference genome, S indicates that the alignment fragment is not aligned but the alignment result is retained, and H indicates that the alignment fragment is not aligned and the alignment result is not retained.

[0132] The specific calculation method is exemplified as follows: For a certain comparison fragment, the CIGAR information of a single nucleotide or a small sequence (or region) composed of multiple nucleotides in the comparison fragment can be confirmed according to the rules shown in Table 1, thereby obtaining a string composed of M, I, D, N, S and H, in which each letter corresponds to a certain nucleotide length. Starting from a specific position of the string, the lengths of single nucleotides and small sequences corresponding to the M region in the string are counted in sequence, and the sum of the lengths corresponding to all M is calculated to obtain the comparison length. When I, D and N appear in the string, the corresponding length is counted as 0, that is, these cases are skipped. When S or H appears in the string, the length of the M region is stopped.

[0133] Screening based on alignment length includes setting a threshold of alignment length and removing alignment fragments whose alignment length is less than the threshold. The threshold here is obtained by statistics or is an empirical value, which can be changed according to different application scenarios and requirements. Exemplary thresholds here can be, for example, 10, 20, 30, 50, 90, 100, 120 bases, etc.

[0134] In certain embodiments, the screening according to the comparison length further includes removing the situation where the coordinates of the first comparison segment and the second comparison segment overlap. Specifically, when the coordinate overlap length of the sequence in the M region in the first comparison segment and the second comparison segment is greater than a certain threshold value, it indicates that such segments belong to noise and need to be excluded. The threshold value here is obtained by statistics or is an empirical value, which can be changed according to different application scenarios and requirements. The exemplary threshold value here can be, for example, 10, 20, 30, 50, 90, 100, 120 bases, etc.

[0135] In certain embodiments, before obtaining a candidate alignment fragment set according to the alignment length of the alignment fragment and the coordinate screening of the alignment fragment, it further includes preliminarily confirming the alignment fragment corresponding to the potential fusion gene according to the SA tag. Wherein, the SA tag is known data, which includes a supplementary alignment record of the alignment fragment. For example, if the genes corresponding to the two alignment fragments are fused, the two alignment fragments are divided into two independent parts, and can be aligned to different positions of the reference genome respectively. The SA tag can be used to record the alignment information on the left side of the alignment fragment and the right side of the alignment fragment. If the alignment fragment contains an SA tag, it indicates that gene fusion may occur. Therefore, two alignment fragments or alignment fragment pairs corresponding to potential fusion genes can be preliminarily screened based on the SA tag.

[0136] Step (3) of the present invention is a step of obtaining initial fusion information, which includes determining the gene information corresponding to each alignment fragment in the candidate alignment fragment set according to the reference genome annotation file, and forming an alignment fragment pair with two alignment fragments from the same read segment and corresponding to different genes, thereby obtaining initial fusion information.

[0137] In certain embodiments, the gene information corresponding to each alignment fragment can be determined by calculating |the minimum value of the ending coordinate position of gene X and the ending coordinate position of the alignment fragment M region|-|the maximum value of the coordinate starting position of gene X and the starting coordinate position of the alignment fragment M region|. Specifically, when the obtained value is greater than a certain threshold, it indicates that the alignment fragment is aligned to the gene X. The threshold here is obtained by statistics or is an empirical value, which can be changed according to different application scenarios and requirements. Exemplary thresholds here can be, for example, 40, 45, 50, 55, 60, 65, 70, 75, 80 bases, etc.

[0138] Step (4) of the present invention is a step of obtaining a set of aligned fragments after filtering the fusion information, which includes performing consistency analysis on the fusion breakpoints in the initial fusion information according to the breakpoint coordinate information, excluding inconsistent alignment fragment pairs, and obtaining a set of aligned fragments after filtering the fusion information. Through consistency analysis of fusion breakpoints, false positives and false negatives are reduced.

[0139] In the present invention, the consistency analysis of the fusion breakpoints in the initial fusion information according to the breakpoint coordinate information includes statistically analyzing the distribution of the breakpoint coordinates of the first comparison fragment and the breakpoint coordinates of the second comparison fragment in the comparison fragment pair supporting the same fusion, and removing the comparison fragments corresponding to the breakpoint coordinates that deviate from the threshold. Specifically, the breakpoint coordinates of the first comparison fragment and the second comparison fragment are converted into two coordinates, which contain the comparison fragment name (first element), the breakpoint coordinates of the first comparison fragment (second element), and the breakpoint coordinates of the second comparison fragment (third element). Sort in descending order by the number of occurrences of the second element. If the difference between the second element and the element with the most occurrences is less than or equal to the set biological distance, the comparison fragment is retained, otherwise the comparison fragment is discarded, and the set biological distance is obtained by statistics or is an empirical value, which can vary according to different application scenarios and needs. The exemplary values ​​here can be, for example, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20 bases, etc. Then, the comparison segments are sorted according to the third element with the most occurrences. If the difference between the third element and the element with the most occurrences is less than the biological distance, the comparison segment is retained, otherwise the comparison segment is discarded. In certain embodiments, if all values ​​of the multiple comparison segment pairs in the list do not meet the above conditions, only the first element is output. As an example, if the breakpoint coordinates in the comparison segment group each appear only once and there is no element with the most occurrences, only the breakpoint coordinates of the first comparison segment or the breakpoint coordinates of the second comparison segment in the segment group can be output. That is, the above-mentioned steps of the breakpoint coordinate consistency analysis filter out the comparison segments whose breakpoint coordinate positions are far away from the coordinate positions with the most occurrences.

[0140] In this step, through a primary filtering and screening, the sequencing fragments whose occurrence number of the second element is close to and equal to the maximum occurrence number according to the proximity criterion are retained. The sequencing fragments after the primary filtering and screening are subjected to a secondary filtering and screening, and the sequencing fragments whose occurrence number of the third element is close to and equal to the maximum occurrence number according to the proximity criterion are retained.

[0141] In certain embodiments, the analysis method of the present invention further includes a step of further determining whether the fusion gene has a complete gene structure in the clustered alignment fragment pair set. This step increases the marker of the coding potential of the fusion transcript, thereby providing an important basis for scientific research or clinical research on the abnormal protein expression caused by gene fusion leading to cancer. This step includes predicting whether the fusion gene has coding potential based on the coverage interval and breakpoint information of the positive and negative chains of the gene and each read segment on the genome, and then determining whether it has a complete gene structure.

[0142] In the present invention, the judgment method is based on the coverage interval of the positive and negative strands of the gene and each read segment on the genome, as well as the information of the breakpoints predicted in the previous step. For the two genes A and B of the fusion gene, the breakpoint coordinates on gene A / B, the positive and negative strand information, and the coverage of the read segments containing the breakpoints of gene A / B are recorded as "gene: (strand information, genome coordinate starting position, genome coordinate ending position, breakpoint position)". Specifically, this step includes assuming that gene A and gene B are two genes that are potentially fused, and judging the positions of gene A and gene B in the fusion gene respectively. If gene A and gene B are located at the 5' end and 3' end of the fusion gene respectively, the fusion gene is considered to have coding potential from the 5' end to the 3' end, otherwise it is considered to have no coding potential.

[0143] In a specific embodiment, for gene A, if |genomic coordinate starting position - breakpoint position|<|genomic coordinate ending position - breakpoint position|, and if the chain information is +, the breakpoint of the fusion gene is at the 5' end, and gene A in the fusion gene is located at the 3' end; if the chain information is -, the breakpoint of the fusion gene is at the 3' end, and gene A in the fusion gene contributes to the 5' end; if |genomic coordinate starting position - breakpoint position|>|genomic coordinate ending position - breakpoint position|, and if the chain information is +, the breakpoint of the fusion gene is at the 3' end, and gene A in the fusion gene is located at the 5' end; if the chain information is -, the breakpoint of the fusion gene is at the 5' end, and gene A in the fusion gene is located at the 3' end.

[0144] Equipment for analyzing samples based on sequence data

[0145] In a third aspect of the present invention, there is provided a device for analyzing a sample based on sequence data, referred to herein as a first device, which comprises at least:

[0146] A data acquisition unit configured to acquire sequencing data;

[0147] A data processing unit is configured to be able to perform the method described in the first aspect of the present invention.

[0148] In the present invention, the data acquisition unit can obtain sequencing data during or after sequencing. Optionally, the data acquisition unit can be connected to a sequencer for communication, thereby obtaining data during the sequencing process. Optionally, the data acquisition unit can be connected to a database for communication, thereby obtaining sequencing data in the database. The database can be known public data.

[0149] The fourth aspect of the present invention provides a device for analyzing samples based on sequence data, referred to herein as a second device, which includes a memory, a processor, and a program stored in the memory and capable of running the method described in the first or second aspect of the present invention on the processor.

[0150] In certain embodiments, the second device of the present invention is a computer device, which may be a server. The computer device includes a processor, a memory and a network interface connected by a system bus. Wherein, the processor of the computer device is used to provide computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system, a computer program and a database. The internal memory provides an environment for the operation of the operating system and the computer program in the non-volatile storage medium. The database of the computer device is used to store sequencing data to be processed. The network interface of the computer device is used to communicate with an external terminal via a network connection. When the computer program is executed by the processor, the method described in the first aspect or the second aspect of the present invention is implemented.

[0151] In the first device or the second device of the present invention, in addition to the above-mentioned units or devices and components, other units or devices may be further included, exemplarily including a display unit or device, etc. Each unit in the first device or the second device of the present invention may be implemented in whole or in part by software, hardware, and a combination thereof. Each unit may be embedded in or independent of any computer device in the form of hardware, or may be stored in a memory in the computer device in the form of software, so that the processor can call and execute the operations corresponding to the above-mentioned units.

[0152] Exemplarily, the first device or the second device of the present invention further comprises one or more output devices, which can be connected to, for example, a computer system. The example of an output device comprises a display, a printer, a communication device such as a modem and an audio output. Alternatively, further comprises one or more input devices, which can also be connected to a computer system. The example of an input device comprises a keyboard, a mouse, a writing pen and a writing board, a communication device and a data input device such as a sensor.

[0153] Computer readable storage medium

[0154] According to a fifth aspect of the present invention, a computer-readable storage medium is provided, wherein the computer-readable storage medium stores a program for analyzing samples based on sequence data, and when the program is executed by a processor, the steps of the method described in the first aspect or the second aspect of the present invention are implemented.

[0155] The storage medium of the present invention includes a readable medium in the form of a magnetic disk, an optical disk, a volatile storage unit, such as a random access memory unit (RAM) and / or a cache memory unit, and may further include a read-only memory unit (ROM).

[0156] Example

[0157] The following Figure 1 Take the analysis process of the present invention as an example:

[0158] 1. Sequencing data alignment

[0159] For example, the offline data Figure 2 Read1-read5 in the example is aligned to the reference genome to obtain a bam file, which is used as the input file of the system. The supported formats are fastq or bam files. If the input file is in fastq format, the software can call the alignment software minimap2 to align and obtain a bam file.

[0160] 2. Filter and remove multiple alignment fragments

[0161] Filter out the alignment fragments whose alignment information on the reference genome is multiple alignments. Figure 2 For example, Figure 2 The alignment fragments 1-5 in the figure are non-multiple alignment sequences obtained by aligning reads 1-5 of the downloaded data with the reference genome.

[0162] 3. Determination of potential fusion genes

[0163] Filter the alignment fragments containing SA tags from the bam file. SA tags are supplementary alignment records for alignment fragments. In the case where the genes corresponding to the two alignment fragments are fused, there are two alignment fragments that are divided into two independent parts. Align the alignment fragments with different positions of the reference genome respectively. If the alignment fragment contains SA tags, gene fusion may occur.

[0164] 4. Obtain the alignment lengths of the first alignment fragment (left) and the second alignment fragment (right) containing the SA tag and the gene names corresponding to the alignment fragments.

[0165] Figure 2As shown in the figure, in a comparison segment pair 1, the black box represents the first comparison segment (left), and the white box represents the second comparison segment (right) in the comparison segment pair 1. The first value from the left in the bracket of the black box or white box is the starting coordinate position of the comparison segment in the reference genome, and the second value from the left in the bracket is the ending coordinate position of the comparison segment in the reference genome.

[0166] The calculation method of the alignment length is as follows:

[0167] Using the CIGAR information of the alignment fragment (the meaning of the CIGAR of the alignment fragment is shown in Table 1), if the CIGAR of the alignment fragment is M, the length of the M region of the alignment fragment on the reference genome (in bp) is calculated. If the starting position of the CIGAR is S or H, it is skipped and not included in the length; if the CIGAR of the alignment fragment is I, D or N, it is skipped and not included in the length; until the end position of the CIGAR of the alignment fragment is H / S, it is terminated, and finally the length of the M interval of the alignment fragment is accumulated. Figure 2 For example, if the start and end coordinates of the right side of alignment fragment pair 1 to the reference genome are (500, 800), the CIGAR of the right side of alignment fragment pair 1 is 10S15M1I10M2D200M25N75M25H. The length of the alignment is calculated by the sum of the lengths of the characters M in CIGAR, that is, 15+10+200+75, which is 300bp.

[0168] Table 1 CIGAR meaning of alignment fragments

[0169]

[0170] Based on this method, the alignment lengths of the left and right fragments in the same alignment fragment pair are calculated. If the alignment length of the left or right fragment of the alignment fragment is less than the threshold of 90 bases, the alignment fragment is filtered out.

[0171] The confirmation of the gene name corresponding to the alignment fragment is based on the reference genome annotation file (gff / gtf). In the annotation file, if |the minimum value of the ending coordinate position of the gene and the ending coordinate position of the alignment fragment M|-|the maximum value of the starting coordinate position of the gene and the starting coordinate position of the alignment fragment M| is greater than the threshold, then it is considered that the gene is aligned. Here, the threshold is an empirical value, and its examples include but are not limited to 30, 50, 60, 70, 80, 90, 100, etc.

[0172] by Figure 2For example, the coordinates of the M region on the left side of the alignment fragment pair 1 on the reference genome are (5, 200), and the start and end coordinates of geneA on the reference genome are (2, 360). The calculation method is: |200-5|=195. When the threshold is 60, since 195 is greater than 60, the gene corresponding to the alignment fragment is annotated as geneA. The same method can be used to calculate that the right side of the alignment fragment pair 1 is annotated to geneB.

[0173] Thus, the initial fusion information of alignment fragment pair 1 is obtained (annotation information (gene name) on the left side of the alignment fragment pair - annotation information (gene name) on the right side of the alignment fragment pair, chromosome number on the left side of the alignment fragment pair on the reference genome, breakpoint coordinates on the left side of the alignment fragment pair on the reference genome, chromosome number on the right side of the alignment fragment pair on the reference genome, breakpoint coordinates on the right side of the alignment fragment pair on the reference genome). For example,

[0174] The initial fusion information of alignment fragment pair 1 is: (geneA-geneB, X, 200, X, 500);

[0175] The initial fusion information of alignment fragment pair 2 is: (geneA-geneB, X, 340, X, 507);

[0176] The initial fusion information of alignment fragment pair 3 is: (geneA-geneB, X, 201, X, 504);

[0177] The initial fusion information of alignment fragment pair 4 is: (geneA--geneB, X, 202, X, 580);

[0178] The initial fusion information of alignment fragment pair 5 is: (geneA--geneB, X, 200, X, 505).

[0179] The set of initial fusion information of aligned fragment pairs 1-5 is named fusion_primer_result=[initial fusion information of aligned fragment pair 1, initial fusion information of aligned fragment pair 2, initial fusion information of aligned fragment pair 3, initial fusion information of aligned fragment pair 4, initial fusion information of aligned fragment pair 5].

[0180] 5. Screening

[0181] If the length of the coordinate overlap of the sequence of the CIGAR M region on the left and right of the alignment fragment containing the SA tag is greater than the threshold, the alignment fragment is filtered. The threshold here is an empirical value, which can vary according to different situations and specific applications, and can be set by those skilled in the art based on experience. Exemplarily, the threshold here is 50, 60, 70, 80, 90, 100, 120 bases, etc.

[0182] Figure 2 The left side of the aligned fragment pair 6 in is aligned to geneA, and the coordinates on the right side span geneA and geneB and are aligned to geneB. Figure 2 The coordinate position of the M region on the left of the alignment fragment pair 6 is (90, 300), and the coordinate position of the M region on the right of the alignment fragment pair 6 is (205, 450). At this time, the coordinate overlap length of the left and right sides of the alignment fragment is 95 bases. When the threshold here is 90, since the coordinate overlap length is 95 bases, which is greater than the threshold, the alignment fragment pair 6 is filtered.

[0183] 6. If the annotated gene information on the left and right sides of the alignment fragment pair corresponds to the same gene, the alignment fragment will be filtered.

[0184] VII. Consistency Analysis of Fusion Breakpoints

[0185] The purpose of fusion breakpoint coordinate consistency analysis is to filter out alignment fragments whose breakpoint coordinate positions are far away from the coordinate positions with the highest number of occurrences, thereby reducing false positives and false negatives.

[0186] First, the left and right breakpoints of the filtered alignment fragments that support the same fusion (the annotation information on the left and right sides of the alignment fragments are the same, i.e., the gene names are the same) are sorted into two-dimensional coordinates. At this time, since the selection of alignment fragments that support the same fusion gene is limited, it is assumed that the chromosome numbers corresponding to the two fusion genes in this type of fragments are the same, and the chromosome information is ignored when the coordinate information is subsequently sorted. Figure 2 Take the example to illustrate, the form of (aligned fragment pair name, aligned fragment left breakpoint coordinates, aligned fragment right breakpoint coordinates) is used as an initial fusion information unit or array, where the first data in the brackets is the name of the aligned fragment pair, and the second and third data are the coordinates of the left and right breakpoints of the aligned fragment in the reference genome (ignoring the chromosome number), that is, Figure 2 The position indicated by the arrow in the alignment fragment 1.

[0187] According to the above format, the fragment pairs 1-5 are sorted and aligned to obtain the following list:

[0188] bp_info_list = [('Alignment fragment pair 1', 200, 500), ('Alignment fragment pair 2', 340, 507), ('Alignment fragment pair 3', 201, 504), ('Alignment fragment pair 4', 202, 800), ('Alignment fragment pair 5', 200, 505)].

[0189] In the list, first sort in descending order according to the frequency or number of occurrences of the second data in each array. If the difference between the second data in the array and the value with the highest frequency or number of occurrences is less than or equal to the set threshold, the array is retained, otherwise the array is deleted. The threshold here is an empirical value, such as biological distance, such as nucleotide length. The specific value can vary according to different situations and different needs, such as 10, 15, 20, 25, 30, 35, 40, 45, 50bp, etc.

[0190] by Figure 1 Taking the alignment fragment pairs 1-5 in the example, when the threshold is 15bp, the value of the left breakpoint coordinate that appears the most times is 200 for the alignment fragment pairs 1 and 5, and the left breakpoint coordinates of the alignment fragment pairs 3 and 4 are 201, 202 and 200, which are all within 15 distances from 200 and are therefore retained. The left breakpoint coordinates of the alignment fragment pairs 2 and 6 are 340 and 300, which are beyond 15 distances from 200, and the arrays corresponding to the alignment fragment pairs 2 and 6 are discarded.

[0191] Next, the third data and the corresponding comparison fragment pairs are screened in the same manner as the second data in the array. Specifically, the third data in the array is sorted according to the number of times it appears the most. At this time, the difference between the third data in the array and the element with the most occurrences is less than a threshold value such as 15, and is retained, otherwise it is removed.

[0192] When processing the second data or the third data in the array, if all the values ​​in the second and / or third data do not have the most frequent or the most frequent values, the breakpoint coordinates in one of the alignment fragment pairs can be randomly selected. For example, the breakpoint coordinates in alignment fragment 1 are taken. For example, in the example of alignment fragment pairs 1-5, the right breakpoint coordinates each appear only once, and there is no element with the most frequent occurrences. At this time, the right breakpoint coordinate 500 of alignment fragment pair 1 can be selected. Thus, a list of the left and right breakpoints of the alignment fragment pair supporting the fusion of geneA--geneB after the consistency analysis of the fusion breakpoints is obtained:

[0193] bp_filterlist = [('Comparison fragment pair 1', 200, 500), ('Comparison fragment pair 3', 201, 504), ('Comparison fragment pair 5', 200, 505)]

[0194] Based on the alignment fragment pairs filtered by bp_filterlist, the initial fusion information obtained in the previous step is filtered to retain only the filtered fusion information of the alignment fragment pairs filtered by bp_filterlist: fusion information of alignment fragment pair 1: (geneA--geneB, X, 200, X, 500); fusion information of alignment fragment pair 3: (geneA--geneB, X, 201, X, 504); initial fusion information of alignment fragment pair 5: (geneA--geneB, X, 200, X, 505). The set of filtered fusion information is named fusion_filter_result = [fusion information of alignment fragment pair 1, fusion information of alignment fragment pair 3, fusion information of alignment fragment pair 5].

[0195] 8. Density Cluster Analysis

[0196] All fusion results obtained in step 7 are clustered based on density to obtain highly reliable alignment fragment pairs that support fusion. This process is a key step to further reduce false positives and false negatives.

[0197] When the fused transcript is long, for example, insertions or deletions caused by sequencing may be introduced, resulting in large fluctuations in the position of the fusion breakpoints, which in turn may produce a large number of false positive and false negative results.

[0198] In view of the above situation, the present invention designs a highly sensitive clustering algorithm to more accurately narrow the breakpoint range, which is a key step in accurately predicting fusion genes. The clustering algorithm designed by the present invention does not require the number of clusters to be specified in advance, and the cluster structure is determined by calculating the density around each data point, thereby automatically discovering different clusters and filtering out noise points.

[0199] The clustering algorithm designed by the present invention is obtained by optimizing the DBSCAN algorithm. Firstly, the relevant concepts and definitions in the algorithm of the present invention are explained.

[0200] Object neighborhood Neps(q): The neighborhood within a radius eps (Euclidean distance) of a given object is called the neighborhood of the object.

[0201] Neps(q) is defined as: Neps(q):{p belongs to D|distance(p,q)<=eps}, p, q is a point in D, and D is all the points in the data set.

[0202] distance(p,q) refers to the Euclidean distance from p to q.

[0203] Neps(q) contains two elements: the domain radius eps and the core object of the point in the domain: if the neighborhood of an object contains at least the minimum number MinPts (MinPts is the minimum number of domain points for a given point to become a core object in the domain), then the object is called a core object, |Neps(p)|>=MinPts.

[0204] Boundary object: An object whose neighborhood is less than MinPts objects, but is in the neighborhood of a core object.

[0205] Noise point: The neighborhood of the object is less than MinPts objects and it is not a boundary point.

[0206] Direct density reachability: Given a set of objects D, if p is in the neighborhood of q and q is a core object, then object p is directly density reachable from object q.

[0207] Density reachable: If there exists a chain of objects p1, p2, ..., pn such that for pi belongs to D(1≤i≤n-1), pi+1 is directly density reachable from pi with respect to eps and MinPts, then object pn is density reachable from p1 with respect to eps and MinPts.

[0208] The input parameters of the clustering algorithm of the present invention are: the dataset of the breakpoint coordinates of the left and right sides of the alignment fragment pair, the dataset after depth weighted reconstruction, the domain radius (epslion) calculated based on the weighted reconstruction dataset, and the minimum domain point number k of each cluster. The clustering process is divided into two steps:

[0209] Step 1: Perform deep weighted reconstruction on the obtained fusion gene dataset. Given the example dataset of {'9:130854065_22:23290411':[['130854063','130854063','130854063','130854063','130854063','130854063','130854063','130854063','130854063','130854063','130854063','130854063','130854063','130854063','130854063','130854063','130854063','130854063','130854063','130854063','1 30854063','130854063','130854063','130854063','130854063','130854063','130854063','130854063','130854063','130854063','130854064','130854067','130854067','130854067','130854067','130854067','130854067','130854067','130854067','130854067','130854067','130854067' 4067','130854068','130854068','130854069','130854074','130854075','130854078','130854079'],['23290405','23290408','23290409','23290409','23290410','23290410','23290411','23290411','23290412','23290412','23290412','23290412' ','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412','23290412']],'19:16328114_5:150053307':[['16328110','16328111','16328111','16328111','16328111','16328111','16328111','16328111','16328111',' 16328111','16328111','16328111','16328113','16328113','16328113','16328113','16328113','16328113','16328113','16328114','16328114','16328114','16328114','16328115','16328115','16328115','16328115','16328118','16328118' ,'16328120','16328123','16328123','16328124'],['150053303','150053303','150053304','150053305','150053305','150053305','150053305','150053305','150053305','150053305','150053305','150053305','150053305','150053305','1 50053305','150053305','150053305','150053305','150053308','150053310','150053310','150053310','150053310','150053310','150053310','150053310','150053311','150053311','150053312','150053315','150053315','150053315','150053316','150053316']]},

[0210] The dataset after processing by taking the last 3 digits of the genomic coordinate value is: A =

[0211] {'9:130854065_22:23290411':([[63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,64,67,67,67,67,67,67,67,67,67,67,68,68,69,74,75,78,79] ],[[405,408,409,409,410,410,411,411,411,412,412,412,412,412,412,412,412,412,412,412,412,412,412,412,412,412,412,412,412,412,412,412,412,412,412,412,412,412,412,412,4 12,412,412,412,412]]),'19:16328114_5::150053307':([[110,111,111,111,111,111,111,111,111,111,111,113,113,113,113,113,114,114,114,114,115,115,115,115,1 18,118,120,123,123,124]],[[303,303,304,305,305,305,305,305,305,305,305,305,305,305,305,305,305,305,305,305,305,305,305,305,308,310,310,310,310,310,310,311,311,312,315,315,315,316,316]])}

[0212] Among them, the average depth of 20bp near the breakpoint of the candidate fusion gene "9:130854065_22:23290411" is 7934.9, the number of breakpoints is 45, and the coefficient calculated according to the formula for calculating the weighted coefficient fn is 1.35. The average depth of 20bp near the breakpoint of the candidate fusion gene '19:16328114_5:150053307' is 13498.65, the number of breakpoints is 34, and the coefficient calculated according to the formula for calculating the weighted coefficient fn is 0.79.

[0213] The weighting coefficient formula is:

[0214] Based on the above data set and the weighting coefficient of each fusion gene, the depth-weighted data set obtained is A'={'9:130854065_22:23290411':[[63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,64,67,67,67,67,67,67,67,67,67,68,68 ,69,74,75,78,79,63 ... 12,412,412,412,412,412,412,412,412,412,412,412,412,412,412,405,408,409,409,410,410,411,411,411,412,412,412,412,412,412,412,412]],'19:16328114_5:150053307':[[110,111,111,111,111,111,111,11 The number of breakpoints in the first fusion gene increased to 61, and the number of breakpoints in the second fusion gene increased to 27.

[0215] Based on the breakpoint information of each fusion gene in A', the Euclidean distances between all points in the data set are calculated, and the median of the Euclidean distance values ​​of all sites is taken as the neighbor distance epslion required by the clustering method. In this example, epslion = 4.

[0216] Euclidean distance calculation formula: Among them, x and y correspond to the coordinate values ​​of the upstream and downstream breakpoints of F'n respectively. For example, the first point (63,405) is calculated as x1 and y1, and the second point (63,408) is calculated as x2 and y2.

[0217] Based on the breakpoint information of each fusion gene in A', the parameter k required for the clustering method is calculated. The number of fusion gene breakpoints from F'1 to F'n is sorted from small to large to form a list B = [27,61]. At this time, the number of values ​​in B is less than or equal to three. k is obtained by multiplying the value in the middle of list B by 0.6. The specific position is calculated by int(len(B) / 2)+1, where len represents the length of B and int represents the value to be rounded down. In this example, the number of breakpoint coordinates of the first fusion gene is 61, and the breakpoint coordinates of the second fusion gene are 27, so B = [27,61], the length of B is 2, and the position = int(2 / 2+1), which is the second position. The value of the second position in B is 61, multiplied by 0.6, and rounded down to 36, that is, k = 36.

[0218] Step 2: Clustering is performed based on the k value obtained in the first step and epslion as input. The specific process flow of clustering is as follows:

[0219] 1. If the number of objects contained in Neps(p) is less than k: mark object p as a boundary point or a noise point, indicating that the density around the point is low. Skip processing the object p.

[0220] 2. If the number of objects contained in Neps(p) is greater than or equal to k: mark object p as a core point, indicating that the density around this point is high. Create a new cluster Cluster and add object p to the cluster Cluster.

[0221] 3. For all unprocessed objects q in Neps(p): check the domain Neps(q) of object q. If the number of objects contained in Neps(q) is greater than or equal to k: add the objects in Neps(q) that are not included in any cluster to cluster Cluster, and expand the members of cluster Cluster. Recursively process the domain of object q and continue to explore more objects that are density-reachable from the core point object p.

[0222] 4. When all objects in the dataset have been processed, the clustering process ends. At the same time, cluster labels are added to each point in the dataset, such as Figure 2As shown in the figure, if the value of the cluster label is -1 (white circle), it means that it is not classified, that is, it is a noise point and is filtered. The clusters whose labels are not -1 (black circles) are the sets of credible fusion gene breakpoints. The clustering-based method can further filter the fusion information of the alignment fragments in the fusion_filter_result of the seventh step.

[0223] 9. Determine whether the fusion gene has a complete gene structure, that is, whether the fusion transcript satisfies the 5'-3' connection

[0224] Markers that increase the coding potential of fusion transcripts can provide important evidence for scientific research or clinical research on abnormal protein expression caused by gene fusion leading to cancer. The judgment method is based on the coverage interval of the positive and negative strands of the gene and each read segment on the genome, as well as the information of the breakpoints predicted in the previous step. For the two genes A and B of the fusion gene, the breakpoint coordinates on gene A / B, the positive and negative strand information, and the coverage of the read segments containing the breakpoints of gene A / B are recorded as "gene: (strand information, genome coordinate starting position, genome coordinate ending position, breakpoint position)". For example, A: (+, 100, 1000, 103), A: (+, 105, 800, 103), A: (+, 104, 960, 103), etc.; B: (-, 2000, 3000, 3001), B: (-, 1988, 3001, 3001), B: (-, 2023, 3003, 3001).

[0225] Based on the above information, the gene link structure is determined as follows:

[0226] First, if "|genomic coordinate starting position - breakpoint position|<|genomic coordinate ending position - breakpoint position|, if the chain information is +, the breakpoint is at the 5' end, and the gene in the fusion gene contributes to the 3' end sequence; if the chain information is -, the breakpoint is at the 3' end, and the gene in the fusion gene contributes to the 5' end sequence"; if "|genomic coordinate starting position - breakpoint position|>|genomic coordinate ending position - breakpoint position|, if the chain information is +, the breakpoint is at the 3' end, and the gene in the fusion gene contributes to the 5' end sequence; if the chain information is -, the breakpoint is at the 5' end, and the gene in the fusion gene contributes to the 3' end sequence";

[0227] Second, if gene A and gene B contribute 5' and 3' sequences respectively, they are considered to have coding potential from 5' to 3' ends; otherwise, they are considered to have no coding potential (for example, if both contribute 5' sequences or both contribute 3' sequences, they are considered to have no coding potential). In the final result, they will be marked as having no coding potential to indicate an unreliable fusion. Figure 3-6 Four common types of fusion genes are shown. Figure 3Take fusion type 1 as an example. The breakpoints of gene A and gene B are indicated by vertical arrows, and the fusion direction of gene A and gene B is indicated by horizontal arrows. The fusion transcript is finally generated: A--B. Fusion type 1 includes the complete promoter and start codon region of gene A, and also includes the stop codon and terminator region of gene B. The fusion transcript formed in this way will be marked as having the characteristics of translating fusion protein, and is a candidate fusion type that needs to be retained; fusion type 2 has no stop codon and terminator region, and is marked as the fusion transcript has no function of translating protein; fusion type 3 has neither promoter nor terminator, and is marked as the fusion transcript has no function of translating protein; fusion type 4 has no promoter region and start codon, and is marked as the fusion transcript has no function of translating protein.

[0228] 10. Output fusion gene detection results

[0229] The sequencing data of the targeted captured cDNA library obtained by the carbon nanopore sequencing platform were collated based on the analysis process from the first step to the ninth step, and the exemplary results are shown in Table 2.

[0230] The meaning of each header is as follows: Gene_pair: the name of the fused gene; Upstream_gene: the name of the upstream fusion gene (chain information of the upstream fusion gene); Downstream_gene: the name of the downstream fusion gene (chain information of the downstream fusion gene); Upstream_gene_breakpoint: the location of the breakpoint of the upstream fusion gene (chromosome information of the fusion gene: coordinates of the breakpoint); Downstream_gene_breakpoint: the location of the breakpoint of the downstream fusion gene (chromosome information of the fusion gene: coordinates of the breakpoint); Support_read_number: the number of aligned fragments that support fusion; Fusion_annots: whether the fused gene has coding potential (protein_coding indicates coding potential, no_protein_coding indicates no coding potential).

[0231] Table 2 Fusion gene detection results

[0232]

[0233] Test Case

[0234] The performance of the method of the present invention was compared with the commonly used fusion gene detection tool Fusionseeker. The data set analyzed was 23 positive samples and 9 negative samples. The performance comparison indicators were: sensitivity, specificity, and detection rate.

[0235] Sensitivity calculation formula: (TP / P)*100, specificity calculation formula: (TN / N)*100. TP represents the number of correctly detected positive samples, P represents the number of all positive samples, TN represents the number of correctly detected negative samples, and N represents the number of all negative samples. The detection rate means that as long as the detection software can detect the gene that actually fused in the sample, the detection rate is 100%, and if the gene that actually fused in the sample is not detected, the detection rate is 0.

[0236] The positive sample data set is a targeted captured cDNA library + full-length transcriptome sequencing constructed based on a domestic nanopore sequencing platform. The results of different detection methods for positive samples are shown in Table 3, and the performance comparison of different detection methods for positive samples is shown in Table 4.

[0237] Table 3 Positive sample test results

[0238]

[0239]

[0240] Note: The actual fusion gene is the actual fusion gene of the real sample.

[0241] Table 4 Performance comparison of different detection methods for positive samples

[0242] Detection Methods Sensitivity The present invention 100％ Fusionseeker 34.79％

[0243] It can be seen from Table 3 and Table 4 that the gene fusion detection sensitivity of all samples using the present invention can reach 100%. The gene fusion detection sensitivity of the Fusionseeker software is only 34.79%; the method of the present invention has a higher sensitivity.

[0244] The negative sample data set is a targeted captured cDNA library + full-length transcriptome sequencing constructed based on a domestic nanopore sequencing platform. The detection results of negative samples are shown in Table 5, and the performance comparison of different detection methods for negative samples is shown in Table 6.

[0245] Table 5 Negative sample detection rate results

[0246]

[0247] Table 6 Performance comparison of different detection methods for negative samples

[0248] Detection Methods Specificity The present invention 100％ Fusionseeker 11.11％

[0249] As shown in Table 5 and Table 6, the gene fusion detection specificity of all samples using the present invention can reach 100%. The gene fusion detection specificity using Fusionseeker software is only 11.11%; the method of the present invention has a higher gene fusion detection specificity.

[0250] In short, the fusion gene detection method of the present invention utilizes the fusion breakpoint of the fusion transcript generated by the fusion gene, which should be relatively fixed in theory. The consistency analysis of the fusion breakpoint is first performed to improve the sensitivity and specificity of the detection. After the consistency analysis of the fusion breakpoint is completed, a clustering algorithm is designed. The advantage of the clustering algorithm is that it does not need to pre-specify the number of clusters. The structure of the clustering is determined by calculating the density around each data point, thereby automatically finding different clusters and eliminating false positive fusion noise. The clustering algorithm eliminates the factors of the distance and number of clusters of the subjective design clustering. The fusion gene results obtained after clustering, combined with biological characteristics, select fusion genes with coding potential that meet the 5'-3' connection rule as the final candidate fusion genes, obtain higher detection rates and accuracy relative to the prior art, and have biologically significant fusion gene detection results, which can assist in the selection of clinical targeted drugs, and can meet the needs of scientific research and the vast majority of clinical analysis.

[0251] In addition, the present invention can simultaneously meet the fusion gene detection requirements of full-length transcript sequencing or targeted capture+full-length transcript sequencing based on DNA / cDNA library based on the sequencing platform of Pacbio, Ont, and QNome, including but not limited to DNA library+full-length transcript sequencing constructed based on Qnome (Qi carbon nanopore sequencing platform); DNA library+full-length transcript sequencing constructed based on Qnome; cDNA library+full-length transcript sequencing constructed based on Qnome; cDNA library+full-length transcript sequencing constructed based on Qnome; DNA library / cDNA library+full-length transcript sequencing constructed based on Qnome; DNA library / cDNA library+full-length transcript sequencing of Ont platform; DNA library / cDNA library+full-length transcript sequencing of Pacbio platform. The present invention can also be adapted to annotation files of different source databases, and in addition to gtf, it is also compatible with gff as an input file, overcoming the requirement of annotation file singularity.

[0252] Although the present invention has been described with reference to exemplary embodiments, it should be understood that the present invention is not limited to the disclosed exemplary embodiments. Various adjustments or changes may be made to the exemplary embodiments of the present specification without departing from the scope or spirit of the present invention. The scope of the claims should be based on the broadest interpretation to cover all modifications and equivalent structures and functions.

Claims

1. A method for analyzing samples based on sequence data, characterized in that: include: Performing density-based spatial cluster analysis on a data set obtained based on sequence data from the sample, judging the true situation of the sample based on the clusters obtained by the analysis, and / or judging the untrue situation or abnormal situation of the sample based on the discrete data relative to the clusters obtained by the analysis, The minimum number k of data points contained in the cluster is obtained by calculation.

2. The method according to claim 1, characterized in that The data set is obtained by weighted reconstruction based on depth; And / or, the cluster-dependent distance parameter e is obtained by calculation, Preferably, the cluster-dependent distance parameter e is the median, mean or percentile of the neighbor distance between any two data points in the data set. More preferably, the neighbor distance is the Euclidean distance. And / or, the minimum number k is the median, mean or percentile of the number of supports for each data point in the data set.

3. The method according to claim 1 or 2, characterized in that: The sample is a biological sample, And / or, the sequence data is selected from the sequencing data of any one or more of polynucleotides, polypeptides, polysaccharides, and lipids, preferably polynucleotides, And / or, the density-based spatial clustering is selected from DBSCAN, OPTICS, or DENCLUE, preferably DBSCAN, and / or, the data set obtained based on the sequence data from the sample is a data set that has been subjected to normalization and / or denoising processing by non-density-based clustering analysis, And / or, each data point in the data set obtained based on the sequence data from the sample includes at least two dimensions, preferably, the at least two dimensions include a time dimension and / or a space dimension.

4. The method according to claim 3, characterized in that The true situation of the sample includes the true situation of gene fusion in the biological sample, And / or, each data point in the data set obtained based on the sequence data from the sample is the breakpoint position coordinates x of two fragments of the same sequencing read of the sample in two different genes corresponding to the same reference sequence or reference genome. i and Y j .

5. The method according to claim 4, characterized in that The method for obtaining a data set based on sequence data from a sample comprises: Aligning the reads in the sequencing data from the sample to a reference sequence or a reference genome, and obtaining aligned fragments, wherein the aligned fragments are portions of the continuous sequence regions in the reads that can be aligned to the reference sequence or the reference genome; Screening according to the alignment lengths of the alignment fragments to obtain a set of candidate alignment fragments; Determine the gene information corresponding to each alignment fragment in the candidate alignment fragment set according to the reference sequence or the reference genome annotation file, and combine two alignment fragments from the same read segment and corresponding to different genes into an alignment fragment pair, thereby obtaining initial fusion information, which includes the gene information and breakpoint coordinate information corresponding to each of the first alignment fragment and the second alignment fragment; Performing consistency analysis on the fusion breakpoints in the initial fusion information according to the breakpoint coordinate information, excluding inconsistent alignment fragment pairs, obtaining an alignment fragment set after filtering the fusion information, and obtaining the data set; The data set is weightedly reconstructed based on the total sequencing depth at the breakpoint position, which includes: calculating a weight coefficient fn based on the larger sequencing depth Dn near the two breakpoints of each aligned fragment pair, and scaling the coordinate matrix of the fusion information of each aligned fragment pair based on fn, wherein: The stretching is to multiply fn by the coordinate matrix of the fusion information of the corresponding aligned fragment pairs; Preferably, the coordinate matrix of the fusion information of each aligned fragment pair is simplified, and more preferably, at least the last three digits of the genomic coordinate values ​​of the coordinates are retained after the simplification.

6. A method for analyzing fusion genes based on sequencing data, comprising: A cluster analysis is performed on a data set obtained based on sequencing data from a sample, and the true and / or false or abnormal gene fusion in the sample is determined based on the analysis, wherein: Each data point in the data set is the breakpoint position coordinates X of two fragments of the same sequencing read of the sample in two different genes corresponding to the same reference sequence or reference genome. i and Y j ; The data set is obtained by weighted reconstruction based on the total sequencing depth at the breakpoint position.

7. The method according to claim 6, characterized in that Aligning the reads in the sequencing data from the sample to a reference sequence or a reference genome, and obtaining aligned fragments, wherein the aligned fragments are portions of the reads that can be aligned to a continuous sequence region of the reference sequence or the reference genome, Screening according to the alignment lengths of the alignment fragments to obtain a set of candidate alignment fragments; Determine the gene information corresponding to each alignment fragment in the candidate alignment fragment set according to the reference sequence or the reference genome annotation file, and combine two alignment fragments from the same read segment and corresponding to different genes into an alignment fragment pair, thereby obtaining initial fusion information, which includes the gene information and breakpoint coordinate information corresponding to each of the first alignment fragment and the second alignment fragment; The fusion breakpoints in the initial fusion information are analyzed for consistency according to the breakpoint coordinate information, and inconsistent alignment fragment pairs are excluded to obtain an alignment fragment set after filtering the fusion information, thereby obtaining the data set.

8. The method according to claim 7, characterized in that The data set is obtained by weighted reconstruction based on the total sequencing depth at the breakpoint position, including: Based on the larger sequencing depth Dn near the two breakpoints of each aligned fragment pair, the weighting coefficient fn is calculated, and the coordinate matrix of the fusion information of each aligned fragment pair is scaled based on fn, where: The stretching is to multiply fn by the coordinate matrix of the fusion information of the corresponding aligned fragment pairs; Preferably, the coordinate matrix of the fusion information of each aligned fragment pair is simplified, and more preferably, at least the last three digits of the genomic coordinate values ​​of the coordinates are retained after the simplification.

9. The method according to claim 6, characterized in that The cluster analysis adopts density-based spatial cluster analysis; Preferably, the true situation of the sample is judged based on the clusters obtained by the analysis, and / or the untrue situation or abnormal situation of the sample is judged based on the discrete data relative to the clusters obtained by the analysis. Preferably, the density-based spatial clustering is selected from DBSCAN, OPTICS, or DENCLUE, more preferably, DBSCAN.

10. The method according to claim 5 or 7, characterized in that: The obtaining of the alignment fragments also includes the step of filtering and removing the alignment fragments with multiple alignments, And / or, the screening of the candidate alignment fragment set according to the alignment length of the alignment fragment includes, after preliminarily confirming the alignment fragment corresponding to the potential fusion gene according to the SA tag, further screening the candidate alignment fragment set according to the alignment length of the alignment fragment, and / or, the calculation of the alignment length in the candidate alignment segment set obtained by screening according to the alignment length of the alignment segment includes analyzing the sequence in the alignment segment according to the CIGAR information to obtain a string containing at least one of M, I, D, N, H and S, counting the number of M in the string in order until H or S appears in the string, and terminating the counting, and calculating the sum of the lengths corresponding to each M as the alignment length, and / or, the screening according to the comparison length of the comparison fragments to obtain the candidate comparison fragment set includes removing the comparison fragments whose comparison length is less than a threshold, and simultaneously removing the comparison fragments whose coordinates of the sequences of the M region in the first comparison fragment and the second comparison fragment overlap by a length greater than a threshold, and / or, when |the minimum value between the ending coordinate position of gene X and the ending coordinate position of the M region of the aligned fragment|-|the maximum value between the starting coordinate position of gene X and the starting coordinate position of the M region of the aligned fragment| is greater than a threshold value, the aligned fragment is determined to correspond to gene X; And / or, performing consistency analysis on the fusion breakpoints in the initial fusion information based on the breakpoint coordinate information includes counting the distribution of the breakpoint coordinates of the first aligned fragment in each aligned fragment pair and the breakpoint coordinates of the second aligned fragment in each aligned fragment pair among multiple aligned fragment pairs supporting the same fusion, and removing the aligned fragments corresponding to the breakpoint coordinates that deviate from the threshold.

11. The method according to claim 4 or 6, characterized in that: It also includes a step of determining whether the fusion gene has a complete gene structure. If so, the result of the actual situation of the gene fusion in the sample is credible; if not, the result of the actual situation of the gene fusion in the sample is not credible.

12. A device for analyzing a sample based on sequence data, characterized in that include: a data acquisition unit configured to acquire sequence data; A data processing unit, configured to be able to perform the method according to any one of claims 1 to 11.

13. A device for analyzing a sample based on sequence data, characterized in that The method comprises a memory, a processor, and a program stored in the memory and capable of running on the processor, wherein the program is configured to implement the method according to any one of claims 1 to 11.

14. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a program for analyzing samples based on sequence data, and when the program is executed by a processor, the method according to any one of claims 1 to 11 is implemented.

Images