Single cell sequencing-based transcriptional start site quantitative analysis system

Abstract

The invention relates to the field of single cell sequencing, in particular to a single cell sequencing-based transcriptional start site quantitative analysis system, which comprises a data acquisition module, a gene determination module, a sequence comparison module, an analysis module and a corrected PSI calculation module. The system provided by the invention can be directly butted with a single cell conventional analysis process, and single cell level quantification, difference analysis and visualization of a selective transcription start site are realized. According to the practical application of the system on physiological and pathological samples, a large number of high-quality selective transcriptional start sites are successfully identified, and the sites not only help cell clustering analysis on the single cell level, but also reveal how the dynamic expression profiles of selective transcriptional start site events affect physiological and pathological transcriptomes.

Description

Technical Field

[0001] The present invention relates to the field of single-cell sequencing, and in particular to a transcription start site quantitative analysis system based on single-cell sequencing. Background Art

[0002] Alternative splicing (AS) is a crucial mechanism of gene transcriptional regulation, encompassing seven main pathways: exon skipping, exon mutual exclusion, exon 5' alternative splicing, exon 3' alternative splicing, alternative transcription start sites (ATS), alternative polyadenylation, and intron retention. Alternative splicing refers to the process by which precursor RNA (pre-RNA) undergoes different splicing pathways (selecting different splice site combinations) to produce mRNA splice isoforms. Among these, alternative transcription start sites, as a form of alternative splicing, regulate mRNA translation efficiency and protein production by altering the open reading frame (ORF), sequence motifs, and structural motifs within the 5' untranslated region (UTR), thereby increasing the complexity of the eukaryotic transcriptome and proteome. It has been reported that 40–50% of mammalian genes possess alternative transcription start sites. Furthermore, alternative transcription start sites play crucial roles in cell differentiation and organ development, representing a crucial gene regulatory mechanism in eukaryotes. Promoter selectivity in many known and potential cancer-related genes is associated with patient survival. At the same time, more and more studies have shown that abnormalities in selective transcription start sites are closely related to the occurrence of tumors such as leukemia, lung cancer and breast cancer.

[0003] With the rapid development of life science research methods, the emergence and advancement of single-cell sequencing technology has become one of the most representative emerging technologies. Single-cell transcriptome sequencing, which quantifies gene expression in individual cells, is widely used in research fields such as immunology, embryonic development, oncology, and COVID-19. Compared with traditional bulk RNA-seq, single-cell transcriptome sequencing can deeply analyze cellular heterogeneity during tissue development and disease progression, thereby better elucidating the relevant biological mechanisms.

[0004] Although some quantitative analysis methods for selective transcription start sites based on single-cell RNA sequencing (such as SCAFE and CamoTSS) have been developed, there are problems such as insufficient accuracy of quantitative results, which limits their practical application. Summary of the Invention

[0005] The present invention provides a transcription start site quantitative analysis system based on single-cell sequencing, comprising: The data acquisition module is configured to acquire single-cell RNA sequencing data of target cells.

[0006] In some embodiments, the single-cell RNA sequencing is 5'-end single-cell sequencing.

[0007] In some embodiments, the single-cell RNA sequencing data includes paired-end sequencing data of multiple sequencing fragments.

[0008] In some embodiments, the double-end sequencing data includes 5'-end complementary sequence data and 3'-end complementary sequence data corresponding to each sequencing fragment.

[0009] In some embodiments, the 5' complementary sequence is the base sequence data obtained by sequencing from the 5' end of the sequencing fragment. In some embodiments, the 3' complementary sequence is the base sequence data obtained by sequencing from the 3' end of the sequencing fragment.

[0010] In some embodiments, the data acquisition module may also be configured to extract high-quality sequencing fragment data from the single-cell RNA sequencing data.

[0011] In some embodiments, the system further comprises: a gene determination module configured to obtain determined genes desired to be analyzed.

[0012] In some embodiments, the expression level of the gene to be analyzed should be greater than a set first expression level threshold. In some embodiments, the first expression level threshold may be ≥200.

[0013] In some embodiments, the system further comprises: a sequence alignment module configured to retrieve the single-cell RNA sequencing data for alignment according to the gene desired to be analyzed.

[0014] In some embodiments, the sequence alignment module may include: A first alignment unit is configured to align the single-cell RNA sequencing data to a reference genome to obtain first data; The second comparison unit is configured to extract uniquely compared data from the first data to obtain second data.

[0015] In some embodiments, the system further comprises: an analysis module configured to identify the transcription start site of the gene to be analyzed.

[0016] In some embodiments, the analysis module specifically includes: a first analyzing unit configured to retrieve the second data and obtain third data based on the 3' complementary sequence located in the annotated region of the gene desired to be analyzed in the second data; A second analysis unit is configured to align the starting position of the 5'-end complementary sequence in the third data to the reference genome to serve as the first transcription start site of the transcript corresponding to the sequencing fragment; a third analysis unit, configured to obtain an empirical cumulative distribution function of the first transcription start site of the gene desired to be analyzed; a fourth analyzing unit configured to sequentially move the first window of the first length W forward from the farthest end of the distribution area L of the first transcription start site based on the third data, each time by moving the window forward by the second length, to obtain L-W+1 first windows; In some embodiments, the first length W comprises 10-100 nt.

[0017] In some embodiments, the first length W is 10 nt.

[0018] In some embodiments, the second length comprises 1-5 nt.

[0019] In some embodiments, the second length is 1-5 nt.

[0020] a fifth analyzing unit, configured to calculate a growth rate of the first start site in each first window based on an empirical cumulative distribution function of the first transcription start site, retaining first windows having a growth rate greater than a first preset value or a growth rate greater than a second preset value and containing an annotated transcription start site, and merging overlapping first windows to obtain a second window; In some embodiments, the first preset value includes 5-10%.

[0021] In some embodiments, the first preset value is 10%.

[0022] In some embodiments, the second preset value includes 1-5%.

[0023] In some embodiments, the second preset value is 5%.

[0024] The sixth analysis unit is configured to obtain the second transcription start site of the gene to be analyzed according to whether the second window contains an annotated transcription start site.

[0025] In some embodiments, when the second window contains an annotated transcription start site, the annotated transcription start site with the highest growth rate is used as the second transcription start site of the second window; when the second window does not contain an annotated transcription start site, the first start site with the highest growth rate is used as the second transcription start site of the second window.

[0026] In some embodiments, the system further comprises: a corrected PSI calculation module configured to calculate the degradation rate and corrected PSI of the second transcription start site of the gene to be analyzed based on a preset first model and an expectation-maximization algorithm.

