A transcript assembly method based on reference genome

Abstract

The present invention discloses a method for assembling transcripts based on a reference genome, which belongs to the field of transcript assembly in bioinformatics, and the specific steps are: (a) using at least two alignment tools to analyze the original sequencing data obtained from the machine; (b) constructing a splicing graph for the reply results of each alignment tool; (c) merging the splicing graphs of each alignment tool to generate a label splicing graph; (d) extracting the double-end paths marked in the label graph; (e) searching for the coverage of the paths in the label splicing graph, and finally obtaining the coverage of the paths representing the transcripts, wherein all the double-end paths are covered by the assembled transcripts. The method of the present invention combines the reply results of different alignment tools, and constructs a label splicing graph based on the reply results of different alignment tools, traverses the graph through a dynamic path extension algorithm, thereby finding the path representing the transcript, taking into account the characteristics and advantages of different alignment tools, and improving the accuracy of the reconstructed transcripts.

Description

Technical Field

[0001] The invention relates to a transcript assembly method based on a reference genome which combines reply results of different alignment tools, and belongs to the field of transcript assembly in bioinformatics. Background Art

[0002] As a powerful transcriptome analysis technology, RNA-seq has been widely used worldwide. Especially in the past five years, this technology has transitioned from research to clinical applications, which has provided clues for studying complex diseases such as cancer associated with abnormal splicing events or differential expression levels. In addition, it provides the opportunity to observe the complexity of the eukaryotic transcriptome, identify expressed transcriptomes, and accurately quantify their expression abundance at the whole transcriptome level. One of the key steps in RNA-seq data analysis is to accurately assemble a large number of sequencing reads into full-length transcripts, which is a rather computationally challenging task.

[0003] The explosive growth of RNA-seq data has driven the development of transcriptome assembly algorithms. A considerable number of algorithms have been developed to assemble RNA-seq reads into full-length transcripts. They are generally divided into two strategies, de novo assembly and reference genome-based assembly. De novo assembly usually attempts to find overlaps between RNA-seq reads and uses extension technology to reconstruct full-length transcripts. This method is usually used in the absence of a reference genome, so this strategy mostly produces highly fragmented and error-prone transcripts. The most advanced de novo assembly software currently includes TransLiG, BinPacker, Bridger, Trinity, ABySS, SOAPdenovo-Trans, and IDBA-Tran. In contrast, if some species have high-quality reference genome templates (such as humans), reference genome-based transcript assembly software such as StringTie2, StringTie, Scallop, TransComb, Cufflinks, CLASS2, iPAC, Traph, CEM, IsoLasso, and Bayesembler can be used. Such methods usually first use alignment tools such as Hisat2, Hisat, Star, Tophat2, Tophat, SpliceMap, MapSplice or GSNAP to post RNA-seq reads back to the reference genome; then, based on the alignment results, a splicing graph model is established for each gene site, and then different computational models are used to traverse the graph to generate paths representing transcripts. Reference genome-based methods are usually more accurate than de novo assembly methods. In addition, many tools have been developed to merge transcripts from multiple RNA-seq samples, such as the merging methods of StringTie2 and TACO, but the reconstructed transcripts still have many shortcomings.

[0004] None of the current assembly tools are specifically designed to merge alignments generated by different alignment tools, which has important practical significance. In addition, through experiments, it was found that assembly software that relies on only one alignment tool cannot completely reconstruct a specific transcript, but it can be reconstructed by the proposed Tiglon software. This shows that merging the alignment results of different alignment tools will be a useful and meaningful strategy for transcriptome assembly. Summary of the invention

[0005] In response to the above-mentioned prior art, the present invention provides a transcript assembly method based on a reference genome that combines the reply results of different alignment tools. This method, on the basis of previous transcript assembly based on the reference genome, combines the reply results of different alignment tools, and constructs a label splicing graph based on the reply results of different alignment tools. The graph is traversed by a dynamic path extension algorithm to find the path representing the transcript, and the characteristics and advantages of different alignment tools are comprehensively considered to improve the accuracy of the reconstructed transcripts.

[0006] The present invention is achieved through the following technical solutions:

[0007] A transcript assembly method based on a reference genome comprises the following steps:

[0008] (i) extracting mRNA from the sample to be tested, amplifying it, obtaining circularized cDNA, and sequencing it on a machine; using at least two alignment tools to analyze the raw sequencing data obtained off the machine, and obtaining the result mapped back to the reference genome.

[0009] The alignment tool is selected from Hisat2 and Star.

[0010] (II) Construct a splicing graph for each comparison tool’s reply results:

[0011] Based on the results of reads being mapped back to the reference genome, they are clustered into different gene loci, and exon-exon splice sites are derived from these spliced ​​reads;

[0012] Specifically, a traditional splicing graph G = (V, E) is constructed for each gene locus, where each node v corresponds to an exon and each edge e corresponds to a splicing site between two exons;

