A fusion gene detection method for long-read RNA-seq sequencing data

Abstract

A fusion gene detection method for long-read RNA-seq sequencing data, the present invention relates to a fusion gene detection method for long-read RNA-seq sequencing data. The purpose of the present invention is to solve the problems of excessive computing resource usage and low accuracy of identification results in the existing method for identifying fusion genes in a variety of human cancer sequencing data. The process is: one: generating recombinant reference transcripts containing only transcript sequences and recombinant transcript annotations; two: obtaining filtered sequencing fragments; three: searching for fusion gene pairs; four: obtaining each fusion gene pair after breakpoint restoration; hierarchical clustering of the sequencing fragments of each fusion gene pair after breakpoint restoration according to the breakpoint position, each clustering cluster is a fusion gene pair; five: obtaining RNA-seq comparison results; six: according to the RNA-seq comparison results, screening the four obtained fusion gene pairs to obtain the final fusion gene pairs. The present invention is used in the field of fusion gene screening.

Description

Technical Field

[0001] The present invention relates to a fusion gene detection method for long-read RNA-seq sequencing data. Background Art

[0002] Gene fusion generally refers to the process in which two or more different genes are connected due to some reason (such as chromosomal structural variation) to form a new gene. This gene is usually called a fusion gene, and the connection point between these different genes is generally called a breakpoint. Generally speaking, fusion may change the function of the protein by encoding it, thereby stimulating tumor development. Past medical research has shown that gene fusion drives 16.5% of human cancers and plays a sole driving role in more than 1% of cancers. At the gene sequence level, fusion can cause cell cancer by changing the expression of tumor suppressors or proto-oncogenes; at the gene expression level, fusion can also cause protein function to change through the expression of fusion proteins, thereby leading to the occurrence of cancer.

[0003] Over the past 20 years, next-generation sequencing (NGS) technology has significantly enhanced the scalability and cost-effectiveness of RNA-seq experiments. However, limited read lengths and systematic biases have hindered accurate alignment of sequence fragments and precise genotyping. However, due to the limited sequencing length, NGS data struggle to cover variants occurring within transcripts, making it difficult to determine breakpoint locations during fusion gene detection.

[0004] In recent years, long-read sequencing technologies such as Pacific Biosciences (PacBio) single-molecule real-time sequencing and Oxford Nanopore Technologies (ONT) sequencing have provided longer sequencing fragment lengths, capable of generating full-length transcripts and addressing the challenges posed by alternative splicing events and complex transcript structures. Compared to second-generation sequencing technologies, sequencing fragments generated by third-generation sequencing technologies have longer read lengths, the longest of which can even exceed 1Mbp, and the average length can reach 10kbp, so they can cover very long genomic regions. However, compared to the high precision of second-generation sequencing, the error rate of third-generation sequencing technology is higher. Based on the characteristics of long-read sequencing data that can cover full-length transcripts, fusion gene detection becomes easier.

[0005] For long-read sequencing data, the reference sequence-based alignment method has more advantages than the direct clustering method. Therefore, the method based on the alignment results on the reference genome is the mainstream design of the current fusion gene detection. The existing fusion gene detection tools based on long-read RNA-seq mainly include LongGF, Genion, FusionSeeker, and JAFFAL. Among them, except for JAFFAL, which is based on the alignment on the transcripts for filtering, the other tools are based on the alignment based on the reference genome for filtering. The general background technology route is as follows: Figure 2 shown.

[0006] Compared to DNA alignments against reference transcripts, RNA alignments against reference sequences do not consider intronic regions, potentially leading to higher false positive rates. Furthermore, compared to whole-genome alignments, transcript alignments consume fewer resources and can precisely identify transcripts corresponding to fusion genes. Therefore, methods based on transcript alignment results filtering are more advantageous. However, JAFFAL can sometimes be overly sensitive, resulting in a high number of false positives. Summary of the Invention

[0007] The purpose of this invention is to solve the problems of excessive computational resources occupied and low accuracy of identification results in the existing methods for identifying fusion genes in various human cancer sequencing data, and to propose a fusion gene detection method for long-read RNA-seq sequencing data.

[0008] A fusion gene detection method for long-read RNA-seq sequencing data. The specific process is as follows:

[0009] Step 1: Recombining the reference sequence based on the gene annotation file to generate a recombinant reference transcript containing only the transcript sequence and the recombinant transcript annotation;

[0010] Step 2: Use the DNA alignment mode of the open-source alignment tool minimap2 to align the sequence fragments to be tested to the recombinant reference transcripts, filter the sequence fragments to be tested on the alignment, and obtain the filtered sequence fragments;

[0011] The genes corresponding to the filtered sequencing fragments constitute the gene set G;

[0012] The filtered sequencing fragments contain a sequence identifier flag;

[0013] Step 3: Search for fusion gene pairs in the filtered sequencing fragments;

[0014] Step 4: Perform breakpoint restoration on each fusion gene pair to obtain each fusion gene pair after breakpoint restoration;

[0015] After breakpoint restoration, the sequencing fragments of each fusion gene pair are hierarchically clustered according to the breakpoint position, and clusters with the number of sequencing fragments greater than or equal to two are found and retained. Each cluster is a fusion gene pair;

[0016] Step 5: Use the open source sequence alignment tool abpoa to align the sequencing fragments in each cluster and output the consensus sequence;

[0017] Use the RNA-seq alignment function of the open source alignment tool minimap2 to align the consensus sequence to the reference sequence in step 1 to obtain the RNA-seq alignment results;

[0018] Step 6: Based on the RNA-seq comparison results, the fusion gene pairs obtained in step 4 are screened to obtain the final fusion gene pairs.

[0019] Preferably, in step 1, the reference sequence is recombined based on the gene annotation file to generate a recombinant reference transcript containing only the transcript sequence and the recombinant transcript annotation;

[0020] The specific process is:

[0021] Step 1: Read the gene annotation file and retain the information of the transcripts and the exons corresponding to the retained transcripts in the gene annotation file as required;

[0022] Step 12: Determine the sequence of each exon based on the start and end positions of each exon information retained in step 11 and the reference sequence, and link the sequences of all exons to obtain the transcript sequence;

[0023] Step 13: Connect the transcript sequences according to chromosomes to obtain recombinant reference transcripts;

[0024] Record the start and end positions on the reference sequence and the start and end positions of the recombinant reference transcript to obtain the recombinant transcript annotation;

[0025] Step 14: Use the index function of minimap2 to create an alignment index between the recombinant reference transcript and the reference sequence.

[0026] Preferably, in step 2, the sequenced fragments to be detected on the alignment are filtered to obtain filtered sequenced fragments;

[0027] The specific process is:

[0028] Step 21: Filter according to the sequence length of the sequencing fragments;

[0029] Step 22: Filter based on the number of genes covered by the sequencing fragments;

[0030] Step 23: Filter according to the boundaries of each sequencing fragment retained in step 22 and the boundaries of the exons of the corresponding recombinant reference transcript to obtain filtered sequencing fragments.

[0031] Preferably, in step 21, filtering is performed according to the sequence length of the sequencing fragments;

[0032] The specific process is;

[0033] The sequence fragments to be detected in the alignment whose sequence length is greater than the length threshold are retained;

[0034] The sequencing fragments to be detected in the alignment whose sequence length is less than or equal to the length threshold are deleted.

[0035] Preferably, in step 22, filtering is performed according to the number of genes covered by the sequencing fragments;

[0036] The specific process is:

[0037] Based on the annotation of the recombinant transcripts, determine the gene corresponding to the sequence information of each sequencing fragment retained in step 2-1;

