Single-cell transcriptome computational analysis method and system integrating deep learning model

Abstract

The present invention provides a single-cell transcriptome computational analysis method and system integrating a deep learning model. Specifically, the present invention provides a polyadenylation site (PAS) detection method and system based on a fusion deep learning model and a polyadenylation site detection and cell typing analysis method and system (SCAPTURE method and system) for single-cell sequencing data. A deep learning model is constructed and integrated in the SCAPTURE system of the present invention to achieve high-accuracy and non-position-dependent polyadenylation site prediction, and is used for the screening of highly reliable polyadenylation sites at the whole genome level; and single-cell transcriptome sequencing data can be used to detect the sequencing read distribution peak of the whole genome region, and the polyadenylation sites at the whole genome level can be identified from scratch; different transcripts on the same gene and their differential expression can also be identified based on the polyadenylation site, and then used for single-cell typing analysis.

Description

Technical Field

[0001] The present invention relates to the fields of biotechnology and bioinformatics, and in particular to a single-cell transcriptome computational analysis method and system based on a fusion deep learning model. Background Art

[0002] With the development of sequencing technology, a series of methods for large-scale determination of single-cell transcriptomes have emerged. Unlike the second-generation sequencing technology that sequences transcripts of cell populations, single-cell transcriptome sequencing technology can distinguish the cell and transcript sources of captured RNA molecules, and improve the accuracy of transcriptome sequencing analysis to the single-cell level. Single-cell transcriptome sequencing is mainly based on differences in gene expression and can reveal cellular heterogeneity within cell populations, tissues and organs, and even entire animals.

[0003] The currently widely used single-cell transcriptome sequencing technologies are mainly divided into 3'tag-based (such as 10x Chromium) and full-length library construction methods (such as Smart-seq2) according to their RNA amplification and library construction methods. By utilizing droplet-based microfluidics technology and cell barcode distribution systems, 3'tag-based single-cell transcriptome sequencing methods such as 10x Chromium can achieve parallel sequencing of thousands of single cells in one experiment. This method can encapsulate barcodes, primer gel beads and single cells in droplets, capture RNA molecules with oligo (dT) primers, obtain cDNA through reverse transcription and amplify, and then break the cDNA and sequence short fragments from the 3' end of the molecule. Compared with the uniform coverage of sequencing reads on genes obtained by the full-length library construction method, the strand-specific sequencing reads obtained by the 3'tag-based single-cell transcriptome sequencing method will be enriched at the 3' end of the gene. Previous studies have found that the enrichment position of 3'tag-based single-cell transcriptome sequencing reads is usually close to the polyadenylation site at the end of the gene. This suggests that we can use 3'tag-based single-cell transcriptome sequencing data to identify polyadenylation sites at the genome-wide level in different single cells / cell types and conduct biological studies.

[0004] The maturation of eukaryotic transcripts usually undergoes polyadenylation processing and modification, and there are widely variable polyadenylation processing and modification forms, resulting in transcript isoforms with different polyadenylation sites. The dynamic regulation of variable polyadenylation is involved in many important biological processes, such as transcript splicing, protein translation localization, tissue-specific expression, disease occurrence and cell development. At present, a variety of experimental and computational methods have been developed to study the polyadenylation modification of the whole genome at the cell population level. However, the widespread use of these methods has limitations. First, the high-throughput sequencing experimental methods developed for the identification of polyadenylation sites usually rely on special compounds, library construction processes and data processing processes, and it is difficult to compare different methods. Secondly, the polyadenylation analysis tools developed based on bulk RNA-seq have limitations in use, such as the need for multi-group experimental design, the length limit of the 3'UTR of the gene itself, or the need to combine additional polyadenylation sequencing data. Finally, there are studies that use machine learning or deep learning models to predict polyadenylation sites on the genome, such as DeepPASTA and APARENT, but these tools still have limitations such as high dependence on the location of polyadenylation sites and insufficient accuracy.

[0005] Recently, research results have shown that single-cell transcriptome sequencing technology has the potential to study polyadenylation. Velten et al. developed BATSeq to analyze single-cell subtypes and found the use of different 3'-end transcripts in mouse embryonic stem cells; Ye et al. developed a histogram difference-based method scDAPA to detect genes with dynamic variable polyadenylation in acute myeloid leukemia. Two recently published single-cell transcriptome sequencing polyadenylation analysis tools, scAPA and Sierra, have achieved genomic polyadenylation site identification and single-cell level polyadenylation quantitative analysis based on peak recognition methods, and verified the dynamic changes of variable polyadenylation and different transcripts in different biological systems, indicating that single-cell transcriptome sequencing can be applied to the study of polyadenylation.

[0006] However, there are still some problems with the above technologies and tools, such as:

[0007] (1) BATSeq analysis of polyadenylation relies on its special library construction and sequencing method, and its applicability is relatively narrow.

[0008] (2) scDAPA determines whether there is variable polyadenylation processing by the histogram changes of sequencing reads on the gene, but it is unable to identify specific genomic polyadenylation sites and conduct single-cell level quantitative analysis.

[0009] (3) scAPA and Sierra use peak identification-based methods, but both perform peak identification based on the alignment results of the entire sequencing reads of the sequencing samples, ignoring the high variability of gene expression signals in single-cell transcriptome sequencing, and lack the sensitivity to analyze relatively low-abundance genes.

[0010] (4) scAPA is limited to using the built-in gene annotations in the software and cannot enable researchers to customize the analysis of research species and other gene annotations. At the same time, scAPA does not support the identification of transcript splicing sites and has limitations in accurate quantification.

[0011] (5) Sierra can identify transcript splicing sites, but like scAPA, both tools lack confidence screening for peaks and have false positive problems in the identification of polyadenylation sites. Since 3'tag-based single-cell transcriptome sequencing requires the capture of RNA molecules through oligo(dT) primers, this process will not only capture RNA with polyadenylic acid tails, but also RNA rich in adenylic acid inside the molecule, introducing false positive interference from internal priming.

[0012] (6) Polyadenylation site prediction tools such as DeepPASTA and APARENT also have limitations such as high dependence on base position and poor accuracy.

[0013] (7) The above tools are also unable to identify transcript expression through polyadenylation sites, and lack single-cell analysis methods that couple polyadenylation sites with transcript expression.

[0014] Therefore, there is an urgent need in this field to develop a method that can use 3'tag-based single-cell transcriptome sequencing to accurately identify and screen high-confidence polyadenylation sites at the single-cell level for polyadenylation site detection and cell typing analysis of single-cell sequencing data. Summary of the invention

[0015] The purpose of the present invention is to provide a single-cell transcriptome sequencing analysis system, which can use single-cell transcriptome sequencing to accurately identify and screen high-confidence polyadenylation sites at the single-cell level, and can use polyadenylation site information to identify different transcript expressions, and further be used for single-cell typing analysis at the transcript level.

[0016] In a first aspect of the present invention, a method for constructing a prediction model of polyadenylation sites is provided, the method comprising the steps of:

[0017] (S1) providing a data set for model training, the data set comprising a positive set and a negative set, wherein the sequences contained in the positive set are sequences containing polyadenylation sites, i.e., positive sequences; and the sequences contained in the negative set are sequences not containing polyadenylation sites, i.e., negative sequences;

[0018] Wherein, the positive sequence includes: a basic positive sequence for the original polyadenylation site and a shifted positive sequence obtained by shift amplification based on the basic positive sequence, and both the basic positive sequence and the shifted positive sequence contain the original polyadenylation site;

[0019] (S2) using the dataset in (S1) to train the model architecture;

[0020] Wherein, the model architecture includes a model architecture based on a convolutional neural network (CNN) module and a recurrent neural network (RNN) module;

[0021] The CNN module includes two convolutional layers and a max pooling layer, wherein the output result of the convolutional layer is compressed by dimensionality reduction through the max pooling layer, and the maximum value in a 4×4 neighborhood is selected and output to the downstream RNN module;

[0022] The RNN module is provided with a bidirectional long short-term memory network structure (BiLSTM), and transmits the output content of the CNN module to the downstream fully connected layer;

[0023] The fully connected layer is provided with a plurality of hidden units and a predetermined dropout ratio (e.g., about 30%), and outputs the final binary classification result through a normalized exponential function (softmax activation function);

[0024] (S3) When the training result reaches a predetermined termination condition, the model training is terminated, and the prediction model of the polyadenylation site (ie, the DeepPASS model) is obtained.

[0025] In another preferred embodiment, the offset amplification includes: random or non-random offset amplification.

[0026] In another preferred embodiment, the offset amplification includes: normal distribution random offset.

[0027] In another preferred embodiment, the offset amplification refers to: placing a basic positive sequence on a coordinate axis x1, and taking the polyadenylated original genomic site of the basic positive sequence as the coordinate zero point, and then making the 5' end and 3' end of the basic sequence independently perform a positive or negative offset of p bp, and the offset sequence is called an offset positive sequence, wherein p is a positive integer (bp).

[0028] In another preferred example, the length of the basic positive sequence and the length of the corresponding offset positive sequence are the same or substantially the same.

[0029] In another preferred embodiment, when the coordinates of the basic positive sequence are (X Left , X Right ) and the lengths of the base positive sequence and the corresponding offset positive sequence are the same, then after a positive offset of p bps, the coordinates of the offset positive sequence are (X Left +p,X Right +p); and after a negative shift of p bp, the coordinates of the shifted positive sequence are (X Left -p,X Right -p).

[0030] In another preferred example, assuming that the length of the basic positive sequence is 2L bp, p is ≤1 / 4L, preferably ≤1 / 5L, and more preferably ≤1 / 10L.

[0031] In another preferred embodiment, the original polyadenylation site of the base positive sequence is set as the coordinate zero point, and the length of the base positive sequence is 2L bp, L is about 100 bp, then the coordinate of the base positive sequence is (X Left , X Right ), X Right -X Left =2L.

[0032] In another preferred embodiment, L is 60-250, preferably 80-150, more preferably 90-110.

[0033] In another preferred embodiment, |X Left |with|X Right | are the same or substantially the same, for example |X Left |with|X Right | is 0.8-1.2, preferably 0.9-1.1, and more preferably about 1.0.

[0034] In another preferred example, in the positive set, the ratio of the total number of basic positive sequences to the total number of offset positive sequences is 1:2 to 1:20, preferably 1:5-1:10, and more preferably 1:10-1:20.

[0035] In another preferred embodiment, the positive sequences contained in the positive set include T0 (about 250,000) polyadenylation original sites.

[0036] In another preferred embodiment, in the positive set, for each polyadenylation original site, it corresponds to 1 basic positive sequence and D offset positive sequences, wherein D is a positive integer ≥ 2, preferably D is 2-5, more preferably 5-10, and most preferably 10-20.

[0037] In another preferred embodiment, in the positive set, the total number T1 of the basic positive sequences is the same or substantially the same as the total number T0 of the original polyadenylation sites.

[0038] In another preferred embodiment, the sequences included in the positive set and the negative set are from the same species or different species.

[0039] In another preferred embodiment, the sequences contained in the positive set and the negative set are from the same species, such as humans.

[0040] In another preferred example, the positive set includes N1 sequences containing polyadenylation sites (positive sequences); and the negative set includes N2 sequences not containing polyadenylation sites (negative sequences).

[0041] In another preferred embodiment, the ratio of the number of sequences N1 in the positive set to the number of sequences N2 in the negative set is not particularly limited, and is generally 1:10 to 10:1, preferably 1:5 to 5:1, more preferably 1:2 to 2:1, such as about 1:1.