[0013] Furthermore, edges and nodes are weighted by the number of reads supporting them; edges or nodes with low weights that may be caused by sequencing errors or unreasonable alignments are removed.

[0014] (III) Merge the splicing graphs of each alignment tool to generate a label splicing graph, as follows:

[0015] (1) Assume that we have obtained the reply results of N comparison tools, and the splicing graph generated by comparison tool i is G i =(E i , V i ), i≤N, for each edge e i ∈E i , support e i The set of reads is R ei ; Label splice graph G L =(E L , VL ) is obtained by merging each splice graph G i To generate, i∈[1,N], where the edges and points represent all the nodes that appear in the graph G i The edges and points in , i∈[1,N];

[0016] (2) Label the splicing graph G L :Through the following method, for each L For an edge e in , if there exists j∈[1,N] such that the edge Then mark the edge e as 1, otherwise mark it as 2; In addition, for every edge e∈E L Marked as 2, if there exist i and j such that R ei ≠R ej , i, j ∈ [1, N], and i ≠ j, label it as 2-, otherwise label it as 2+; then, the weight of the edge of the label splicing graph is assigned to each edge that appears in G i The average weight of the edges in , i≤N;

[0017] (IV) Extract the two-terminal circuit marked in the label diagram, as follows:

[0018] (1) For each gene locus, generate a set of two-terminal paths for each alignment tool i, denoted as PP i ;

[0019] (1.1) For each pair of double-end reads r1 and r2, based on the response of alignment tool i, if L In the example, r1 spans a path p1=n j1 →n j2 →…→n jk , r2 spans a path p2 = n j’1 →n j’2 →…→n j’q , and p1 and p2 share a compatible subpath p in =n m1 →n m2 →…→n ms Satisfy k+qs≥3;

[0020] (1.2) Through the shared compatible subpath p in , connect the two paths p1 and p2 to generate a two-terminal path p;

[0021] (2) After processing all replies to G L After the paired-end reads are obtained, the set PP of all paired-end reads that depend on alignment tool i is obtained. i ;

[0022] (3) Along each path belonging to p∈PP i Paired-end reads set R ei The generated paired-end pathway set PP from each alignment tool i i Then, let PP L =∪ i∈[1,N] PP i The double-ended path of the label splicing graph is then labeled as PP L , marking method and marking edge e∈E L The method is the same as: for each L A two-terminal path p in a Then mark the two-terminal path p as 1, otherwise mark it as 2; In addition, for every two-terminal path p∈PP i Marked as 2, if there exist i and j such that R ei ≠R ej , i, j∈[1, N], and i≠j, mark it as 2-, otherwise mark it as 2+;

[0023] (V) Find the coverage of the path in the label splicing graph, and use the label-based dynamic programming path extraction algorithm, as follows:

[0024] (1) Select a seed and generate a subgraph using a seed:

[0025] A seed is an edge or a paired path in the labeled graph, which can be further extended to a path representing the full-length transcript. An unused paired path with the largest weight and labelled 2+ is selected as a seed. If such paired paths do not exist or they are all included in the assembled transcript, the seeds are selected in the following order: paired paths with a 2-label, edges with a 2+ label, edges with a 2-label, paired paths with a label of 1, or edges with a label of 1.

[0026] Paired-end paths or edges with label 2 represent that they are captured by multiple alignment tools, that is, they are most likely part of positive transcripts, and these seeds will be extended to full-length transcripts with a higher probability;

[0027] With a two-ended path S = v k →v k +1 →… →v l Or side S = v k →v l is the seed, extending all its right / left adjacent points v k (v l ), and extends until its adjacent point encounters a point without an outgoing / incoming edge, labeling the spliced ​​graph G L The subgraph GS is generated by the selected seed; then, the source / sink is added to the subgraph G S In the example, if all the nodes connected to it have no inbound / outbound edges, the weight of the newly added edge is set to ∞;

[0028] (2) Expanding the path through a label-based dynamic programming algorithm:

[0029] Traversing the graph G through a dynamic programming algorithm S The edges of the label clipping graph G S The idea of ​​dynamic programming is to find a path with a higher weight from the seed to each other adjacent edge;

[0030] Assume that graph G S The edges in are topologically sorted. Without loss of generality, they are numbered 1, 2...N. For each edge i=v→w, where v and w are two points in the label splicing graph, predecessor(i) is set to be the set of all incoming edges of point v, successor(i) is set to be the set of all outgoing edges of point w, and weight(i) is the coverage of edge i.

[0031] Assume that the edge v of seed S k →v k+1 corresponds to edge m, and edge v of seed S l-1 →v l Corresponding to edge n, for edge i from n to N (from m to l) there are two variables Exp i and Pre i (Suc i ), where Exp i represents the upper bound of the highest expression value of the path from seed to edge i, Pre i (Suc i ) is the predecessor / successor of edge i, if the seed is an edge then m=n; the dynamic programming algorithm is derived from Exp i =∞,Exp m =cov(S),Exp n = cov(S), and Pre i =nullFor i≥n, Suc i = null for each edge i≤m start;