[0038] If the number of genes corresponding to the sequence information of the sequencing fragment is greater than or equal to two, the sequencing fragment is retained;

[0039] If the number of genes corresponding to the sequence information of the sequencing fragment is less than two, the sequencing fragment is deleted.

[0040] Preferably, in step 23, filtering is performed based on the boundaries of each sequencing fragment retained in step 22 and the boundaries of the exons of the corresponding recombinant reference transcript to obtain filtered sequencing fragments;

[0041] The specific process is:

[0042] 1) Calculate the distance between the start position of a sequencing fragment A retained in step 22 and the start position of the exon of the corresponding recombinant reference transcript;

[0043] 2) Calculate the distance between the starting position of a sequencing fragment A retained in step 22 and the ending position of the exon of the corresponding recombinant reference transcript;

[0044] 3) Calculate the distance between the end position of a sequencing fragment A retained in step 22 and the start position of the exon of the corresponding recombinant reference transcript;

[0045] 4) Calculate the distance between the end position of a sequencing fragment A retained in step 22 and the end position of the exon of the corresponding recombinant reference transcript;

[0046] If any of the distances in 1), 2), 3), and 4) is less than the distance threshold, sequencing fragment A is retained; otherwise, sequencing fragment A is deleted;

[0047] This process continues until all the sequenced fragments retained in step 22 are determined to obtain the final sequenced fragments retained.

[0048] Preferably, in step 3, searching for fusion gene pairs in the filtered sequencing fragments obtained in step 2;

[0049] The specific process is:

[0050] Step 31: Assume that the number of filtered sequencing fragments obtained in step 2 is m;

[0051] After determining the gene set G, each gene in the gene set is represented as an m-dimensional vector;

[0052] gene α =(q α1 ,q α2 ,…,q αm ),where q αβ =0or 1,α∈[1,n],β∈[1,m]

[0053] Among them, gene α represents the αth gene;

[0054] q α1 represents the correlation coefficient of the αth gene relative to the filtered sequencing fragment 1 obtained in step 2. The filtered sequencing fragment 1 obtained in step 2 is aligned to the αth gene, q α1 The value is 1; otherwise, q α1 The value is 0;

[0055] q α2 represents the correlation coefficient of the αth gene relative to the filtered sequencing fragment 2 obtained in step 2. The filtered sequencing fragment 2 obtained in step 2 is aligned to the αth gene, q α2 The value is 1; otherwise, q α2 The value is 0;

[0056] q αm represents the correlation coefficient of the αth gene relative to the filtered sequencing fragment m obtained in step 2. The filtered sequencing fragment m obtained in step 2 is aligned to the αth gene, q αm The value is 1; otherwise, q αm The value is 0;

[0057] β represents the sequence number of the filtered sequencing fragment obtained in step 2;

[0058] n represents the total number of genes in the gene set G;

[0059] m represents the total number of sequencing fragments after filtering obtained in step 2;

[0060] Step 32: Get an n×m dimensional matrix Q based on n genes n×m ;

[0061]

[0062] Among them, Q n×m Represents a matrix of n×m dimensions;

[0063] gene1 indicates the first gene, gene2 indicates the second gene, gene n represents the nth gene;

[0064] q 11 It represents the correlation coefficient of the first gene relative to the filtered sequencing fragment 1 obtained in step 2. The filtered sequencing fragment 1 obtained in step 2 is aligned to the first gene, q 11 The value is 1; otherwise, q 11 The value is 0;

[0065] q 12 represents the correlation coefficient of the first gene relative to the filtered sequencing fragment 2 obtained in step 2. The filtered sequencing fragment 2 obtained in step 2 is aligned to the first gene, q 12 The value is 1; otherwise, q 12 The value is 0;

[0066] q 1m represents the correlation coefficient of the first gene relative to the filtered sequencing fragment m obtained in step 2. The filtered sequencing fragment m obtained in step 2 is aligned to the first gene, q 1m The value is 1; otherwise, q 1m The value is 0;

[0067] q 21 It represents the correlation coefficient of the second gene relative to the filtered sequencing fragment 1 obtained in step 2. The filtered sequencing fragment 1 obtained in step 2 is aligned to the second gene, q 21 The value is 1; otherwise, q 21 The value is 0;

[0068] q 22 represents the correlation coefficient of the second gene relative to the filtered sequencing fragment 2 obtained in step 2. The filtered sequencing fragment 2 obtained in step 2 is aligned to the second gene, q 22 The value is 1; otherwise, q 22 The value is 0;

[0069] q 2mrepresents the correlation coefficient of the second gene relative to the filtered sequencing fragment m obtained in step 2. The filtered sequencing fragment m obtained in step 2 is aligned to the second gene, q 2m The value is 1; otherwise, q 2m The value is 0;

[0070] q n1 represents the correlation coefficient of the nth gene relative to the filtered sequencing fragment 1 obtained in step 2. The filtered sequencing fragment 1 obtained in step 2 is aligned to the nth gene, q n1 The value is 1; otherwise, q n1 The value is 0;

[0071] q n2 represents the correlation coefficient of the nth gene relative to the filtered sequencing fragment 2 obtained in step 2. The filtered sequencing fragment 2 obtained in step 2 is aligned to the nth gene, q n2 The value is 1; otherwise, q n2 The value is 0;

[0072] q nm represents the correlation coefficient of the nth gene relative to the filtered sequencing fragment m obtained in step 2. The filtered sequencing fragment m obtained in step 2 is aligned to the nth gene, q nm The value is 1; otherwise, q nm The value is 0;

[0073] Step 3: Convert the n×m dimensional matrix Q n×m With Q n×m Transpose Multiply them together to get an n×n dimensional symmetric matrix R n×n ;

[0074]

[0075]

[0076] Among them, R n×n represents an n×n dimensional symmetric matrix; r αα′ Represents the symmetric matrix R n×n Items in

[0077] q αβ Denotes the matrix Q n×m The correlation coefficient of the αth gene in step 2 relative to the filtered sequencing fragment β obtained in step 2, the filtered sequencing fragment β obtained in step 2 is aligned to the αth gene, q αβ The value is 1; otherwise, q αβ The value is 0;

[0078] q α′β Represents Q n×mThe correlation coefficient of the α′th gene in step 2 relative to the filtered sequencing fragment β obtained in step 2, the filtered sequencing fragment β obtained in step 2 is aligned to the α′th gene, q α′β The value is 1; otherwise, q α′β The value is 0; α′∈[1,n];

[0079] q′ βα′ express The correlation coefficient of the α′th gene in step 2 relative to the filtered sequencing fragment β obtained in step 2, the filtered sequencing fragment β obtained in step 2 is aligned to the α′th gene, q′ βα′ The value is 1; otherwise q′ βα′ The value is 0; α′∈[1,n]; β∈[1,m];

[0080] Steps 3 and 4: Symmetric matrix R n×n Each non-diagonal term r in αα′ Represents the number of sequencing fragments supporting the existence of fusion between gene α and gene α′;

[0081] Traverse the symmetric matrix R n×n All fusion gene pairs with fusion can be found.

[0082] Preferably, in step 4, each fusion gene pair is subjected to breakpoint restoration processing to obtain each fusion gene pair after breakpoint restoration;

[0083] After breakpoint restoration, the sequencing fragments of each fusion gene pair are hierarchically clustered according to the breakpoint position, and clusters with the number of sequencing fragments greater than or equal to two are found and retained. Each cluster is a fusion gene pair;

[0084] The specific process is:

[0085] Step 41: Perform breakpoint restoration on each fusion gene pair to obtain each fusion gene pair after breakpoint restoration; the specific process is as follows:

[0086] There are two types of breakpoints: back-end breakpoints and front-end breakpoints.

[0087] Assume that gene g1 is in front and gene g2 is in the back;

[0088] The flag = 0 of the sequencing fragment corresponding to gene g1 is the back-end breakpoint;

[0089] Flag = 16 of the sequencing fragment corresponding to gene g2 is the back-end breakpoint;

[0090] Flag = 16 of the sequencing fragment corresponding to gene g1 is the front breakpoint;

[0091] The flag = 0 of the sequencing fragment corresponding to gene g2 is the front breakpoint;

[0092] Restore the backend breakpoint; the specific process is:

[0093] By subtracting the breakpoint position of the sequenced fragment filtered in step 2 from the end position of the exon of the reconstructed reference transcript in step 1;

[0094] If the absolute value of the difference is less than the threshold, the breakpoint is considered to belong to the exon of the corresponding recombinant reference transcript in step 1, and the breakpoint position is restored to the reference sequence; the expression is:

[0095] POS ref =pos trans -exon trans +exon ref

[0096] Among them, pos trans Indicates the breakpoint position of the sequencing fragment after filtering in step 2; exon trans Indicates the end position of the exon corresponding to the recombinant reference transcript in step 1, exon ref Indicates the end position of the exon in the reference sequence, pos ref Indicates the breakpoint position on the reference sequence;

[0097] If the absolute value of the difference is greater than or equal to the threshold, the breakpoint is considered not to belong to the exon of the corresponding recombinant reference transcript in step 1, and the corresponding fusion gene pair is filtered;

[0098] The specific process of restoring the front-end breakpoint is as follows:

[0099] By subtracting the breakpoint position of the sequenced fragment filtered in step 2 from the start position of the exon of the reconstructed reference transcript in step 1;

[0100] If the absolute value of the difference is less than the threshold, the breakpoint is considered to belong to the exon of the corresponding recombinant reference transcript in step 1, and the breakpoint position is restored to the reference sequence; the expression is:

[0101] POS ref =pos trans -exon′ trans +exon′ ref

[0102] Among them, pos trans Indicates the breakpoint position of the sequencing fragment after filtering in step 2; exon′ trans Indicates the starting position of the exon corresponding to the recombinant reference transcript in step 1, exon′ ref Indicates the starting position of the exon in the reference sequence, pos ref Indicates the breakpoint position on the reference sequence;

[0103] If the absolute value of the difference is greater than or equal to the threshold, the breakpoint is considered not to belong to the exon of the corresponding recombinant reference transcript in step 1, and the corresponding fusion gene pair is filtered;

[0104] Step 42: After breakpoint restoration, the sequencing fragments of each fusion gene pair are hierarchically clustered according to the breakpoint position, and clusters with the number of sequencing fragments greater than or equal to two are found and retained. Each cluster is a fusion gene pair.

[0105] Preferably, in step 42, hierarchical clustering is performed on the sequencing fragments of each fusion gene pair after breakpoint restoration according to the breakpoint positions, and clusters with the number of sequencing fragments greater than or equal to two are found and retained, each cluster being a fusion gene pair;

[0106] The specific process is:

[0107] Calculate the Manhattan distance between the breakpoint positions of the sequenced fragment genes, cluster the sequenced fragments according to the Manhattan distance, obtain clusters, and find and retain clusters with more than or equal to two sequenced fragments;

[0108] The sequencing fragments are clustered according to the Manhattan distance to obtain clusters; the process is:

[0109] Add the sequencing fragments whose distance from the cluster center is less than N to the cluster and recalculate the cluster center until all sequencing fragments are clustered;

[0110] Manhattan distance = |x c -x d |+|y c -y d |

[0111] in,

[0112] x c Indicates the position of the breakpoint of gene g1 of the first sequenced fragment; y c Indicates the position of the breakpoint of gene g2 of the first sequenced fragment;

[0113] x d Indicates the position of the breakpoint of gene g1 in the second sequenced fragment, y d Indicates the position of the breakpoint of gene g2 in the second sequenced fragment.

[0114] Preferably, in step 6, the fusion gene pairs obtained in step 4 are screened according to the RNA-seq comparison results to obtain the final fusion gene pairs;

[0115] The specific process is:

[0116] If the RNA-seq alignment result obtained in step 5 fails to map to the reference sequence in step 1, or maps to the reference sequence in step 1 but the alignment record contains only one gene, there is no split fragment signal and it is considered a false fusion gene. The fusion gene pair obtained in step 4 is filtered;

[0117] If the RNA-seq alignment result obtained in step 5 is mapped to the reference sequence in step 1 and the alignment record contains two or more genes, there is a split fragment signal. However, if the position difference between the breakpoint corresponding to the RNA-seq alignment result obtained in step 5 and the center of the cluster corresponding to step 4 exceeds the threshold, it is considered to be a possible false fusion gene, and the fusion gene pair obtained in step 4 is retained.

[0118] If the RNA-seq alignment result obtained in step 5 is mapped to the reference sequence in step 1 and the alignment record contains two or more genes, there is a split fragment signal, and the position difference between the breakpoint corresponding to the RNA-seq alignment result obtained in step 5 and the center of the cluster corresponding to step 4 is less than or equal to the threshold, then it is determined to be a true fusion gene, and the fusion gene pair obtained in step 4 is retained.

[0119] The beneficial effects of the present invention are:

[0120] By aligning sequenced fragments to pre-prepared recombinant reference transcripts, this method significantly reduces computational resource consumption for fusion gene detection, improves runtime speed, and reduces memory usage. Furthermore, through clustering and multiple sequence alignment, the method successfully characterizes fusion gene pairs, reducing computational resource consumption while improving detection accuracy.

[0121] The present invention determines the confidence level of fusion gene detection results through a retesting method. This ensures that while maintaining high accuracy for high-confidence fusion gene detection results, the present invention can also guarantee sensitivity for lower-confidence detection results. By analyzing the possible fusion gene scenarios, the present invention can accurately calculate the fusion gene breakpoint locations, avoiding some false positives. The present invention introduces multi-threaded parallel computing during the detection process, reducing detection time and overhead. BRIEF DESCRIPTION OF THE DRAWINGS

[0122] Figure 1 Flowchart of the method of the present invention, SA represents split fragment;

[0123] Figure 2 The figure is a flow chart of the existing method. DETAILED DESCRIPTION

[0124] Specific implementation method 1: Combination Figure 1 This embodiment describes a method for detecting fusion genes based on long-read RNA-seq sequencing data. The specific process is as follows:

[0125] Step 1: Reconstruct the reference sequence based on the gene annotation file (gene annotation file and reference sequence https: / / www.gencodegenes.org / human / ) to generate a reconstructed reference transcript containing only the transcript sequence and the reconstructed transcript annotation;

[0126] Step 2: Use the DNA alignment mode of the open-source alignment tool minimap2 to align the sequence fragments to be tested to the recombinant reference transcripts, filter the sequence fragments to be tested on the alignment, and obtain the filtered sequence fragments;

[0127] The genes corresponding to the filtered sequencing fragments constitute the gene set G;

[0128] The filtered sequencing fragments contain a sequence identifier flag;

[0129] Step 3: Search for fusion gene pairs in the filtered sequencing fragments;

[0130] Step 4: Perform breakpoint restoration on each fusion gene pair to obtain each fusion gene pair after breakpoint restoration;