[0027] In some embodiments, the first model comprises: Among them, g represents the gene to be analyzed, t represents the effective region of the transcription start site, and c represents a constant. represents the degradation rate, represents the corrected PSI, j represents the effective region of the transcription start site using the third window t Cut into equal length blocks (each block has length 1), Representative Block j length, Represents the slave block j The center of the transcription start site is located in the effective region t The normalized distance to the 3' end of For blocks j The number of first transcription start sites within the region.

[0028] In some embodiments, the transcription start site effective region refers to the position distribution region of the first transcription start site in the reference genome.

[0029] In some embodiments, the third length of the third window is 1 nt.

[0030] In some embodiments, the system may further include: The database is configured to store one or more of the single-cell RNA sequencing data of the target cell, the first data, the second data, the third data, the first transcription start site information, the second transcription start site information, the degradation rate information, and the corrected PSI information.

[0031] In some embodiments, the system may further include: The visualization module is configured to perform graphical processing on the data related to the second transcription start site and generate a graphical display result.

[0032] In some embodiments, the graphical display result may display the relative position and / or density distribution of the second transcription start site in the gene to be analyzed.

[0033] In some embodiments, the graphical display result can also display the expression distribution of the second transcription start site in different cells.

[0034] In some embodiments, the system may further include: The differential analysis module is configured to further analyze, based on the second transcription start site information, at least one of the number of genes to which the second transcription start site belongs, the number of cells associated with the second transcription start site, the expression level of the second transcription start site at the single-cell level, and the expression level of the second transcription start site in different cells.

[0035] Compared with the prior art, the beneficial effects of the present invention include at least the following aspects: With the application of 5'-enriched single-cell transcriptome sequencing technology, sequencing costs have been significantly reduced, but existing methods for identifying and analyzing selective transcription start sites have difficulty in effectively utilizing this technology. The present invention found that due to RNA degradation, short reads tend to be significantly more generated at the 3' end and decrease exponentially towards the 5' end. RNA degradation is mainly caused by the natural degradation of RNA molecules and the incomplete extension of RNA molecules during PCR amplification during library preparation. In selective splicing, especially in the quantitative analysis of selective transcription start sites, existing quantitative analysis methods for selective transcription start sites based on single-cell RNA sequencing (such as SCAFE and CamoTSS) actually ignore the uneven distribution of short reads, so the accuracy of transcript abundance estimation is insufficient, and the calculated PSI value is also biased. Therefore, the existing methods have the following problems: (1) They mainly focus on the analysis of selective transcription start site regions or clusters and cannot accurately locate the position of transcription start sites; (2) They do not consider the bias caused by RNA degradation in the quantitative analysis of transcription start sites.

[0036] Different transcription start sites may lead to changes in the splicing pattern of mRNA, thereby affecting the structure and function of the protein product. The PSI value (Percent Spliced ​​In) can reveal the impact of different transcription start sites on the splicing structure. As mentioned above, it is difficult for existing methods to accurately analyze selective transcription start sites, and it is difficult to ensure the accuracy of the quantitative results of selective transcription start sites. Therefore, the PSI calculated based on the existing methods also has the problem of insufficient accuracy. Based on the discovery that there is an uneven coverage of transcripts in single-cell RNA sequencing data, the present invention provides a transcription start site quantitative analysis system based on single-cell sequencing, which can not only realize the accurate identification of transcription start sites, but also calculate the corrected PSI at the transcription start site level. Experiments of the present invention have shown that the corrected PSI calculated by the system of the present invention can rescue the expression level of the transcription start site. Compared with the PSI calculated by traditional methods, it is more in line with the actual situation of the target cells, and the corrected PSI is also difficult to calculate by other methods.

[0037] In addition, it should be emphasized that the system provided by the present invention can be directly connected to the conventional single-cell analysis process to achieve single-cell level quantification, differential analysis and visualization of selective transcription start sites. Verification in real and simulated samples shows that the present invention is superior to existing methods in terms of accuracy and recall rate. The practical application of the present invention on physiological and pathological samples has successfully identified a large number of high-quality selective transcription start sites. These sites not only help to perform cell clustering analysis at the single-cell level, but also reveal how the dynamic expression profile of selective transcription start site events affects the physiological and pathological transcriptome.

[0038] In summary, the present invention can more realistically reflect biological differences, better reveal the role of transcription start sites in regulating splicing variations, more accurately identify alternative splicing events, and facilitate the analysis of specific regulatory patterns at the single cell, cell type, and even tissue levels, providing new research directions and opportunities for further exploring cellular dynamics and diversity in health and cancer. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following is a brief introduction to the drawings required for the embodiments or the description of the prior art. In all drawings, similar elements or parts are generally identified by similar reference numerals. In the drawings, the various elements or parts are not necessarily drawn according to the actual scale. Obviously, the drawings described below are some embodiments of the present invention. For those of ordinary skill in the art, other drawings can also be obtained based on these drawings without inventive work.

[0040] Figure 1 is a flow chart of the method of embodiment 1; Figure 2 This is a module diagram of the system of Example 2; Figure 3 A schematic diagram showing the visualization of alternative transcription start sites; Figure 4 A graph showing the uniformity of transcript coverage across different single-cell sequencing datasets. Figure 5 A result diagram showing the performance of the present invention on a simulated data set; Figure 6 The results are shown in the figure showing the performance of the present invention on a real data set; Figure 7 A diagram showing the results of the present invention applied to COVID-19 samples; Figure 8 A diagram showing the results of the differential analysis of COVID-19 samples performed by the present invention; Figure 9 A diagram showing the results of the present invention's analysis of COVID-19 sample cell clustering; Figure 10 This is a diagram showing the results of differential analysis of tumor samples performed by the present invention. DETAILED DESCRIPTION

[0041] To make the purpose, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0042] Herein, suffixes such as "module," "component," or "unit" used to represent elements are only used to facilitate description of the present invention and have no specific meaning. Therefore, "module," "component," or "unit" may be used interchangeably.

[0043] As used herein, terms such as "upper," "lower," "inner," "outer," "front," "back," "one end," and "the other end" indicate positions or locations based on those shown in the accompanying drawings. These terms are intended solely to facilitate and simplify the description of the present invention and are not intended to indicate or imply that the devices or components referred to must have, be constructed, or operate in a specific orientation. Therefore, they should not be construed as limitations on the present invention. Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance.

[0044] As used herein, unless otherwise expressly specified or limited, the terms "installed," "provided with," and "connected" should be understood broadly. For example, "connected" may refer to a fixed connection, a detachable connection, or an integral connection; it may refer to a mechanical connection, a direct connection, an indirect connection via an intermediate medium, or internal communication between two components. Those skilled in the art will understand the specific meanings of the above terms in the present invention on a case-by-case basis.

[0045] As used herein, "and / or" includes any and all combinations of one or more of the associated listed items.

[0046] Herein, "plurality" means two or more than two, ie, it includes two, three, four, five, etc.