[0032] For i from n+1 to N, if predecessor(i) is not empty, select label 2+Exp from predecessor(i) i’The largest edge i'. Note that if there is no edge with label 2+ or the coverage of the edge with label 2+ is 0, then select edge i' in the following order: edge with label 2-, edge with label 1; then set Exp i is min{weight(i), Exp i’} and let Prei be i'. For i from m-1 to 1, use the same method to set Exp i and Suc i , starting from the edge with the largest Exp value connecting the sink / source, backward / forward based on the label Pre(Suc) until reaching the seed S, then graph G S The most highly expressed transcript in h will be covered;

[0033] (3) Update the label splice map and repeat the process:

[0034] Define f min is the extended path p h The minimum edge weight on p h For each edge e in, update the weight weight(e) to weight(n)-f min Repeat (1) and (2) in the path extraction process until graph G L All edges in are covered by the predicted transcripts, and finally, a path coverage representing the transcripts is obtained, in which all double-end paths are covered by the assembled transcripts.

[0035] The method of the present invention considers merging the reply results of different alignment tools on the basis of the transcript assembly method based on the reference genome in the past, and based on the dynamic programming algorithm based on the label to find the coverage of the path in the label splicing graph, it makes a great contribution to reconstructing a more accurate transcript. The existing methods all only consider the reply results of a kind of alignment tool, and this assembly method that only considers a single reply result will inevitably lead to the loss of some reads information. Therefore, on this basis, the present invention merges the reply results of different alignment tools, reconstructs a more accurate transcript, and lays the foundation for downstream analysis.

[0036] The method of the present invention has the following beneficial effects:

[0037] (1) The present invention combines the reply results of different comparison tools and constructs a label splicing graph, thereby being able to extract more reliable two-terminal circuits from the splicing graph.

[0038] (2) The present invention uses a dynamic programming algorithm to strategically use label information to reconstruct transcripts.

[0039] (3) The present invention constructs a label splicing graph by combining the reply results of different alignment tools, and finds the coverage of the paths representing transcripts in the splicing graph through a label-based dynamic programming algorithm, thereby improving the accuracy of the reconstructed transcripts. It is expected to play a key role in the discovery of transcriptome research using RNA-seq, especially in complex human diseases related to abnormal splicing events and expression levels, such as cancer.

[0040] Various terms and phrases used herein have the general meanings that are well known to those skilled in the art. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] Figure 1 : Extract the schematic diagram of the two-terminal circuit marked in the label diagram.

[0042] Figure 2 : Schematic diagram of the performance evaluation of various related software on the simulated data set, where A) represents the accuracy and correct output number evaluation of each software on the simulated data set; B) represents the F-Score of each program on the simulated data set; C) represents the evaluation of the correct output number of correctly assembled transcripts detected on the simulated data set at low, medium and high expression levels. The abbreviation ST stands for StringTie2, SC stands for Scallope, and MA stands for MergedAlignments.

[0043] Figure 3 : Schematic diagram of performance evaluation on eight Homo sapiens sample H1-H8 datasets, where A) represents the evaluation of the accuracy and correct output number of each software on the eight Homo sapiens sample datasets; B) represents the F-Score of each software on the eight Homo sapiens sample datasets; C) represents the average number evaluation of the correctly assembled transcripts detected on the eight Homo sapiens sample datasets at low, medium and high expression levels. The abbreviation ST stands for StringTie2, SC stands for Scallope, and MA stands for MergedAlignments.

[0044] Figure 4 : Schematic diagram of performance evaluation on four mouse muscle sample datasets, where A) represents the evaluation of the accuracy and correct output quantity of each software on four mouse muscle sample datasets; B) represents the F-Score of each software on four mouse muscle sample datasets; C) represents the average number evaluation of the correctly assembled transcripts detected on four mouse muscle sample datasets at low, medium and high expression levels. The abbreviation ST stands for StringTie2, SC stands for Scallope, and MA stands for MergedAlignments. . DETAILED DESCRIPTION

[0045] The present invention will be further described below in conjunction with the examples. However, the scope of the present invention is not limited to the following examples. Those skilled in the art will appreciate that various changes and modifications may be made to the present invention without departing from the spirit and scope of the present invention.

[0046] The instruments, reagents, materials, etc. involved in the following embodiments, unless otherwise specified, are all conventional instruments, reagents, materials, etc. in the prior art and can be obtained through regular commercial channels. The experimental methods, detection methods, etc. involved in the following embodiments, unless otherwise specified, are all conventional experimental methods, detection methods, etc. in the prior art.

[0047] Example 1

[0048] A transcript assembly method based on a reference genome, the steps are as follows:

[0049] (i) Extract the mRNA of the sample to be tested, amplify it, obtain circularized cDNA, and sequence it on a machine; use two alignment tools (Hisat2 and Star) to analyze the raw sequencing data obtained off the machine and obtain the results mapped back to the reference genome.