[0131] After breakpoint restoration, the sequencing fragments of each fusion gene pair are hierarchically clustered from bottom to top according to the breakpoint position, and clusters with the number of sequencing fragments greater than or equal to two are found and retained. Each cluster is a fusion gene pair;

[0132] Step 5: Use the open source multiple sequence alignment tool abpoa to perform multiple sequence alignment on the sequencing fragments in each cluster and output the consensus sequence (a sequence that reflects the common characteristics of the sequencing fragment sequences in the cluster by combining the most identical parts of the sequencing fragment sequences in the cluster);

[0133] Use the RNA-seq alignment function of the open source alignment tool minimap2 to align the consensus sequence to the reference sequence in step 1 to obtain the RNA-seq alignment results;

[0134] Step 6: Based on the RNA-seq comparison results, the fusion gene pairs obtained in step 4 are screened to obtain the final fusion gene pairs.

[0135] Specific embodiment 2: This embodiment differs from specific embodiment 1 in that in step 1, the reference sequence is recombined based on the gene annotation file to generate a recombinant reference transcript containing only the transcript sequence and the recombinant transcript annotation;

[0136] The specific process is:

[0137] Step 1: Read the gene annotation file and retain the specific types of transcripts in the gene annotation file (by default, only protein-coding transcripts are retained) and the exon information corresponding to the retained transcripts as required;

[0138] Step 12: Determine the sequence of each exon based on the start and end positions of each exon information retained in step 11 and the reference sequence, and link the sequences of all exons to obtain the transcript sequence;

[0139] Step 13: Connect the transcript sequences according to chromosomes to obtain recombinant reference transcripts;

[0140] Record the start and end positions on the reference sequence and the start and end positions of the recombinant reference transcript to obtain the recombinant transcript annotation;

[0141] Step 14: Use the index function of minimap2 to create an alignment index between the recombinant reference transcript and the reference sequence.

[0142] Other steps and parameters are the same as those in the first embodiment.

[0143] Specific embodiment three: This embodiment differs from specific embodiment one or two in that, in step two, the sequenced fragments to be detected on the alignment are filtered to obtain filtered sequenced fragments;

[0144] The specific process is:

[0145] Step 21: Filter according to the sequence length of the sequencing fragments;

[0146] Step 22: Filter based on the number of genes covered by the sequencing fragments;

[0147] Step 23: Filter according to the boundaries of each sequencing fragment retained in step 22 and the boundaries of the exons of the corresponding recombinant reference transcript to obtain filtered sequencing fragments.

[0148] Other steps and parameters are the same as those in the first or second embodiment.

[0149] Specific embodiment 4: This embodiment differs from specific embodiments 1 to 3 in that in step 21, filtering is performed based on the sequence length of the sequencing fragments; the specific process is as follows;

[0150] The sequence fragments to be detected in the alignment whose sequence length is greater than the length threshold are retained;

[0151] The sequencing fragments to be detected in the alignment whose sequence length is less than or equal to the length threshold are deleted.

[0152] The default length threshold is 50 bp.

[0153] The other steps and parameters are the same as those in the first to third embodiments.

[0154] Specific embodiment 5: This embodiment differs from specific embodiments 1 to 4 in that in step 22, filtering is performed based on the number of genes covered by the sequencing fragments; the specific process is as follows:

[0155] Based on the annotation of the recombinant transcripts, determine the gene corresponding to the sequence information of each sequencing fragment retained in step 2-1;

[0156] If the number of genes corresponding to the sequence information of the sequencing fragment is greater than or equal to two, the sequencing fragment is retained;

[0157] If the number of genes corresponding to the sequence information of the sequencing fragment is less than two, the sequencing fragment is deleted.

[0158] The sequence information of the sequencing fragment includes the starting position, ending position, insertion position, deletion position, and alignment position.

[0159] The other steps and parameters are the same as those in the first to fourth embodiments.

[0160] Specific embodiment 6: This embodiment differs from any one of specific embodiments 1 to 5 in that, in step 23, filtering is performed based on the boundaries of each sequencing fragment retained in step 22 and the boundaries of the exons of the corresponding recombinant reference transcript to obtain filtered sequencing fragments; the specific process is as follows:

[0161] 1) Calculate the distance between the start position of a sequencing fragment A retained in step 22 and the start position of the exon of the corresponding recombinant reference transcript;

[0162] 2) Calculate the distance between the starting position of a sequencing fragment A retained in step 22 and the ending position of the exon of the corresponding recombinant reference transcript;

[0163] 3) Calculate the distance between the end position of a sequencing fragment A retained in step 22 and the start position of the exon of the corresponding recombinant reference transcript;

[0164] 4) Calculate the distance between the end position of a sequencing fragment A retained in step 22 and the end position of the exon of the corresponding recombinant reference transcript;

[0165] If any of the distances in 1), 2), 3), and 4) is less than the distance threshold, sequencing fragment A is retained; otherwise, sequencing fragment A is deleted;

[0166] This process continues until all the sequenced fragments retained in step 22 are determined to obtain the final sequenced fragments retained.

[0167] The distance threshold is 30 bp.

[0168] After these filters, useless sequencing fragments in subsequent detection can be minimized.

[0169] The filtered sequencing fragments all correspond to at least two genes.

[0170] The other steps and parameters are the same as those in the first to fifth embodiments.

[0171] Specific embodiment seven: This embodiment differs from any one of specific embodiments one to six in that, in step three, a fusion gene pair is searched for in the filtered sequencing fragments obtained in step two; the specific process is as follows:

[0172] Step 31: Assume that the number of filtered sequencing fragments obtained in step 2 is m;

[0173] After determining the gene set G, each gene in the gene set is represented as an m-dimensional vector;

[0174] gene α =(q α1 ,q α2 ,…,q αm ),where q αβ =0or 1,α∈[1,n],β∈[1,m]

[0175] Among them, gene α represents the αth gene;

[0176] q α1 represents the correlation coefficient of the αth gene relative to the filtered sequencing fragment 1 obtained in step 2. The filtered sequencing fragment 1 obtained in step 2 is aligned to the αth gene, q α1 The value is 1; otherwise, q α1 The value is 0;

[0177] q α2 represents the correlation coefficient of the αth gene relative to the filtered sequencing fragment 2 obtained in step 2. The filtered sequencing fragment 2 obtained in step 2 is aligned to the αth gene, q α2 The value is 1; otherwise, q α2 The value is 0;

[0178] q αm represents the correlation coefficient of the αth gene relative to the filtered sequencing fragment m obtained in step 2. The filtered sequencing fragment m obtained in step 2 is aligned to the αth gene, q αm The value is 1; otherwise, q αm The value is 0;

[0179] β represents the sequence number of the filtered sequencing fragment obtained in step 2;

[0180] n represents the total number of genes in the gene set G;

[0181] m represents the total number of sequencing fragments after filtering obtained in step 2;

[0182] Step 32: Get an n×m dimensional matrix Q based on n genes n×m ;

[0183]

[0184] Among them, Q n×m Represents a matrix of n×m dimensions;

[0185] gene1 indicates the first gene, gene2 indicates the second gene, gene n represents the nth gene;

[0186] q 11 It represents the correlation coefficient of the first gene relative to the filtered sequencing fragment 1 obtained in step 2. The filtered sequencing fragment 1 obtained in step 2 is aligned to the first gene, q 11 The value is 1; otherwise, q 11 The value is 0;

[0187] q 12 represents the correlation coefficient of the first gene relative to the filtered sequencing fragment 2 obtained in step 2. The filtered sequencing fragment 2 obtained in step 2 is aligned to the first gene, q 12 The value is 1; otherwise, q 12 The value is 0;