[0047] It should be noted that, in this document, the terms "comprises," "includes," or any other variations thereof are intended to encompass non-exclusive inclusion, such that a process, method, article, or apparatus comprising a series of elements includes not only those elements but also other elements not explicitly listed, or elements inherent to such process, method, article, or apparatus. In the absence of further limitations, an element defined by the phrase "comprising a ..." does not exclude the presence of other identical elements in the process, method, article, or apparatus comprising the element.

[0048] As used in this specification, the term "about" typically means + / - 5% of the stated value, more typically + / - 4% of the stated value, more typically + / - 3% of the stated value, more typically + / - 2% of the stated value, even more typically + / - 1% of the stated value, and even more typically + / - 0.5% of the stated value.

[0049] In this specification, certain embodiments may be disclosed in a format that is within a range. It should be understood that this description of "within a range" is merely for convenience and brevity and should not be interpreted as a rigid limitation on the disclosed range. Therefore, the description of a range should be considered to have specifically disclosed all possible subranges and individual numerical values ​​within this range. For example, the description of a range of 1-6 should be considered to have specifically disclosed subranges such as from 1 to 3, from 1 to 4, from 1 to 5, from 2 to 4, from 2 to 6, from 3 to 6, etc., as well as individual numbers within this range, such as 1, 2, 3, 4, 5, and 6. Regardless of the breadth of the range, the above rules apply.

[0050] Example 1 like Figure 1 As shown, this embodiment provides a method for quantitative analysis of transcription start sites based on single-cell sequencing, comprising the following steps: S101 obtains single-cell RNA sequencing data of target cells; Single-cell RNA sequencing (scRNA-seq) is a technique used to analyze the transcriptome expression profile of individual cells. This method obtains target cell-specific gene expression information by reverse transcribing, amplifying, and constructing libraries of target cell mRNA, then sequencing the libraries using a high-throughput sequencing platform.

[0051] It should be understood that the present invention is not limited to the specific type of target cells, as long as single-cell RNA sequencing can be performed.

[0052] After constructing a single-cell RNA sequencing library, paired-end sequencing can be used to read the library. Paired-end sequencing involves simultaneously sequencing both ends of the same sequencing fragment in a DNA library, generating two complementary but non-overlapping reads: read1 (also referred to as the 5' complementary sequence, forward read sequence, or R1 in this context) and read2 (also referred to as the 3' complementary sequence, reverse read sequence, or R2 in this context).

[0053] In some embodiments, the single-cell RNA sequencing is 5'-end single-cell sequencing.

[0054] 5' end single-cell sequencing can obtain the transcription start site (TSS) of each transcript in the target cell at the single-cell level by reverse transcribing, amplifying and sequencing the 5' end of mRNA.

[0055] The experimental results of the present invention show that in single-cell RNA sequencing data (especially 5'-end single-cell sequencing), there is a common position-dependent sequencing fragment coverage pattern, that is, the distribution of sequencing fragments along the entire length of the transcript is not uniform, but may be over-enriched at the 3' and 5' ends.

[0056] In some embodiments, the single-cell RNA sequencing data includes paired-end sequencing data of multiple sequencing fragments.

[0057] In some embodiments, the double-end sequencing data includes 5'-end complementary sequence data and 3'-end complementary sequence data corresponding to each sequencing fragment.

[0058] In some embodiments, the 5'-complementary sequence is the base sequence data obtained by sequencing starting from the 5' end of the sequencing fragment. In some embodiments, the 3'-complementary sequence is the base sequence data obtained by sequencing starting from the 3' end of the sequencing fragment. In some embodiments, the 5'-complementary sequence and the 3'-complementary sequence are generally complementary and may not overlap.

[0059] In some embodiments, the sequencing fragments are high-quality sequencing fragments, that is, sequencing data units that meet certain quality control standards and can be accurately used for subsequent alignment, annotation, and analysis, such as high base quality values ​​(Phred quality value (Q score) of each site is usually ≥20 or ≥30), complete length, no adapter contamination, low-complexity sequences, rRNA, bacterial or host contamination sequences, high alignment rate, no redundancy / no duplication, etc.

[0060] S102: Determine the genes to be analyzed; In some embodiments, the expression level of the gene to be analyzed should be greater than a set first expression level threshold. In some embodiments, the first expression level threshold may be ≥200.

[0061] S103 aligning the single-cell RNA sequencing data to a reference genome to obtain first data; A reference genome refers to a standard set of nucleic acid sequences at the genomic level for a species (e.g., human, mouse, zebrafish, Arabidopsis, etc.). It typically consists of high-quality, continuous DNA sequences and represents a typical reference version of the genomic structure of that species. In some embodiments, the reference genome can be derived from a public database. In some embodiments, the obtained sequencing data (e.g., paired-end sequencing data) can be aligned with the reference genome using a suitable alignment algorithm to obtain one or more of its location information, annotation information, and variation information within the genome.

[0062] S104 extracting unique matched data from the first data to obtain second data; Unique alignment means that the sequencing fragment obtained by sequencing (for example, the paired sequence composed of the 5' complementary sequence and the 3' complementary sequence) can only match (mapping) one position on the reference genome.

[0063] S105 obtaining third data based on the 3'-end complementary sequence located in the annotated region of the gene to be analyzed in the second data; An annotated region refers to the entire region of the gene you wish to analyze in a reference genome, containing sequence regions known and annotated to have specific biological functions. In some embodiments, the annotated region typically includes the coding sequence (CDS), exons, introns, promoters, terminators, non-coding RNA (e.g., miRNA, lncRNA) regions, and UTRs (5' / 3' untranslated regions).

[0064] S106 aligning the starting position of the 5'-end complementary sequence in the third data to the reference genome as the first transcription start site of the transcript corresponding to the sequenced fragment; The transcription start site of a transcript is the genomic location corresponding to the first nucleotide on the DNA template strand that is recognized by RNA polymerase and initiates transcription. This site defines the 5' end of the transcript and is the starting point of transcription, typically located downstream of the promoter region. In eukaryotes, the same gene can have multiple transcription start sites, resulting in multiple transcripts with different 5' untranslated regions (5'UTRs) or regulatory structures.

[0065] As mentioned above, single-cell RNA sequencing data is generally subject to position-dependent sequencing fragment coverage patterns. The present invention's analysis reveals that this problem is primarily caused by RNA degradation. Specifically, due to RNA degradation, short reads tend to be generated more frequently at the 3' end, with the number decreasing exponentially towards the 5' end. In the context of alternative splicing, especially the quantitative analysis of alternative transcription start sites, existing single-cell RNA sequencing-based quantitative analysis methods for alternative transcription start sites (e.g., SCAFE and CamoTSS) effectively ignore this uneven distribution of short reads. Consequently, transcript abundance estimates are inaccurate, and the calculated PSI value (i.e., the proportion of sequence fragments that include the target exon and those that skip the target exon) is also biased.