[0050] (II) Construct a splicing graph for each comparison tool’s reply results:

[0051] Based on the results of reads being mapped back to the reference genome, they are clustered into different gene loci, and exon-exon splice sites are derived from these spliced ​​reads;

[0052] Specifically, a traditional splicing graph G = (V, E) is constructed for each gene locus, where each node v corresponds to an exon and each edge e corresponds to a splicing site between two exons;

[0053] Furthermore, edges and nodes are weighted by the number of reads supporting them; edges or nodes with low weights that may be caused by sequencing errors or unreasonable alignments are removed.

[0054] (III) Merge the splicing graphs of each alignment tool to generate a label splicing graph, as follows:

[0055] (1) Assume that we have obtained the reply results of N comparison tools, and the splicing graph generated by comparison tool i is G i =(E i , V i ), i≤N, for each edge e i ∈E i , support e i The set of reads is R ei ; Label splice map G L =(E L , VL ) is obtained by merging each splice graph G i To generate, i∈[1,N], where the edges and points represent all the nodes that appear in the graph G i The edges and points in , i∈[1,N];

[0056] (2) Label the splicing graph G L :Through the following method, for each L For an edge e in , if there exists j∈[1,N] such that the edge Then mark the edge e as 1, otherwise mark it as 2; In addition, for every edge e∈E L Marked as 2, if there exist i and j such that R ei ≠R ej , i, j∈[1, N], and i≠N, label it as 2-, otherwise label it as 2+; then, the weight of the edge of the label splicing graph is assigned to each edge that appears in G i The average weight of the edges in , i≤N;

[0057] (IV) Extract the two-terminal circuit marked in the label diagram, such as Figure 1 As shown, the details are as follows:

[0058] (1) For each gene locus, generate a set of two-terminal paths for each alignment tool i, denoted as PP i ;

[0059] (1.1) For each pair of double-end reads r1 and r2, based on the response of alignment tool i, if L In the example, r1 spans a path p1=n j1 →n j2 →…→n jk , r2 spans a path p2 = n j’1 →n j’2 →…→n j’q , and p1 and p2 share a compatible subpath p in =n m1 →n m2 →…→n ms Satisfy k+qs≥3;

[0060] (1.2) Through the shared compatible subpath p in , connect the two paths p1 and p2 to generate a two-terminal path p;

[0061] (2) After processing all replies to G L After the paired-end reads are obtained, the set PP of all paired-end reads that depend on alignment tool i is obtained. i ;

[0062] (3) Along each path belonging to p∈PP i Paired-end reads set R ei The generated paired-end pathway set PP from each alignment tool i i Then, let PP L =U i∈[1,N] PP i The double-ended path of the label splicing graph is then labeled as PP L , marking method and marking edge e∈E L The method is the same as: for each L A two-terminal path p in a Then mark the two-terminal path p as 1, otherwise mark it as 2; In addition, for every two-terminal path p∈PP i Marked as 2, if there exist i and j such that R ei ≠R ej , i, j∈[1, N], and i≠j, mark it as 2-, otherwise mark it as 2+;

[0063] (V) Find the coverage of the path in the label splicing graph, and use the label-based dynamic programming path extraction algorithm, as follows:

[0064] (1) Select a seed and generate a subgraph using a seed:

[0065] A seed is an edge or a paired path in the labeled graph, which can be further extended to a path representing the full-length transcript. An unused paired path with the largest weight and labelled 2+ is selected as a seed. If such paired paths do not exist or they are all included in the assembled transcript, the seeds are selected in the following order: paired paths with a 2-label, edges with a 2+ label, edges with a 2-label, paired paths with a label of 1, or edges with a label of 1.

[0066] Paired-end paths or edges with label 2 represent that they are captured by multiple alignment tools, that is, they are most likely part of positive transcripts, and these seeds will be extended to full-length transcripts with a higher probability;

[0067] With a two-ended path S = v k →v k +1 →… →v l Or side S = v k →v l is the seed, extending all its right / left adjacent points v k (v l ), and extends until its adjacent point encounters a point without an outgoing / incoming edge, labeling the spliced ​​graph GL The subgraph G S is generated by the selected seed; then, the source / sink is added to the subgraph G S In the example, if all the nodes connected to it have no inbound / outbound edges, the weight of the newly added edge is set to ∞;

[0068] (2) Expanding the path through a label-based dynamic programming algorithm:

[0069] Traversing the graph G through a dynamic programming algorithm S The edges of the label clipping graph G S The idea of ​​dynamic programming is to find a path with a higher weight from the seed to each other adjacent edge;

[0070] Assume that graph G S The edges in are topologically sorted. Without loss of generality, they are numbered 1, 2...N. For each edge i=v→w, where v and w are two points in the label splicing graph, predecessor(i) is set to be the set of all incoming edges of point v, successor(i) is set to be the set of all outgoing edges of point w, and weight(i) is the coverage of edge i.