[0188] q 1m represents the correlation coefficient of the first gene relative to the filtered sequencing fragment m obtained in step 2. The filtered sequencing fragment m obtained in step 2 is aligned to the first gene, q 1m The value is 1; otherwise, q 1m The value is 0;

[0189] q 21 It represents the correlation coefficient of the second gene relative to the filtered sequencing fragment 1 obtained in step 2. The filtered sequencing fragment 1 obtained in step 2 is aligned to the second gene, q 21 The value is 1; otherwise, q 21 The value is 0;

[0190] q 22 represents the correlation coefficient of the second gene relative to the filtered sequencing fragment 2 obtained in step 2. The filtered sequencing fragment 2 obtained in step 2 is aligned to the second gene, q 22 The value is 1; otherwise, q 22 The value is 0;

[0191] q 2m represents the correlation coefficient of the second gene relative to the filtered sequencing fragment m obtained in step 2. The filtered sequencing fragment m obtained in step 2 is aligned to the second gene, q 2m The value is 1; otherwise, q 2m The value is 0;

[0192] q n1 represents the correlation coefficient of the nth gene relative to the filtered sequencing fragment 1 obtained in step 2. The filtered sequencing fragment 1 obtained in step 2 is aligned to the nth gene, q n1 The value is 1; otherwise, q n1 The value is 0;

[0193] q n2 represents the correlation coefficient of the nth gene relative to the filtered sequencing fragment 2 obtained in step 2. The filtered sequencing fragment 2 obtained in step 2 is aligned to the nth gene, q n2 The value is 1; otherwise, q n2 The value is 0;

[0194] q nm represents the correlation coefficient of the nth gene relative to the filtered sequencing fragment m obtained in step 2. The filtered sequencing fragment m obtained in step 2 is aligned to the nth gene, q nm The value is 1; otherwise, q nm The value is 0;

[0195] Since there are few non-zero entries in the matrix Q, it can be regarded as a sparse matrix and the storage and algorithms of sparse matrices can be used;

[0196] Step 3: Convert the n×m dimensional matrix Q n×m With Q n×m Transpose Multiply them together to get an n×n dimensional symmetric matrix R n×n ;

[0197]

[0198]

[0199] Among them, R n×n represents an n×n dimensional symmetric matrix; r αα′ Represents the symmetric matrix R n×n Items in

[0200] q αβ Denotes the matrix Q n×mThe correlation coefficient of the αth gene in step 2 relative to the filtered sequencing fragment β obtained in step 2, the filtered sequencing fragment β obtained in step 2 is aligned to the αth gene, q αβ The value is 1; otherwise, q αβ The value is 0;

[0201] q α′β Represents Q n×m The correlation coefficient of the α′th gene in step 2 relative to the filtered sequencing fragment β obtained in step 2, the filtered sequencing fragment β obtained in step 2 is aligned to the α′th gene, q α′β The value is 1; otherwise, q α′β The value is 0; α′∈[1,n];

[0202] q′ βα′ express The correlation coefficient of the α′th gene in step 2 relative to the filtered sequencing fragment β obtained in step 2, the filtered sequencing fragment β obtained in step 2 is aligned to the α′th gene, q′ βα′ The value is 1; otherwise q′ βα′ The value is 0; α′∈1,n]; β∈[1,m];

[0203] Steps 3 and 4: Symmetric matrix R n×n Each non-diagonal term r in αα′ Represents the number of sequencing fragments supporting the existence of fusion between gene α and gene α′;

[0204] Traverse the symmetric matrix R n×n All fusion gene pairs with fusion can be found.

[0205] The other steps and parameters are the same as those in the first to sixth embodiments.

[0206] Specific embodiment eight: This embodiment differs from any one of specific embodiments one to seven in that, in step four, each fusion gene pair is subjected to breakpoint restoration processing to obtain each fusion gene pair after breakpoint restoration;

[0207] After breakpoint restoration, the sequencing fragments of each fusion gene pair are hierarchically clustered from bottom to top according to the breakpoint position, and clusters with the number of sequencing fragments greater than or equal to two are found and retained. Each cluster is a fusion gene pair;

[0208] The specific process is:

[0209] Step 41: Perform breakpoint restoration on each fusion gene pair to obtain each fusion gene pair after breakpoint restoration; the specific process is as follows:

[0210] For each candidate fusion gene pair, there are four fusion methods;

[0211] Specify the order of two genes, labeled as g1 and g2, and detect the supported sequencing fragments. The possible situations are shown in Table 1. Therefore, the exons to which the breakpoints belong are restored according to the situation of each sequencing fragment.

[0212] There are two types of breakpoints: back-end breakpoints and front-end breakpoints (the start position of the gene is complete, which is the front-end breakpoint; the end position of the gene is complete, which is the back-end breakpoint);

[0213] Assume that gene g1 is in front and gene g2 is in the back;

[0214] The flag = 0 of the sequencing fragment corresponding to gene g1 is the back-end breakpoint;

[0215] Flag = 16 of the sequencing fragment corresponding to gene g2 is the back-end breakpoint;

[0216] Flag = 16 of the sequencing fragment corresponding to gene g1 is the front breakpoint;

[0217] The flag = 0 of the sequencing fragment corresponding to gene g2 is the front breakpoint;

[0218] Restore the backend breakpoint; the specific process is:

[0219] By subtracting the breakpoint position of the sequenced fragment filtered in step 2 from the end position of the exon of the reconstructed reference transcript in step 1 (both positions are known);

[0220] If the absolute value of the difference is less than the threshold, the breakpoint is considered to belong to the exon of the corresponding recombinant reference transcript in step 1, and the breakpoint position is restored to the reference sequence; the expression is:

[0221] POS ref =pos trans -exon trans +exon ref

[0222] Among them, pos trans Indicates the breakpoint position of the sequencing fragment after filtering in step 2; exon trans Indicates the end position of the exon corresponding to the recombinant reference transcript in step 1, exon ref Indicates the end position of the exon in the reference sequence, pos ref Indicates the breakpoint position on the reference sequence;

[0223] The threshold is 30 bp;

[0224] If the absolute value of the difference is greater than or equal to the threshold, the breakpoint is considered not to belong to the exon of the corresponding recombinant reference transcript in step 1, and the corresponding fusion gene pair is filtered;

[0225] The specific process of restoring the front-end breakpoint is as follows:

[0226] By subtracting the breakpoint position of the sequenced fragment filtered in step 2 from the start position of the exon of the reconstructed reference transcript in step 1;

[0227] If the absolute value of the difference is less than the threshold, the breakpoint is considered to belong to the exon of the corresponding recombinant reference transcript in step 1, and the breakpoint position is restored to the reference sequence; the expression is:

[0228] POS ref =pos trans -exon′ trans +exon′ ref

[0229] Among them, pos trans Indicates the breakpoint position of the sequencing fragment after filtering in step 2; exon′ trans Indicates the starting position of the exon corresponding to the recombinant reference transcript in step 1, exon′ ref Indicates the starting position of the exon in the reference sequence, pos ref Indicates the breakpoint position on the reference sequence;

[0230] If the absolute value of the difference is greater than or equal to the threshold, the breakpoint is considered not to belong to the exon of the corresponding recombinant reference transcript in step 1, and the corresponding fusion gene pair is filtered;

[0231] Step 42: After breakpoint restoration, perform bottom-up hierarchical clustering on the sequencing fragments of each fusion gene pair according to the breakpoint position, and find and retain clusters with the number of sequencing fragments greater than or equal to two. Each cluster is a fusion gene pair.