[0042] In another preferred embodiment, N1 is ≥ 500,000, more preferably ≥ 1,000,000, and most preferably ≥ 2,000,000.

[0043] In another preferred embodiment, N2 is ≥ 500,000, more preferably ≥ 1,000,000, and most preferably ≥ 2,000,000.

[0044] In another preferred embodiment, the prediction model is a prediction model for constructing polyadenylation sites in the present invention (hereinafter referred to as DeepPASS).

[0045] In another preferred embodiment, the prediction model is used to determine whether the input sequence has a high-confidence polyadenylation site.

[0046] In another preferred example, the sequence may be a sequence in RNA or DNA format, or the sequencing data may be RNA sequencing data.

[0047] In another preferred embodiment, the high confidence is the prediction that the DeepPASS output result is 1 (containing a polyadenylation site).

[0048] In another preferred example, in step (S1), the offset amplification includes: taking the original polyadenylation genomic site as the coordinate zero point, randomly offsetting the site using a normal distribution, offsetting a polyadenylation site 4-20 times, and extracting the genomic sequence from 100±20bp upstream to 100±20bp downstream of the site after the offset amplification as the offset positive sequence in the positive set of model training.

[0049] In another preferred example, the offset amplification includes: taking the original polyadenylation site as the coordinate zero point, randomly offsetting the site using a normal distribution (μ=0, σ=10), performing about 10 offsets for each polyadenylation site (ten times the number of amplification sites), and extracting the genomic sequence from 100bp upstream to 100bp downstream of the amplified site as the offset positive sequence in the positive set of DeepPASS model training.

[0050] In another preferred embodiment, four polyadenylation site databases from different sources are collected: PolyA_DB3, PolyA-seq, Poly(A)Site2.0 and GENCODE. The common sites in the three published databases (PolyA_DB3, PolyA-seq and Poly(A)Site2.0) are screened and merged with GENCODE as a known polyadenylation site library for DeepPASS model training.

[0051] In another preferred embodiment, four polyadenylation site databases from different sources are collected: PolyA_DB3, PolyA-seq, Poly(A)Site2.0 and GENCODE. Highly reliable polyadenylation sites in the three published databases (PolyA_DB3, PolyA-seq and Poly(A)Site2.0) are screened, and the screening condition is: the polyadenylation sites in one database overlap with those in the other two databases within 12bp upstream and downstream of its genomic position. The highly reliable sites in these three databases are merged with GENCODE as a known polyadenylation site library for polyadenylation site annotation.

[0052] In another preferred embodiment, the negative sequences in the negative set are screened in the following manner: by taking the union of four polyadenylation site databases, the intergenic regions of the genome are reversely screened, and only the intergenic region sequences outside the upstream and downstream distance of 100 bp of the known polyadenylation sites are selected.

[0053] In another preferred example, the DeepPASS model selects the classification with a higher probability in the two classifications as the prediction result for output, that is, the classification result of 1 (containing polyadenylation site) or 0 (not containing polyadenylation site) is used to judge whether the input sequence contains a high-confidence polyadenylation site.

[0054] In another preferred embodiment, in step (S2), the model training sets the data batch size to 3000-7000 (such as 5000), and the model is trained for 50-150 periods (such as 100 periods).

[0055] In another preferred embodiment, the predetermined termination condition includes: there is no improvement in the training results within 5-20 periods.

[0056] In another preferred embodiment, in step (S2), it also includes: evaluating the training results of the model in training, preferably the evaluation is performed using a receiver operating characteristic curve (ROC) and an area under the curve (AUC).

[0057] In another preferred example, the chimeric deep learning model DeepPASS based on CNN and RNN can theoretically identify sequence features such as motifs near polyadenylation sites, such as AAUAAA and upstream and downstream GU repeats, U enrichment, etc., and establish potential correlations between motif features to achieve the effect of overall recognition of polyadenylation signals.

[0058] In another preferred embodiment, the DeepPASS model prediction only requires 200 bp of base sequence as input, and the output is a binary (1 or 0) polyadenylation site judgment.

[0059] In another preferred example, the evaluation of the DeepPASS model and other models uses the sequences from 100bp upstream to 100bp downstream of the original polyadenylation site in the GENCODE database as the positive test set, and the sequences of 200bp in the intergenomic region beyond 100bp upstream and downstream of the known polyadenylation site as the negative test set.

[0060] In a second aspect of the present invention, a system for predicting polyadenylation sites is provided, the system comprising:

[0061] An input unit, wherein the input unit is configured to input data, wherein the data includes an input sequence to be predicted, wherein it is necessary to determine whether a polyadenylation site exists in the input sequence to be predicted;

[0062] a polyadenylation site prediction unit, the prediction unit being configured to execute a polyadenylation site prediction model, thereby obtaining a prediction result of whether a polyadenylation site exists in the input sequence to be predicted; wherein the prediction model is constructed using the method described in the first aspect of the present invention;

[0063] An output unit is configured to output the prediction result of the polyadenylation site prediction unit.

[0064] In another preferred embodiment, the prediction model is a DeepPASS model, which selects the classification with a larger probability in the two classifications as the prediction result, that is, it predicts whether a given input sequence contains a real polyadenylation site through the classification results of "yes" (or 1) or "no" (or 0).

[0065] In a third aspect of the present invention, a method for analyzing polyadenylation sequencing read signal peaks based on single-cell transcriptome sequencing data is provided, characterized in that it comprises the steps of:

[0066] (Y1) Aligning the single-cell transcriptome sequencing data with the genome sequence to obtain a sequencing read alignment result file stored in the BAM format. The BAM file should contain: sequencing read alignment result information required by the BAM format, sequencing read cell barcode (CB) and sequencing read molecular identifier (UMI) information;

[0067] (Y2) converting the gene-level annotation file into a gene-transcript-level annotation file labeled “gene-transcript” according to its gene and transcript information;

[0068] (Y3) Based on the sequencing read alignment result file and the gene transcript level annotation file, HOMERfindPeaks is used to identify the sequencing read distribution peaks at the transcript level, identify all peaks corresponding to each transcript, and mark the interval position of the peak based on the coordinate information of the corresponding transcript;

[0069] (Y4) combining the gene transcript level annotation file obtained in the above steps (Y2)-(Y3) and the peak interval position based on the transcript coordinates, mapping the transcript position system to the genome coordinate position system, and calculating the genome coordinate system annotation file of the peak;

[0070] (Y5) combining the gene transcript level annotation file obtained in the above steps (Y2)-(Y4) and the peak file based on the genome coordinate system, further calculating the peak splicing site information of the corresponding transcript by identifying the transcript exon splicing site in the gene transcript level annotation file;

[0071] (Y6) performing a normality test on the sequencing read signal distribution curve of the identified peaks according to the peak annotation results obtained in the above step (Y5), removing peaks that may be caused by abnormal amplification, multiple alignments or low signal abundance; and performing a unimodality test to retain the overlapping peaks formed by adjacent polyadenylation signals;

[0072] (Y7) The above steps (Y2)-(Y6) can calculate peaks from different transcript sources of the same gene, and the peaks can be deduplicated according to their gene information; preferably, the deduplication method is: for peaks with interval overlap ratio exceeding 50%, the abundance of their sequencing reads is sorted, and the peak with the highest expression level is retained;

[0073] (Y8) recalculating the abundance of sequencing reads for the peaks after deduplication in the above step (Y7), filtering out the peaks whose sequencing reads account for ≤1% in the same gene, and obtaining the preliminary screened peaks (raw peaks);

[0074] (Y9) for each raw peak obtained from the steps (Y3) to (Y8) in the plurality of samples, by marking the sample information of the peak source, and marking the peaks of different samples with an overlap ratio of more than 60% as the same group, a unique value is selected from the peaks in the same group as the peak of the polyadenylation site;

[0075] (Y10) The end of the raw peak annotation interval identified in the above step is approximated as the polyadenylation site, and the sequence from 100bp upstream to 100bp downstream of the site is extracted and input into the DeepPASS model for prediction, to determine whether each raw peak is a high-confidence polyadenylation peak (PAS), thereby obtaining the polyadenylation sequencing read signal peak analysis results of the single-cell transcriptome sequencing data.

[0076] In another preferred embodiment, in step (Y10), the prediction model is a prediction model constructed using the method described in the first aspect of the present invention, namely, the DeepPASS model.

[0077] In another preferred example, in step (Y10), the base sequence from 100 bp upstream to 100 bp downstream of the polyadenylation site is input into the prediction model for prediction.

[0078] In another preferred embodiment, the DeepPASS model is used for prediction to determine whether each raw peak is a high-confidence polyadenylation peak (PAS), and the selected high-confidence peaks are used for downstream analysis.

[0079] In another preferred embodiment, the analysis method is used to analyze polyadenylation sequencing read signal peaks at the whole genome level.

[0080] In another preferred embodiment, in step (Y2), the gene annotation is preprocessed, and an annotation file based on transcripts is constructed using the correspondence between genes and transcripts to replace the original annotation file based on genes.

[0081] In another preferred embodiment, in step (Y3), the peak is identified in units of transcripts and the annotation interval of the transcript coordinates is obtained. At the same time, the peak width is set to a fixed length (such as 400bp) according to a priori values ​​to make it suitable for single-cell transcriptome sequencing data of 10xChromium technology. It can also be set to other length values ​​according to different sequencing methods.

[0082] In another preferred example, in steps (Y4)-(Y5), the annotation interval markers of the peak are mapped from the transcript position to the genomic position, and the splicing sites of the exons in the transcript annotations are parsed, and the splicing sites of the transcripts corresponding to the peaks are used to further annotate the peaks to obtain genomic coordinate peak annotations that identify the splicing information of the sequencing reads.

[0083] In another preferred embodiment, in step (Y6), the normal distribution test and the unimodal test are introduced for the quality control of the peak. Each base position in the peak interval is taken as the independent variable x, and the sequencing read depth of each base is taken as the dependent variable y, and a distribution curve is made according to the values ​​of x and y; the normal distribution test is used to evaluate whether the curve is close to the normal distribution, and the abnormally amplified transcript fragments in the library construction, the abnormal alignment in the genome alignment process, and the peaks with low abundance and coverage are removed and optimized; the unimodal test is used to further evaluate whether the curve conforms to the multimodal distribution, and the overlapping peaks formed by adjacent polyadenylation signals are retained.

[0084] In another preferred embodiment, in step (Y7), for the sequencing read distribution peak generated by the same polyadenylation site, multiple transcript annotations covering the tailing site in the same gene annotation will calculate multiple peaks, and the method for de-redundanting these peaks is: by calculating the peak interval overlap ratio, the peaks with an overlap ratio greater than 50% are classified and grouped, and for all overlapping peaks in the same group, only the peak annotation with the highest sequencing read coverage depth is retained.

[0085] In another preferred embodiment, in step (Y8), the abundance of peak sequencing reads within the same gene is estimated using stringtie, and peaks with a read abundance ratio ≥ 1% are retained.

[0086] In another preferred example, in step (Y9), the end position of the identified raw peak genomic annotation interval can be used as an approximate polyadenylation site, and the base sequence from 100bp upstream to 100bp downstream extracted from the site can be input into the DeepPASS model for polyadenylation prediction, and the peak with a prediction result of 1 is used as a highly reliable polyadenylation site.

[0087] In another preferred embodiment, the single-cell transcriptome data is 3'tag-based single-cell sequencing data, which is a high-throughput single-cell sequencing technology that relies on oligo(dT) to capture RNA molecules and enrich the 3' ends of RNA molecules.

[0088] In a fourth aspect of the present invention, a method for quantitatively analyzing transcripts and / or performing single-cell typing analysis at a single-cell level is provided, the method comprising:

[0089] (Z1) providing single-cell transcriptome sequencing data of a single-cell sample to be analyzed;

[0090] (Z2) analyzing the single-cell transcriptome sequencing data using the analysis method described in the third aspect of the present invention to obtain a high-confidence peak;

[0091] (Z3) Based on the high-confidence peaks obtained in the previous step, perform quantitative analysis of transcripts and / or perform single-cell typing analysis.

[0092] In another preferred embodiment, in step (Z3), comprising:

[0093] (W1) For the peak integration analysis of multiple samples, the integration method is: mark the overlapping peaks between different samples as the same group, sort the peaks in the same group by confidence, and select the peak with the highest confidence as the only peak of the polyadenylation site;

[0094] (W2) According to the peaks selected in the above step (W1), the peaks are mapped to the transcripts of the corresponding sources through the splicing sites and genomic interval information of the peaks, so as to identify the transcripts generated by variable polyadenylation (APA) in the single-cell transcriptome sequencing samples, thereby realizing the identification of transcript expression;

[0095] (W3) Construct a GTF format annotation file with the PAS annotations of the mapped transcript information, and recalculate the distribution of PAS annotations in the aligned BAM reads through featureCounts to obtain a re-assigned BAM file marked with PAS annotations; UMI-tools can be used to count the unique molecular identifier (UMI) information of all PAS annotations corresponding to each cell barcode in the re-assigned BAM file to obtain a transcript expression matrix at the single-cell level. The matrix includes the number of unique molecular identifiers of each transcript in each cell, which can be used for single-cell transcriptome analysis;

[0096] (W4) For single-cell analysis at the transcript level, the peak-annotated transcript level expression matrix obtained in the above step (W3) is used as the input of the single-cell analysis tool Seurat for single-cell analysis; preferably, the single-cell analysis system based on the peak-annotated transcript level includes the following main parts: expression quality control, feature extraction, dimensionality reduction analysis, unsupervised cell clustering, and transcript molecular marker expression analysis;

[0097] (W5) performing expression quality control on the input transcript expression matrix, and retaining transcripts with at least ≥0.5% of the expressed cells; preferably, using the built-in function of Seurat to standardize the transcript expression and calculate the expression variation of the screened expression matrix, and using the 2000 highly variable transcripts with the highest expression variation as feature transcripts for downstream dimensionality reduction analysis;

[0098] (W6) The characteristic transcripts obtained in the above step (W5) can be subjected to dimensionality reduction analysis and unsupervised cell clustering through the built-in functions of Seurat; preferably, the principal component analysis method (PCA) is used for dimensionality reduction, and the main (such as the first 50) principal component results are retained; a shared nearest neighbor (SNN) similarity graph is constructed for the overall single cells based on the 50 principal components, and cell clustering is performed in an unsupervised manner at a specific resolution (such as 0.3), and finally a cell classification result based on the expression of PAS annotated transcripts is obtained;

[0099] (W7) Based on the peak annotation information obtained in steps (W2) and (W3), the transcripts expressed by the given molecular marker genes are calculated and screened. The expression levels of the molecular marker transcripts of different cell classifications can be calculated and visualized using the built-in functions of Seurat to identify the cell types in the unsupervised cell classification results.

[0100] In another preferred embodiment, in step (W1), for the peak integration analysis of multiple samples, the integration method is: (1) for the peaks of different samples originating from the same polyadenylation site, the number of known polyadenylation sites covered by each peak end interval (50 bp upstream to 25 bp downstream) in the interval is calculated, and the peaks are sorted in descending order according to the numerical value of the number of covered sites; (2) the interval sequence from 100 bp upstream to 100 bp downstream of the peak end is extracted, input into the DeepPASS model for prediction, and sorted from 1 to 0 according to the prediction results; (3) in the same group of peaks, the peak with the best comprehensive ranking is selected as the only peak of the polyadenylation site in the multi-sample analysis.

[0101] In another preferred embodiment, in step (W2), the peak splicing site and genome annotation information are used to identify transcripts that match the peak source; preferably, the identification method is: the transcript exon splicing site must completely match the peak splicing site information, and the transcript with the closest end annotation distance is selected as the transcript of the peak source, and finally the expression of different transcripts in the single-cell transcriptome sequencing sample at the whole genome level is obtained.

[0102] In another preferred example, in step (W3), the PAS annotation information and transcript source information are used to construct a PAS transcription annotation file in Genetransfer format (GTF), and the sequencing reads of the BAM file are redistributed in the PAS transcription annotation file in the GTF format, marking the PAS annotation to which each sequencing read belongs.

[0103] In another preferred example, in step (W3), the sequencing reads of each PAS interval of each single cell are calculated according to the PAS transcription annotation file, and the cumulative number of unique molecular identifiers on the sequencing reads is further counted to obtain a single-cell transcript expression matrix based on PAS annotation.

[0104] In another preferred embodiment, in steps (W4)-(W6), the transcript level expression matrix annotated by the peak in step (W3) is used instead of the traditional gene level expression matrix to carry out cell typing analysis of single-cell transcriptome data based on transcript quantitative data; preferably, the analysis process includes expression quality control, feature extraction, dimensionality reduction analysis and unsupervised cell clustering and transcript molecular marker expression analysis based on transcript quantitative data.

[0105] In another preferred embodiment, the single-cell transcriptome data is 3'tag-based single-cell sequencing data, including a high-throughput single-cell sequencing technology that relies on oligo(dT) to capture RNA molecules and enrich the 3' ends of RNA molecules.

[0106] In a fifth aspect of the present invention, a single-cell transcriptome sequencing data analysis system is provided, the system comprising:

[0107] (a) an input unit, wherein the input unit is configured to input sequence data, wherein the input sequence data is selected from the following group: (i) a DeepPASS prediction model and / or an input sequence to be predicted, wherein it is necessary to determine whether a polyadenylation site exists in the input sequence to be predicted; (ii) single-cell transcriptome sequencing data; (iii) a combination of (i) and (ii);

[0108] (b) an analysis module, the analysis module being configured to perform a predetermined analysis on the input data to obtain an analysis result, and the analysis module is selected from the following group:

[0109] (M1) a prediction unit for polyadenylation sites, wherein the prediction unit (M1) is configured to execute a prediction model for polyadenylation sites, thereby obtaining a prediction result of whether a polyadenylation site exists in the input sequence to be predicted; wherein the prediction model is constructed using the method described in the first aspect of the present invention;

[0110] (M2) a polyadenylation sequencing read signal peak analysis unit, wherein the polyadenylation sequencing read signal peak analysis unit (M1) is configured to perform the analysis method described in the third aspect of the present invention, thereby obtaining a polyadenylation sequencing read signal peak analysis result of the single-cell transcriptome sequencing data;

[0111] (M3) single-cell level transcript quantitative analysis and / or single-cell typing analysis unit, wherein the (M3) single-cell level transcript quantitative analysis and / or single-cell typing analysis unit is configured to perform the analysis method described in the fifth aspect of the present invention, thereby obtaining single-cell level transcript quantitative analysis results and / or single-cell typing analysis results;

[0112] (M4) A combination of the above (M1), (M2) and (M3).

[0113] (c) an output unit, wherein the output unit is configured to output the analysis result of the analysis module.

[0114] It should be understood that within the scope of the present invention, the above-mentioned technical features of the present invention and the technical features specifically described below (such as embodiments) can be combined with each other to form a new or preferred technical solution. Due to space limitations, they will not be described one by one here. BRIEF DESCRIPTION OF THE DRAWINGS

[0115] Figure 1The data collection, preprocessing, model building and training evaluation of the DeepPASS model in Example 1 are shown.

[0116] Figure 2 The identification and screening of polyadenylation site sequencing read peaks of single-cell transcriptome data in Example 2 are shown.

[0117] Figure 3 The identification of transcript expression based on polyadenylation sites and application to single cell typing analysis in Example 3 are shown. DETAILED DESCRIPTION

[0118] After extensive and in-depth research, the inventors have developed a single-cell transcriptome computational analysis method and system based on a fusion deep learning model for the first time. Specifically, the inventors have designed a polyadenylation site (PAS) detection method that integrates a deep learning model for polyadenylation site detection and cell typing analysis of single-cell sequencing data (hereinafter referred to as the computational analysis system SCAPTURE). First, a deep learning model is constructed and integrated in the SCAPTURE system of the present invention, achieving high accuracy and non-position-dependent polyadenylation site prediction, and used for the screening of highly reliable polyadenylation sites at the whole genome level; second, the SCAPTURE system of the present invention can also use single-cell transcriptome sequencing data to detect the sequencing read distribution peak of the whole genome region and identify the polyadenylation sites at the whole genome level from scratch; third, the SCAPTURE system of the present invention also identifies different transcripts and their differential expression on the same gene based on variable polyadenylation sites, and is applied to single-cell typing analysis. The present invention was completed on this basis.

[0119] the term

[0120] As used herein, the term "the method of the present invention" refers to a single-cell transcriptome computational analysis method based on a fusion deep learning model. It should be understood that the method of the present invention includes any one of the following three aspects or a combination thereof: (a) constructing a DeepPASS model to predict and screen highly reliable polyadenylation sites; (b) identifying and screening polyadenylation sequencing read sites at the whole genome level from single-cell transcriptome data; (c) identifying different transcripts and their differential expression in single-cell transcriptome sequencing data based on polyadenylation sites, and performing transcript quantification and single-cell typing analysis at the single-cell level.

[0121] As used herein, the term "DeepPASS Deep neural network for p oly a denylationprediction withs equences s The term "hyper-elevation model" refers to a model for predicting polyadenylation sites constructed using a deep learning method in the method of the present invention.

[0122] As used herein, the term "SCAPTURE sc RNA-seq a analysis of P ASs and their corresponding t ranscript expression u sed to re "fine cell identities)" specifically refers to the analysis system established by the present invention for identifying and predicting polyadenylation tailing sites using single-cell transcriptomes and quantifying selectively expressed transcripts at the single-cell level, including the following three aspects: (a) constructing a DeepPASS model to predict and screen highly reliable polyadenylation sites; (b) de novo identifying and screening polyadenylation sequencing read sites at the whole genome level from single-cell transcriptome data; (c) identifying different transcripts and their differential expression from single-cell transcriptome sequencing data based on polyadenylation sites, and performing transcript quantification and single-cell typing analysis at the single-cell level.

[0123] Single-cell RNA sequencing (scRNA-seq) is a new technology that amplifies and performs high-throughput sequencing on the whole transcriptome at the level of a single cell. It studies the gene expression at the overall level of a single cell, aiming to solve the common problems in the study of cell molecular mechanisms, such as cell heterogeneity and the inability to perform conventional high-throughput sequencing due to small cell numbers.

[0124] 3'tag-based single-cell transcriptome sequencing (3'tag-based scRNA-seq): A single-cell sequencing method that uses oligo(dT) primers with unique molecular identifiers and cell barcodes to capture and reverse transcribe single-cell RNA molecules, amplify and build cDNA libraries, and perform 3'-end specific sequencing of cDNA.

[0125] Sequence Read: A single sequenced fragment generated by sequencing technology. In this specification, it refers to the cDNA sequenced fragment in 3'tag-based single-cell transcriptome sequencing.

[0126] Cell barcode (CB): A short deoxyribonucleotide sequence, usually 10-20 nt in length. During single-cell transcriptome sequencing, RNA molecules in the same cell are labeled. Different cells have different cell labels to distinguish the cell sources of sequencing reads.