[0071] Assume that the edge v of seed S k →v k+1 corresponds to edge m, and edge v of seed S l-1 →v l Corresponding to edge n, for edge i from n to N (from m to l) there are two variables Exp i and Pre i (Suc i ), where Exp i represents the upper bound of the highest expression value of the path from seed to edge i, Pre i (Suc i ) is the predecessor / successor of edge i, if the seed is an edge then m=n; the dynamic programming algorithm is derived from Exp i =∞,Exp m =cov(S),Exp n = cov(S), and Pre i =nullFor i≥n, Suc i = null for each edge i≤m start;

[0072] For i from n+1 to N, if predecessor(i) is not empty, select label 2+Exp from predecessor(i) i’The largest edge i'. Note that if there is no edge with label 2+ or the coverage of the edge with label 2+ is 0, then select edge i' in the following order: edge with label 2-, edge with label 1; then set Exp i is min{weight(i), Exp i’} and assume Pre i For i from m-1 to 1, use the same method to set Exp i and Suc i , starting from the edge with the largest Exp value connecting the sink / source, backward / forward based on the label Pre(Suc) until reaching the seed S, then graph G S The most highly expressed transcript in h will be covered;

[0073] (3) Update the label splice map and repeat the process:

[0074] Define f min is the extended path p h The minimum edge weight on p h For each edge e in, update the weight weight(e) to weight(n)-f min Repeat (1) and (2) in the path extraction process until graph G L All edges in are covered by the predicted transcripts, and finally, a path coverage representing the transcripts is obtained, in which all double-end paths are covered by the assembled transcripts.

[0075] Application Example 1

[0076] We tested one simulated dataset and twelve real datasets (eight from Homo sapiens samples and four from mouse muscle samples). In this study, the reference genome and transcriptome versions of Homo sapiens and Mus musculus samples used were GRCh38 / hg38 and GRCm38 / mm10, respectively, downloaded from UCSC genomeBrowser. In this study, we focused the alignment task on two of the latest alignment tools, Hisat2 and Star, to generate alignments. In addition, we also used the merge function of samtools to generate a set of merged data based on the alignment results generated by the above two alignment tools. Then, the alignment results generated by the two alignment tools were used as input to run Tiglon. To evaluate its performance, we ran two of the best and representative assembly software, StringTie2 and scallop, based on the alignment results of Hisat2, the alignment results of Star, and the merged alignment results, respectively. We further ran the merge option in StringTie2, taking the transcripts generated by StringTie2 based on each alignment tool as input. The comparison standard commonly used in this study is that an assembled transcript is considered correct if and only if the intron chain of a reference transcript matches that of an assembled transcript exactly. In this process, we use the Cuffcompare tool in the Cufflinks software to detect correctly assembled transcripts. The accuracy of the assembly tool is measured by the number and accuracy of correctly assembled transcripts (the percentage of correctly assembled transcripts among all predicted transcripts).

[0077] We tested different samples separately. The basic process of the test is as follows:

[0078] 1. First, we evaluated the performance of the software on a simulated dataset. We used the RNA-seq data simulator RSEM to generate a simulated dataset containing about 52 million 100bp long paired-end reads. The parameters of the simulation model were learned from the real human RNA-seq dataset numbered SRR7536920 in NCBI SRA. Next, we used Hisat2 and Star to post the simulated reads to the reference genome. Subsequently, we used the merge function in samtools to generate a merged alignment based on the posting results of the two alignment tools mentioned above. We ran Tiglon using the alignment results generated by the above alignment tools as input. Similarly, we ran StringTie2 and Scallop using the alignment results generated by the above alignment tools as input. The results of each run are represented as: StringTie+Hisat2, StringTie2+Star, StringTie2+MergedAlignments, Scallop+Hisat2, Scallop+Star and Scallop+MergedAlignments. In addition, we ran the merge option in StringTie2, using the transcripts generated by StringTie2+Hisat2 and StringTie2+Star as input to generate a set of consistent transcripts, denoted as String2-Merge. The results after the test are as follows Figure 2 shown.

[0079] The test results of the simulated data set show that Tiglon performs best and can correctly reconstruct more transcripts while maintaining the highest accuracy. Specifically, Tiglon correctly assembled 22,928 transcripts, while StringTie2-merge had 20,630, StringTie2+MergedAlignments had 20,654, Scallop+MergedAlignments had 20,277, StringTie2+Hisat2 had 18,800, StringTie2+Star had 19,267, Scallop+Hisat2 had 18,786, and Scallop+Star had 19,463. Broadly speaking, using these two alignment tools strategically, Tiglon reconstructed 11.1% more correct transcripts than StringTie2-Merge, 11.01% more than StringTie2+MergedAlignments, 13.1% more than Scallop+MergedAlignments, and between 17.8% and 22% more correct transcripts than StringTie2 and Scallop using different alignment tools (e.g. Figure 2 A). The comparison results show that Tiglon also achieved the highest accuracy of 64.01%, while the accuracy of StingTie2-Merge was 56.40%, the accuracy of StringTie2+MergedAlignments was 56.87%, the accuracy of Scallop+MergedAlignments was 36.81%, the accuracy of StringTie2+Hisat2 was 58.94%, the accuracy of StringTie2+Star was 58.82%, the accuracy of Scallop+Hisat2 was 52.86%, and the accuracy of Scallop+Star was 41.60%. Therefore, whether it is the number of correctly assembled transcripts or the accuracy, Tiglon has shown superior performance among all compared assembly tools.

