A method and system for filtering full-length transcripts

Abstract

The application discloses a method and system for filtering full-length transcripts, which comprises the following steps: aligning sequencing obtained transcript reads with a reference genome, clustering and deduplicating according to transcript splicing patterns, filtering according to average sequencing depth, exon support number and transcript read support number, re-aligning with the reference genome, and clustering and deduplicating again. The system is used for executing the filtering method. The method and system provided by the application can be used for processing transcripts with unknown splicing patterns in the reference genome, effectively retaining low-abundance reliable transcripts, and the filtering result has higher sensitivity and accuracy compared with prior art.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics, specifically to a method and system for filtering full-length transcripts. Background Technology

[0002] Traditional gene expression studies are conducted at the individual or tissue level, treating the millions or more cells within an individual or tissue sample as a whole. However, it is clear that different types of cells within the same sample, and cells of the same type at different stages, exhibit different gene expression patterns. These differences primarily lie in which genes are expressed, how the same gene is expressed (transcript differences), and the amount of each gene expressed (expression level differences). Traditional, averaged gene expression studies cannot adequately address this technical challenge, thus leading to the development of single-cell technology.

[0003] In 2009, the first single-cell transcriptome sequencing technology was successfully developed (Tang F, Barbacioru C, Wang Y, et al. mRNA-Seq whole-transcriptome analysis of a single cell[J]. NATUREMETHODS, 2009, 6(5):377-382.), followed by Smart-Seq technology ( Daniel, Luo, Shujun, Wang, Yu-Chieh, et al. Full-length mRNA-Seq from single-cell levels of RNA and individual circulating tumor cells[J]. Nature Biotechnology, 2012, 30(8):777-782.), the development of Smart-Seq2 technology (Serra, Lorrayne, Chang, Dennis, Macchietto, Marissa, et al. Adapting the Smart-seq2 Protocol for Robust Single Worm RNA-seq[J]. BioProtocol, 2018, 8(4).), and the launch of commercial platforms such as 10×Genomics Chromium (https: / / www.10xgenomics.com / ) and Illumina Bio-Rad Single-Cell Sequencing Solution (https: / / www.bio-rad.com / ) have led to the rapid development of single-cell research.

[0004] Currently, most mature single-cell research is conducted in conjunction with next-generation sequencing (NGS) technology. NGS is characterized by short read lengths, high accuracy, high throughput, and low cost. NGS data (short reads) can be compared at the gene level, determining the relationship between each short read and a specific gene. However, because each gene exhibits multiple splicing patterns during transcription, producing different transcripts, the short read length limitation prevents direct determination of the relationship between each individual short read and a transcript. Therefore, NGS technology can only utilize bioinformatics algorithms to assemble data and obtain transcript information, resulting in very limited applicability in alternative splicing patterns and transcript classification.

[0005] Third-generation sequencing platforms, represented by Pacific Biosciences and Oxford Nanopore Technologies, are characterized by their long read lengths, bringing significant changes to the sequencing field. As an emerging application, the combination of third-generation sequencing and single-cell technology still faces many technical challenges. Long-read sequencing has a higher error rate and lower sequencing depth than short-read sequencing. In data analysis, using traditional "one-size-fits-all" methods will inevitably fail to obtain low-depth data that accurately reflects the true transcriptional processes of organisms. This data is crucial for cellular-level research and is likely related to single-cell heterogeneity. Summary of the Invention

[0006] In view of this, the present invention provides a method and system for filtering full-length transcripts, which achieves accurate filtering of full-length transcript data and has high accuracy and sensitivity.

[0007] To achieve the above objectives, the technical solution of the present invention is as follows:

[0008] A method for filtering full-length transcripts includes the following steps:

[0009] S1. The transcript reads obtained from sequencing are compared with the reference genome. Based on the transcript splicing pattern, clustering and redundancy removal are performed to obtain the first transcript. The first transcript includes a set A of transcripts with the same splicing pattern as the reference genome and a set B of transcripts with different splicing patterns than the reference genome.

[0010] S2. Based on the average sequencing depth, exon support number, and transcript read support number, filter transcript set B to obtain transcript set C;

[0011] S3. Take the union of transcript set A and transcript set C, re-align it with the reference genome, and cluster again to remove redundancy, to obtain the second transcript, which is the filtered full-length transcript.

[0012] In the above technical solution, it can be understood that the first transcript is the union of transcript set A and transcript set B; the reference genome includes sequence information and structural annotation information; after the transcript reads are aligned with the reference genome, the transcript splicing pattern of each read can be obtained.