[0127] Unique molecular identifier (UMI): A short deoxyribonucleotide sequence, usually 8-15 nt in length. It is used to track amplification products from the same transcript molecule source, to exclude duplication and bias introduced by PCR amplification and sequencing, and to obtain accurate readings for quantitative analysis.

[0128] Gene: In biology, it refers to the nucleotide sequence in DNA or RNA that encodes the synthesis of gene products (RNA or protein). Gene can also be regarded as the basic genetic unit, that is, a functional DNA or RNA sequence.

[0129] Transcript: One or more mature RNA molecules formed by transcription of a gene. The RNA formed by gene transcription can form different transcripts through different processing processes.

[0130] Exon: A region of DNA in a eukaryotic gene that can be transcribed and preserved after splicing. In this manual, it refers to the name of the exon used to analyze exon splicing sites.

[0131] Intron: A non-coding DNA segment in a gene that separates adjacent exons. Introns are sequences that block the linear expression of a gene. Introns on DNA are transcribed into pre-mRNA, but introns on RNA are spliced ​​out before the RNA leaves the nucleus for translation.

[0132] Gene annotation file: In this manual, it refers to the annotation file format that records the name of the gene and its transcript on the genome, the gene location and exon splicing information, which can be in genepred, GTF or BED format.

[0133] Peak annotation file: In this manual, it refers to the annotation file format that records the gene and source transcript name on the genome, the location of the sequencing read peak of the gene polyadenylation site in single-cell transcriptome sequencing, and the exon splicing information. It can be in genepred, GTF or BED format.

[0134] Polyadenylation: The process of adding a poly(A) tail to an RNA transcript. The poly(A) tail is composed of multiple adenosine monophosphates and is a stretch of RNA that contains only adenine bases. In eukaryotes, polyadenylation is part of the process of producing mature mRNA for translation.

[0135] Polyadenylation site (PAS): The site where the poly(A) tail is added during RNA polyadenylation processing. It is located downstream of the polyadenylation signal. The cleavage site is usually but not limited to the CA base.

[0136] Polyadenylation signal: The classic core sequence is AAUAAA, which is usually located 20-30bp upstream of the polyadenylation signal and can be recognized by cleavage and polyadenylation specific factors.

[0137] Alternative polyadenylation (APA): Alternative polyadenylation is a widespread gene regulatory mechanism that can generate different 3' ends in transcripts produced by RNA polymerase II transcription.

[0138] Constructing the DeepPASS model to predict and screen highly reliable polyadenylation sites

[0139] The present invention provides a method and system for predicting and screening highly reliable polyadenylation sites based on a deep learning model. The model (method or system) for predicting or determining highly reliable polyadenylation sites of the present invention is called the "DeepPASS model".

[0140] Typically, in the present invention, the biological characteristics of polyadenylation are extracted by constructing a deep learning model DeepPASS, and then the polyadenylation sites with high confidence are screened out.

[0141] In the present invention, DeepPASS model construction includes four parts: data collection, data preprocessing, model building and model training evaluation.

[0142] In a preferred embodiment of the present invention, the construction of the DeepPASS model includes the following steps (see Figure 1 A, 1B and 1C):

[0143] (1) For DeepPASS model training data, the present invention collects polyadenylation site databases from four different sources: PolyA_DB3, PolyA-seq, Poly(A)Site2.0 and GENCODE. Highly reliable polyadenylation sites in the three published databases (PolyA_DB3, PolyA-seq and Poly(A)Site2.0) are screened, and the screening conditions are: the polyadenylation sites in one database overlap with those in the other two databases within 12bp upstream and downstream of its genomic position. The highly reliable sites in these three databases are merged with GENCODE as a known polyadenylation site library. The polyadenylation site library is as follows: Figure 1 As shown in A.

[0144] (2) Based on the polyadenylation site library of step (1) above, the present invention uses a method based on normal distribution random offset to perform data amplification. The original polyadenylation site is used as the coordinate zero point, and the site is randomly offset using normal distribution (μ=0, σ=10). Each polyadenylation site is offset about 10 times (ten times the number of amplified sites), and the genomic sequence from 100bp upstream to 100bp downstream of the amplified site is extracted as the positive set for DeepPASS model training. The schematic diagram of this step can be seen in Figure 1 B.

[0145] (3) All polyadenylation sites in the four databases PolyA_DB3, PolyA-seq, Poly(A)Site2.0 and GENCODE (data amplification was performed using a normal distribution random offset method in step (2)) were taken as the union as the positive set (about 1,907,290 sequences). In addition, sequences (200 bp in length) in the intergenic regions of the human genome that did not contain polyadenylation sites were extracted as the negative set (about 1,907,000 sequences) for DeepPASS model training.

[0146] Preferably, in the present invention, the number of sequences N1 in the positive set and the number of sequences N2 in the negative set are each independently ≥ 500,000, preferably ≥ 1,000,000, and more preferably ≥ 2,000,000. In addition, the ratio of the number of sequences N1 in the positive set to the number of sequences N2 in the negative set is not particularly limited, and is generally 1:10 to 10:1, preferably 1:5 to 5:1, and more preferably 1:2 to 2:1, such as about 1:1.

[0147] (4) For the training and construction of the DeepPASS model, the present invention mainly adopts the convolutional neural network (CNN) module and the recurrent neural network (RNN) module as the model architecture. CNN includes two convolutional layers and a maximum pooling layer (max pooling layer). The first convolutional layer includes 12-mer width, 4 channels, and a total of 128 filters; the second convolutional layer includes 6-mer, 64 filters can cover the 128 channels output by the first convolutional layer. The output result of the convolutional layer is compressed by the maximum pooling layer (Maxpooling layer), and the maximum value in the 4×4 neighborhood is selected and output to the downstream RNN layer. The RNN adopts a bidirectional long short-term memory network structure (BiLSTM), sets 128 units to process the output content of the CNN module, and passes the result to the downstream fully connected layer. The fully connected layer sets 1024 hidden units and a dropout ratio of 30%, and outputs the final binary classification probability through a normalized exponential function (softmax activation function). The DeepPASS model selects the classification with a higher probability in the binary classification as the prediction result for output, that is, it predicts whether the given input sequence contains a true polyadenylation site through the classification result of 1 or 0.

[0148] (5) The positive and negative set sequences obtained in the above steps (2)-(3) are used as DeepPASS model training data, and the model architecture of the above step (4) is used for training. The model training data batch size is set to 5000, and the model training is performed for 100 periods. If the training results do not improve within 10 periods, the model training is terminated early. The model training results are evaluated using the receiver operating characteristic curve (ROC) and the area under the curve (AUC).

[0149] The effect of the DeepPASS model of the present invention can be further evaluated. In the present invention, the inventors compared DeepPASS with similar polyadenylation site prediction tools (DeepPASTA and APARENT) that have been published so far and have good performance. That is, according to the highly reliable polyadenylation site (site 49942 of GENCODE) obtained in the above step (1), the sequences 100bp upstream and 100bp downstream of the original site were extracted as the positive test set. The 200bp sequence of the intergenomic region that does not contain the polyadenylation site obtained in the above step (3) was used as the negative test set. After merging the positive test set and the negative test set, approximately 50,000 sequences were randomly sampled for model evaluation of DeepPASS, DeepPASTA and APARENT, and the receiver operating characteristic curve (ROC) and area under the curve (AUC) were used to evaluate the results of different models. The evaluation results are shown in Figure 2. Figure 1D. The results show that the performance of the DeepPASS model of the present invention is much better than DeepPASTA and APARENT.

[0150] Analysis method of polyadenylation sequencing read signal peaks based on single-cell transcriptome sequencing data

[0151] The present invention also provides an analysis method for polyadenylation sequencing read signal peaks based on single-cell transcriptome sequencing data, which includes: peak identification at the transcript level, peak coordinates and splicing annotation, peak quality inspection and redundancy removal, and peak merging of multiple sequencing samples based on the signal peaks of the single-cell transcriptome sequencing data.

[0152] In a preferred embodiment, the analysis method of the present invention comprises:

[0153] (Y1) Aligning the single-cell transcriptome sequencing data with the genome sequence to obtain a sequencing read alignment result file stored in the BAM format. The BAM file should contain: sequencing read alignment result information and cell barcode (CB) information required by the BAM format;

[0154] (Y2) converting the gene-level annotation file from GENCODE into a gene-transcript-level annotation file labeled “gene-transcript” based on its gene and transcript information;

[0155] (Y3) Based on the sequencing read alignment result file and the gene transcript level annotation file, HOMERfindPeaks is used to identify the sequencing read distribution peaks at the transcript level, identify all peaks corresponding to each transcript, and mark the interval position of the peak based on the coordinate information of the corresponding transcript;

[0156] (Y4) combining the gene transcript level annotation file obtained in the above steps (Y2)-(Y3) and the peak interval position based on the transcript coordinates, mapping the transcript position system to the genome coordinate position system, and calculating the peak genome coordinate system annotation file corresponding to the transcript identification;

[0157] (Y5) combining the gene transcript level annotation file obtained in the above steps (Y2)-(Y4) and the peak file based on the genome coordinate system, further calculating the peak splicing site information of the corresponding transcript by identifying the transcript exon splicing site in the gene transcript level annotation file;

[0158] (Y6) Based on the peak annotation results obtained in the above step (Y5), a normality test is performed on the sequencing read signal distribution curve of the identified peak to remove peaks that may be caused by abnormal amplification, multiple alignments or low signal abundance; a unimodality test is performed at the same time to retain the overlapping peaks formed by adjacent polyadenylation signals;

[0159] (Y7) The above steps (Y2)-(Y6) can calculate peaks from different transcript sources of the same gene, and the peaks can be deduplicated based on their gene information. Specifically, for peaks with an overlap ratio of more than 50%, the signal abundance of their sequencing reads is sorted, and the peak with the highest expression level among the overlapping peaks is retained;

[0160] (Y8) recalculating the abundance of sequencing reads for the peaks in the same gene after deduplication in the above step (Y7), filtering out the peaks whose total sequencing reads account for ≤1%, and obtaining the peaks (raw peaks) for preliminary identification and screening for the evaluation of polyadenylation sites;

[0161] (Y9) For the raw peaks obtained from each of the multiple samples through the above steps (Y3)-(Y8), the present invention marks the sample information of the peak source, and marks the peaks with an overlap ratio of more than 60% between different samples as the same group, and selects a unique value from the peaks in the same group as the peak of the polyadenylation site.

[0162] (Y10) The raw peak annotation interval identified in the above step (Y9) is used as the polyadenylation site to be predicted and screened, and the genomic sequence from 100bp upstream to 100bp downstream of the annotation end is input into the DeepPASS model for prediction. It can be determined whether each raw peak is a high-confidence polyadenylation peak (PAS), and the screened high-confidence peak is used for downstream analysis.

[0163] Methods for quantifying transcripts at the single-cell level and / or performing single-cell phenotyping

[0164] The present invention also provides a method for quantitatively analyzing transcripts at the single-cell level and / or performing single-cell typing analysis.

[0165] Typically, the method involves screening the identified raw peaks using the deep learning model DeepPASS to select high-confidence polyadenylated peaks. The peaks obtained by computational screening can then be used to identify transcripts that are selectively expressed in single-cell transcriptome sequencing data, and to perform single-cell level quantification and single-cell typing analysis.

[0166] In a preferred embodiment of the present invention, the method for quantitatively analyzing transcripts and / or performing single cell typing analysis at the single cell level comprises:

[0167] (W1) For the peak integration analysis of multiple samples, the integration method is as follows: (1) For peaks from different samples at the same polyadenylation site, calculate the number of known polyadenylation sites covered by each peak end interval (50 bp upstream to 25 bp downstream) within this interval, and sort the peaks in descending order according to the number of covered sites; (2) extract the interval sequence from 100 bp upstream to 100 bp downstream of the peak end, input it into the DeepPASS model for prediction, and sort it from 1 to 0 according to the prediction results; (3) in the same group of peaks, select the peak with the best comprehensive ranking as the only peak of the polyadenylation site in the multi-sample analysis.

[0168] (W2) According to the peaks selected in the above step (W1), the peaks are mapped to the transcripts of the corresponding sources through the splicing sites and interval information of the peaks, so as to identify the transcripts generated by variable polyadenylation (APA) in the single-cell transcriptome sequencing samples and realize the identification of transcript expression.

[0169] (W3) The PAS annotations of the mapped transcript information are used to construct a GTF format annotation file, and the PAS annotation distribution of the aligned BAM reads is recalculated through featureCounts to obtain a re-assigned BAM file marked with PAS annotations. The UMI-tools can be used to count the unique molecular identifier (UMI) information of all PAS annotations corresponding to each cell barcode in the re-assigned BAM file, and finally obtain a transcript expression matrix at the single-cell level. The matrix includes the number of unique molecular identifiers of each transcript in each cell, which can be used for single-cell transcriptome analysis.

[0170] (W4) For single-cell analysis at the transcript level, the peak-annotated transcript level expression matrix obtained in the above step (W3) is used as the input of the single-cell analysis tool Seurat for single-cell analysis. The single-cell analysis system based on peak-annotated transcript levels established in the present invention includes the following main parts: expression quality control, feature extraction, dimensionality reduction analysis, unsupervised cell clustering, and transcript molecular marker expression analysis.

[0171] (W5) Perform quality control on the input transcript expression matrix and retain transcripts with at least ≥0.5% of the expressed cells. Use the built-in functions of Seurat to normalize the transcript expression and calculate the expression variation of the screened expression matrix, and use the 2000 most highly variable transcripts as feature transcripts for downstream dimensionality reduction analysis.

[0172] (W6) The characteristic transcripts obtained in the above step (W5) can be used for dimensionality reduction analysis and unsupervised cell clustering through the built-in functions of Seurat. The principal component analysis method (PCA) is used for dimensionality reduction, and the first 50 principal component results are retained; based on the 50 principal components, a shared nearest neighbor (SNN) similarity graph is constructed for the overall single cells, and cell clusters are divided at a specific resolution (0.3) in an unsupervised manner, and finally the cell classification results based on the expression of PAS annotated transcripts are obtained.

[0173] (W7) Based on the peak annotation information obtained in steps (W2) and (W3), the transcripts expressed by the given molecular marker genes are calculated and screened. The expression levels of the molecular marker transcripts of different cell classifications can be calculated and visualized using the built-in functions of Seurat to identify the cell types in the unsupervised cell classification results.

[0174] Prediction and analysis systems

[0175] The present invention also provides a prediction system for polyadenylation sites based on the method of the present invention, the system comprising:

[0176] An input unit, wherein the input unit is configured to input data, wherein the data includes an input sequence to be predicted, wherein it is necessary to determine whether a polyadenylation site exists in the input sequence to be predicted;

[0177] a polyadenylation site prediction unit, the prediction unit being configured to execute a polyadenylation site prediction model, thereby obtaining a prediction result of whether a polyadenylation site exists in the input sequence to be predicted; wherein the prediction model is constructed using the method described in the first aspect of the present invention;

[0178] An output unit is configured to output the prediction result of the polyadenylation site prediction unit.

[0179] The present invention also provides a single-cell transcriptome sequencing data analysis system (SCAPTURE system), the system comprising:

[0180] (a) an input unit, wherein the input unit is configured to input sequence data, wherein the input sequence data is selected from the following group: (i) a DeepPASS prediction model and / or an input sequence to be predicted, wherein it is necessary to determine whether a polyadenylation site exists in the input sequence to be predicted; (ii) single-cell transcriptome sequencing data; (iii) a combination of (i) and (ii);

[0181] (b) an analysis module, the analysis module being configured to perform a predetermined analysis on the input data to obtain an analysis result, and the analysis module is selected from the following group:

[0182] (M1) Prediction unit of polyadenylation site;

[0183] (M2) analysis unit of polyadenylation sequencing read signal peak;

[0184] (M3) Single-cell transcript quantification analysis and / or single-cell typing analysis unit;

[0185] (M4) a combination of (M1), (M2) and (M3) above;

[0186] (c) an output unit, wherein the output unit is configured to output the analysis result of the analysis module.

[0187] In another preferred embodiment, the output system includes a display screen, a pad, a printer, a mobile phone, etc.

[0188] The main advantages of the present invention include:

[0189] (1) According to step (2) of invention content 1, during the preprocessing of the training data of the polyadenylation site prediction model DeepPASS, a data offset and data augmentation strategy based on normal distribution random perturbation was adopted; data offset can enable DeepPASS to achieve a non-position-dependent polyadenylation site prediction effect; data augmentation can avoid overfitting of DeepPASS training and improve the training and prediction effects.

[0190] (2) According to steps (Y2)-(Y4) of invention content 2, in the peak identification of polyadenylated sequencing reads, the present invention integrates the gene and transcript information in a given annotation, and uses the sequencing read signal distribution at the transcript level for peak calculation, which can improve the sensitivity of peak identification.

[0191] (3) According to step (Y6) of invention content 2, in the screening of polyadenylated sequencing read peaks, the present invention uses normal distribution test multi-peak test and single peak test screening conditions to remove noise signals in the peak identification process.

[0192] (4) According to step (Y9) of invention content 2, the present invention adopts a strategy of ranking the coverage number of known polyadenylation sites and the results predicted by DeepPASS during the integration of multi-sample sequencing read peaks, and selects the polyadenylation site peak with the highest confidence for consistency analysis in multiple samples.

[0193] (5) According to step (Y10) of invention content 2, the present invention introduces a non-position-dependent deep learning model prediction in the screening of polyadenylation sequencing read peaks in single-cell transcriptome sequencing data, which can specifically screen out polyadenylation peaks with high confidence.

[0194] (6) According to steps (Y4)-(Y5) of content 2 and step (W2) of content 3, the present invention utilizes the read segment cross-exon alignment position in the sequencing read segment, combined with the splice site position and gene end position of the transcript annotation, to perform computational inference on the transcript from which the sequencing read segment peak is derived, thereby achieving the identification of selectively expressed transcripts.

[0195] (7) According to steps (W4) to (W7) of Invention Content 3, the present invention can quantify transcripts at the single-cell level and apply the quantification at the transcript level to the single-cell typing analysis system, replacing the traditional gene-level quantitative analysis.

[0196] The present invention will be further described below in conjunction with specific examples. It should be understood that these examples are intended to illustrate the present invention only and are not intended to limit the scope of the present invention. The experimental methods in the following examples where specific conditions are not specified are usually performed under conventional conditions, such as those described in Sambrook et al., Molecular Cloning: A Laboratory Manual (New York: Cold Spring Harbor Laboratory Press, 1989), or under conditions recommended by the manufacturer. Unless otherwise indicated, percentages and parts are weight percentages and weight parts.

[0197] Example 1

[0198] Construction of DeepPASS, a polyadenylation site prediction model

[0199] In this embodiment, data collection, preprocessing, model building and training evaluation of the polyadenylation site prediction model DeepPASS are included.

[0200] Materials: The published polyadenylation site database was collected, and a total of 290,051 polyadenylation site annotations were downloaded; a total of 514,262 polyadenylation site annotations were downloaded from PolyA-seq (accession number GSE30198); a total of 569,005 polyadenylation site annotations were downloaded; a total of 49,942 polyadenylation site annotations were downloaded from GENCODE polyadenylation site annotations; human genome and gene annotation files were downloaded, with a total of 228,048 gene annotations; DeepPASTA prediction tool and APARENT prediction tool were downloaded and installed from the GitHub distributed version control software website.

[0201] Step 01: Download and install the open source cross-package management and environment management system Conda bedtools annotation file analysis tool (version number: 2.28.0), programming language Python application (version number: 3.7), TensorFlow deep learning framework (version number: 2.0.0), programming language R application (version number: 4.0.0)

[0202] Step 02: Convert the polyadenylation site database files from PolyA_DB3, PolyA-seq and Poly(A)Site2.0 into BED format files, and extend the distance 6bp upstream and downstream in the genomic coordinate interval

[0203] Step 03: Use bedtools intersect to calculate whether there are any genomic interval intersections between the expanded polyadenylation sites, and mark the polyadenylation sites with intersections with the same number.

[0204] Step 04: Count the number of databases with polyadenylation sites with the same number, retain the high-confidence sites shared by the three databases, and merge them with the polyadenylation sites of GENCODE as a high-confidence polyadenylation site library, totaling 251,072 entries.

[0205] Step 05: Use the rnorm function (μ=0, σ=10) of the R language to generate normally distributed data points, add the genomic coordinate values ​​of the high-confidence polyadenylation sites to the normally distributed data points, and generate the polyadenylation sites after normal distribution offset amplification.

[0206] Step 06: Use bedtools getfasta to obtain the 200bp genomic sequence from 100bp upstream to 100bp downstream of the polyadenylation site after normal distribution offset amplification, and use these sequences as the positive set for model training.

[0207] Step 07: Based on gene annotation, calculate the intergenic regions in the genome without gene annotation sites, and use PolyA_DB3, PolyA-seq, Poly(A)Site2.0 and all polyadenylation sites of GENCODE for reverse screening to obtain intergenic regions that contain neither gene annotations nor polyadenylation site annotations. Use bedtools getfasta to obtain 200bp of genomic sequences in random regions within these intergenic regions, and use these sequences as negative sets for model training.

[0208] Step 08: Use TensorFlow to build the DeepPASS deep learning model framework. The basic structure is: (1) The convolutional neural network layer includes two convolutional layers and a max pooling layer. The first convolutional layer includes 12-mer width, 4 channels, and a total of 128 filters; the second convolutional layer includes 6-mer, 64 filters can cover the 128 channels output by the first convolutional layer. The output of the convolutional layer is compressed by the max pooling layer, and the maximum value in the 4×4 neighborhood is selected and output to the downstream RNN layer. (2) The recurrent neural network layer adopts a bidirectional long short-term memory network structure (BiLSTM), sets 128 units to process the output content of the convolutional neural network module, and passes the result to the downstream fully connected layer. The fully connected layer sets 1024 hidden units and a dropout ratio of 30%, and outputs the final binary classification probability through a normalized exponential function (softmax activation function). The overall structure can be summarized as follows:

[0209]

[0210] Step 09: Use the positive and negative set sequences obtained in steps 06 and 07 to train DeepPASS. The DeepPASS model training sets the data batch size to 5000, and the model is trained for 100 periods. If the training results do not improve within 10 periods, the model training is terminated early.

[0211] Step 10: Compare DeepPASS with similar polyadenylation site prediction tools (DeepPASTA and APARENT) published so far. After merging the positive and negative sequences obtained in steps 06 and 07, 50,000 sequences were randomly sampled as the test set for model evaluation of DeepPASS, DeepPASTA, and APARENT.

[0212] Step 11: Use the test set to input DeepPASS for prediction, calculate the original value and predicted probability of the test set sequence, and calculate the true positive rate and false positive rate under different probability thresholds.

[0213] Step 12: Format the test set file according to the DeepPASTA software requirements, input DeepPASTA for prediction, calculate the difference between the maximum and minimum predicted probabilities, take different probability thresholds in a continuous and incremental manner according to the difference interval, and calculate the true positive rate and false positive rate under each threshold.