[0080] In addition, we further calculated the F-score, which is the harmonic mean of recall and precision (calculated as 2*precision*recall / (precision+recall)) to evaluate the overall performance of each assembly tool, where recall refers to the proportion of transcripts that the assembly tool can correctly reconstruct among known transcripts. In the simulated data set, Tiglon obtained an F-score of 0.480, which is a significant increase from StringTie2-Merge's 0.428, StringTie2+MergedAlignments' 0.430, Scallop+MergedAlignments' 0.353, StringTie2+Hisat2's 0.410, StringTie2+Star's 0.412, Scallop+Hisat2's 0.394, and Scallop+Star's 0.365, which shows that Tiglon has a better ability to balance recall and precision (such as Figure 2 B).

[0081] Especially for low-expressed transcripts, the correct reconstruction rate of Tiglon is 30.11% higher than that of StringTie2-merge, 19.08% higher than that of StringTie2+MergedAlignments, 58.06% higher than that of Scallop+MergedAlignments, and 40%-108% higher than that of StringTie2 and Scallop, which rely on the response results of different alignment tools for assembly.

[0082] 2. We also tested on real data sets, and the process was similar to testing the software performance on simulated data. Therefore, we only discuss the analysis of the results after the test. RNA-seq data from eight Homo sapiens samples and four mouse muscle samples were tested separately. The data obtained used all reference transcripts downloaded from the UCSC gene annotation GRCh38 as the true situation of the transcripts contained in the Homo sapiens samples, and GRCm38 as the true situation of the transcripts contained in the mouse muscle samples. Eight Homo sapiens RNA-seq samples H1-H8 and four mouse muscle RNA-seq samples M1-M4 were collected to evaluate the performance of the assembly software. All datasets were downloaded from NCBI Sequence Read Archive (SRA). The numbers of Homo sapiens samples are SRR307911, SRR387662, SRR10517380, ERR2403203, SRR307903, SRR315323, SRR315334 and SRR7536920, and the numbers of mouse muscle samples are DRR205674, DRR205677, ERR3320855 and ERR3320871.

[0083] 2.1. Testing on eight Homo sapiens sample datasets

[0084] The test results are as follows Figure 3 The results show that on the eight Homo sapiens sample data sets tested, Tiglon's number of correctly assembled transcripts and accuracy are significantly higher than all other assembly software. Specifically, the number of correctly reconstructed transcripts by Tiglon on the eight samples ranges from 13244 to 20489, with an average of about 17377, while the range of StringTie2-Merge is 12278-18112, with an average of 15831, and the range of StringTie2+MergedAlignments is 12429-17857, with an average of 15584. The range of Scallop+MergedAlignments is 101. The range is 11397-17371, with an average of 15251, the range of StringTie2+Hisat2 is 10369-16671, with an average of 13834, the range of StringTie2+Star is 11116-16469, with an average of 14303, the range of Scallop+Hisat2 is 9993-16531, with an average of 13802, and the range of Scallop+Star is 10653-16750, with an average of 14470. In general, on the eight tested datasets, Tiglon reconstructed 9.51% more correct transcripts than StringTie2-Merge, 11.25% more correct transcripts than StringTie2+MergedAlignments, 13.68% more correct transcripts than Scallop+MergedAlignments, 25.32% more correct transcripts than StringTie2+Hisat2, 21.21% more correct transcripts than StringTie2+Star, 25.61% more correct transcripts than Scallop+Hisat2, and 19.82% more correct transcripts than Scallop+Star (e.g. Figure 3 A). In terms of accuracy, Tiglon still maintains the highest accuracy on all test sets, reaching 41.27%, among which StringTie2-Merge is 35.89%, StringTie2+MergedAlignment is 32.03%, Scallop+MergedAlignment is 26.90%, StringTie2+Hisat2 is 36.31%, StringTie2+Hisat2 is 36.31%, Scallop+Hisat2 is 30.47%, and Scallop+Star is 29.77% (e.g. Figure 3 A). We further calculated the F-score, and Tiglon still maintained the highest performance. On average, on eight datasets, Tiglon's F-score reached 0.1701, which was about 10.67%-28.08% higher than other assembly methods (e.g. Figure 3 B). After that, we evaluated the ability of the assembly software to reconstruct transcripts on samples with different expression levels. Although we could not know the true expression values ​​of these samples, we used the famous abundance estimation software Kallistor to evaluate the abundance of these eight RNA-seq samples. As we did on the simulated data set, we divided these eight samples into four expression levels: low, medium, and high according to the estimated abundance. After testing, we found that Tiglon still achieved the highest value of the number of correctly reconstructed transcripts on samples with different expression levels, such as Figure 3 As shown in C.