[0232] Table 1: Possible fusion gene situations

[0233]

[0234] The other steps and parameters are the same as those in the first to sixth embodiments.

[0235] Specific embodiment nine: This embodiment differs from any one of specific embodiments one to eight in that, in step four-two, the sequencing fragments of each fusion gene pair after breakpoint restoration are hierarchically clustered from bottom to top according to the breakpoint positions, and clusters with two or more sequencing fragments are found and retained, and each cluster is a fusion gene pair;

[0236] The specific process is:

[0237] The two breakpoints of a sequencing fragment are regarded as a point in a two-dimensional plane space;

[0238] Calculate the Manhattan distance between the breakpoint positions of the sequenced fragment genes, cluster the sequenced fragments according to the Manhattan distance, obtain clusters, and find and retain clusters with more than or equal to two sequenced fragments;

[0239] The sequencing fragments are clustered according to the Manhattan distance to obtain clusters; the process is:

[0240] Add the sequencing fragments whose distance from the cluster center is less than N to the cluster, and recalculate the cluster center (the cluster center is calculated based on the center of the sequencing fragment position in the cluster. The first cluster is a randomly selected sequencing fragment, which is itself a cluster, and the center of this cluster is the position of this sequencing fragment) until all sequencing fragments are clustered;

[0241] Manhattan distance = |x c -x d |+|y c -y d |

[0242] in,

[0243] x c Indicates the position of the breakpoint of gene g1 of the first sequenced fragment; y c Indicates the position of the breakpoint of gene g2 of the first sequenced fragment;

[0244] x d Indicates the position of the breakpoint of gene g1 in the second sequenced fragment, y d Indicates the position of the breakpoint of gene g2 in the second sequenced fragment.

[0245] The other steps and parameters are the same as those in Specific Embodiments 1 to 8.

[0246] Specific embodiment ten: This embodiment differs from any one of specific embodiments one to nine in that, in step six, the fusion gene pairs obtained in step four are screened based on the RNA-seq comparison results to obtain the final fusion gene pairs;

[0247] The specific process is:

[0248] If the RNA-seq alignment result obtained in step 5 fails to map to the reference sequence in step 1, or if it maps to the reference sequence in step 1 but the alignment record contains only one gene (the mapping result may contain many genes), then there is no split fragment signal and it is considered a false fusion gene. The fusion gene pair obtained in step 4 is filtered;

[0249] If the RNA-seq alignment result obtained in step 5 is mapped to the reference sequence in step 1 and the alignment record contains two or more genes (the mapping result may contain many genes), there is a split fragment signal. However, if the position difference between the breakpoint corresponding to the RNA-seq alignment result obtained in step 5 and the center of the cluster corresponding to step 4 (the same gene) exceeds the threshold, it is considered to be a possible false fusion gene, and the fusion gene pair obtained in step 4 is retained.

[0250] If the RNA-seq alignment result obtained in step 5 is mapped to the reference sequence in step 1 and the alignment record contains two or more genes (the mapping result may contain many genes), there is a split fragment signal, and the position difference between the breakpoint corresponding to the RNA-seq alignment result obtained in step 5 and the center of the cluster corresponding to step 4 (the same gene) is less than or equal to the threshold, then it is determined to be a true fusion gene, and the fusion gene pair obtained in step 4 is retained.

[0251] The other steps and parameters are the same as those in Specific Embodiments 1 to 9.

[0252] Example:

[0253] In order to test the performance of the fusion gene detection method proposed in the present invention, the following test method is proposed:

[0254] Compare with existing known results. Existing research provides access to fusion gene pairs that have been verified to exist in relevant samples. This data allows for sensitivity validation of experimental results from six real-world datasets.

[0255] Results from real-world data were cross-validated. By analyzing experimental results in which all methods successfully ran, we identified the correct answer for each dataset that could be found by at least three fusion gene detection methods. The correct answer was then used to determine the precision, recall, and F1-score of each method. The calculation formula is as follows:

[0256]

[0257]

[0258]

[0259] The temporal and spatial performance of the tools. This paper tests the temporal and spatial performance of five fusion gene detection tools in the SK-BR-3PacBio dataset in 1, 2, 4, and 8 threads and 170GB environments.

[0260] The real data used for testing include various data from samples such as MCF-7, SK-BR-3, and HCT-116. In order to explore the advantages and disadvantages of various fusion gene detection methods, the experimental results on these eight long-read datasets were combined with multiple evaluation methods for the above analysis.

[0261] Comparison with Existing Results: The experimental results are shown in Table 2. LongGF and JAFFAL failed to run correctly due to the large size of the HCT-116 cDNA dataset. As can be seen from the table, the proposed method found the most known fusion genes in the HCT-116 and SK-BR-3 datasets. It also achieved superior results in the MCF-7 dataset, falling behind the top-performing JAFFAL method by only one or two genes.

[0262] Results from real-world data were cross-validated. The experimental results are shown in Table 3. Under relatively stringent criteria, the proposed method achieved the best results on most datasets, achieving the highest precision and recall. On the MCF-7PacBio and MCF-7dRNA datasets, despite some shortcomings in precision, it achieved the highest F1-scores of 0.4 and 0.81, respectively.

[0263] Table 4 shows the experimental results of spatial and temporal performance. Due to their different architectures, LongGF, Genion, and FusionSeeker all require data to be aligned to the reference genome first. This alignment step consumes a significant amount of time, resulting in significantly longer total execution times for these tools than the proposed method and JAFFAL. Across all four thread counts, the proposed method consumes approximately eight times as much time as the proposed method. However, because the proposed method supports multithreaded processing for fusion gene detection in addition to minimap2 alignment, increasing the number of threads significantly reduces execution time compared to JAFFAL. In this scenario, with a thread count of 1, the proposed method is significantly faster than JAFFAL under multithreaded conditions, while maintaining a relatively small difference in execution time compared to JAFFAL. With 8 threads, the proposed method is 70% of JAFFAL's time. Therefore, the proposed method achieves relatively ideal execution time performance. Regarding memory, since the memory overhead for aligning data to the reference genome is relatively fixed (approximately 20 GB) and is significantly higher than that for the fusion gene detection step, all other tools, except the proposed method, consume approximately 20 GB of memory. The method of the present invention, however, significantly reduces the number of sequence fragments in the dataset during the pre-filtering step. This results in fewer sequence fragments remaining during the subsequent re-alignment, thus reducing the memory overhead of the second alignment. As a result, the memory usage of the method of the present invention is significantly lower than that of other tools, ranging from 32% to 56%. Therefore, the method of the present invention outperforms other tools in terms of space performance.

[0264] Table 2 Comparison with existing known results

[0265]

[0266] Table 3 Real data results mutual verification

[0267]

[0268] Table 4. Spatial and temporal performance of the tools

[0269]

[0270]

[0271] The present invention may have many other embodiments. Without departing from the spirit and essence of the present invention, those skilled in the art may make various corresponding changes and modifications based on the present invention, but these corresponding changes and modifications should all fall within the scope of protection of the claims attached to the present invention.

Claims