[0013] Preferably, the sequencing in step S1 is third-generation sequencing, and the sample type of the transcript read is single cell.

[0014] Preferably, the clustering and redundancy removal method described in step S1 and step S3 is the same, specifically: based on the complete matching of the 3' exons and splicing sites, transcript reads with the same splicing pattern are clustered into one class, and only the longest transcript with a single splicing pattern is retained. In the above technical solution, since there are many similar transcripts under a gene model, similar transcripts are clustered based on the complete matching of the 3' exons and splicing sites; moreover, since RNA is prone to degradation at the 5' end, transcripts with only 5' differences are grouped into one class in this invention.

[0015] Preferably, the method for counting the number of exons in step S2 is as follows: based on the similarity between exons with overlapping relationships in different transcripts within set B, determine whether they are the same exons and count them.

[0016] More preferably, exons are considered to be the same exons when the similarity between them is greater than or equal to 70%.

[0017] The preferred method for calculating the similarity between exons is as follows: a window region is set and segmented to obtain segmentation windows; the segmentation windows are assigned values ​​and normalized; the sum of the assigned values ​​of the segmentation windows into which each exon falls is counted as the score of that exon; and the similarity between exons is obtained by comparing the scores of each exon.

[0018] More preferably, the starting point of the window region is the upstream 5' end of the exon with the overlapping alignment relationship, and the ending point of the window region is the downstream 3' end of the exon with the overlapping alignment relationship; and the logarithmic value assigned to the segmentation window is uniformized. The window region, segmentation window, and exon are all represented in the direction conventionally used in the art, from 5' to 3'.

[0019] Preferably, the transcript set B filtering method in step S2 specifically includes:

[0020] 1) For transcripts of known genes in the reference genome but with different splicing patterns, each exon in the transcript should be retained if it meets the following condition:

[0021] If an exon overlaps with an exon in the reference genome, the number of exon supports is less than the first threshold and the number of transcript read supports is greater than or equal to the second threshold, or the number of exon supports is greater than or equal to the first threshold.

[0022] If the exon does not overlap with the exon in the reference genome, the exon length is less than the third threshold and the number of exon supports is greater than or equal to the fourth threshold, or the exon length is greater than or equal to the third threshold and the number of exon supports is greater than or equal to the fifth threshold.

[0023] 2) For transcripts of unknown genes in the reference genome, transcripts with an average sequencing depth ≥ the sixth threshold are retained.

[0024] More preferably, the first threshold is 2, the second threshold is 2, the third threshold is 200bp, the fourth threshold is 10, the fifth threshold is 5, and the sixth threshold is 1000.

[0025] In the above technical solution, the average sequencing depth in the transcript set B is the sum of the sequencing depths of each base contained in the transcript divided by the total number of bases in the transcript, the exon support number is the total number of times the exon appears in all transcripts in the transcript set B, and the transcript read support number is the number of transcript reads under each transcript in the transcript set B.

[0026] Preferably, the method further includes step S4, which involves classifying the second transcript based on the splicing site.

[0027] More preferably, the classification types are: fully matched transcripts, partially matched transcripts, transcripts containing recombination of known splice sites, transcripts containing new splice sites, antisense transcripts, fusion transcripts, transcripts of intron regions, transcripts of intergenic regions, and transcripts spanning intron and exon regions. Wherein, a fully matched transcript refers to a transcript whose splice site is completely consistent with the known annotation information in the reference genome; a partially matched transcript refers to a transcript whose splice site is only partially identical to the known annotation information in the reference genome; a transcript containing recombination of known splice sites refers to a transcript where all splice sites are known, but new permutations and combinations are formed between the splice sites, resulting in different transcripts; a transcript containing a new splice site refers to a transcript with a new splice site; an antisense transcript refers to a transcript whose sequence direction is opposite to that of the reference genome; a fusion transcript refers to a transcript that spans two different genes; an intron region transcript refers to a transcript that is completely aligned to an intron region of the genome; an intergenetic region transcript refers to a transcript located in a region between different genes; and a transcript spanning both intron and exon regions refers to a transcript that is half located in an intron region and half in an exon region.

[0028] The present invention also provides a system for filtering full-length transcripts, comprising:

[0029] Alignment module: used to align transcript reads with the reference genome to obtain the transcript splicing pattern of each read;

[0030] Clustering and redundancy removal module: Based on transcript splicing patterns, reads with the same splicing pattern are clustered into one class and only the longest unique transcript in that class is retained;