[0085] 2.2. Testing on four mouse muscle sample datasets

[0086] like Figure 4 The results show that Tiglon correctly assembled 14450-17683 transcripts, with an average of 16153, which is 4.74% higher than StingTie2-merge, 4.28% higher than StringTie2+MergedAlignments, 14.01% higher than Scallop+MergedAlignments, and 12.72%-17.51% higher than StingTie2 and Scallop. In terms of accuracy, Tiglon performed well again. Specifically, Tiglon achieves an average accuracy of 51.71%, ranging from 50.44% to 52.73% in the four tested samples, StringTie2-Merge has an average accuracy of 44.42%, StringTie2+MergedAlignments is 37.75%, Scallop+MergedAlignments is 31.11%, StringTie2+Hisat2 is 44.31%, StringTie2+Star is 43.64%, Scallop+Hisat2 is 35.27%, and Scallop+Star is 34.16%. Overall, Tiglon's accuracy is improved by an average of 16.40%-66.19% compared to other methods (such as Figure 4A). In terms of F-score, Tiglon achieves 0.2341, which is 7.45% higher than StringTie-Merge, 12.19% higher than StringTie2+MergedAlignments, 21.85% higher than Scallop+MergedAlignments, and 14.17%-24.86% higher than StringTie2 and Scallop (e.g. Figure 4 B). Similar to the method used for testing Homo sapiens samples, we evaluated the performance of reconstructing transcripts at different expression levels. After comparison, we found that Tiglon outperformed other assembly software in all four test samples regardless of the expression level, especially for transcripts at low expression levels. Tiglon improved by 25.36%-39.91% over StringTie-Merge, 19.65%-45.05% over StringTie2+MergedAlignments, 55.40%-96.11% over Scallop+MergedAlignments, 67.74%-98.69% over StringTie2+Hisat2, 56.12%-79.09% over StringTie2+Star, 126.85%-154.34% over Scallop+Hisat2, and 108.64%-139.58% over Scallop+Star (e.g. Figure 4 C).

[0087] Through the above tests, we found that Tiglon's results are significantly better than the best assembly tools currently available, including StringTie2 and Scallop, which rely on Hisat2 or Star or merge the alignment results of these two alignment tools, and StringTie2-Merge (the merge option in StringTie2). For example, on eight tested Homo sapiens datasets, Tiglon reconstructed 9.51% more correct transcripts on average than StringTie2-Merge, 11.25%-13.68% more correct transcripts than StringTie2 and Scallop when relying on merged alignment results, and 19.82%-25.61% more correct transcripts than StringTie2 and Scallop when relying on Hisat2 or Star alignment results. More importantly, Tiglon showed a clear advantage in recovering low-expression transcripts. In eight Homo sapiens samples, Tiglon's correct reconstruction rate for low-expressed transcripts was 22.57%-53.48% higher than StirngTie2-Merge, 29.73%-89.05% higher than Scallop and StringTie2 that relied on merged alignment results, and 51.93%-161.18% higher than StringTie2 and Scallop that relied on the results of only one alignment tool.

[0088] The above examples are provided to those skilled in the art to fully disclose and describe how to make and use the claimed embodiments, and are not intended to limit the scope of the disclosure herein. Modifications that are obvious to those skilled in the art will fall within the scope of the appended claims.

Claims