[0066] To "rescue" the expression level of the transcription start site caused by RNA degradation and obtain a more accurate PSI value (which can also be understood as "correcting PSI" in the present invention), the present invention first captures the adjacent potential transcription start site (i.e., the first transcription start site) based on the cumulative distribution of the start site of the 5' complementary sequence in the third data: S107: Obtain an empirical cumulative distribution function of the first transcription start site of the gene to be analyzed, and sequentially move a first window of a first length W forward from the farthest end of the distribution region L of the first transcription start site, each time moving forward by a second length, to obtain L-W+1 first windows; In some embodiments, the first length W comprises 10-100 nt.

[0067] In some embodiments, the first length W may be 10 nt.

[0068] In some embodiments, the second length comprises 1-5 nt.

[0069] In some embodiments, the second length may be 1-5 nt.

[0070] In some embodiments, the distribution region L of the first transcription start site refers to the position distribution range of the gene to which the first transcription start site belongs in the reference genome.

[0071] In some embodiments, the most distal end refers to the first transcription start site that is closest to the 3' end of the corresponding gene among the first transcription start sites.

[0072] In some embodiments, a window refers to a continuous region on a genome or transcript, ie, an interval of fixed length.

[0073] S108: calculating the growth rate of the first start site in each first window based on the empirical cumulative distribution function of the first transcription start site, retaining first windows having a growth rate greater than a first preset value or a growth rate greater than a second preset value and containing an annotated transcription start site, and merging overlapping first windows to obtain a second window; The present invention can effectively screen the first transcription start site through step S108 and setting specific growth rates for different situations, thereby retaining a second window with a greater possibility and thus retaining a more authentic second transcription start site.

[0074] In some embodiments, the first preset value includes 5-10%.

[0075] In some embodiments, the first preset value is 10%.

[0076] In some embodiments, the second preset value includes 1-5%.

[0077] In some embodiments, the second preset value is 5%.

[0078] S109 obtaining a second transcription start site of the gene to be analyzed according to whether the second window contains an annotated transcription start site; In some embodiments, step S109 specifically includes: when the second window contains an annotated transcription start site, the annotated transcription start site with the highest growth rate is used as the second transcription start site of the second window; when the second window does not contain an annotated transcription start site, the first start site with the highest growth rate is used as the second transcription start site of the second window; In some embodiments, the step S109 further includes: obtaining the expression level of the second transcription start site.

[0079] S110 Obtaining the degradation rate and corrected PSI of the second transcription start site of the gene to be analyzed according to the preset first model and the expectation maximization algorithm: In some embodiments, the first model comprises: Among them, g represents the gene to be analyzed, t represents the effective region of the transcription start site, and c represents a constant. represents the degradation rate, represents the corrected PSI, j represents the effective region of the transcription start site using the third window t Cut into equal length blocks (each block has length 1), Representative Block j length, Represents the slave blockj The center of the transcription start site is located in the effective region t The normalized distance to the 3' end of For blocks j The number of first transcription start sites within the region.

[0080] In some embodiments, the transcription start site effective region refers to the position distribution region of the first transcription start site in the reference genome.