[0214] Step 13: Format the test set file according to the requirements of the APATENT software, input APATENT for prediction, select the maximum value of the return value of a single sequence as the prediction result of the test sequence, calculate the difference between the maximum and minimum prediction scores, and take different score thresholds in a continuously increasing manner according to the difference interval, and calculate the true positive rate and false positive rate under each threshold.

[0215] Step 14: The prediction performance of DeepPASS, DeepPASTA, and APATENT in steps 11-13 was evaluated using the receiver operating characteristic curve (ROC) and area under the curve (AUC).

[0216] Step 15: Use the polyadenylation site of GENCODE to shift the genomic coordinates, shifting 20bp upstream, 15bp upstream, 10bp upstream, 5bp upstream, 0bp, 5bp downstream, 10bp downstream, 15bp downstream, and 20bp downstream, respectively, to construct nine different data sets.

[0217] Step 16: Test DeepPASS and DeepPASTA using the nine datasets constructed in Step 15, and calculate the accuracy and stability of DeepPASS and DeepPASTA in predicting polyadenylation sites under different position offset conditions.

[0218] result

[0219] like Figure 1 As shown in A and 1B, a total of 251,072 high-confidence sites were screened from four polyadenylation site databases. After normal distribution random shift and amplification, 1,907,294 positive sequences were extracted. A similar number of negative sequences were extracted from the intergenic regions. A training set of about four million sequences was constructed for DeepPASS model training ( Figure 1 C).

[0220] The DeepPASS model training results are as follows Figure 1D shows that the prediction accuracy of the test set reaches 98%. Comparing the overall prediction performance of DeepPASS, DeepPASTA and APARENT, the area under the curve of DeepPASS prediction is calculated to be 0.99, while that of DeepPASTA and APARENT is 0.92 and 0.72 respectively.

[0221] The results show that DeepPASS performs better than the existing DeepPASTA and APARENT in polyadenylation site prediction.

[0222] In order to verify the non-position dependence of the DeepPASS model, the present invention also constructed nine data sets based on the GENCODE database, with offsets of 20bp upstream, 15bp upstream, 10bp upstream, 5bp upstream, 0bp, 5bp downstream, 10bp downstream, 15bp downstream, and 20bp downstream, and conducted comparative tests with DeepPASTA.

[0223] The results are as follows Figure 1 As shown in Figures E and 1F, the performance accuracy of DeepPASS on the nine datasets are 92.6%, 95.5%, 96.9%, 97.5%, 97.6%, 97.2%, 96.5%, 94.8%, and 91.8%, respectively, and the performance accuracy of DeepPASTA on the five datasets are 21.3%, 29.2%, 43.4%, 71.4%, 81.1%, 59.1%, 29.0%, 20.4%, and 12.2%, respectively. This shows that DeepPASS has the advantages of significant position independence and higher accuracy in predicting polyadenylation sites.

[0224] The above results show that the deep learning model DeepPASS constructed by the present invention has the advantages of being position-independent and highly accurate in the prediction of polyadenylation sites. The DeepPASS model can help screen out truly reliable sites from the polyadenylation sites identified by single-cell transcriptome sequencing data, reduce false positives introduced by sequencing library construction technology and sequencing data processing, and promote the development of polyadenylation biology research at the single-cell level.

[0225] Example 2

[0226] Identification and screening of polyadenylation site sequencing read peaks in single-cell transcriptome data

[0227] Materials: Six datasets of human PBMC single-cell transcriptome sequencing were downloaded from the 10x Genomics website, and named 3k, 4k, 5k, 6k, 8k, and 10k according to the number of cells detected in each dataset. The sequencing read numbers of these six datasets were 179147840, 379462522, 383941607, 143988649, 784064148, and 638901019, respectively; human genome and gene annotation files were downloaded from GENCODE, with a total of 228048 gene annotations;

[0228] Step 01: Download and install Cell Ranger single cell data preprocessing software from the 10x Genomics website

[0229] (https: / / support.10xGenomics.com / single-cell-gene-expression / software / pipelines / latest / installation) (version number: 4.0); download the annotation format processing suite UCSCUtilities; download the HOMER peak recognition tool (version number: 4.11.1); download the BAM file processing tool samtools from GitHub (version number: 1.9); download the single-cell sequencing tool Sierra from GitHub (version number: 0.99.24); download the single-cell sequencing tool scAPA from GitHub (version number: 0.1.0);

[0230] Step 02: Use Cell Ranger, a single-cell transcriptome data analysis tool released by 10x Genomics, to create index files for the GENCODE genome and gene annotation files.

[0231] Step 03: Use Cell Ranger to align the 6 PBMC single-cell transcript sequencing data to the index file to obtain the sequencing read segment alignment result BAM file and cell × gene expression expression matrix file for each sample.

[0232] Step 04: Convert the gene-level annotation file from GENCODE into a gene-transcript-labeled gene-transcript-level annotation file based on its gene and transcript information, and use the UCSC Utilities tool to construct genepred, GTF, and BED format files respectively.

[0233] Step 05: Extract the genomic coordinate range and chain direction information of each gene according to the gene annotation file, and use bedtools and samtools to extract the sequencing read segment alignment results within the corresponding genome range of each gene from the BAM file.

[0234] Step 06: Based on the constructed transcript-level annotation file and the extracted sequencing read alignment results within each gene, for each transcript of the same gene, use makeTagDirectory and findPeaks in the HOMER tool to identify the peak region where the transcript sequencing reads are enriched, i.e., the sequencing read peak.

[0235] Step 07: For the peak interval calculated based on the transcript horizontal coordinates, the transcript position system is mapped to the genome coordinate position system, and the peak genome coordinate system annotation file corresponding to the transcript identification is calculated.

[0236] Step 08: For the peak annotation after the coordinate conversion, identify the exon splicing sites in the transcript annotation according to its gene and transcript level annotation information, and calculate the peak splicing site information of the corresponding transcript.

[0237] Step 09: Execute steps 05-08 for each expressed gene in the sequencing sample to calculate the peak annotation of the sequencing reads in the whole genome of the sample.

[0238] Step 10: Perform preliminary screening on the obtained peaks to remove noise and low-abundance peaks. The screening conditions include: using nortest in R to perform a normality test on the sequencing read distribution curve of the peaks, and removing peaks that do not conform to the normal distribution; using diptest in R to perform a unimodal test on the sequencing read distribution curve of the peaks, and removing peaks that do not conform to the multimodal distribution and the normal distribution.

[0239] Step 11: For the sequencing read distribution peaks generated by the same polyadenylation site, multiple transcripts covering the polyadenylation site in the gene annotation will calculate multiple peaks. The method for de-redundanting these peaks is as follows: by calculating the peak interval overlap ratio, the peaks with an overlap ratio greater than 50% are classified and grouped. For all overlapping peaks in the same group, only the peak annotation with the highest sequencing read coverage depth is retained.

[0240] Step 12: Recalculate the abundance of sequencing reads for the peaks after deduplication within the same gene, and remove the peaks whose read count is less than 1% of the total read count on the same gene.

[0241] Step 13: For all peaks found in the sequencing sample, perform steps 10-12 to obtain peaks that are preliminarily identified and screened at the whole genome level for evaluation of polyadenylation sites.

[0242] Step 14: For the existing tool Sierra, use the same six PBMC single-cell sequencing data sets and gene annotations for analysis, and obtain the peak annotation file of the results according to the instructions for use of Sierra.

[0243] Step 15: For the existing tool scAPA, use the same six PBMC single-cell sequencing datasets and scAPA built-in gene annotations for analysis, and obtain the peak annotation files of the results according to the instructions for use of scAPA.

[0244] Step 16: Count the peak annotations obtained from SCAPTURE, Sierra, and scAPA analyses of the six PBMC single-cell sequencing datasets, count the number of peaks in the intronic and exonic regions identified in each sample, and compare between the tools.

[0245] Step 17: The DeepPASS deep learning model constructed by the present invention is used for sequencing read peak screening. The specific method is: the end position of the sequencing read peak annotation interval is used as a potential polyadenylation site, and the genome sequence from 100bp upstream to 100bp downstream is input into the DeepPASS model for prediction. It can be determined whether each sequencing read peak is a high-confidence polyadenylation peak, and the high-confidence polyadenylation peak is screened out.

[0246] result

[0247] Six single-cell sequencing data sets of human PBMC from 10x Genomics have good alignment results, with genome alignment efficiencies of 92.2%, 89.8%, 90.2%, 92.0%, 88.8% and 89.0%, respectively. Based on the genome alignment results, the SCAPTURE system established by the present invention can use single-cell sequencing data sets to identify and screen polyadenylated sequencing read peaks (such as Figure 2 A).

[0248] See also Figure 2 A, taking NADK gene as an example, shows the steps of identifying and screening the peak of polyadenylation sequencing reads in the single-cell transcriptome sequencing data set of human PBMC 10k by the SCAPTURE computing system. Figure 2SCAPTURE first extracts the sequencing read alignment results within the NADK gene range, and identifies the sequencing read peaks for each transcript based on the constructed transcript-level annotation information. Then, the peaks from different transcript sources within NADK are de-redundant and expression quality controlled, and the polyadenylation sequencing read peak information of NADK is identified, and the results are Peak1 and Peak2 (see Figure 2 A).

[0249] like Figure 2 As shown in Figure B, the present invention compares SCAPTURE with the existing Sierra and scAPA processes. For the peak sites of sequencing reads in the exon regions of six PBMC single-cell sequencing samples, SCAPTURE can identify 32569, 45135, 50970, 33578, 48906 and 54236 peak sites respectively; Sierra's identification results are 11684, 18635, 19956, 10025, 24106 and 24590 peak sites respectively; scAPA's identification results are 11684, 18635, 19956, 10025, 24106 and 24590 peak sites respectively. For the sequencing read peak sites in the intron region and the exon region, SCAPTURE was able to identify 54303, 75816, 92108, 57772, 82056 and 98328 peak sites respectively; Sierra identified 10392, 36095, 49495, 8659, 59925 and 86912 peak sites respectively; scAPA identified 7067, 16947, 19747, 5631, 18987 and 25423 peak sites respectively. This result shows that for the same single-cell sequencing sample, SCAPTURE can identify more sequencing read peak sites in the exon region and intron, and the sensitivity of identifying sequencing read peaks is significantly better than existing tools.

[0250] like Figure 2 As shown in Figure C, the present invention uses the constructed DeepPASS model to screen high-confidence sequencing read peaks to reduce false positives in polyadenylation site identification. After screening by the DeepPASS model, the enrichment multiples of the proportion of classic polyadenylation tailing signals in the exon peak and intron peak are: exon peak enrichment 6.2 to 7.0 times, intron peak enrichment 2.9 to 7.4 times.

[0251] The results show that the DeepPASS model constructed in the present invention can be used for polyadenylation peak confidence screening in single-cell transcriptome sequencing, thereby improving the accuracy of polyadenylation research at the single-cell level.

[0252] like Figure 2 As shown in C, after the constructed DeepPASS model was used to screen and remove false positive polyadenylation sites, the number of exon and intron peaks retained by DeepPASS constructed in the present invention was 25247 and 12286, which were significantly more than Sierra (14039 and 9830) and scAPA (8378 and 2429).

[0253] Based on the above results, SCAPTURE has significantly higher sensitivity and accuracy compared with existing single-cell transcriptome polyadenylation research tools.

[0254] Example 3

[0255] Identify transcript expression based on polyadenylation sites and apply it to single-cell phenotyping analysis.