[0031] Transcript filtering module: Using average sequencing depth, exon support number, and transcript read support number as parameters, transcripts with new splicing patterns are filtered; the transcripts with new splicing patterns are those whose splicing patterns are not reported in the reference genome after clustering and redundancy removal, and the "not reported" means that at least one splicing pattern is different from the splicing patterns known in the reference genome.

[0032] Filtered Transcript Output Module: Outputs the filtered transcripts.

[0033] The beneficial effects of this invention are as follows: The full-length transcript filtering method provided by this invention can specifically process transcripts with unknown splicing patterns in the reference genome, effectively retaining low-abundance reliable transcripts; using average sequencing depth, exon support number, and transcript read support number as parameters, it can filter transcripts with new splicing patterns, and the filtering results have higher sensitivity and accuracy compared with existing technologies; in practical applications, the filtering results can be optimized by controlling the similarity between exons. Attached Figure Description

[0034] Figure 1 This is a schematic diagram of the structure of the segmented window in Example 1 (in the figure, the double vertical lines indicate the omitted parts);

[0035] Figure 2 A comparison diagram showing the technical effects of the filtering method provided by this invention with existing technologies. Detailed Implementation

[0036] To better understand the present invention, the following detailed description is provided in conjunction with embodiments.

[0037] It should be noted that the terms "A", "B", "first", "second", etc. in this invention are used to distinguish different objects, rather than to describe a specific order; "set A" and "set B" in this invention contain at least one full-length transcript.

[0038] Example 1

[0039] Single-cell cDNA sequencing of human cell line samples was performed using the ONT sequencing platform. The sequencing data volume is shown in the table below.

[0040]

[0041] The sequencing data is filtered using the method and system of this invention, and the specific steps are as follows:

[0042] (1) The sequencing data above were compared with the reference genome to obtain the transcript splicing pattern of each transcript sequencing read.

[0043] In the specific implementation process, alignment software such as Minimap2 and Gmap can be used. The alignment with the highest score is selected as the final alignment result, while alignments containing chimeric sequences are filtered out. Chimeric sequences refer to non-naturally occurring DNA molecules formed when two DNA molecules from different locations in the genome are joined together due to external factors (such as library preparation) during sequencing. In the alignment results, this is reflected as a single read aligning to two or more locations in the genome.

[0044] (2) Based on the transcript splicing pattern, the sequencing data are clustered to remove redundancy, and the first transcript is obtained. The specific process of clustering to remove redundancy is as follows: based on the complete match between the 3' exon and the splicing site, transcript reads with the same splicing pattern are clustered into one class, and only the longest transcript with a single splicing pattern is retained.

[0045] Based on whether the transcript splicing pattern is known, the first transcript is divided into transcript set A and transcript set B. Transcripts in transcript set A have splicing patterns identical to those reported in the reference genome, while transcripts in transcript set B have splicing patterns different from those reported in the reference genome.

[0046] (3) Based on the average sequencing depth, exon support number and transcript read support number, filter transcript set B to obtain transcript set C.

[0047] First, the overlapping exons in transcript set B are identified. The specific process is as follows:

[0048] ① Set a window region, where the starting point of the window region is the uppermost 5' end of the exon with overlapping alignment relationship, and the ending point of the window region is the lowermost 3' end of the exon with overlapping alignment relationship;

[0049] ② Divide the window area into segmented windows; in this embodiment, the window area is divided into 200 segmented windows, such as... Figure 1 As shown;

[0050] ③ Assign values ​​to the segmentation windows and normalize them. The sum of the values ​​assigned to the segmentation windows into which each exon falls is the score for that exon. In this embodiment, the segmentation windows are assigned values ​​starting from 2, and the logarithmic values ​​of the assigned values ​​are normalized to achieve the following: Figure 1 Taking the four transcripts a, b, c, and d, which have overlapping relationships, as examples, the scores of each exon are as follows:

[0051] score(a)A=log23+log24+log25…+log2201;

[0052] score(b)=log22+log23+log24…+log2157;

[0053] score(c)=log23+log24+log25…+log279;

[0054] score(d)=log22+log23+log24…+log278;

[0055] The scores above are compared pairwise, and those with a difference of less than 30% (i.e., similarity greater than or equal to 70%) are recorded as the same exon.

[0056] Secondly, record the number of times each exon appears in transcript set B, denoted as the support number of that exon. The higher the support number of an exon, the more reliable that exon is in the transcript.