1. A method for detecting fusion genes in long-read RNA-seq sequencing data, characterized by: The specific process of the method is: Step 1: Recombining the reference sequence based on the gene annotation file to generate a recombinant reference transcript containing only the transcript sequence and the recombinant transcript annotation; Step 2: Use the DNA alignment mode of the open-source alignment tool minimap2 to align the sequence fragments to be tested to the recombinant reference transcripts, filter the sequence fragments to be tested on the alignment, and obtain the filtered sequence fragments; The genes corresponding to the filtered sequencing fragments constitute the gene set G; The filtered sequencing fragments contain a sequence identifier flag; Step 3: Search for fusion gene pairs in the filtered sequencing fragments; Step 4: Perform breakpoint restoration on each fusion gene pair to obtain each fusion gene pair after breakpoint restoration; After breakpoint restoration, the sequencing fragments of each fusion gene pair are hierarchically clustered according to the breakpoint position, and clusters with the number of sequencing fragments greater than or equal to two are found and retained. Each cluster is a fusion gene pair; Step 5: Use the open source sequence alignment tool abpoa to align the sequencing fragments in each cluster and output the consensus sequence; Use the RNA-seq alignment function of the open source alignment tool minimap2 to align the consensus sequence to the reference sequence in step 1 to obtain the RNA-seq alignment results; Step 6: Based on the RNA-seq comparison results, the fusion gene pairs obtained in step 4 are screened to obtain the final fusion gene pairs; In the step 1, the reference sequence is recombined based on the gene annotation file to generate a recombined reference transcript containing only the transcript sequence and the recombinant transcript annotation; The specific process is: Step 1: Read the gene annotation file and retain the information of the transcripts and the exons corresponding to the retained transcripts in the gene annotation file as required; Step 12: Determine the sequence of each exon based on the start and end positions of each exon information retained in step 11 and the reference sequence, and link the sequences of all exons to obtain the transcript sequence; Step 13: Connect the transcript sequences according to chromosomes to obtain recombinant reference transcripts; Record the start and end positions on the reference sequence and the start and end positions of the recombinant reference transcript to obtain the recombinant transcript annotation; Step 14: Use the index function of minimap2 to create an alignment index between the recombinant reference transcript and the reference sequence; In the step 2, the sequenced fragments to be detected are filtered to obtain filtered sequenced fragments; The specific process is: Step 21: Filter according to the sequence length of the sequencing fragments; Step 22: Filter based on the number of genes covered by the sequencing fragments; Step 23: Filtering is performed based on the boundaries of each sequencing fragment retained in step 22 and the boundaries of the exons of the corresponding recombinant reference transcripts to obtain filtered sequencing fragments; In the step 3, searching for fusion gene pairs in the filtered sequencing fragments obtained in the step 2; The specific process is: Step 31: Assume that the number of filtered sequencing fragments obtained in step 2 is m; After determining the gene set G, each gene in the gene set is represented as an m-dimensional vector; gene α =(q α1 ,q α2 ,…,q αm ),where q αβ =0 or 1,α∈[1,n],β∈[1,m] Among them, gene α represents the αth gene; q α1 represents the correlation coefficient of the αth gene relative to the filtered sequencing fragment 1 obtained in step 2. The filtered sequencing fragment 1 obtained in step 2 is aligned to the αth gene, q α1 The value is 1; otherwise, q α1 The value is 0; q α2 represents the correlation coefficient of the αth gene relative to the filtered sequencing fragment 2 obtained in step 2. The filtered sequencing fragment 2 obtained in step 2 is aligned to the αth gene, q α2 The value is 1; otherwise, q α2 The value is 0; q αm represents the correlation coefficient of the αth gene relative to the filtered sequencing fragment m obtained in step 2. The filtered sequencing fragment m obtained in step 2 is aligned to the αth gene, q αm The value is 1; otherwise, q αm The value is 0; β represents the sequence number of the filtered sequencing fragment obtained in step 2; n represents the total number of genes in the gene set G; m represents the total number of sequencing fragments after filtering obtained in step 2; Step 32: Get an n×m dimensional matrix Q based on n genes n×m ; Among them, Q n×m Represents a matrix of n×m dimensions; gene1 indicates the first gene, gene2 indicates the second gene, gene n represents the nth gene; q 11 It represents the correlation coefficient of the first gene relative to the filtered sequencing fragment 1 obtained in step 2. The filtered sequencing fragment 1 obtained in step 2 is aligned to the first gene, q 11 The value is 1; otherwise, q 11 The value is 0; q 12 represents the correlation coefficient of the first gene relative to the filtered sequencing fragment 2 obtained in step 2. The filtered sequencing fragment 2 obtained in step 2 is aligned to the first gene, q 12 The value is 1; otherwise, q 12 The value is 0; q 1m represents the correlation coefficient of the first gene relative to the filtered sequencing fragment m obtained in step 2. The filtered sequencing fragment m obtained in step 2 is aligned to the first gene, q 1m The value is 1; otherwise, q 1m The value is 0; q 21 It represents the correlation coefficient of the second gene relative to the filtered sequencing fragment 1 obtained in step 2. The filtered sequencing fragment 1 obtained in step 2 is aligned to the second gene, q 21 The value is 1; otherwise, q 21 The value is 0; q 22 represents the correlation coefficient of the second gene relative to the filtered sequencing fragment 2 obtained in step 2. The filtered sequencing fragment 2 obtained in step 2 is aligned to the second gene, q 22 The value is 1; otherwise, q 22 The value is 0; q 2m represents the correlation coefficient of the second gene relative to the filtered sequencing fragment m obtained in step 2. The filtered sequencing fragment m obtained in step 2 is aligned to the second gene, q 2m The value is 1; otherwise, q 2m The value is 0; q n1 represents the correlation coefficient of the nth gene relative to the filtered sequencing fragment 1 obtained in step 2. The filtered sequencing fragment 1 obtained in step 2 is aligned to the nth gene, q n1 The value is 1; otherwise, q n1 The value is 0; q n2 represents the correlation coefficient of the nth gene relative to the filtered sequencing fragment 2 obtained in step 2. The filtered sequencing fragment 2 obtained in step 2 is aligned to the nth gene, q n2 The value is 1; otherwise, q n2 The value is 0; q nm represents the correlation coefficient of the nth gene relative to the filtered sequencing fragment m obtained in step 2. The filtered sequencing fragment m obtained in step 2 is aligned to the nth gene, q nm The value is 1; otherwise, q nm The value is 0; Step 3: Convert the n×m dimensional matrix Q n×m With Q n×m Transpose Multiply them together to get an n×n dimensional symmetric matrix R n×n ; Among them, R n×n represents an n×n dimensional symmetric matrix; r αα′ Represents the symmetric matrix R n×n Items in q αβ Denotes the matrix Q n×m The correlation coefficient of the αth gene in step 2 relative to the filtered sequencing fragment β obtained in step 2, the filtered sequencing fragment β obtained in step 2 is aligned to the αth gene, q αβ The value is 1; otherwise, q αβ The value is 0; q α′β Indicates Q n×m The correlation coefficient of the α′th gene in step 2 relative to the filtered sequencing fragment β obtained in step 2, the filtered sequencing fragment β obtained in step 2 is aligned to the α′th gene, q α′β The value is 1; otherwise, q α′β The value is 0; α′∈[1,n]; q′ βα′ express The correlation coefficient of the α′th gene in step 2 relative to the filtered sequencing fragment β obtained in step 2, the filtered sequencing fragment β obtained in step 2 is aligned to the α′th gene, q′ βα′ The value is 1; otherwise q′ βα′ The value is 0; α′∈1,n]; β∈[1,m]; Steps 3 and 4: Symmetric matrix R n×n Each non-diagonal term r in αα′ Represents the number of sequencing fragments supporting the existence of fusion between gene α and gene α′; Traverse the symmetric matrix R n×n You can find all fusion gene pairs that have fusion; In the step 4, each fusion gene pair is subjected to breakpoint restoration processing to obtain each fusion gene pair after breakpoint restoration; After breakpoint restoration, the sequencing fragments of each fusion gene pair are hierarchically clustered according to the breakpoint position, and clusters with the number of sequencing fragments greater than or equal to two are found and retained. Each cluster is a fusion gene pair; The specific process is: Step 41: Perform breakpoint restoration on each fusion gene pair to obtain each fusion gene pair after breakpoint restoration; the specific process is as follows: There are two types of breakpoints: back-end breakpoints and front-end breakpoints. Assume that gene g1 is in front and gene g2 is in the back; The flag = 0 of the sequencing fragment corresponding to gene g1 is the back-end breakpoint; Flag = 16 of the sequencing fragment corresponding to gene g2 is the back-end breakpoint; Flag = 16 of the sequencing fragment corresponding to gene g1 is the front breakpoint; The flag = 0 of the sequencing fragment corresponding to gene g2 is the front breakpoint; Restore the backend breakpoint; the specific process is: By subtracting the breakpoint position of the sequenced fragment filtered in step 2 from the end position of the exon of the reconstructed reference transcript in step 1; If the absolute value of the difference is less than the threshold, the breakpoint is considered to belong to the exon of the corresponding recombinant reference transcript in step 1, and the breakpoint position is restored to the reference sequence; the expression is: pos ref =pos trans -exon trans +exon ref Among them, pos trans Indicates the breakpoint position of the sequencing fragment after filtering in step 2; exon trans Indicates the end position of the exon corresponding to the recombinant reference transcript in step 1, exon ref Indicates the end position of the exon in the reference sequence, pos ref Indicates the breakpoint position on the reference sequence; If the absolute value of the difference is greater than or equal to the threshold, the breakpoint is considered not to belong to the exon of the corresponding recombinant reference transcript in step 1, and the corresponding fusion gene pair is filtered; The specific process of restoring the front-end breakpoint is as follows: By subtracting the breakpoint position of the sequenced fragment filtered in step 2 from the start position of the exon of the reconstructed reference transcript in step 1; If the absolute value of the difference is less than the threshold, the breakpoint is considered to belong to the exon of the corresponding recombinant reference transcript in step 1, and the breakpoint position is restored to the reference sequence; the expression is: pos ref =pos trans -exon′ trans +exon′ ref Among them, pos trans Indicates the breakpoint position of the sequencing fragment after filtering in step 2; exon′ trans Indicates the starting position of the exon corresponding to the recombinant reference transcript in step 1, exon′ ref Indicates the starting position of the exon in the reference sequence, pos ref Indicates the breakpoint position on the reference sequence; If the absolute value of the difference is greater than or equal to the threshold, the breakpoint is considered not to belong to the exon of the corresponding recombinant reference transcript in step 1, and the corresponding fusion gene pair is filtered; Step 42: After breakpoint restoration, the sequencing fragments of each fusion gene pair are hierarchically clustered according to the breakpoint position, and clusters with the number of sequencing fragments greater than or equal to two are found and retained. Each cluster is a fusion gene pair.