[0256] Materials: Six data sets of human PBMC single-cell transcriptome sequencing downloaded from the 10x Genomics company website used in Example 2; gene annotations downloaded from GENCODE used in Example 2.

[0257] Step 01: Download and install Seurat single-cell analysis tool (version number: 3.2.2) from GitHub; download featureCounts sequencing read counting tool (version number: 1.6.2) from SourceForge; download molecular identification counting software UMI-tools (version number: 1.0.1) from GitHub.

[0258] Step 02: According to the peak annotation information and GENCODE transcript annotation information obtained in Example 2, the specific transcript from which the peak comes is identified. The identification method is as follows: the transcript exon splicing site must completely match the peak splicing site information, and the transcript with the peak closest to the annotation end is selected as the transcript from which the peak comes, and finally the transcript selectively expressed by the gene in the single-cell transcriptome sequencing sample is obtained.

[0259] Step 03: Convert the peak annotation format after mapping the transcript information in step 02, and use the UCSCUtilities tool to construct a GTF format mapped transcript information peak annotation file for downstream single-cell transcript quantitative analysis.

[0260] Step 04: Use Cell Ranger to align the 6 PBMC single-cell transcript sequencing data to the index file to obtain the sequencing read segment alignment result BAM file and cell × gene expression expression matrix file for each sample.

[0261] Step 05: Input the cell × gene expression matrix of the six samples into Seurat for cell quality control. The cell screening criteria are: the proportion of unique analysis identifiers of the mitochondrial genome ≤ 20%, the number of unique analysis identifiers of the cells ≥ 800, and the number of genes detected in the cells ≥ 500. The cell barcodes of the six samples after screening are retained.

[0262] Step 06: Map the peak annotation file of the transcript information, and reallocate the sequencing reads in the BAM files of the six samples after alignment through featureCounts to obtain a re-assigned BAM file that marks the peak annotation of the transcript information.

[0263] Step 07: UMI-tools can be used to quantify the transcript expression of the cell barcodes after screening in the re-assigned BAM file, and calculate the cell × transcript expression matrix. The specific method is: the unique molecular identifier information on the sequencing reads assigned to each transcript peak annotation is counted, and the number of unique molecular identifiers of each transcript peak in each cell is calculated. This value can be used as the expression level of the transcript in the cell.

[0264] Step 08: Use the gene-level expression matrix and transcript-level expression matrix obtained in step 05 and step 08 respectively to perform single-cell transcriptome analysis, including quality control, standardization, dimensionality reduction analysis, unsupervised classification and cell type identification.

[0265] Step 09: Perform quality control on the single-cell expression data, remove transcripts with transcript expression in ≤ 0.5% of cells, and standardize gene expression and transcript expression.

[0266] Step 10: Calculate the inter-cell expression variability of genes and transcripts, and select the top 2000 highly variable features for downstream dimensionality reduction analysis.

[0267] Step 11: Use Seurat's built-in function RunPCA to perform principal component analysis on the highly variable features of genes and transcripts to reduce their dimensionality, and retain the first 50 principal component results.

[0268] Step 12: Use the built-in function FindNeighbors of Seurat to construct a shared nearest neighbor similarity graph of the overall single cell using the 50 principal components obtained from gene and transcript analysis.

[0269] Step 13: Use the built-in function FindClusters of Seurat to divide the gene and transcript analysis into cell clusters with a resolution of 0.3 in an unsupervised manner, and obtain the unsupervised classification results of cells.

[0270] Step 14: Find the marker genes of different cell types in human PBMC reported in the literature, perform marker gene expression analysis on the gene and transcript expression levels to identify the cell clusters, and annotate the cell clusters as corresponding cell types based on the marker gene expression of different cell types.

[0271] result:

[0272] like Figure 3 As shown in A and 3B, the calculation process of the present invention can identify the selective transcript expression of the gene based on the peak information of the polyadenylation sequencing read segment. Taking RASGRP2 as an example, the gene expression transcript 1 (ENST00000377486.7) and transcript 2 (ENST00000354024.7) can be identified by Peak1 and Peak2, and transcript 1 and transcript 2 can characterize the gene expression amount at the single cell level. This result shows that the present invention can identify and quantify transcripts based on the polyadenylation sequencing read segment peak.

[0273] like Figure 3 As shown in C, for the polyadenylated sequencing read peak after screening, its selectively expressed transcript can be identified based on its splicing site information. In six PBMC single-cell sequencing samples, the number of expressed genes obtained by the traditional method was 13,626, and the number of expressed transcripts obtained by the SCAPTURE system constructed by the present invention was 34,316, which expanded the single-cell analysis from the gene level to the transcript level.

[0274] like Figure 3 As shown in D, 3E, 3F, 3G, and 3H, in the case of single-cell sample integration and unsupervised clustering analysis, the analysis power based on transcript expression is close to that of gene expression. Figure 3 G and 3H, some cells of the unclassified cell populations identified at the gene level can be reclassified into known cell populations based on transcript expression, such as reclassified dendritic cells and reclassified natural killer cells. This result shows that the transcript expression analysis established by the present invention can replace gene expression for conventional single-cell computational analysis; in addition, the advantage of transcript expression is that it can more accurately distinguish different cell types in a single-cell population.

[0275] The results show that the SCAPTURE computing system established by the present invention can identify selectively expressed transcripts based on polyadenylation sites, and perform transcript quantification and single-cell typing analysis at the single-cell level.

[0276] discuss

[0277] The present invention provides a polyadenylation site (PAS) detection method integrating a deep learning model, which is used to detect polyadenylation sites and differential expression of transcripts in single-cell sequencing data, and is applied to single-cell typing analysis (i.e., the SCAPTURE system).

[0278] The significance of the SCAPTURE computing system established by the present invention lies in that the constructed DeepPASS model realizes high-accuracy and position-independent prediction of polyadenylation sites, and the prediction tool only needs the genome sequence as input, so it can also be widely used in other research systems and biological problems related to polyadenylation; SCAPTURE can highly sensitively identify single-cell polyadenylation sites at the whole genome level from scratch, overcoming the shortcomings of low-abundance signal loss and high false positives in single-cell transcriptome sequencing; SCAPTURE can also identify different transcripts of selective gene expression based on polyadenylation sites, improve single-cell analysis from the gene level to the transcript level, and has advantages in single-cell analysis; SCAPTURE can also be used for the identification, screening and quantification of polyadenylation sites in single-cell transcriptome sequencing, and can also be applied to polyadenylation research in systems such as tissue specificity, cell development and disease occurrence, so that polyadenylation research is improved to the single-cell level.

[0279] In the present invention, the known polyadenylation sites that can be used are collected from PolyA_DB3, PolyA-seq, Poly(A)Site2.0 and GENCODE. These polyadenylation sites can also come from other databases or data sets.

[0280] When constructing the gene-transcript level annotation file, human GENCODE annotation information is downloaded. The annotation information can also be gene annotations from other research species, database sources, or manually established.

[0281] In the present invention, the 3'tag-based single-cell transcriptome sequencing can obtain sequencing data using the 10x Chromium technology of 10x Genomics, and is also applicable to other similar sequencing data based on the 3'tag-based single-cell transcriptome sequencing principle.

[0282] In the present invention, sequencing data can use Cell Ranger or other suitable alignment processes to perform genome alignment positions of sequencing reads of single-cell transcriptome data, labeling of cell barcodes and unique molecular identifiers, and the obtained results are stored in a BAM file.

[0283] In the present invention, similar peak identification tools based on high-throughput sequencing read distribution can also be used when identifying sequencing read peaks formed by polyadenylation.

[0284] In the present invention, when calculating the transcript expression based on peak annotation, other similar tools that can process cell barcodes and unique molecular identifiers of sequencing reads can also be used.

[0285] It should be understood that the above embodiments of the present invention are described based on single-cell sequencing data of human PBMC. Obviously, the present invention can also be used for PAS identification and transcript quantification of 3'tag-based single-cell transcriptome sequencing data of other species.

[0286] All documents mentioned in the present invention are cited as references in this application, just as each document is cited as reference individually. In addition, it should be understood that after reading the above teachings of the present invention, those skilled in the art can make various changes or modifications to the present invention, and these equivalent forms also fall within the scope defined by the claims attached to this application.

Claims

1. A method for constructing a prediction model for polyadenylation sites, characterized in that: The method comprises the steps of: (S1) providing a data set for model training, the data set comprising a positive set and a negative set, wherein the sequences contained in the positive set are sequences containing polyadenylation sites, i.e., positive sequences; and the sequences contained in the negative set are sequences not containing polyadenylation sites, i.e., negative sequences; Wherein, the positive sequence includes: a basic positive sequence for the original polyadenylation site and a shifted positive sequence obtained by shift amplification based on the basic positive sequence, and both the basic positive sequence and the shifted positive sequence contain the original polyadenylation site; (S2) using the dataset described in (S1) for training the model architecture; Wherein, the model architecture includes a model architecture based on a convolutional neural network (CNN) module and a recurrent neural network (RNN) module; The CNN module includes two convolutional layers and a max pooling layer, wherein the output result of the convolutional layer is compressed by dimensionality reduction through the max pooling layer, and the maximum value in a 4×4 neighborhood is selected and output to the downstream RNN module; The RNN module is provided with a bidirectional long short-term memory network structure (BiLSTM), and transmits the output content of the CNN module to the downstream fully connected layer; The fully connected layer is provided with a plurality of hidden units and a predetermined dropout ratio, and outputs the final binary classification result through a normalized exponential function (softmax activation function); (S3) When the training result reaches a predetermined termination condition, the model training is terminated, and the prediction model of the polyadenylation site, namely the DeepPASS model, is obtained.

2. The method according to claim 1, characterized in that The offset amplification refers to: placing a basic positive sequence on a coordinate axis x1, and taking the polyadenylated original genomic site of the basic positive sequence as the coordinate zero point, and then making the 5' end and 3' end of the basic sequence independently perform a positive or negative offset of p bp, and the offset sequence is called the offset positive sequence, where p is a positive integer (bp).

3. The method according to claim 1, characterized in that In step (S1), the offset amplification includes: taking the original polyadenylation genomic site as the coordinate zero point, randomly offsetting the site using a normal distribution, offsetting a polyadenylation site 4-20 times, and extracting the genomic sequence from 100±20 bp upstream to 100±20 bp downstream of the site after offset amplification as the offset positive sequence in the positive set of model training.

4. The method according to claim 1, wherein: The offset amplification includes: taking the original polyadenylation site as the coordinate zero point, randomly offsetting the site using a normal distribution, performing about 10 offsets on each polyadenylation site, and extracting the genomic sequence from 100 bp upstream to 100 bp downstream of the amplified site as the offset positive sequence in the positive set of DeepPASS model training.

5. The method according to claim 1, characterized in that The offset amplification includes: random or non-random offset amplification.

6. The method according to claim 1, characterized in that The length of the basic positive sequence and the length of the corresponding offset positive sequence are the same or substantially the same.

7. The method according to claim 1, characterized in that When the coordinates of the basic positive sequence are (X Left , X Right ) and the lengths of the base positive sequence and the corresponding offset positive sequence are the same, then after a positive offset of p bps, the coordinates of the offset positive sequence are (X Left +p,X Right +p); and after a negative shift of p bp, the coordinates of the shifted positive sequence are (X Left -p,X Right -p).

8. The method according to claim 1, characterized in that In the positive set, the ratio of the total number of basic positive sequences to the total number of offset positive sequences is 1:2 to 1:

20.

9. The method according to claim 1, characterized in that In the positive set, for each polyadenylation original site, it corresponds to 1 basic positive sequence and D offset positive sequences, wherein D is a positive integer ≥2.