[0057] The process of filtering transcript set B is as follows:

[0058] ① For transcripts of known genes in the reference genome but with different splicing patterns, each exon in the transcript should satisfy the following condition:

[0059] If an exon overlaps with an exon in the reference genome, the exon support number is <2 and the transcript read support number is ≥2, or the exon support number is ≥2.

[0060] If the exon does not overlap with the exon in the reference genome, the exon length is <200bp and the exon support number is ≥10, or the exon length is ≥200bp and the exon support number is ≥5;

[0061] ② For transcripts of unknown genes in the reference genome, transcripts with an average sequencing depth ≥1000 are retained.

[0062] The set of transcripts obtained after filtering is denoted as transcript set C.

[0063] (4) Merge set A and set C, and use the same method as in steps (1) and (2) on the union to obtain the second transcript, which is the filtered full-length transcript.

[0064] It should be noted that the number of window area segments, the method of assigning values ​​to the segmented windows, the specific method of normalization, the base of the logarithmic value, and the specific values ​​of each indicator in the filtering conditions in this embodiment are all specific examples and do not limit the scope of protection of this invention.

[0065] Example 2

[0066] The second transcript obtained in Example 1 was classified, specifically using a method based on splicing sites. For a detailed implementation process of this example, please refer to patent 201910388054.2.

[0067] Using this classification method, the second transcript in Example 1 can be divided into the following nine types, specifically:

[0068] A perfectly matched transcript: a transcript whose splicing site is exactly the same as the known annotation information in the reference genome;

[0069] Partially matched transcripts: The splicing sites only share partial sequence information with the known annotation information in the reference genome;

[0070] Transcripts containing recombination of known splice sites: Each splice site is known, but different transcripts are generated due to the new permutations and combinations formed between the splice sites;

[0071] Transcripts containing new splicing sites: Transcripts containing new splicing sites;

[0072] Antisense transcripts: Transcripts whose sequence orientation is opposite to that of the reference genome;

[0073] Fusion transcripts: transcripts that span two different genes;

[0074] Transcripts from intron regions: Transcripts that have been fully aligned to intron regions of the genome;

[0075] Intergenic transcripts: Transcripts located in regions between different genes;

[0076] Transcripts spanning both intron and exon regions: transcripts that are half located in the intron region and half in the exon region.

[0077] Transcripts “containing known splicing sites recombined” and “containing new splicing sites” in the classification results are classified as new transcripts.

[0078] Comparative Example 1

[0079] Using the same data as in Example 1, filtered full-length transcript data were first obtained using the software Flair (https: / / github.com / BrooksLabUCSC / flair), and then the transcripts were classified using the method in Example 2.

[0080] Comparative Example 2

[0081] Using the same data as in Example 1, the filtered full-length transcript data were first obtained using the ONT official software pinfish (https: / / github.com / nanoporetech / pipeline-pinfish-analysis), and then the transcripts were classified using the method in Example 2.

[0082] Comparative Example 3

[0083] Unlike Example 1, in step (3), when comparing exon scores, those with a difference of less than 40% (i.e., similarity greater than or equal to 60%) are recorded as the same exon; the other processes are completely consistent with Example 1; then, the transcripts are classified using the method of Example 2.

[0084] The accuracy and sensitivity of the classification results from Example 2, Comparative Example 1, Comparative Example 2, and Comparative Example 3 were evaluated to verify the filtering effect of the full-length transcript data. Sensitivity refers to the proportion of reported transcripts identified out of the total reported transcripts in the reference genome, reflecting the filtering method's ability to discover transcripts; accuracy refers to the proportion of reported transcripts identified out of the total identified transcripts, reflecting the filtering method's ability to correctly identify transcripts.

[0085] The classification results of Example 2 are compared with those of Comparative Examples 1 and 2 in the table below:

[0086]

[0087] As shown in the table above, Example 2 exhibits better accuracy and sensitivity than Comparative Example 1, while its accuracy is slightly lower than Comparative Example 2, but its sensitivity is significantly higher. A screenshot of the gene analysis results is provided for visual visualization, as shown below. Figure 2As shown in the figure (where filled squares represent exons, lines represent introns, and arrows represent transcription directions), 1 represents the reported transcripts in the human reference genome hg38; 2 represents the transcripts obtained in Comparative Example 2, which only yielded one transcript, indicating poor sensitivity; 3 represents the transcripts obtained in Comparative Example 1, which, although numerous, clearly shows many erroneous transcripts (with a different splicing pattern compared to 1); 4 represents the transcripts obtained in Example 2, which not only have a large number of transcripts but also lack erroneous transcripts, demonstrating excellent performance in both accuracy and sensitivity. Therefore, the filtering method described in this invention can obtain more accurate, reliable, and comprehensive full-length transcript data.