1. A transcript assembly method based on a reference genome, characterized in that: The following steps are involved: (1) extracting mRNA of the sample to be tested, amplifying it, obtaining circularized cDNA, and sequencing it; using at least two alignment tools to analyze the raw sequencing data obtained from the machine, and obtaining the result mapped back to the reference genome; (II) Construct a splicing graph for each comparison tool’s reply results: Based on the results of reads being mapped back to the reference genome, they are clustered into different gene loci, and exon-exon splice sites are derived from these spliced ​​reads; (III) Merge the splicing graphs of each alignment tool to generate a label splicing graph, as follows: (1) Assume that we have obtained the reply results of N comparison tools, and the splicing graph generated by comparison tool i is G i =(E i , V i ), i≤N, for each edge e i ∈E i , support e i The set of reads is R ei ; Label splice graph G L =(E L , V L ) is obtained by merging each splice graph G i To generate, i∈[1,N], where the edges and points represent all the nodes that appear in the graph G i The edges and points in , i∈[1,N]; (2) Label the splicing graph G L :Through the following method, for each L For an edge e in the equation, if there exists j∈[1,N] such that the edge Then mark the edge e as 1, otherwise mark it as 2; In addition, for every edge e∈E L Marked as 2, if there exist i and j such that R ei ≠R ej , i, j ∈ [1, N], and i ≠ j, label it as 2-, otherwise label it as 2+; then, the weight of the edge of the label splicing graph is assigned to each edge that appears in G i The average weight of the edges in , i≤N; (IV) Extract the two-terminal circuit marked in the label diagram, as follows: (1) For each gene locus, generate a set of two-terminal paths for each alignment tool i, denoted as PP i ; (1.1) For each pair of double-end reads r1 and r2, based on the response of alignment tool i, if L In the example, r1 spans a path p1=n j1 →n j2 →…→n jk , r2 spans a path p2 = n j′1 →n j′2 →…→n j′q , and p1 and p2 share a compatible subpath p in =n m1 →n m2 →…→n ms Satisfy k+qs≥3; (1.2) Through the shared compatible subpath p in , connect the two paths p1 and p2 to generate a two-terminal path p; (2) After processing all replies to G L After the paired-end reads are obtained, the set PP of all paired-end reads that depend on alignment tool i is obtained. i ; (3) Along each path belonging to p∈PP i Paired-end reads set R ei The generated paired-end pathway set PP from each alignment tool i i Then, let PP L =∪ i∈[1，N] PP i The double-ended path of the label splicing graph is then labeled as PP L , marking method and marking edge e∈E L The method is the same as: for each L A two-terminal path p in a Then mark the two-terminal path p as 1, otherwise mark it as 2; In addition, for every two-terminal path p∈PP i Marked as 2, if there exist i and j such that R ei ≠R ej , i, j∈[1, N], and i≠j, mark it as 2-, otherwise mark it as 2+; (V) Find the coverage of the path in the label splicing graph, and use the label-based dynamic programming path extraction algorithm, as follows: (1) Select a seed and generate a subgraph using a seed: Select an unused, weighted paired path with label 2+ as a seed; if such paired paths do not exist or they are all already included in the assembled transcripts, then select seeds in the following order: paired paths with 2-label, edges with 2+ label, edges with 2-label, paired paths with label 1, or edges with label 1; With a two-ended path S = v k →v k +1 →… →v l Or side S = v k →v l is the seed, extending all its right / left adjacent points v k (v l ), and extends until its adjacent point encounters a point without an outgoing / incoming edge, labeling the spliced ​​graph G L The subgraph G S is generated by the selected seed; then, the source / sink is added to the subgraph G S In the example, if all the nodes connected to it have no inbound / outbound edges, the weight of the newly added edge is set to ∞; (2) Expanding the path through a label-based dynamic programming algorithm: Traversing the graph G through a dynamic programming algorithm S The edges of the label clipping graph G S Find the transcript with the highest expression level in Assume that graph G S The edges in are topologically sorted. Without loss of generality, they are numbered 1, 2...N. For each edge i=v→w, where v and w are two points in the label splicing graph, predecessor(i) is set to be the set of all incoming edges of point v, successor(i) is set to be the set of all outgoing edges of point w, and weight(i) is the coverage of edge i. Assume that the edge v of seed S k →v k+1 corresponds to edge m, and edge v of seed S l-1 →v l Corresponding to edge n, for edge i from n to N (from m to l) there are two variables Exp i and Pre i (Suc i ), where Exp i represents the upper bound of the highest expression value of the path from seed to edge i, Pre i (Suc i ) is the predecessor / successor of edge i, if the seed is an edge then m=n; the dynamic programming algorithm is derived from Exp i =∞,Exp m =cov(S),Exp n = cov(S), and Pre i =nullFor i≥n, Suc i = null for each edge i≤m start; For i from n+1 to N, if predecessor(i) is not empty, select label 2+Exp from predecessor(i) i′ The largest edge i'. Note that if there is no edge with label 2+ or the coverage of the edge with label 2+ is 0, then select edge i' in the following order: edge with label 2-, edge with label 1; then set Exp i is min{weight(i), Exp i’ } and assume Pre i For i from m-1 to 1, use the same method to set Exp i and Suc i , starting from the edge with the largest Exp value connecting the sink / source, backward / forward based on the label Pre(Suc) until reaching the seed S, then graph G s The most highly expressed transcript in h will be covered; (3) Update the label splice map and repeat the process: Define f min is the extended path p h The minimum edge weight on p h For each edge e in, update the weight weight(e) to weight(n)-f min Repeat (1) and (2) in the path extraction process until graph G L All edges in are covered by the predicted transcripts, and finally, a pathway coverage representing the transcripts is obtained, in which all double-end pathways are covered by the assembled transcripts.

2. The method for assembling transcripts based on a reference genome according to claim 1, characterized in that: In the step (1), the comparison tool is selected from Hisat2 and Star.

3. The method for assembling transcripts based on a reference genome according to claim 1, characterized in that: In step (2), a traditional splicing graph G=(V, E) is constructed for each gene locus, wherein each node v corresponds to an exon, and each edge e corresponds to a splicing site between two exons.

4. The method for assembling transcripts based on a reference genome according to claim 1, characterized in that: In step (2), edges and nodes are weighted by the number of reads supporting them; low-weighted edges or nodes that may be caused by sequencing errors or unreasonable alignments are deleted.

Images