[0081] It should be noted that in order to further improve the accuracy of the quantitative determination of transcription start sites, the present invention assigns a degradation rate ( ) and corrected PSI ( ), see the first model: (First model) Specifically, for a gene that is expected to be analyzed g , a fixed-length third window (1 bp) was used to separate the transcription start site into the effective region t Cut into pieces of equal length j (j=1), Representative Block j length. For the slave block j The center of the transcription start site is located in the effective region t (This can be understood as the normalized distance from the 3' end of the transcript to the transcription start site cluster, which represents the effective area on the transcript used to find the transcription start site) ), For blocks j The number of the first transcription start site in the region. g All the first transcription start sites of the align to the reference genome can be obtained belonging to each block j of .

[0082] The transcription start site effective region t The effective length This can be obtained through the second model: ≡ (Second model) Among them, g represents the gene to be analyzed, t represents the effective region of the transcription start site, represents the degradation rate, j represents the effective region of the transcription start site using the third window t Cut into equal length blocks (each block has length 1), Representative Blockj length, Represents the effective length of the effective region t of the transcription start site, is the index matrix of the effective region t and block j of the transcription start site, Represents the slave block j The center of the transcription start site is located in the effective region t The normalized distance to the 3' end of the

[0083] The index is 1 when the transcription start site valid region t contains this block j, otherwise it is 0.

[0084] The present invention can use the EM (maximum expectation) algorithm to obtain the degradation rate of all TSS (second transcription start site) of the gene g to be analyzed ( ) and corrected PSI ( ).

[0085] As an example, the present invention constructs a method for estimating the degradation rate of the second transcription start site of the gene g that is expected to be analyzed ( ) and corrected PSI ( ), that is, the degradation rate ( ) and corrected PSI ( ) is the optimal estimate of .

[0086] The likelihood function of the gene g to be analyzed can be obtained using the third model: (Third model) Among them, g represents the gene to be analyzed. represents the joint likelihood function constructed for the gene g expected to be analyzed; t represents the effective region of the transcription start site, and j represents the effective region of the transcription start site using the third window. t Cut into pieces of equal length, Represents the effective length of the effective region t of the transcription start site; represents the degradation rate, Represents the slave block j The center of the transcription start site is located in the effective region t The normalized distance to the 3' end of the gene; h represents the index variable, indicating all the effective regions of the transcription start site of the gene g that is expected to be analyzed, while t is only one of them; Represents each transcription start site region h of the gene g to be analyzed Calculate them all and add them all up. In some embodiments, the constraints Next, As The function is maximized to obtain the optimal .

[0087] The EM process can be: Step E: Given and An initial value and calculate The expected value of Among them, h represents the index variable, which indicates all the effective regions of the transcription start sites of the gene g that is expected to be analyzed, and t is just one of them; x represents the iteration round of the EM algorithm, is the intermediate expected value calculated in the xth iteration.

[0088] M step: maximize , and order Here, x and x+1 represent the iterations of the EM algorithm. The whole process repeats the E and M steps, substituting the results of the x round into the calculation to obtain the parameters of the x+1 round. The parameters of the x+1 round are then used to calculate the x+2 round, and this cycle continues until the parameter values ​​converge (no longer change). ^ indicates that it is an "estimator" or "estimated value." In other words, β is the true parameter that we want to know, and is the calculated estimate of β; Represents the updated estimated value of the contribution ratio parameter β of the effective region t of the transcription start site of the gene g expected to be analyzed in the x+1th iteration.

[0089] in and is the number of all second transcription start sites on the gene to be analyzed. The E step and M step are cycled multiple times (default 100 times) until the model converges, that is, the degradation rate and corrected PSI The value of no longer changes, and the gene that is expected to be analyzed is fitted. g Optimal degradation rates at different transcription start sites and corrected PSI .

[0090] Because of the degradation, the PSI calculated based on the traditional method will be relatively low. However, the corrected PSI obtained by the first model and the expected maximum algorithm in the present invention "rescues" the expression amount that should belong to this transcription start site, so the degradation rate ( ) Larger transcription start sites tend to have higher corrected PSI. Chchd2 For example, the degradation degree of its transcription start site 1 is low ( = 0.14), while transcription start site 2 showed a higher degradation level ( = 2.76). In terms of PSI calculation, the traditional PSI of the transcription start site 2, which is more severely degraded, is only 0.1, while the corrected PSI calculated based on the present invention increases to 0.27.

[0091] When the transcription start site is identified, for avoiding the randomly generated fragment error identification as the transcription start site, the existing method for identifying the transcription start site based on cumulative distribution assumes that the distribution of the sequencing fragments based on random distribution is different from the distribution of the transcription start site, so it needs to fit a random distribution, and then the distribution that may be the transcription start site is compared with the random distribution. The present invention finds that the result obtained by this method is affected by the reliability of random distribution to a great extent. If the sequencing depth is not enough or the transcription start site is larger at intervals, it is possible to be inseparable from the random distribution, and then it is difficult to accurately identify the transcription start site. The present invention does not adopt the random distribution fitting of above-mentioned complexity, and from single cell RNA sequencing data, according to the cumulative distribution of the first transcription start site and by the threshold value of the retention rate provided, the transcription start site can be accurately identified on the contrary, and then the degradation rate and more real correction PSI of the transcription start site level can be calculated.

[0092] Example 2 like Figure 2 As shown, this embodiment provides a transcription start site quantitative analysis system 100 based on single-cell sequencing corresponding to the method of embodiment 1, including: The data acquisition module 102 is configured to acquire single-cell RNA sequencing data of target cells.

[0093] In some embodiments, the single-cell RNA sequencing is 5'-end single-cell sequencing.

[0094] In some embodiments, the single-cell RNA sequencing data includes paired-end sequencing data of multiple sequencing fragments.

[0095] In some embodiments, the double-end sequencing data includes 5'-end complementary sequence data and 3'-end complementary sequence data corresponding to each sequencing fragment.

[0096] In some embodiments, the 5' complementary sequence is the base sequence data obtained by sequencing from the 5' end of the sequencing fragment. In some embodiments, the 3' complementary sequence is the base sequence data obtained by sequencing from the 3' end of the sequencing fragment.

[0097] In some embodiments, the data acquisition module 102 may also be configured to extract high-quality sequencing fragment data from the single-cell RNA sequencing data.

[0098] In some embodiments, the system 100 further includes a gene determination module 104 configured to obtain genes that are desired to be analyzed.

[0099] In some embodiments, the expression level of the gene to be analyzed should be greater than a set first expression level threshold. In some embodiments, the first expression level threshold may be ≥200.

[0100] In some embodiments, the system 100 further includes: a sequence alignment module 106 configured to retrieve the single-cell RNA sequencing data for alignment according to the gene desired to be analyzed.

[0101] In some embodiments, the sequence alignment module 106 may include: The first comparison unit 202 is configured to compare the single-cell RNA sequencing data to a reference genome to obtain first data.

[0102] The second comparing unit 204 is configured to extract uniquely matched data from the first data to obtain second data.

[0103] In some embodiments, the system further includes: an analysis module 108 configured to identify the transcription start site of the gene to be analyzed.

[0104] In some embodiments, the analysis module 108 specifically includes: The first analysis unit 302 is configured to retrieve the second data and obtain third data based on the 3' complementary sequence located in the annotated region of the gene to be analyzed in the second data; The second analysis unit 304 is configured to align the starting position of the 5'-end complementary sequence in the third data to the reference genome to serve as the first transcription start site of the transcript corresponding to the sequencing fragment; The third analysis unit 306 is configured to obtain an empirical cumulative distribution function of the first transcription start site of the gene desired to be analyzed; The fourth analyzing unit 308 is configured to sequentially move the first window of the first length W forward from the farthest end of the distribution area L of the first transcription start site based on the third data, each time by moving the window forward by the second length, to obtain L-W+1 first windows; In some embodiments, the first length W comprises 10-100 nt.

[0105] In some embodiments, the first length W may be 10 nt.

[0106] In some embodiments, the second length comprises 1-5 nt.

[0107] In some embodiments, the second length may be 1-5 nt.

[0108] a fifth analyzing unit 310 configured to calculate a growth rate of the first start site in each first window based on an empirical cumulative distribution function of the first transcription start site, retaining first windows having a growth rate greater than a first preset value or a growth rate greater than a second preset value and containing an annotated transcription start site, and merging overlapping first windows to obtain a second window; In some embodiments, the first preset value includes 5-10%.

[0109] In some embodiments, the first preset value is 10%.

[0110] In some embodiments, the second preset value includes 1-5%.

[0111] In some embodiments, the second preset value is 5%.

[0112] The sixth analyzing unit 312 is configured to obtain the second transcription start site of the gene to be analyzed according to whether the second window contains an annotated transcription start site.

[0113] In some embodiments, when the second window contains an annotated transcription start site, the annotated transcription start site with the highest growth rate is used as the second transcription start site of the second window; when the second window does not contain an annotated transcription start site, the first start site with the highest growth rate is used as the second transcription start site of the second window.

[0114] In some embodiments, the system 100 further includes: a corrected PSI calculation module 110 configured to calculate the degradation rate and corrected PSI of the second transcription start site of the gene to be analyzed based on a preset first model and an expectation-maximization algorithm.

[0115] In some embodiments, the first model comprises: Among them, g represents the gene to be analyzed, t represents the effective region of the transcription start site, and c represents a constant. represents the degradation rate, represents the corrected PSI, j represents the effective region of the transcription start site using the third window t Cut into equal length blocks (each block has length 1), Representative Block j length, Represents the slave block j The center of the transcription start site is located in the effective region t The normalized distance to the 3' end of For blocks j The number of first transcription start sites within the region.

[0116] In some embodiments, the system 100 may further include: The database 112 is configured to store one or more of the single-cell RNA sequencing data, the first data, the second data, the third data, the first transcription start site information, the second transcription start site information, the degradation rate information, and the corrected PSI information of the target cell.

[0117] In some embodiments, the system 100 may further include: The visualization module 114 is configured to perform graphical processing on the data related to the second transcription start site and generate a graphical display result.

[0118] In some embodiments, the graphical display result may display the relative position and / or density distribution of the second transcription start site in the gene to be analyzed (eg, a distribution map of the second transcription start site on the sequencing fragments).

[0119] In some embodiments, the graphical display result can also display the expression distribution of the second transcription start site in different cells.

[0120] In some embodiments, the system 100 may further include: The difference analysis module 116 is configured to further analyze at least one of the number of genes to which the second transcription start site belongs, the number of cells associated with the second transcription start site, the expression level of the second transcription start site at the single-cell level, and the expression level of the second transcription start site in different cells based on the second transcription start site information.

[0121] Example 3 For transcription start site quantification, this example only considers uniquely aligned, high-quality sequencing reads. To successfully quantify novel alternative transcription start sites (TSSs) within annotated regions, the reverse read sequence (R2) of each pair of fragments must be completely within the annotated region of a given gene. The start position of the forward read sequence (R1) of all fragments is then mapped to the genomic position (or the end position for negative-strand genes) to determine the start site of the transcript corresponding to that fragment. The empirical cumulative distribution function (eCDF) of all forward read sequence start sites for that gene is then calculated.

[0122] The growth rate of the R1 start site within each fixed-length window (100 nt by default) is then calculated, starting from the farthest end, and advancing 1 nt at a time, for a total of L-w+1 windows (L is the genomic distribution of the forward read start site). Windows with a growth rate greater than 10% or greater than 5% but containing a known transcription start site are retained, and all overlapping windows are merged. For the merged window W, the expression level of all forward read start sites within the window is calculated. If the window contains an annotated transcription start site, the known transcription start site with the highest growth rate is used as the final transcription start site for the window. Otherwise, the growth rate is calculated for each site, and the site with the highest growth rate is used as the final transcription start site for the window.

[0123] This study demonstrates that RNA degradation at the 5' end of transcripts and the uneven distribution of short reads in 5'-enriched single-cell transcriptome sequencing significantly reduce the accuracy of quantification based on the uniformity assumption. The methods and systems provided by this study can accurately infer the relative usage of alternative transcription start sites from 5'-enriched single-cell transcriptome sequencing.

[0124] The specific steps are as follows: The input is as follows: Prepare the R result object, sequencing result BAM file, and reference GTF file obtained from the regular Seurat analysis process as quantitative input files and input them into the object, bam, and gtf parameters respectively.

[0125] Select the genes to be quantified and enter them into the genes parameter. If all expressed genes are quantified by default, this parameter does not need to be set.

[0126] Based on sequencing depth, set a gene expression threshold. Genes below this threshold will not be included in the analysis of selective transcription start sites. Select a value greater than 0 as the threshold and enter it into the MinGeneReads parameters. The default value is 500.

[0127] Select whether the quantification region of alternative transcription start sites is limited to the untranslated region (UTR) by inputting it into the UTROnly parameter, which defaults to 0.

[0128] Select the sequencing fragment quality for quantification and input it into the mapqFilter, isSecondaryAlignment, isSupplementaryAlignment, and isDuplicate parameters. The default values ​​are 255, FALSE, FALSE, and NA.

[0129] Set the minimum percentage of transcription start sites as input to the min.TSS.percent parameter, which defaults to 5.

[0130] Set the p-value threshold for transcription start sites, input to the p.cutoff parameter, the default is 0.01.

[0131] Set the parameter used to evaluate the degree of local variation in transcription start sites. Input to min.local.percent, the default value is 1.

[0132] Set the window size used to estimate the local mean. Enter the window parameter. The default value is 10.

[0133] Sets whether to use a greedy algorithm to analyze the first transcription start site. If this parameter is set to TRUE, the first transcription start site will be printed even if it does not meet the specified conditions. Input to the greedy parameter, the default is TRUE.

[0134] Set the number of CPUs to run, not exceeding the maximum number of cores on the computer, and enter it into the cores parameter. The default value is 1.

[0135] Set whether to visualize quantitative progress, input to the verbose parameter, the default is FALSE.

[0136] Result output: The output of selective transcription start site quantification contains detailed information about each transcription start site, stored in `rowRange` objects and organized into a `Granges` object, as shown in Table 1.

[0137] Table 1 Schematic diagram of the quantitative results of selective transcription start sites Note: Columns 1 through 12 in Table 1 are related to quantitative analysis. Columns 1 through 5 contain the genomic location of the transcription start site and the gene to which it belongs. The TSS column contains the transcription start site ID, and the PSI (traditional PSI) and Percent columns contain the quantitative expression level of that transcription start site. The Annotated column indicates whether the site is a known annotated transcription start site, and the Greedy column indicates whether the transcription start site was derived using a greedy algorithm. Columns 13 through 19 are related to RNA degradation, measuring the degradation level of the alternative transcription start site. Higher RNA_decay values ​​indicate greater degradation. Columns 30 through 21 provide other sequencing read-related information.

[0138] Based on the results obtained in the above steps, difference analysis can also be performed, as shown in Table 2: Table 2 Difference analysis results Note: TSS column: This column indicates the transcription start sites involved in the analysis. Each transcription start site represents a potential transcription start position. G1 / G2 column: This column indicates different groups (G1 and G2). These groups may represent different experimental conditions, treatment methods or different time points. n / N / Cells column: This column contains quantitative information related to specific transcription start sites under different groups. Specifically including: n: the number of the transcription start site itself in different groups; N: the number of genes to which the transcription start site belongs; Cells: the number of cells related to the transcription start site; PSI / Corrected PSI column: This column shows the expression level of the transcription start site at different levels: PSI: represents the traditional PSI of the transcription start site at the single-cell level, which can reflect the transcription activity of this position in a single cell. Corrected PSI: represents the corrected PSI. wald.test, wilcox.test, and prop.test columns: These three columns correspond to the p-values ​​of different hypothesis testing methods, including: wald.test: the p-value based on the Wald test, used to evaluate the significance of parameter estimation; wilcox.test: the p-value based on the Wilcoxon rank sum test, suitable for comparing non-normally distributed data; prop.test: the p-value based on the proportion test, mainly used to evaluate the proportion differences in categorical data.

[0139] Based on the results obtained in the above steps, selective transcription start site visualization can also be performed. For details, see Figure 3 The results are divided into two parts: The upper part: shows the distribution of 5'-enriched single-cell transcriptome sequencing fragments. The y-axis shows the detailed expression level, and the indicator line indicates the inferred transcription start site. The lower part: shows the location of the transcripts and transcription start sites of the reference transcriptome. The indicator line indicates the inferred transcription start site, and different transcripts are marked.

[0140] Figure 3 It shows that the positions of the transcription start sites in the upper and lower parts are consistent, indicating that the results of the quantitative analysis of the present invention are very accurate.

[0141] Example 4 First, the present invention was used to compare the transcript coverage uniformity between different single-cell sequencing datasets (including SMART-seq2, 3' end and 5' end single-cell sequencing data). Figure 4 As shown, the reads in SMART-seq2 sequencing are evenly distributed, while the 3' end single-cell sequencing method shows over-enrichment in the 3' end untranslated region, and the 5' end single-cell sequencing method is enriched in the 5' end untranslated region. In addition, amplification of the first exon region also found that both 5' end single-cell sequencing and SMART-seq2 showed over-enrichment in the 3' end.

[0142] By analyzing the uneven distribution of short reads, the present invention found that it was mainly caused by RNA degradation, resulting in a position-dependent sequencing fragment coverage pattern. In order to further improve the accuracy of transcription start site quantification, the present invention assigns a degradation rate ( ) and corrected PSI ( Because of degradation, the traditional PSI will be relatively low. The corrected PSI proposed in the present invention "rescues" the expression level that should belong to this transcription start site, so the degradation rate ( ) Larger transcription start sites tend to have higher corrected PSI. Chchd2 For example, the degradation degree of its transcription start site 1 is low ( = 0.14), while transcription start site 2 showed a higher degradation level ( = 2.76), a result supported by the scONT and Smart-seq2 datasets. In the calculation of PSI, the traditional PSI of the transcription start site 2 with more severe degradation was only 0.1, while the corrected PSI increased to 0.27 ( Figure 4 ).

[0143] In summary, the present invention has found that the transcript degradation rate has a significant impact on the accuracy of selective transcription start site quantification. Existing single-cell sequencing methods, such as SMART-seq2, 3'-end and 5'-end single-cell sequencing, have the problem of uneven transcript coverage, which may lead to deviations in the quantification of transcription start sites. The present invention deduces the RNA degradation results of the transcription start site and two different PSIs (traditional PSI and corrected PSI) by modeling the forward read sequence start site. This can more accurately and robustly quantify the abundance of transcription start sites at the single-cell level, and has important application prospects.

[0144] Example 5 Comparison of the robustness of the proposed method in simulated and real data sets In order to verify the robustness of the present invention, the present invention uses a simulated single-cell sequencing data set ( Figure 5 a) Comparative analysis with two existing quantitative methods, SCAFE and CamoTSS. The simulated dataset used covers 920 genes and 1,078 cells, and the predicted transcription start sites are evaluated using high-confidence transcription start sites identified in the reference genome as the "true value."

[0145] The experimental results show that the estimated transcription start sites are highly consistent with the true values, and the difference between the estimated transcription start sites and the annotated transcription start sites does not exceed 85 base pairs ( Figure 5 b and Figure 5 c). Further comparative analysis showed that the present invention significantly outperformed SCAFE and CamoTSS in terms of recall rate and F score for transcription start site identification ( Figure 5 d).

[0146] To further evaluate the performance of the present invention in real datasets, the present invention conducted cross-validation based on single-cell sequencing data of mouse bone marrow samples. This dataset covers 10x single-cell sequencing, third-generation sequencing ONT, and PacBio data. In the 10x single-cell sequencing data, the present invention successfully detected 11,396 transcription start sites involving 8,332 genes, of which approximately 31% of the genes had multiple transcription start sites ( Figure 6 a). In addition, 92% of the transcription start sites were annotated in both reference files ( Figure 6 b).

[0147] Overall, the present invention demonstrated higher precision, recall, and F-score ( ) than SCAFE and CamoTSS in identifying transcription start sites within 50 bp. Figure 6 c). To demonstrate the ability of the present invention to identify annotated, novel, and multiple transcription start sites, we cross-validated our predicted transcription start sites using data from different sequencing platforms. For example, all three annotated transcription start sites for the Emg1 gene were identified by the present invention in the 10x data set and further validated using ONT and PacBio data, whereas SCAFE and CamoTSS failed to accurately identify these sites. Importantly, the present invention identified both annotated and novel sites for the Idh3g gene, whereas SCAFE and CamoTSS failed to identify any sites in this gene ( Figure 6 d).

[0148] Example 6 The present invention constructs the expression profile of the new crown transcription start site In the study of COVID-19 infection, although population RNA sequencing has revealed that selective promoter usage is associated with disease specificity and patient prognosis, the potential of precise transcription start site usage at single-cell resolution to explain the heterogeneity of COVID-19 infection remains underexplored. To address this issue, the present invention analyzed a human COVID-19 dataset, which included samples from COVID-19 patients (Cov; n=6) and healthy controls (NHC; n=3), totaling 90,127 cells. Using UMAP, the present invention clustered the immune cell populations between Cov and NHC samples based on gene expression profiles and identified 18 immune cell clusters ( Figure 7 a. Figure 7 b), the results showed that the proportion of lymphocytes was higher in the Cov group ( Figure 7 c), indicating a significant immune response.

[0149] Next, we analyzed the entire dataset and identified 12,487 transcription start sites, of which 3,216 (from 956 genes) showed significant differential usage across all 18 cell types ( Figure 7 d. Figure 7 e), identified by the “FindMarkers” function (p value < 0.05, proportion test).

[0150] It is worth noting that among the host genes with differential transcription start sites in COVID-19 patients, several known immune-related regulatory genes were found, such as CDC25B, which plays a key role in a variety of human tumors. The distal transcription start site of CDC25B (chr20:3795760:+) is highly expressed in NHC, while the proximal TSS (chr20:3796310:+) is highly expressed in Cov ( Figure 8 Interestingly, when analyzing from the perspective of cell type, the present invention found that the most significant differences occurred in the initial CD8 cells ( Figure 8 b). Further investigation of the biological functional relevance of Cov-specific transcription start sites in different cell types revealed that host genes with differential transcription start sites in naive CD8 T cells were enriched in immune-related GO pathways, such as "production of immune response mediators" and "production of interleukin-2." This finding highlights the diversity of transcription start sites required for these immune-related genes, an important result observed in naive CD8 T cells ( Figure 8 c).

[0151] Example 7 Transcription start site expression profiling promotes cell clustering To elucidate the contribution of gene expression and transcription start site expression profiles to cell clustering, the present invention revealed partial overlap between gene expression and transcription start site expression by comparing differences in cell type-specific marker gene and transcription start site levels ( Figure 9 a), indicating that the transcription start site can be used as an additional predictive feature. Next, the present invention combines the gene expression profile and the transcription start site expression profile and analyzes them using the weighted nearest neighbor (WNN) UMAP method ( Figure 9 b). Independent analyses of gene and transcription start site expression broadly confirmed the consistency of the cell classifications, although some details differed. In particular, naive CD8 cells were indistinguishable from naive CD4 cells in terms of gene expression but were clearly separated in the transcription start site data ( Figure 9 c). In addition, when comparing the unsupervised weighted WNN in Seurat ( Figure 9 d), naive CD8 cells showed higher transcription start site modal weighting, while plasma cell B cells showed higher gene expression modal weighting. The present invention identified 959 genes with differential transcription start sites in naive CD8 cells (p value < 0.05, Wilcoxon test). For example, Figure 9 e shows the transcription start site of the LDLRAP1 gene, which plays a regulatory role in the metabolic process of adipose tissue. Although the gene expression of LDLRAP1 is not significantly different between different cell types, the transcription start site 2 (chr1:25562794:+) detected by the present invention shows significant expression in naive CD8 cells ( Figure 9 f). These results indicate that transcription start site expression profiling significantly enhances the cell clustering effect across different cell types.

[0152] Example 8 Application of tumor-specific transcription start sites in lung cancer To characterize tumor-specific transcription start site events, the present invention compared the expression of transcription start sites in tumor cells and normal cells in human lung cancer. Clustering results showed that there were immune and non-immune cells ( Figure 10 a). This example specifically focuses on epithelial cells (n=4,299) and further analyzes them by dividing them into malignant and non-malignant cell types ( Figure 10 b). For each cell type, the present invention uses the inferCNV method ( Figure 10 c) Copy number variations (CNVs) were inferred based on the average expression of 100 genes within each chromosomal region, and large-scale amplification regions were identified in malignant cells, including those on chromosomes 7, 11, and 20. At the same time, large-scale deletions were found, such as the EZR gene deletion on chromosome 6, and deletions on chromosomes 14 and 21, which are consistent with previous studies ( Figure 10 d). In addition, the present invention also identified 2,760 differential transcription start sites, the usage ratios of which were significantly changed between malignant cells and non-malignant cells ( Figure 10 e). Overall, the present invention provides new tools and methods for studying the expression of tumor-specific transcription start sites and their role in immune response and tumor progression, and has broad application potential.

[0153] Through the description of the above embodiments, those skilled in the art can clearly understand that the above-mentioned embodiment methods can be implemented by means of software plus the necessary general hardware platform, and of course can also be implemented by hardware, but in many cases the former is a better embodiment. Based on this understanding, the technical solution of the present invention is essentially or the part that contributes to the prior art can be embodied in the form of a software product, which is stored in a storage medium (such as ROM / RAM, magnetic disk, optical disk), including a number of instructions for enabling a computer terminal (which can be a mobile phone, computer, server, or network device, etc.) to execute the methods described in each embodiment of the present invention.

[0154] The embodiments of the present invention are described above in conjunction with the accompanying drawings, but the present invention is not limited to the above-mentioned specific implementation methods. The above-mentioned specific implementation methods are merely illustrative and not restrictive. Under the guidance of the present invention, ordinary technicians in this field can also make many forms without departing from the scope of protection of the present invention and the claims, all of which are protected by the present invention.

Claims

1. A transcription start site quantitative analysis system based on single-cell sequencing, characterized in that: include: A data acquisition module is configured to acquire single-cell RNA sequencing data of target cells; the single-cell RNA sequencing data includes paired-end sequencing data of multiple sequencing fragments, and the paired-end sequencing data includes 5' end complementary sequence data and 3' end complementary sequence data corresponding to each sequencing fragment; A gene determination module is configured to obtain genes that are determined to be analyzed; A sequence alignment module is configured to retrieve the single-cell RNA sequencing data for alignment based on the desired gene to be analyzed; the sequence alignment module includes: a first alignment unit configured to align the single-cell RNA sequencing data to a reference genome to obtain first data; a second alignment unit configured to extract unique aligned data from the first data to obtain second data; The analysis module is configured to identify the transcription start site of the gene to be analyzed; the analysis module includes: a first analysis unit, configured to retrieve the second data and obtain third data based on the 3' complementary sequence located in the annotation region of the gene to be analyzed in the second data; a second analysis unit, configured to align the starting position of the 5' complementary sequence in the third data to the reference genome as the first transcription start site of the transcript corresponding to the sequencing fragment; a third analysis unit, configured to obtain an empirical cumulative distribution function of the first transcription start site of the gene to be analyzed; a fourth analysis unit, configured to align the starting position of the 5' complementary sequence in the third data to the reference genome as the first transcription start site of the transcript corresponding to the sequencing fragment; A first window of length W is sequentially moved forward from the farthest end of the distribution area L of the first transcription start site, each time moving forward by a second length, to obtain L-W+1 first windows; a fifth analysis unit is configured to calculate the growth rate of the first start site in each first window based on the empirical cumulative distribution function of the first transcription start site, retain the first window with a growth rate greater than a first preset value, or a growth rate greater than a second preset value and containing an annotated transcription start site, and merge overlapping first windows to obtain a second window; a sixth analysis unit is configured to obtain the second transcription start site of the gene desired to be analyzed based on whether the second window contains an annotated transcription start site; The corrected PSI calculation module is configured to calculate the degradation rate and corrected PSI of the second transcription start site of the gene to be analyzed according to a preset first model and an expectation maximization algorithm, wherein the first model includes: Among them, g represents the gene to be analyzed, t represents the effective region of the transcription start site, and c represents a constant. represents the degradation rate, represents the corrected PSI, j represents the effective region of the transcription start site using the third window t Cut into pieces of equal length, Representative Block j length, Represents the slave block j The center of the transcription start site is located in the effective region t The normalized distance to the 3' end of For blocks j The number of first transcription start sites within the region.

2. The system according to claim 1, wherein The sixth analysis unit is configured to, when the second window contains an annotated transcription start site, use the annotated transcription start site with the highest growth rate as the second transcription start site of the second window; When the second window does not contain an annotated transcription start site, the first start site with the highest growth rate is used as the second transcription start site in the second window.

3. The system according to claim 1, wherein: The system may further include: a visualization module configured to perform graphical processing on the data related to the second transcription start site to generate a graphical display result.

4. The system according to claim 1, wherein: The system may further include: a differential analysis module configured to further analyze, based on the second transcription start site information, at least one of the number of genes to which the second transcription start site belongs, the number of cells associated with the second transcription start site, the expression level of the second transcription start site at the single-cell level, and the expression level of the second transcription start site in different cells.

5. The system according to claim 1, wherein: The expression level of the gene to be analyzed is greater than a set first expression level threshold, and the first expression level threshold is ≥200.

6. The system according to claim 1, wherein: The single-cell sequencing is 5'-end single-cell sequencing.

7. The system according to claim 1, wherein The first length comprises 10-100 nt.

8. The system according to claim 1, wherein: The second length comprises 1-5 nt.

9. The system according to claim 1, wherein: The first preset value includes 5-10%.

10. The system according to claim 1, wherein: The second preset value includes 1-5%.

Images