[0088] The classification results of Example 2 and Comparative Example 3 are compared in the table below:

[0089]

[0090]

[0091] Therefore, adjusting the similarity of exons with overlapping relationships from 70% to 60% will decrease both sensitivity and accuracy.

[0092] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method of filtering full-length transcripts, characterized by, Includes the following steps: S1. The transcript reads obtained from sequencing are compared with the reference genome. Based on the transcript splicing pattern, clustering and redundancy removal are performed to obtain the first transcript. The first transcript includes a set A of transcripts with the same splicing pattern as the reference genome and a set B of transcripts with different splicing patterns than the reference genome. S2. Based on the average sequencing depth, exon support number, and transcript read support number, filter set B to obtain set C; S3. Take the union of set A and set C, re-align it with the reference genome, cluster again to remove redundancy, and obtain the second transcript. The filtering method for set B in step S2 is as follows: 1) For transcripts of known genes in the reference genome but with different splicing patterns, each exon in the transcript should be retained if it meets the following condition: If an exon overlaps with an exon in the reference genome, the exon support number is <2 and the transcript read support number is ≥2, or the exon support number is ≥2. If the exon does not overlap with the exon in the reference genome, the exon length is <200bp and the exon support number is ≥10, or the exon length is ≥200bp and the exon support number is ≥5; 2) For transcripts of unknown genes in the reference genome, transcripts with an average sequencing depth ≥1000 are retained; The method for counting the number of exons in step S2 is as follows: based on the similarity between exons with overlapping relationships in different transcripts within set B, it is determined whether they are the same exons and counted. When the similarity between the exons is greater than or equal to 70%, they are determined to be the same exons.

2. The method of filtering full-length transcripts according to claim 1, wherein, The clustering and redundancy removal method described in step S1 is the same as that described in step S3. Specifically, based on the complete matching of the 3' end exons and splicing sites, transcript reads with the same splicing pattern are clustered into one class, and only the longest transcript with a single splicing pattern is retained.

3. The method of filtering full-length transcripts according to claim 1, wherein, The method for calculating the similarity between exons is as follows: a window region is set and segmented to obtain segmentation windows. The segmentation windows are assigned values ​​and normalized. The sum of the assigned values ​​of the segmentation windows into which each exon falls is counted as the score of that exon. The similarity between exons is obtained by comparing the scores of each exon.

4. The method of filtering full-length transcripts according to claim 3, wherein, The starting point of the window region is the upstream 5' end of the exon with overlapping alignment, and the ending point of the window region is the downstream 3' end of the exon with overlapping alignment.

5. The method of filtering full-length transcripts according to claim 3, wherein, The logarithmic values ​​assigned to the segmentation window are then normalized.

6. The method of filtering full-length transcripts of claim 1, wherein, The method further includes the following steps: S4. Classify the second transcript based on the splicing site.

7. A system for filtering full-length transcripts, characterized in that, include: Alignment module: used to align transcript reads with the reference genome to obtain the transcript splicing pattern of each read; Clustering and redundancy removal module: Based on transcript splicing patterns, reads with the same splicing pattern are clustered into one class and only the longest unique transcript in that class is retained; Transcript filtering module: Filters transcripts with new splicing patterns using parameters such as average sequencing depth, exon support number, and transcript read support number; Filtered Transcript Output Module: Outputs the filtered transcripts; The method for filtering transcripts with new splicing patterns is as follows: 1) For transcripts of known genes in the reference genome but with different splicing patterns, each exon in the transcript should be retained if it meets the following condition: If an exon overlaps with an exon in the reference genome, the exon support number is <2 and the transcript read support number is ≥2, or the exon support number is ≥2. If the exon does not overlap with the exon in the reference genome, the exon length is <200bp and the exon support number is ≥10, or the exon length is ≥200bp and the exon support number is ≥5; 2) For transcripts of unknown genes in the reference genome, transcripts with an average sequencing depth ≥1000 are retained; The method for counting the number of exons is as follows: based on the similarity between exons with overlapping relationships in different transcripts within set B, it is determined whether they are the same exons and the count is performed. When the similarity between the exons is greater than or equal to 70%, they are determined to be the same exons.

Images