2. A fusion gene detection method for long-read RNA-seq sequencing data according to claim 1, characterized in that: In the step 21, filtering is performed according to the sequence length of the sequencing fragments; The specific process is; The sequence fragments to be detected in the alignment whose sequence length is greater than the length threshold are retained; The sequencing fragments to be detected in the alignment whose sequence length is less than or equal to the length threshold are deleted.

3. A fusion gene detection method for long-read RNA-seq sequencing data according to claim 2, characterized in that: In step 22, filtering is performed based on the number of genes covered by the sequencing fragments; The specific process is: Based on the annotation of the recombinant transcripts, determine the gene corresponding to the sequence information of each sequencing fragment retained in step 2-1; If the number of genes corresponding to the sequence information of the sequencing fragment is greater than or equal to two, the sequencing fragment is retained; If the number of genes corresponding to the sequence information of the sequencing fragment is less than two, the sequencing fragment is deleted.

4. A method for detecting fusion genes in long-read RNA-seq sequencing data according to claim 3, characterized in that: In step 23, filtering is performed based on the boundaries of each sequencing fragment retained in step 22 and the boundaries of the exons of the corresponding recombinant reference transcript to obtain filtered sequencing fragments; The specific process is: 1) Calculate the distance between the start position of a sequencing fragment A retained in step 22 and the start position of the exon of the corresponding recombinant reference transcript; 2) Calculate the distance between the starting position of a sequencing fragment A retained in step 22 and the ending position of the exon of the corresponding recombinant reference transcript; 3) Calculate the distance between the end position of a sequencing fragment A retained in step 22 and the start position of the exon of the corresponding recombinant reference transcript; 4) Calculate the distance between the end position of a sequencing fragment A retained in step 22 and the end position of the exon of the corresponding recombinant reference transcript; If any of the distances in 1), 2), 3), and 4) is less than the distance threshold, sequencing fragment A is retained; otherwise, sequencing fragment A is deleted; This process continues until all the sequenced fragments retained in step 22 are determined to obtain the final sequenced fragments retained.

5. A method for detecting fusion genes based on long-read RNA-seq sequencing data according to claim 4, characterized in that: In the step 42, hierarchical clustering is performed on the sequencing fragments of each fusion gene pair after breakpoint restoration according to the breakpoint positions, and clusters with the number of sequencing fragments greater than or equal to two are found and retained, each cluster being a fusion gene pair; The specific process is: Calculate the Manhattan distance between the breakpoint positions of the sequenced fragment genes, cluster the sequenced fragments according to the Manhattan distance, obtain clusters, and find and retain clusters with more than or equal to two sequenced fragments; The sequencing fragments are clustered according to the Manhattan distance to obtain clusters; the process is: Add the sequencing fragments whose distance from the cluster center is less than N to the cluster and recalculate the cluster center until all sequencing fragments are clustered; Manhattan distance = |x c -x d |+|y c -y d | in, x c Indicates the position of the breakpoint of gene g1 of the first sequenced fragment; y c Indicates the position of the breakpoint of gene g2 of the first sequenced fragment; x d Indicates the position of the breakpoint of gene g1 in the second sequenced fragment, y d Indicates the position of the breakpoint of gene g2 of the second sequenced fragment.

6. A method for detecting fusion genes based on long-read RNA-seq sequencing data according to claim 5, characterized in that: In step 6, the fusion gene pairs obtained in step 4 are screened according to the RNA-seq comparison results to obtain the final fusion gene pairs; The specific process is: If the RNA-seq alignment result obtained in step 5 fails to map to the reference sequence in step 1, or maps to the reference sequence in step 1 but the alignment record contains only one gene, there is no split fragment signal and it is considered a false fusion gene. The fusion gene pair obtained in step 4 is filtered; If the RNA-seq alignment result obtained in step 5 is mapped to the reference sequence in step 1 and the alignment record contains two or more genes, there is a split fragment signal. However, if the position difference between the breakpoint corresponding to the RNA-seq alignment result obtained in step 5 and the center of the cluster corresponding to step 4 exceeds the threshold, it is considered to be a possible false fusion gene, and the fusion gene pair obtained in step 4 is retained. If the RNA-seq alignment result obtained in step 5 is mapped to the reference sequence in step 1 and the alignment record contains two or more genes, there is a split fragment signal, and the position difference between the breakpoint corresponding to the RNA-seq alignment result obtained in step 5 and the center of the cluster corresponding to step 4 is less than or equal to the threshold, then it is determined to be a true fusion gene, and the fusion gene pair obtained in step 4 is retained.

Images