10. The method according to claim 1, characterized in that In the positive set, the total number T1 of the basic positive sequences is the same or substantially the same as the total number T0 of the polyadenylation original sites.

11. The method according to claim 1, characterized in that The sequences contained in the positive set and the negative set are from the same species or different species.

12. The method according to claim 1, characterized in that Four polyadenylation site databases from different sources were collected: PolyA_DB3, PolyA-seq, Poly(A)Site2.0 and GENCODE. The common sites in the three published databases PolyA_DB3, PolyA-seq and Poly(A)Site2.0 were screened and merged with GENCODE as a library of known polyadenylation sites for DeepPASS model training.

13. The method according to claim 1, wherein: Four polyadenylation site databases from different sources were collected: PolyA_DB3, PolyA-seq, Poly(A)Site2.0 and GENCODE. High-confidence polyadenylation sites in the three published databases, PolyA_DB3, PolyA-seq and Poly(A)Site2.0, were screened. The screening criteria were: the polyadenylation sites in one database overlapped with those in the other two databases within 12 bp upstream and downstream of its genomic position. The high-confidence sites in these three databases were merged with GENCODE as a known polyadenylation site library for polyadenylation site annotation.

14. The method according to claim 1, wherein: The negative sequences in the negative set are screened in the following manner: by taking the union of four polyadenylation site databases, the intergenic regions of the genome are reversely screened, and only the intergenic region sequences outside the upstream and downstream distance of 100 bp of the known polyadenylation sites are selected.

15. A prediction system for polyadenylation sites, characterized in that: The system comprises: An input unit, wherein the input unit is configured to input data, wherein the data includes an input sequence to be predicted, wherein it is necessary to determine whether a polyadenylation site exists in the input sequence to be predicted; a polyadenylation site prediction unit, the prediction unit being configured to execute a polyadenylation site prediction model, thereby obtaining a prediction result of whether a polyadenylation site exists in the input sequence to be predicted; wherein the prediction model is constructed using the method of claim 1; An output unit is configured to output the prediction result of the polyadenylation site prediction unit.

16. An analysis method for identifying polyadenylation sequencing read signal peaks based on single-cell transcriptome sequencing data, characterized in that: Includes steps: (Y1) Align the single-cell transcriptome sequencing data with the genome sequence to obtain a sequencing read alignment result file stored in the BAM format. The BAM file should contain: sequencing read alignment result information required by the BAM format, sequencing read cell barcode (CB) and sequencing read molecular identifier (UMI) information; (Y2) converting the gene-level annotation file into a gene-transcript-level annotation file labeled "gene-transcript" according to its gene and transcript information; (Y3) Based on the sequencing read alignment result file and gene transcript level annotation file, HOMER was used to findPeaks Identify sequencing read distribution peaks at the transcript level, identify all peaks corresponding to each transcript, and mark the interval positions of the peaks based on the coordinate information of the corresponding transcripts; (Y4) combining the gene transcript level annotation file obtained in the above steps (Y2)-(Y3) and the peak interval position based on the transcript coordinates, mapping the transcript position system to the genome coordinate position system, and calculating the genome coordinate system annotation file of the peak; (Y5) combining the gene transcript level annotation file obtained in the above steps (Y2)-(Y4) and the peak file based on the genome coordinate system, further calculating the peak splicing site information of the corresponding transcript by identifying the transcript exon splicing site in the gene transcript level annotation file; (Y6) Based on the peak annotation results obtained in the above step (Y5), a normality test is performed on the sequencing read signal distribution curve of the identified peak to remove peaks that may be caused by abnormal amplification, multiple alignments or low signal abundance; and a unimodality test is performed to retain the overlapping peaks formed by adjacent polyadenylation signals; (Y7) The above steps (Y2)-(Y6) can calculate peaks from different transcript sources of the same gene, and the peaks can be deduplicated according to their gene information; the deduplication method is: for peaks with interval overlap ratio exceeding 50%, the abundance of their sequencing reads is sorted, and the peak with the highest expression level is retained; (Y8) Recalculating the abundance of sequencing reads for the peaks after deduplication in the above step (Y7), filtering out the peaks whose sequencing reads account for ≤1% in the same gene, and obtaining the preliminary screened peaks (raw peaks); (Y9) for each raw peak obtained from steps (Y3) to (Y8) of the plurality of samples, the sample information of the peak source is marked, and the peaks with an overlap ratio of more than 60% between different samples are marked as the same group, and a unique value is selected from the peaks in the same group as the peak of the polyadenylation site; (Y10) The end of the raw peak annotation interval identified in the above step is approximated as the polyadenylation site, and the sequence from 100 bp upstream to 100 bp downstream of the site is extracted and input into the DeepPASS model for prediction, and it is determined whether each raw peak is a high-confidence polyadenylation peak (PAS), thereby obtaining the polyadenylation sequencing read signal peak analysis results of the single-cell transcriptome sequencing data; Wherein, in step (Y10), the prediction model is a prediction model constructed using the method described in claim 1, namely, the DeepPASS model.

17. The analysis method according to claim 16, characterized in that The DeepPASS model was used for prediction to determine whether each raw peak was a high-confidence polyadenylation peak (PAS), and the selected high-confidence peaks were used for downstream analysis.

18. A method for quantitative analysis of transcripts and / or single cell typing analysis at the single cell level, characterized in that: The method comprises: (Z1) Provide single-cell transcriptome sequencing data of single-cell samples to be analyzed; (Z2) analyzing the single-cell transcriptome sequencing data using the analysis method of claim 16 to obtain a high-confidence peak; (Z3) Based on the high-confidence peaks obtained in the previous step, perform transcript quantitative analysis and / or single-cell typing analysis.

19. The method according to claim 18, characterized in that In step (Z3), comprising: (W1) For the peak integration analysis of multiple samples, the integration method is: mark the overlapping peaks between different samples as the same group, sort the peaks in the same group by confidence, and select the peak with the highest confidence as the only peak of the polyadenylation site; (W2) According to the peaks selected in the above step (W1), the peaks are mapped to the transcripts of the corresponding sources through the splicing sites and genomic interval information of the peaks, so as to identify the transcripts generated by alternative polyadenylation (APA) in the single-cell transcriptome sequencing samples and realize the identification of transcript expression; (W3) Construct a GTF format annotation file with the PAS annotations of the mapped transcript information, and recalculate the distribution of PAS annotations in the aligned BAM reads through featureCounts to obtain a re-assigned BAM file marked with PAS annotations; UMI-tools can be used to count the unique molecular identifier (UMI) information of all PAS annotations corresponding to each cell barcode in the re-assigned BAM file to obtain a transcript expression matrix at the single-cell level. The matrix includes the number of unique molecular identifiers of each transcript in each cell, which can be used for single-cell transcriptome analysis; (W4) For single-cell analysis at the transcript level, the peak-annotated transcript level expression matrix obtained in the above step (W3) is used as the input of the single-cell analysis tool Seurat for single-cell analysis; the single-cell analysis system based on the peak-annotated transcript level includes the following parts: expression quality control, feature extraction, dimensionality reduction analysis, unsupervised cell clustering, and transcript molecular marker expression analysis; (W5) Perform quality control on the input transcript expression matrix and retain transcripts with expression cell numbers ≥ 0.5%; use the built-in functions of Seurat to standardize transcript expression and calculate expression variation on the screened expression matrix, and use the 2000 most highly variable transcripts as feature transcripts for downstream dimensionality reduction analysis; (W6) The characteristic transcripts obtained in the above step (W5) can be subjected to dimensionality reduction analysis and unsupervised cell clustering through the built-in functions of Seurat; the principal component analysis method (PCA) is used for dimensionality reduction, and the first 50 principal component results are retained; based on the 50 principal components, a shared nearest neighbor (SNN) similarity graph is constructed for the overall single cells, and cell clusters are divided at a specific resolution in an unsupervised manner, and finally the cell classification results based on the expression of PAS annotated transcripts are obtained; (W7) Based on the peak annotation information obtained in steps (W2) and (W3), the transcripts expressed by the given molecular marker genes are calculated and screened. The expression levels of the molecular marker transcripts of different cell classifications can be calculated and visualized using the built-in functions of Seurat to identify the cell types in the unsupervised cell classification results.

20. The method according to claim 19, wherein In step (W1), for the peak integration analysis of multiple samples, the integration method is as follows: (1) For the peaks of different samples originating from the same polyadenylation site, the number of known polyadenylation sites covered by the interval from 50 bp upstream to 25 bp downstream of the end of each peak in the interval is calculated, and the peaks are sorted in descending order according to the number of covered sites; (2) The interval sequence from 100 bp upstream to 100 bp downstream of the end of the peak is extracted, input into the DeepPASS model for prediction, and sorted from 1 to 0 according to the prediction results; (3) In the same group of peaks, the peak with the best comprehensive ranking is selected as the only peak of the polyadenylation site in the multi-sample analysis.

21. The method of claim 19, wherein: In step (W2), the peak splicing site and genome annotation information are used to identify the transcripts that match the peak source; the identification method is: the transcript exon splicing site must completely match the peak splicing site information, and the transcript with the closest end annotation distance is selected as the peak source transcript, and finally the expression of different transcripts in the single-cell transcriptome sequencing sample at the whole genome level is obtained.

22. The method of claim 19, wherein: In step (W3), the PAS annotation information and transcript source information are used to construct a PAS transcription annotation file in the Gene transfer format (GTF) format, and the sequencing reads of the BAM file are redistributed in the PAS transcription annotation file in the GTF format, marking the PAS annotation to which each sequencing read belongs.

23. The method of claim 19, wherein: In steps (W4)-(W6), the transcript level expression matrix annotated by the peak in step (W3) is used instead of the traditional gene level expression matrix to carry out cell typing analysis of single-cell transcriptome data based on transcript quantitative data; the analysis process includes expression quality control, feature extraction, dimensionality reduction analysis and unsupervised cell clustering and transcript molecular marker expression analysis based on transcript quantitative data.

24. A single-cell transcriptome sequencing data analysis system, characterized in that: The system comprises: (a) an input unit, wherein the input unit is configured to input sequence data, wherein the input sequence data is selected from the following group: (i) a DeepPASS prediction model and / or an input sequence to be predicted, wherein it is necessary to determine whether a polyadenylation site exists in the input sequence to be predicted; (ii) single-cell transcriptome sequencing data; (iii) a combination of (i) and (ii); (b) an analysis module, wherein the analysis module is configured to perform a predetermined analysis on the input data to obtain an analysis result, and the analysis module is selected from the following group: (M1) a polyadenylation site prediction unit, wherein the (M1) prediction unit is configured to execute a polyadenylation site prediction model to obtain a prediction result of whether a polyadenylation site exists in the input sequence to be predicted; wherein the prediction model is constructed using the method of claim 1; (M2) a polyadenylation sequencing read signal peak analysis unit, wherein the polyadenylation sequencing read signal peak analysis unit (M1) is configured to perform the analysis method of claim 16, thereby obtaining a polyadenylation sequencing read signal peak analysis result of the single-cell transcriptome sequencing data; (M3) a single-cell level transcript quantitative analysis and / or single-cell typing analysis unit, wherein the (M3) single-cell level transcript quantitative analysis and / or single-cell typing analysis unit is configured to perform the analysis method of claim 18, thereby obtaining single-cell level transcript quantitative analysis results and / or single-cell typing analysis results; (M4) a combination of (M1), (M2) and (M3) above; (c) an output unit, wherein the output unit is configured to output the analysis result of the analysis module.

Images