WO2022188785A1 - 融合深度学习模型的单细胞转录组计算分析方法和系统 - Google Patents
融合深度学习模型的单细胞转录组计算分析方法和系统 Download PDFInfo
- Publication number
- WO2022188785A1 WO2022188785A1 PCT/CN2022/079788 CN2022079788W WO2022188785A1 WO 2022188785 A1 WO2022188785 A1 WO 2022188785A1 CN 2022079788 W CN2022079788 W CN 2022079788W WO 2022188785 A1 WO2022188785 A1 WO 2022188785A1
- Authority
- WO
- WIPO (PCT)
- Prior art keywords
- cell
- polyadenylation
- analysis
- peak
- transcript
- Prior art date
Links
- 238000004458 analytical method Methods 0.000 title claims abstract description 159
- 238000013136 deep learning model Methods 0.000 title abstract description 18
- 230000008488 polyadenylation Effects 0.000 claims abstract description 324
- 238000012163 sequencing technique Methods 0.000 claims abstract description 192
- 108090000623 proteins and genes Proteins 0.000 claims abstract description 139
- 230000014509 gene expression Effects 0.000 claims abstract description 102
- 238000000034 method Methods 0.000 claims abstract description 84
- 238000012549 training Methods 0.000 claims description 39
- 238000011144 upstream manufacturing Methods 0.000 claims description 33
- 239000011159 matrix material Substances 0.000 claims description 30
- 230000006870 function Effects 0.000 claims description 25
- 238000012216 screening Methods 0.000 claims description 25
- 238000012360 testing method Methods 0.000 claims description 25
- 230000003321 amplification Effects 0.000 claims description 23
- 238000003199 nucleic acid amplification method Methods 0.000 claims description 23
- 238000013527 convolutional neural network Methods 0.000 claims description 15
- 238000004445 quantitative analysis Methods 0.000 claims description 15
- 230000009467 reduction Effects 0.000 claims description 15
- 238000004364 calculation method Methods 0.000 claims description 12
- 238000003908 quality control method Methods 0.000 claims description 12
- 230000010354 integration Effects 0.000 claims description 10
- 239000003147 molecular marker Substances 0.000 claims description 10
- 239000000284 extract Substances 0.000 claims description 9
- 230000000717 retained effect Effects 0.000 claims description 9
- 238000011176 pooling Methods 0.000 claims description 8
- 102000010029 Homer Scaffolding Proteins Human genes 0.000 claims description 6
- 108010077223 Homer Scaffolding Proteins Proteins 0.000 claims description 6
- 238000007405 data analysis Methods 0.000 claims description 6
- 238000010195 expression analysis Methods 0.000 claims description 6
- 238000001422 normality test Methods 0.000 claims description 6
- 238000013528 artificial neural network Methods 0.000 claims description 5
- 239000003550 marker Substances 0.000 claims description 5
- 230000002159 abnormal effect Effects 0.000 claims description 4
- 230000004913 activation Effects 0.000 claims description 4
- 230000002457 bidirectional effect Effects 0.000 claims description 4
- 238000000605 extraction Methods 0.000 claims description 4
- 230000000306 recurrent effect Effects 0.000 claims description 4
- 230000006403 short-term memory Effects 0.000 claims description 4
- 230000013595 glycosylation Effects 0.000 claims description 3
- 238000006206 glycosylation reaction Methods 0.000 claims description 3
- 230000020477 pH reduction Effects 0.000 claims description 3
- 238000012800 visualization Methods 0.000 claims description 3
- 230000006835 compression Effects 0.000 claims description 2
- 238000007906 compression Methods 0.000 claims description 2
- 238000001514 detection method Methods 0.000 abstract description 6
- 230000001419 dependent effect Effects 0.000 abstract description 2
- 238000001914 filtration Methods 0.000 abstract 1
- 210000004027 cell Anatomy 0.000 description 143
- 108091032973 (ribonucleotides)n+m Proteins 0.000 description 21
- 108091034057 RNA (poly(A)) Proteins 0.000 description 15
- 210000003819 peripheral blood mononuclear cell Anatomy 0.000 description 14
- 230000008569 process Effects 0.000 description 14
- 238000010276 construction Methods 0.000 description 13
- 238000005516 engineering process Methods 0.000 description 12
- 238000011002 quantification Methods 0.000 description 11
- 238000011156 evaluation Methods 0.000 description 9
- 238000012545 processing Methods 0.000 description 9
- 230000001413 cellular effect Effects 0.000 description 8
- 238000011160 research Methods 0.000 description 7
- 239000002299 complementary DNA Substances 0.000 description 6
- 241000894007 species Species 0.000 description 6
- 108091092195 Intron Proteins 0.000 description 5
- 230000008901 benefit Effects 0.000 description 5
- 238000010205 computational analysis Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 5
- 239000012634 fragment Substances 0.000 description 5
- 230000035945 sensitivity Effects 0.000 description 5
- VYZAMTAEIAYCRO-UHFFFAOYSA-N Chromium Chemical compound [Cr] VYZAMTAEIAYCRO-UHFFFAOYSA-N 0.000 description 4
- 108020004414 DNA Proteins 0.000 description 4
- 229910052804 chromium Inorganic materials 0.000 description 4
- 239000011651 chromium Substances 0.000 description 4
- 238000002474 experimental method Methods 0.000 description 4
- 238000012165 high-throughput sequencing Methods 0.000 description 4
- 239000000463 material Substances 0.000 description 4
- 230000004048 modification Effects 0.000 description 4
- 238000012986 modification Methods 0.000 description 4
- 239000002773 nucleotide Substances 0.000 description 4
- 125000003729 nucleotide group Chemical group 0.000 description 4
- 238000007781 pre-processing Methods 0.000 description 4
- 238000013518 transcription Methods 0.000 description 4
- 230000035897 transcription Effects 0.000 description 4
- 208000028399 Critical Illness Diseases 0.000 description 3
- 108700024394 Exon Proteins 0.000 description 3
- 108091029795 Intergenic region Proteins 0.000 description 3
- 238000003559 RNA-seq method Methods 0.000 description 3
- 230000033228 biological regulation Effects 0.000 description 3
- 238000013480 data collection Methods 0.000 description 3
- 238000013135 deep learning Methods 0.000 description 3
- 201000010099 disease Diseases 0.000 description 3
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 3
- 230000004927 fusion Effects 0.000 description 3
- 210000002865 immune cell Anatomy 0.000 description 3
- 102000004169 proteins and genes Human genes 0.000 description 3
- 230000001105 regulatory effect Effects 0.000 description 3
- 238000012174 single-cell RNA sequencing Methods 0.000 description 3
- 230000014616 translation Effects 0.000 description 3
- 230000009385 viral infection Effects 0.000 description 3
- 101001112714 Homo sapiens NAD kinase Proteins 0.000 description 2
- 108091027974 Mature messenger RNA Proteins 0.000 description 2
- 102100023515 NAD kinase Human genes 0.000 description 2
- 108091028043 Nucleic acid sequence Proteins 0.000 description 2
- 230000000840 anti-viral effect Effects 0.000 description 2
- 210000003719 b-lymphocyte Anatomy 0.000 description 2
- 230000011712 cell development Effects 0.000 description 2
- 238000003776 cleavage reaction Methods 0.000 description 2
- 230000001186 cumulative effect Effects 0.000 description 2
- 210000004443 dendritic cell Anatomy 0.000 description 2
- 239000005547 deoxyribonucleotide Substances 0.000 description 2
- 125000002637 deoxyribonucleotide group Chemical group 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000018109 developmental process Effects 0.000 description 2
- 238000013507 mapping Methods 0.000 description 2
- 210000001616 monocyte Anatomy 0.000 description 2
- 101150084853 nadK gene Proteins 0.000 description 2
- 210000004180 plasmocyte Anatomy 0.000 description 2
- 238000000513 principal component analysis Methods 0.000 description 2
- 238000010839 reverse transcription Methods 0.000 description 2
- 230000002441 reversible effect Effects 0.000 description 2
- 230000007017 scission Effects 0.000 description 2
- 238000011222 transcriptome analysis Methods 0.000 description 2
- 238000013519 translation Methods 0.000 description 2
- 102100036657 26S proteasome non-ATPase regulatory subunit 7 Human genes 0.000 description 1
- 108020005345 3' Untranslated Regions Proteins 0.000 description 1
- 208000031261 Acute myeloid leukaemia Diseases 0.000 description 1
- 229930024421 Adenine Natural products 0.000 description 1
- GFFGJBXGBJISGV-UHFFFAOYSA-N Adenine Chemical compound NC1=NC=NC2=C1N=CN2 GFFGJBXGBJISGV-UHFFFAOYSA-N 0.000 description 1
- OIRDTQYFTABQOQ-KQYNXXCUSA-N Adenosine Natural products C1=NC=2C(N)=NC=NC=2N1[C@@H]1O[C@H](CO)[C@@H](O)[C@H]1O OIRDTQYFTABQOQ-KQYNXXCUSA-N 0.000 description 1
- 101100339431 Arabidopsis thaliana HMGB2 gene Proteins 0.000 description 1
- 239000002126 C01EB10 - Adenosine Substances 0.000 description 1
- 241000206602 Eukaryota Species 0.000 description 1
- 102100036241 HLA class II histocompatibility antigen, DQ beta 1 chain Human genes 0.000 description 1
- 108010065026 HLA-DQB1 antigen Proteins 0.000 description 1
- 102000055207 HMGB1 Human genes 0.000 description 1
- 108700010013 HMGB1 Proteins 0.000 description 1
- 101150021904 HMGB1 gene Proteins 0.000 description 1
- 101001136696 Homo sapiens 26S proteasome non-ATPase regulatory subunit 7 Proteins 0.000 description 1
- 101000868279 Homo sapiens Leukocyte surface antigen CD47 Proteins 0.000 description 1
- 101001129610 Homo sapiens Prohibitin 1 Proteins 0.000 description 1
- 101001080429 Homo sapiens Proteasome inhibitor PI31 subunit Proteins 0.000 description 1
- 101000905936 Homo sapiens RAS guanyl-releasing protein 2 Proteins 0.000 description 1
- 108700005091 Immunoglobulin Genes Proteins 0.000 description 1
- 102100032913 Leukocyte surface antigen CD47 Human genes 0.000 description 1
- 241001465754 Metazoa Species 0.000 description 1
- 208000033776 Myeloid Acute Leukemia Diseases 0.000 description 1
- 108091092724 Noncoding DNA Proteins 0.000 description 1
- 238000012408 PCR amplification Methods 0.000 description 1
- 102100031169 Prohibitin 1 Human genes 0.000 description 1
- 102100027565 Proteasome inhibitor PI31 subunit Human genes 0.000 description 1
- 108010029485 Protein Isoforms Proteins 0.000 description 1
- 102000001708 Protein Isoforms Human genes 0.000 description 1
- 102100023488 RAS guanyl-releasing protein 2 Human genes 0.000 description 1
- 102000009572 RNA Polymerase II Human genes 0.000 description 1
- 108010009460 RNA Polymerase II Proteins 0.000 description 1
- 239000002253 acid Substances 0.000 description 1
- 229960000643 adenine Drugs 0.000 description 1
- 229960005305 adenosine Drugs 0.000 description 1
- UDMBCSSLTHHNCD-KQYNXXCUSA-N adenosine 5'-monophosphate Chemical compound C1=NC=2C(N)=NC=NC=2N1[C@@H]1O[C@H](COP(O)(O)=O)[C@@H](O)[C@H]1O UDMBCSSLTHHNCD-KQYNXXCUSA-N 0.000 description 1
- -1 adenosine monophosphates Chemical class 0.000 description 1
- 230000006154 adenylylation Effects 0.000 description 1
- 239000011324 bead Substances 0.000 description 1
- 230000031018 biological processes and functions Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000010835 comparative analysis Methods 0.000 description 1
- 150000001875 compounds Chemical class 0.000 description 1
- 238000000205 computational method Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000004069 differentiation Effects 0.000 description 1
- 210000001671 embryonic stem cell Anatomy 0.000 description 1
- 238000013401 experimental design Methods 0.000 description 1
- 230000002068 genetic effect Effects 0.000 description 1
- 230000036737 immune function Effects 0.000 description 1
- 230000036039 immunity Effects 0.000 description 1
- 230000004957 immunoregulator effect Effects 0.000 description 1
- 238000009434 installation Methods 0.000 description 1
- 230000004807 localization Effects 0.000 description 1
- 238000010801 machine learning Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000035800 maturation Effects 0.000 description 1
- 239000012528 membrane Substances 0.000 description 1
- 108020004999 messenger RNA Proteins 0.000 description 1
- 230000002438 mitochondrial effect Effects 0.000 description 1
- 238000010369 molecular cloning Methods 0.000 description 1
- 230000009456 molecular mechanism Effects 0.000 description 1
- 210000000822 natural killer cell Anatomy 0.000 description 1
- 238000007481 next generation sequencing Methods 0.000 description 1
- 210000000056 organ Anatomy 0.000 description 1
- 235000015927 pasta Nutrition 0.000 description 1
- 230000035790 physiological processes and functions Effects 0.000 description 1
- 230000037452 priming Effects 0.000 description 1
- 238000012847 principal component analysis method Methods 0.000 description 1
- 230000008844 regulatory mechanism Effects 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B5/00—ICT specially adapted for modelling or simulations in systems biology, e.g. gene-regulatory networks, protein interaction networks or metabolic networks
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/02—Neural networks
- G06N3/04—Architecture, e.g. interconnection topology
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/02—Neural networks
- G06N3/08—Learning methods
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
- G16B30/10—Sequence alignment; Homology search
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Health & Medical Sciences (AREA)
- Theoretical Computer Science (AREA)
- Biophysics (AREA)
- General Health & Medical Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Medical Informatics (AREA)
- Evolutionary Biology (AREA)
- Biotechnology (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Data Mining & Analysis (AREA)
- Molecular Biology (AREA)
- Software Systems (AREA)
- Artificial Intelligence (AREA)
- Evolutionary Computation (AREA)
- Mathematical Physics (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Biomedical Technology (AREA)
- Computational Linguistics (AREA)
- Computing Systems (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Public Health (AREA)
- Epidemiology (AREA)
- Databases & Information Systems (AREA)
- Bioethics (AREA)
- Chemical & Material Sciences (AREA)
- Analytical Chemistry (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Physiology (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本发明提供了融合深度学习模型的单细胞转录组计算分析方法和系统。具体地，本发明提供了基于融合深度学习模型的多聚腺苷酸化位点(PAS)检测方法和系统以及用于单细胞测序数据的多聚腺苷酸化位点检测及细胞分型分析方法和系统(SCAPTURE方法和系统)。本发明SCAPTURE系统中构建和融合了一个深度学习模型，实现了高准确性和非位置依赖的多聚腺苷酸化位点预测，并用于全基因组水平高可信多聚腺苷酸化位点的筛选；并可利用单细胞转录组测序数据，通过检测全基因组区域的测序读段分布峰值，从头鉴定全基因组水平的多聚腺苷酸化位点；且可基于多聚腺苷酸化位点鉴定出同一基因上不同转录本及其差异表达，进而应用于单细胞分型分析。
Description
本发明涉及生物技术和生物信息学领域，具体地涉及基于融合深度学习模型的单细胞转录组计算分析方法和系统。
随着测序技术的发展，目前涌现了一系列大规模测定单细胞转录组的方法。不同于对细胞群体转录本进行测序的二代测序技术，单细胞转录组测序技术可以区分捕获的RNA分子的细胞和转录本来源，将转录组测序分析精度提高到单细胞水平。单细胞转录组测序主要基于基因表达差异，可揭示细胞群体、组织器官、甚至整个动物内的细胞异质性。
目前广泛使用的单细胞转录组测序技术根据其RNA扩增建库方法，主要分为3'tag-based(如10x Chromium)和full-length建库方法(如Smart-seq2)。通过利用基于微滴的微流控技术和细胞条形码分配系统，10x Chromium等3'tag-based单细胞转录组测序方法可以在一次实验中实现数千个单细胞平行测序。该方法可将条形码、引物凝胶珠与单个细胞包裹在微滴中，以oligo(dT)引物进行RNA分子的捕获，经过反转录获得cDNA并进行扩增，随后对cDNA进行打断，从分子3'端进行短片段测序。与full-length建库方法得到的测序读段在基因上的均匀覆盖相比，3'tag-based单细胞转录组测序方法获得的链特异性测序读段会在基因的3'末端富集。通过前期的研究发现，3'tag-based单细胞转录组测序读段富集位置通常接近基因末端的多聚腺苷酸化位点。这提示了我们可利用3'tag-based单细胞转录组测序数据，在不同单细胞/细胞类型中开展全基因组水平的多聚腺苷酸化位点的鉴定并进行生物学研究。
真核生物转录本成熟过程通常经历多聚腺苷酸化加工修饰，且广泛存在可变多聚腺苷酸化加工修饰形式，从而产生具有不同多聚腺苷酸化位点的转录本亚型。可变多聚腺苷酸化的动态调节参与许多重要生物学过程，如转录本剪接、蛋白翻译定位、组织特异性表达、疾病发生和细胞发育过程。目前已经开发了多种实验和计算方法，在细胞群体水平对全基因组多聚腺苷酸化修饰情况开展研究。但是，这些方法的广泛使用存在局限性。首先，针对多聚腺苷酸化位点鉴定开发的高通量测序实验方法通常依赖于特殊的化合物、建库流程和数据处理流程，不同方法之间的比较存在困难。其次，基于Bulk RNA-seq开发的多聚腺苷酸化分析工具存在使用局限性，例如：需要多组别实验设计、基因本身3' UTR长度限制或需要结合额外多聚腺苷酸化测序数据等。最后，还有研究通过机器学习或深度学习模型对基因组上的多聚腺苷酸化位点进行预测，如DeepPASTA和APARENT等，但这些工具仍存在多聚腺苷酸化位点位置高度依赖、准确度不足等局限。
近期，有研究结果表明单细胞转录组测序技术具有研究多聚腺苷酸化的潜力。Velten等开发了BATSeq分析单细胞亚型，发现小鼠胚胎干细胞中存在不同3'端转录本的使用；Ye等开发了基于直方图差异的方法scDAPA来检测急性髓系白血病中具有动态可变多聚腺苷酸化的基因。最近新发表的两个单细胞转录组测序的多聚腺苷酸化分析工具scAPA和Sierra，基于peak识别的方法实现了基因组多聚腺苷酸化位点识别和单细胞水平的多聚腺苷酸化定量分析，并验证了不同生物体系中可变多聚腺苷酸化和不同转录本的动态变化，表明单细胞转录组测序能够应用于多聚腺苷酸化的研究。
但是，上述技术和工具仍存在一些问题，例如：
(1)BATSeq分析多聚腺苷酸化依赖于其特殊的建库测序方法，适用性较窄。
(2)scDAPA通过测序读段在基因上的直方图变化来判断是否存在可变多聚腺苷酸化加工，无法识别具体的基因组多聚腺苷酸化位点和单细胞水平定量等分析。
(3)scAPA与Sierra采用基于peak识别的方法，但都基于测序样品整体测序读段比对结果进行peak识别，忽略了单细胞转录组测序中基因表达信号的高变异度，分析相对低丰度基因的灵敏度不足。
(4)scAPA仅限于使用软件内建的基因注释，无法实现研究人员自定义研究物种和其它基因注释的分析；同时scAPA不支持转录本剪接位点的识别，在准确定量上存在局限性。
(5)Sierra能够识别转录本剪接位点，但与scAPA相同，这两个工具都缺乏对peak的置信度筛选，在多聚腺苷酸化位点的识别上存在假阳性问题。由于3'tag-based单细胞转录组测序需通过oligo(dT)引物进行RNA分子的捕获，该过程不仅会捕获到多聚腺苷酸尾的RNA，也会捕获到分子内部富含腺苷酸的RNA，引入internal priming的假阳性干扰。
(6)DeepPASTA和APARENT等多聚腺苷酸化位点预测工具也存在碱基位置高度依赖、准确度不佳等局限。
(7)以上工具也无法通过多聚腺苷酸化位点识别转录本表达，缺少将多聚腺苷酸化位点与转录本表达耦联的单细胞分析方法。
因此，本领域迫切需要开发一种能够利用3'tag-based单细胞转录组测序在 单细胞水平精确鉴定和筛选高可信多聚腺苷酸化位点的方法，用于单细胞测序数据的多聚腺苷酸化位点检测及细胞分型分析。
发明内容
本发明目的是提供一套基于单细胞转录组测序分析系统，该系统能够利用单细胞转录组测序在单细胞水平精确鉴定和筛选高可信多聚腺苷酸化位点，并且能够利用多聚腺苷酸化位点信息识别不同转录本表达，进一步用于转录本水平的单细胞分型分析。
在本发明的第一方面，提供了一种构建多聚腺苷酸化位点的预测模型的方法，所述方法包括步骤：
(S1)提供用于模型训练的数据集，所述数据集包括正集和负集，其中，所述的正集中包含的序列为含多聚腺苷酸化位点的序列，即正序列；而所述的负集中包含的序列为不含多聚腺苷酸化位点的序列，即负序列；
其中，所述的正序列包括：对于多聚腺苷酸化原始位点而言的基础正序列和基于所述基础正序列通过偏移扩增所获得的偏移正序列，并且所述的基础正序列和偏移正序列均含有所述的多聚腺苷酸化原始位点；
(S2)将(S1)中的数据集用于模型架构进行训练；
其中，所述的模型架构包括基于卷积神经网络(CNN)模块和循环神经网络(RNN)模块的模型架构；
所述CNN模块包括两个卷积层和最大池化层(max pooling layer)，其中，所述卷积层的输出结果通过最大池化层(max pooling layer)进行降维压缩，并选择4×4邻域内的最大值输出到下游RNN模块；
所述RNN模块设有双向长短期记忆网络结构(BiLSTM)，并将CNN模块的输出内容传递到下游全连接层；
所述的全连接层设置多个隐藏单元(hidden units)和预定的丢码率(dropout ratio)(如约30％)，并通过归一化指数函数(softmax activation function)对最终二分类的分类结果进行输出；
(S3)当训练结果达到预定终止条件，则终止模型训练，并获得所述的多聚腺苷酸化位点的预测模型(即DeepPASS模型)。
在另一优选例中，所述的偏移扩增包括：随机的或非随机的偏移扩增。
在另一优选例中，所述的偏移扩增包括：正态分布随机偏移。
在另一优选例中，所述的偏移扩增指：将一个基础正序列置于一坐标轴x1上，并以基础正序列的多聚腺苷酸化原始基因组位点作为坐标零点，然后使得 该基础序列的5'端和3'端各自独立地进行p个bp的正偏移或负偏移，并将偏移后序列称为偏移正序列，其中p为正整数(bp)。
在另一优选例中，所述的基础正序列的长度和对应的偏移正序列的长度是相同或基本相同的。
在另一优选例中，当基础正序列的坐标为(X
Left，X
Right)且基础正序列和对应的偏移正序列的长度是相同时，则进行p个bp的正偏移后，偏移正序列的坐标为(X
Left+p，X
Right+p)；而进行p个bp的负偏移后，偏移正序列的坐标为(X
Left-p，X
Right-p)。
在另一优选例中，设基础正序列的长度为2L个bp，则p为≤1/4L，较佳地≤1/5L，更佳地≤1/10L。
在另一优选例中，设础正序列的多聚腺苷酸化原始位点作为坐标零点，且基础正序列的长度为2L个bp，L为约100bp，则所述的基础正序列的坐标为(X
Left，X
Right)，X
Right-X
Left＝2L。
在另一优选例中，L为60-250，较佳地80-150，更佳地90-110。
在另一优选例中，|X
Left|与|X
Right|是相同或基本相同的，例如|X
Left|与|X
Right|之比为0.8-1.2，较佳地0.9-1.1，更佳地约1.0。
在另一优选例中，在所述的正集中，基础正序列的总数与偏移正序列的总数之比为1:2至1:20，较佳地为1:5-1:10，更佳地为1:10-1:20。
在另一优选例中，所述的正集中所含的正序列包括针对T0(约25万)个多聚腺苷酸化原始位点。
在另一优选例中，在所述的正集中，对于每一个多聚腺苷酸化原始位点，其对应具有1个基础正序列与D个偏移正序列，其中D为≥2的正整数，较佳地D为2-5，更佳地5-10，最佳地10-20。
在另一优选例中，在所述的正集中，所述的基础正序列的总数T1与多聚腺苷酸化原始位点的总数T0是相同或基本相同的。
在另一优选例中，所述的正集和负集所包含的序列来自同一物种或不同物种。
在另一优选例中，所述的正集和负集所包含的序列来自同一物种，如人。
在另一优选例中，正集包括N1条含多聚腺苷酸化位点的序列(正序列)；而负集包括N2条不含多聚腺苷酸化位点的序列(负序列)。
在另一优选例中，正集的序列数量N1与负集的序列数量N2之比没有特别限制，通常为1:10至10:1，较佳地1:5至5:1，更佳地1:2至2:1，如约1:1。
在另一优选例中，N1≥50万，更佳地≥100万，最佳地≥200万条。
在另一优选例中，N2≥50万，更佳地≥100万，最佳地≥200万条。
在另一优选例中，所述的预测模型为本发明构建多聚腺苷酸化位点的预测模型(以下简称该预测模型为DeepPASS)。
在另一优选例中，所述的预测模型用于判断输入的序列是否存在高置信度的多聚腺苷酸化位点。
在另一优选例中，所述的序列可为RNA或DNA格式的序列，或所述的测序数据为RNA测序数据。
在另一优选例中，所述的高置信度为DeepPASS输出结果为1(含有多聚腺苷酸化位点)的预测。
在另一优选例中，在步骤(S1)中，所述的偏移扩增包括：以多聚腺苷酸化原始基因组位点作为坐标零点，用正态分布对位点进行随机偏移，对一个多聚腺苷酸化位点进行4-20次偏移，提取偏移扩增后位点的上游100±20bp到下游100±20bp的基因组序列作为模型训练的正集中的偏移正序列。
在另一优选例中，所述的偏移扩增包括：以多聚腺苷酸化原始位点作为坐标零点，用正态分布(μ＝0，σ＝10)对位点进行随机偏移，每个多聚腺苷酸化位点进行约10次偏移(十倍扩增位点数量)，提取扩增后位点的上游100bp到下游100bp的基因组序列作为DeepPASS模型训练的正集中的偏移正序列。
在另一优选例中，收集了四个不同来源的多聚腺苷酸化位点数据库：PolyA_DB3、PolyA-seq、Poly(A)Site2.0和GENCODE。筛选三个已发表数据库(PolyA_DB3、PolyA-seq和Poly(A)Site2.0)中共有的位点，与GENCODE合并作为已知多聚腺苷酸化位点库，用于DeepPASS模型训练。
在另一优选例中，收集了四个不同来源的多聚腺苷酸化位点数据库：PolyA_DB3、PolyA-seq、Poly(A)Site2.0和GENCODE。筛选三个已发表数据库(PolyA_DB3、PolyA-seq和Poly(A)Site2.0)中高可信的多聚腺苷酸化位点，筛选条件为：一数据库中的多聚腺苷酸化位点在其基因组位置上下游12bp内与其它两个数据库中都有重合。将这三个数据库中的高可信位点，与GENCODE合并作为已知多聚腺苷酸化位点库，用于多聚腺苷酸化位点注释。
在另一优选例中，所述的负集中的负序列是通过以下方式进行筛选：通过取四个多聚腺苷酸化位点数据库的并集，对基因组的基因间区进行反向筛选，只选取已知多聚腺苷酸化位点的上下游距离100bp以外的基因间区序列。
在另一优选例中，所述DeepPASS模型选取二分类中概率较大的分类作为预测结果进行输出，即通过1(含有多聚腺苷酸化位点)或0(不含有多聚腺苷酸化位点)的分类结果判断输入序列中是否包含高可信多聚腺苷酸化位点。
在另一优选例中，在步骤(S2)中，模型训练设置数据量批次大小为3000-7000(如5000)，模型训练50-150期(如100期)。
在另一优选例中，所述的预定终止条件包括：训练结果在5-20期内无提升。
在另一优选例中，在步骤(S2)中，还包括：训练中的模型训练结果进行评估，较佳地所述评估是用接受者操作特征曲线(ROC)和曲线下面积(AUC)进行评估。
在另一优选例中，所述的基于CNN与RNN建立的嵌合深度学习模型DeepPASS，理论上可识别多聚腺苷酸化位点附近的基序(motif)如AAUAAA和上下游GU重复、U富集等的序列特征，并建立基序特征之间潜在的关联性，达到整体识别的多聚腺苷酸化信号的效果。
在另一优选例中，所述的DeepPASS的模型预测仅需200bp的碱基序列作为输入，且输出为二分类(1或0)的多聚腺苷酸化位点判断。
在另一优选例中，DeepPASS模型与其它模型的评估使用了GENCODE数据库的多聚腺苷酸化原始位点上游100bp至下游100bp的序列作为正测试集，已知多聚腺苷酸化位点的上下游距离100bp以外区域的基因组间区200bp的序列作为负测试集。
在本发明的第二方面，提供了一种多聚腺苷酸化位点的预测系统，所述的系统包括：
输入单元，所述输入单元被配置为输入数据，所述数据包括待预测的输入序列，其中，需要判断在所述待预测的输入序列中是否存在多聚腺苷酸化位点；
多聚腺苷酸化位点的预测单元，所述的预测单元被配置为执行一多聚腺苷酸化位点的预测模型，从而获得所述待预测的输入序列中是否存在多聚腺苷酸化位点的预测结果；其中，所述的预测模型是用本发明第一方面所述的方法构建的；
输出单元，所述的输出单元被配置为输出所述多聚腺苷酸化位点的预测单元的预测结果。
在另一优选例中，所述预测模型为DeepPASS模型，其选取二分类中概率较大的分类作为预测结果，即通过“是”(或1)或“否”(或0)的分类结果预测给定的输入序列中是否包含真实的多聚腺苷酸化位点。
在本发明的第三方面，提供了一种基于单细胞转录组测序数据识别多聚腺苷酸化测序读段信号峰(peak)的分析方法，其特征在于，包括步骤：
(Y1)将单细胞转录组测序数据与基因组序列进行比对，得到BAM格式存储的测序读段比对结果文件，BAM文件中应包含：BAM格式要求的测序读段比对结果信息、测序读段细胞条形码(CB)和测序读段分子标识(UMI)信息；
(Y2)将基因水平注释文件，根据其基因和转录本信息，转换成“基因-转录本”标记的基因转录本水平注释文件；
(Y3)根据测序读段比对结果文件和基因转录本水平注释文件，利用HOMER findPeaks在转录本水平识别测序读段分布peak，识别每个转录本对应的所有peak，并基于对应转录本的坐标信息标记peak的区间位置；
(Y4)结合上述步骤(Y2)-(Y3)得到的基因转录本水平注释文件和基于转录本坐标的peak区间位置，将转录本位置系统映射到基因组坐标位置系统，计算出peak的基因组坐标系统注释文件；
(Y5)结合上述步骤(Y2)-(Y4)得到的基因转录本水平注释文件和基于基因组坐标系统的peak文件，通过识别基因转录本水平注释文件中的转录本外显子剪接位点，进一步计算出对应转录本的peak剪接位点信息；
(Y6)根据上述步骤(Y5)中得到的peak注释结果，对识别的peak的测序读段信号分布曲线进行正态分布检验(normality test)，去除可能源于异常扩增、多重比对或信号丰度较低的peak；并进行单峰检验(unimodality test)，保留来自邻近多聚腺苷酸化信号形成的相互重叠的peak；
(Y7)上述步骤(Y2)-(Y6)计算可得到来自于同一基因的不同转录本来源的peaks，根据其基因信息可对peak进行去重；优选地，去重方法为：对于区间重合比例超过50％的peak，对其测序读段丰度进行排序，保留表达量最高的peak；
(Y8)对上述步骤(Y7)去重后的peak重新计算测序读段丰度，过滤掉同一基因内测序读段数量占比≤1％的peak，得到初步筛选的peak(raw peak)；
(Y9)对于多个样品中每个样品经过上述步骤(Y3)-(Y8)得到的raw peak，通过标记peak来源的样品信息，并将不同样品之间重合比例超过60％的peak标记为同一组，对同一组的peak中选取唯一值作为该多聚腺苷酸化位点的peak；
(Y10)将上述步骤识别的raw peak注释区间的末端近似为多聚腺苷酸化位点，提取位点上游100bp到下游100bp序列并输入DeepPASS模型进行预测，判断每一个raw peak是否为高可信多聚腺苷酸化的peak(PAS)，从而获得单细胞转录组测序数据的多聚腺苷酸化测序读段信号峰分析结果。
在另一优选例中，在步骤(Y10)中，所述的预测模型为用本发明第一方面所述的方法构建的预测模型，即DeepPASS模型。
在另一优选例中，在步骤(Y10)中，取所述多聚腺苷酸化位点上游100bp到下游100bp的碱基序列输入预测模型进行预测。
在另一优选例中，用DeepPASS模型进行预测，以判断每一个raw peak是 否为高可信多聚腺苷酸化的peak(PAS)，并筛选出的高可信peak用于下游分析。
在另一优选例中，所述的分析方法用于对全基因组水平的多聚腺苷酸化测序读段信号峰进行分析。
在另一优选例中，在步骤(Y2)中，对基因注释进行了预处理，用基因和转录本的对应关系构建转录本为单位的注释文件，代替原有以基因为单位的注释文件。
在另一优选例中，在步骤(Y3)中，以转录本为单位对peak进行识别并得到转录本坐标的注释区间，同时根据先验值设定peak宽度为固定长度(如400bp)使其适用于10x Chromium技术的单细胞转录组测序数据，也可根据不同测序方法设定为其它长度数值。
在另一优选例中，在步骤(Y4)-(Y5)中，对peak的注释区间标进行了转录本位置到基因组位置的映射，并解析转录本注释中外显子的剪接位点，将peak对应转录本的剪接位点用于进一步注释该peak，得到识别测序读段剪接信息的基因组坐标peak注释。
在另一优选例中，在步骤(Y6)中，对peak的质控引入了正态分布检验及单峰检验。将peak区间每个碱基位置作为自变量x，每个碱基的测序读段深度作为因变量y，根据x与y的值做分布曲线；用正态分布检验评估该曲线是否接近正态分布，对建库中异常扩增的转录本片段、基因组比对过程中异常比对、丰度和覆盖度较低的peak进行去除优化；用单峰检验进一步评估该曲线是否符合多峰分布，对来自邻近多聚腺苷酸化信号形成的相互重叠的peak进行保留。
在另一优选例中，在步骤(Y7)中，对于同一个多聚腺苷酸化位点产生的测序读段分布峰，在同一基因注释中覆盖该加尾位点的多个转录本注释会计算出多个peak，对这些peak的去冗余方法为：通过计算peak区间重合比例，将重合比例大于50％的peak进行归类分组，对于同一组重合的所有peak仅保留测序读段覆盖深度最高的peak注释。
在另一优选例中，在步骤(Y8)中，对同一基因内的peak测序读段丰度用stringtie进行估计，保留读段丰度占比≥1％的peak。
在另一优选例中，在步骤(Y9)中，识别的raw peak基因组注释区间的末端位置可作为近似多聚腺苷酸化的位点，该位点提取的上游100bp到下游100bp的碱基序列可输入DeepPASS模型进行多聚腺苷酸化预测，预测结果为1的peak作为高可信的多聚腺苷酸化位点。
在另一优选例中，所述单细胞转录组数据为3'tag-based单细胞测序数据，依靠oligo(dT)捕获RNA分子并富集RNA分子3'末端的高通量单细胞测序技术。
在本发明的第四方面，提供了一种在单细胞水平对转录本进行定量分析和/或进行单细胞分型分析的方法，所述方法包括：
(Z1)提供待分析的单细胞样本的单细胞转录组测序数据；
(Z2)用本发明第三方面所述的分析方法，对单细胞转录组测序数据进行分析，从而获得高可信peak；
(Z3)基于上一步骤中获得的高可信peak，进行转录本进行定量分析和/或进行单细胞分型分析。
在另一优选例中，在步骤(Z3)中，包括：
(W1)对于多样品的peak整合分析，整合方法为：将不同样品之间重合的peak标记为同一组，对同一组内的peak进行置信度排序，选取置信度最高的peak作为该多聚腺苷酸化位点的唯一peak；
(W2)根据上述步骤(W1)筛选的peak，通过peak的剪接位点和基因组区间信息将其映射到对应来源的转录本上，可鉴定单细胞转录组测序样品中由可变多聚腺苷酸化(APA)产生的转录本，实现转录本表达鉴定；
(W3)将映射转录本信息的PAS注释构建GTF格式注释文件，通过featureCounts将比对过的BAM中reads重新计算PAS注释的分配情况，得到标记PAS注释的re-assigned BAM文件；利用UMI-tools可对re-assigned BAM文件中每一个细胞条形码对应的所有PAS注释的唯一分子标识(UMI)信息进行统计，获得单细胞水平的转录本表达矩阵，矩阵包括每个转录本在每个细胞的唯一分子标识数目，可用于单细胞转录组的分析；
(W4)对于转录本水平的单细胞分析，将上述步骤(W3)得到基于peak注释转录本水平表达矩阵作为单细胞分析工具Seurat的输入进行单细胞分析；优选地，基于peak注释转录本水平的单细胞分析体系包含如下主要部分：表达质控，特征提取，降维分析，无监督细胞聚类，转录本分子标记物表达分析；
(W5)对于输入的转录本表达矩阵进行表达质控，保留至少表达细胞数占比≥0.5％的转录本；优选地，利用Seurat的内置函数对筛选后的表达矩阵进行转录本表达的标准化和表达变异计算，将表达变异最高的2000个高变异转录本作为特征转录本用于下游降维分析；
(W6)上述步骤(W5)得到的特征转录本可通过Seurat的内置函数进行降维分析和无监督细胞聚类；优选地，利用主成分分析方法(PCA)进行降维，并保留主要(如前50个)主成分结果；根据该50个主成分对总体单细胞构建共享最近邻(SNN)相似度图，以无监督方式进行特定分辨率(如0.3)进行细胞集群划分，最终得到基于PAS注释转录本表达的细胞分类结果；
(W7)根据步骤(W2)和(W3)中得到peak注释信息，计算筛选给定分子标记 物基因表达的转录本，利用Seurat的内置函数可对不同细胞分类的分子标记物转录本的表达量进行计算和可视化，鉴定无监督细胞分类结果中的细胞类型。
在另一优选例中，在步骤(W1)中，对于多样品的peak整合分析，整合方法为：(1)源于同一多聚腺苷酸化位点的不同样品的peak，计算该区间内每个peak末端区间(上游50bp到下游25bp)覆盖的已知多聚腺苷酸位点的数量，根据覆盖位点数量的数值对peak进行递减排序；(2)提取peak末端上游100bp到下游100bp的区间序列，输入DeepPASS模型进行预测，根据预测结果从1到0排序；(3)在同一组peak中，选取综合排序最优的peak作为多样品分析中该多聚腺苷酸化位点的唯一peak。
在另一优选例中，在步骤(W2)中，利用peak剪接位点和基因组注释信息，对符合peak来源的转录本进行鉴定；优选地，鉴定方法为：转录本外显子剪接位点需与peak剪接位点信息完全匹配，且选取末端注释距离最近的转录本作为peak来源的转录本，最终得到全基因组水平的单细胞转录组测序样品中的不同转录本表达。
在另一优选例中，在步骤(W3)中，利用PAS注释信息和转录本来源信息，构建Gene transfer format(GTF)格式的PAS转录注释本文件文件，并将BAM文件的测序读段在GTF格式的PAS转录注释本文件中进行读段重新分配，标记每一个测序读段所属的PAS注释。
在另一优选例中，在步骤(W3)中，根据PAS转录注释本文件，计算出每个单细胞每个PAS区间的测序读段，并进一步统计测序读段上唯一分子标识的累计数，得到基于PAS注释的单细胞转录本表达矩阵。
在另一优选例中，在步骤(W4)-(W6)中，使用步骤(W3)中peak注释的转录本水平表达矩阵，替代传统的基因水平表达矩阵，开展基于转录本定量数据的单细胞转录组数据细胞分型分析；优选地，该分析流程包括基于转录本定量数据的表达质控、特征提取、降维分析和无监督细胞聚类和转录本分子标记物表达分析。
在另一优选例中，所述单细胞转录组数据为3'tag-based单细胞测序数据，包括依靠oligo(dT)捕获RNA分子并富集RNA分子3'末端的高通量单细胞测序技术。
在另一优选例中，步骤Z3包括多聚腺苷酸化分析。
在另一优选例中，所述的多聚腺苷酸化分析还包括：
(W8)根据单细胞分型和多聚腺苷酸化位点定量结果，进行不同类型不同组别的细胞间差异化比较和统计，计算筛选得到单细胞水平全转录组多聚腺苷酸化趋势以及多聚腺苷酸化差异变化的基因，从而获得单细胞水平多聚腺苷酸 化差异变化的分析结果。
在另一优选例中，所述步骤W8，包括以下子步骤：
(L1)利用所述的注释文件，对多聚腺苷酸化位点进行分类；
优选地，子步骤L1包括：选取含有2个以上多聚腺苷酸化位点的基因，并根据基因组注释坐标的远近，分为“近端”、“远端”和“中间”三种不同多聚腺苷酸化位点；
(L2)利用所述的注释文件和定量结果，计算位点相对使用率表达矩阵；
优选地，子步骤L2包括：该基因的“近端”(或远端)位点表达值，除以该基因其余位点表达值的加和，结果再取2为底的对数；并对每个细胞的每个基因进行计算，即可得到位点相对使用率×细胞的数值矩阵；
(L3)利用(L2)得到的位点相对使用率，计算细胞整体的多聚腺苷酸趋势；
优选地，子步骤L3包括：对每个细胞的所有基因的计算位点相对使用率的平均值，再根据需要将细胞划分为不同组别，对不同组别细胞的位点相对使用率进行累计分布统计，评估细胞整体的多聚腺苷酸化倾向；
(L4)利用(L2)得到的位点相对使用率，筛选多聚腺苷酸化存在差异的基因；
优选地，子步骤L4包括：将位点相对使用率×细胞的数值矩阵构建Seurat对象，通过Seurat的内置函数FindMarkers对任意所需细胞组别间进行差异分析，得到不同组别下位点相对使用率存在差异的基因。
在本发明的第五方面，提供了一种单细胞转录组测序数据分析系统，所述的系统包括：
(a)输入单元，所述输入单元被配置为输入序列数据，所述的输入序列数据选自下组：(i)DeepPASS预测模型和/或待预测的输入序列，其中，需要判断在所述待预测的输入序列中是否存在多聚腺苷酸化位点；(ii)单细胞转录组测序数据；(iii)所述(i)和(ii)的组合；
(b)分析模块，所述分析模块被配置为对输入的数据进行预定的分析，从而获得分析结果，并且所述分析模块选自下组：
(M1)多聚腺苷酸化位点的预测单元，所述的(M1)预测单元被配置为执行一多聚腺苷酸化位点的预测模型，从而获得所述待预测的输入序列中是否存在多聚腺苷酸化位点的预测结果；其中，所述的预测模型是用本发明第一方面所述的方法构建的；
(M2)多聚腺苷酸化测序读段信号峰(peak)的分析单元，其中，所述的(M1)多聚腺苷酸化测序读段信号峰(peak)的分析单元被配置为执行本发明第 三方面所述的分析方法，从而获得单细胞转录组测序数据的多聚腺苷酸化测序读段信号峰分析结果；
(M3)单细胞水平转录本定量分析和/或单细胞分型分析单元，所述的(M3)单细胞水平转录本定量分析和/或单细胞分型分析单元被配置为执行本发明第五方面所述的分析方法，从而获得单细胞水平转录本定量分析结果和/或单细胞分型分析结果；
(Mx)上述(M1)、(M2)和(M3)的组合。
(c)输出单元，所述的输出单元被配置为输出所述分析模块的分析结果。
在另一优选例中，所述的分析系统包括单细胞多聚腺苷酸化分析系统。
在另一优选例中，所述的转录组测序数据分析包括多聚腺苷酸化分析。
在另一优选例中，所述的分析模块(b)还包括：
(M4)多聚腺苷酸化差异变化的分析单元，所述的(M4)多聚腺苷酸化差异变化的分析单元被配置为执行以下步骤W8：根据单细胞分型和多聚腺苷酸化位点定量结果，进行不同类型不同组别的细胞间差异化比较和统计，计算筛选得到单细胞水平全转录组多聚腺苷酸化趋势以及多聚腺苷酸化差异变化的基因，从而获得单细胞水平多聚腺苷酸化差异变化的分析结果；
(Mx)上述(M1)、(M2)、(M3)和(M4)的任意组合。
在另一优选例中，所述步骤W8包括以下子步骤：
(L1)利用所述的注释文件，对多聚腺苷酸化位点进行分类；
优选地，子步骤L1包括：选取含有2个以上多聚腺苷酸化位点的基因，并根据基因组注释坐标的远近，分为“近端”、“远端”和“中间”三种不同多聚腺苷酸化位点；
(L2)利用所述的注释文件和定量结果，计算位点相对使用率表达矩阵；
优选地，子步骤L2包括：该基因的“近端”(或远端)位点表达值，除以该基因其余位点表达值的加和，结果再取2为底的对数；并对每个细胞的每个基因进行计算，即可得到位点相对使用率×细胞的数值矩阵；
(L3)利用(L2)得到的位点相对使用率，计算细胞整体的多聚腺苷酸趋势；
优选地，子步骤L3包括：对每个细胞的所有基因的计算位点相对使用率的平均值，再根据需要将细胞划分为不同组别，对不同组别细胞的位点相对使用率进行累计分布统计，评估细胞整体的多聚腺苷酸化倾向；
(L4)利用(L2)得到的位点相对使用率，筛选多聚腺苷酸化存在差异的基因；
优选地，子步骤L4包括：将位点相对使用率×细胞的数值矩阵构建Seurat对象，通过Seurat的内置函数FindMarkers对任意所需细胞组别间进行差异分 析，得到不同组别下位点相对使用率存在差异的基因。
应理解，在本发明范围内中，本发明的上述各技术特征和在下文(如实施例)中具体描述的各技术特征之间都可以互相组合，从而构成新的或优选的技术方案。限于篇幅，在此不再一一累述。
图1显示了实施例1中DeepPASS模型的数据收集、预处理、模型搭建和训练评估。
图2显示了实施例2中单细胞转录组数据的多聚腺苷酸化位点测序读段peak的识别和筛选。
图3显示了实施例3中基于多聚腺苷酸化位点鉴定转录本表达和应用于单细胞分型分析。
图4显示了实施例4中不同条件下多种细胞类群的多聚腺苷酸化表达差异分析。
本发明人经过广泛而深入的研究，首次开发了一种基于融合深度学习模型的单细胞转录组计算分析方法和系统。具体地，本发明人设计了一种融合深度学习模型的多聚腺苷酸化位点(PAS)检测方法，用于单细胞测序数据的多聚腺苷酸化位点检测及细胞分型分析(以下简称该计算分析系统为SCAPTURE)。第一，本发明SCAPTURE系统中构建和融合了一个深度学习模型，实现了高准确性和非位置依赖的多聚腺苷酸化位点预测，并用于全基因组水平高可信多聚腺苷酸化位点的筛选；第二，本发明SCAPTURE系统还可以利用单细胞转录组测序数据，通过检测全基因组区域的测序读段分布峰值，从头鉴定全基因组水平的多聚腺苷酸化位点；第三，本发明SCAPTURE系统还基于可变多聚腺苷酸化位点鉴定出同一基因上不同转录本及其差异表达，并应用于单细胞分型分析。在此基础上完成了本发明。
术语
如本文所用，术语“本发明方法”指出基于融合深度学习模型的单细胞转录组计算分析方法。应理解，本发明方法包括了以下三方面中任一方法或其组合：(a)构建DeepPASS模型预测和筛选高可信的多聚腺苷酸化位点；(b)从单细胞转录组数据从头识别和筛选出全基因组水平的多聚腺苷酸化测序读段位点； (c)根据多聚腺苷酸化位点鉴定单细胞转录组测序数据中鉴定出不同转录本及其差异表达，并进行单细胞水平的转录本定量和单细胞分型分析。
如本文所用，术语“DeepPASS(
Deep neural network for
poly
adenylation prediction with
sequences
shifting)模型”指本发明方法利用深度学习方法搭建的用于多聚腺苷酸化位点预测的模型。
如本文所用，术语“SCAPTURE(
scRNA-seq
analysis of
PASs and their corresponding
transcript expression
used to
refine cell identities)”特指本发明建立的利用单细胞转录组鉴定和预测评估多聚腺苷酸加尾位点并在单细胞水平对选择性表达转录本定量的分析系统，包括了以下三方面：(a)构建DeepPASS模型预测和筛选高可信的多聚腺苷酸化位点；(b)从单细胞转录组数据从头识别和筛选出全基因组水平的多聚腺苷酸化测序读段位点；(c)根据多聚腺苷酸化位点鉴定单细胞转录组测序数据中鉴定出不同转录本及其差异表达，并进行单细胞水平的转录本定量和单细胞分型分析。
单细胞转录组测序(Single-cell RNA sequencing,scRNA-seq)：是在单个细胞水平对全转录组进行扩增与高通量测序的一项新技术，针对单个细胞研究其整体水平的基因表达情况，旨在解决细胞分子机制研究中常见的细胞异质性、细胞量少而无法进行常规高通量测序等难题。
3'tag-based单细胞转录组测序(3'tag-based scRNA-seq)：通过使用带有唯一分子标识和细胞条形码的oligo(dT)引物对单细胞RNA分子进行捕获和反转录，cDNA扩增和建库，及基于cDNA的3'端特异性测序的单细胞测序方法。
测序读段(Sequence Read)：测序技术所产生的单个测序片段，在本说明书中是指3'tag-based单细胞转录组测序中的cDNA测序片段。
细胞条形码(Cell barcode,CB)：一段短脱氧核糖核苷酸序列，长度通常在10-20nt。在单细胞转录组测序过程中标记同一细胞内的RNA分子，不同的细胞具有不同的细胞标签，用以区分测序读段的细胞来源。
唯一分子标识(Unique molecular identifier,UMI)：一段短脱氧核糖核苷酸序列，长度通常在8-15nt。用于对同一转录本分子来源扩增产物进行追踪，用于排除PCR扩增和测序引入的重复和偏好性，以获得准确的读数进行定量分析。
基因(Gene)：在生物学中指DNA或RNA内编码基因产物(RNA或蛋白质)的合成的核苷酸序列。基因也可视作基本遗传单位，亦即一段具有功能性的DNA或RNA序列。
转录本(Transcript)：由一个基因通过转录形成的一种或多种成熟的RNA分子。基因转录形成的RNA可通过不同加工过程形成不同的转录本。
外显子(Exon)：是真核生物基因中能转录，且在剪接后会被保存下来相应 DNA区域。在本说明书中是指用来分析外显子剪接位点时的外显子名称。
内含子(Intron)：一个基因中非编码DNA片段，它分开相邻的外显子。内含子是阻断基因线性表达的序列。DNA上的内含子会被转录到前体mRNA中，但RNA上的内含子会在RNA离开细胞核进行翻译前被剪除。
基因注释文件：在本说明书中是指记录基因及其转录本在基因组上的名称、基因位置和外显子剪接信息的注释文件格式，可以是genepred、GTF或BED格式。
Peak注释文件：在本说明书中是指记录单细胞转录组测序中基因多聚腺苷酸化位点测序读段峰值在基因组上的基因及来源转录本名称、测序读段峰值区间的位置和外显子剪接信息的注释文件格式，可以是genepred、GTF或BED格式。
多聚腺苷酸化(Polyadenylation)：将poly(A)尾部添加到RNA转录本上的过程。Poly(A)尾巴由多个单磷酸腺苷组成，是只有腺嘌呤碱基的RNA片段。在真核生物中，多聚腺苷酸化是产生成熟mRNA进行翻译的过程的一部分。
多聚腺苷酸化位点(Polyadenylation site,PAS)：在RNA多聚腺苷酸化加工过程中添加poly(A)尾的位点，位于多聚腺苷酸化信号的下游，切割位点通常是但不限于CA碱基。
多聚腺苷酸化信号(Polyadenylation signal)：经典核心序列为AAUAAA，通常位于多聚腺苷酸化信号上游20-30bp，可以被切割和聚腺苷酸化特异性因子识别。
可变多聚腺苷酸化(Alternative polyadenylation,APA)：可变多聚腺苷酸化是一种广泛存在的基因调节机制，可在RNA聚合酶II转录产生的转录本中产生不同的3'末端。
构建DeepPASS模型预测筛选高可信的多聚腺苷酸化位点
本发明提供了一种基于深度学习模型构建可预测筛选高可信的多聚腺苷酸化位点的方法和系统。本发明的用于预测或确定高可信的多聚腺苷酸化位点的模型(方法或系统)被称为“DeepPASS模型”。
典型地，在本发明中，通过构建深度学习模型DeepPASS提取多聚腺苷酸化的生物学特征，进而筛选出高置信度的来源于多聚腺苷酸化位点。
在本发明中，DeepPASS模型构建包括数据收集、数据预处理、模型搭建和模型训练评估等四个部分。
在本发明的一个优选例中，DeepPASS模型的构建包括以下步骤(参见图1A、1B和1C)：
(1)对于DeepPASS模型训练数据，本发明收集了四个不同来源的多聚腺苷酸化位点数据库：PolyA_DB3、PolyA-seq、Poly(A)Site2.0和GENCODE。筛选三个已发表数据库(PolyA_DB3、PolyA-seq和Poly(A)Site2.0)中高可信的多聚腺苷酸化位点，筛选条件为：一数据库中的多聚腺苷酸化位点在其基因组位置上下游12bp内与其它两个数据库中都有重合。将这三个数据库中的高可信位点，与GENCODE合并作为已知多聚腺苷酸化位点库。该多聚腺苷酸化位点库如图1A所示。
(2)根据上述步骤(1)的多聚腺苷酸化位点库，本发明采用基于正态分布随机偏移的方法进行数据扩增。以多聚腺苷酸化原始位点作为坐标零点，用正态分布(μ＝0,σ＝10)对位点进行随机偏移，每个多聚腺苷酸化位点约进行10次偏移(十倍扩增位点数量)，提取扩增后位点的上游100bp到下游100bp的基因组序列作为DeepPASS模型训练的正集。该步骤的示意图可参见图1B。
(3)将四个数据库PolyA_DB3、PolyA-seq、Poly(A)Site2.0和GENCODE的所有多聚腺苷酸化位点(经步骤(2)用基于正态分布随机偏移的方法进行数据扩增)取并集，作为正集(约1907290条序列)，另外，提取人基因组的基因间区中不含多聚腺苷酸化位点的序列(长度为200bp)作为DeepPASS模型训练的负集(约1907000条序列)。
优选地，在本发明中，正集的序列数量N1和负集的序列数量N2各自独立地≥50万，较佳地≥100万，更佳地≥200万。此外，正集的序列数量N1与负集的序列数量N2之比没有特别限制，通常为1:10至10:1，较佳地1:5至5:1，更佳地1:2至2:1，如约1:1。
(4)对于DeepPASS模型训练构建，本发明主要采用卷积神经网络(CNN)模块和循环神经网络(RNN)模块作为模型架构。CNN包括两个卷积层和最大池化层(max pooling layer)。第一个卷积层包括12-mer宽，4通道，共128filter；第二个卷积层包括6-mer，64filter可以覆盖第一个卷积层输出的128通道。卷积层输出结果通过最大池化层(Max pooling layer)进行降维压缩，选择4×4邻域内的最大值输出到下游RNN层。RNN采用双向长短期记忆网络结构(BiLSTM)，设置128个单元处理CNN模块输出内容，并将结果传递到下游全连接层。全连接层设置1024个隐藏单元(hidden units)和30％的丢码率(dropout ratio)，并通过归一化指数函数(softmax activation function)对最终二分类的分类概率进行输出。DeepPASS模型选取二分类中概率较大的分类作为预测结果进行输出，即通过1或0的分类结果预测给定的输入序列中是否包含真实的多聚腺苷酸化位点。
(5)将上述步骤(2)-(3)得到的正集和负集序列作为DeepPASS模型训练数据，采用上述步骤(4)的模型架构进行训练。模型训练设置数据量批次大小为5000， 模型训练100期，训练结果在10期内无提升则提前终止模型训练。模型训练结果用接受者操作特征曲线(ROC)和曲线下面积(AUC)作为评估。
对于本发明的DeepPASS模型，可以进一步评估其效果。在本发明中，本发明人将DeepPASS与目前发表的性能较好的同类多聚腺苷酸化位点预测工具(DeepPASTA和APARENT)进行比较。即根据上述步骤(1)得到的高可信多聚腺苷酸化位点(GENCODE的49942位点)，提取出原始位点上游100bp和下游100bp的序列作为正测试集。根据上述步骤(3)得到的不含多聚腺苷酸化位点的基因组间区200bp的序列作为负测试集。将正测试集和负测试集合并后随机取样约50000条序列用于DeepPASS、DeepPASTA和APARENT的模型评估，并用接受者操作特征曲线(ROC)和曲线下面积(AUC)评估不同模型结果。该评价结果如图1D所示。结果表明，本发明DeepPASS模型的性能远远优于DeepPASTA和APARENT。
基于单细胞转录组测序数据的多聚腺苷酸化测序读段信号峰的分析方法
本发明还提供了一种基于单细胞转录组测序数据的多聚腺苷酸化测序读段信号峰的分析方法，它包括：基于单细胞转录组测序数据的所述信号峰(peak)进行转录本水平的peak识别，peak坐标和剪接注释，peak质量检验和去冗余，及多测序样品peak合并。
在一优选例中，本发明的分析方法包括：
(Y1)将单细胞转录组测序数据与基因组序列进行比对，得到BAM格式存储的测序读段比对结果文件，BAM文件中应包含：BAM格式要求的测序读段比对结果信息、细胞条形码(CB)信息；
(Y2)将来自GENCODE的基因水平注释文件，根据其基因和转录本信息，转换成“基因-转录本”标记的基因转录本水平注释文件；
(Y3)根据测序读段比对结果文件和基因转录本水平注释文件，利用HOMERfindPeaks在转录本水平识别测序读段分布peak，识别每个转录本对应的所有peak，并基于对应转录本的坐标信息标记peak的区间位置；
(Y4)结合上述步骤(Y2)-(Y3)得到的基因转录本水平注释文件和基于转录本坐标的peak区间位置，将转录本位置系统映射到基因组坐标位置系统，计算出对应转录本识别的peak基因组坐标系统注释文件；
(Y5)结合上述步骤(Y2)-(Y4)得到的基因转录本水平注释文件和基于基因组坐标系统的peak文件，通过识别基因转录本水平注释文件中的转录本外显子剪接位点，进一步计算出对应转录本的peak剪接位点信息；
(Y6)根据上述步骤(Y5)中得到的peak注释结果，对识别的peak的测序读 段信号分布曲线进行正态分布检验(normality test)，去除可能源于异常扩增、多重比对或信号丰度较低的peak；通时进行单峰检验(unimodality test)，保留来自邻近多聚腺苷酸化信号形成的相互重叠的peak；
(Y7)上述步骤(Y2)-(Y6)计算可得到来自于同一基因的不同转录本来源的peaks，根据其基因信息可对peak进行去重。具体来说，对于区间重合比例超过50％的peak，通过对其测序读段信号丰度进行排序，保留重合peak中表达量最高的peak；
(Y8)对上述步骤(Y7)同一基因内去重后的peak重新进行测序读段丰度计算，过滤掉总体测序读段数量占比≤1％的peak，得到初步识别筛选的peak(raw peak)用于多聚腺苷酸化位点的评估；
(Y9)对于多个样品中每个样品经过上述步骤(Y3)-(Y8)得到的raw peak，本发明通过标记peak来源的样品信息，并将不同样品之间重合比例超过60％的peak标记为同一组，对同一组的peak中选取唯一值作为该多聚腺苷酸化位点的peak。
(Y10)将上述步骤(Y9)识别的raw peak注释区间作为待预测筛选的多聚腺苷酸化位点，取其注释末端上游100bp到下游100bp的基因组序列输入DeepPASS模型进行预测，可判断每一个raw peak是否为高可信多聚腺苷酸化的peak(PAS)，筛选出的高可信peak用于下游分析。
在单细胞水平对转录本进行定量分析和/或进行单细胞分型分析的方法
本发明还提供了一种单细胞水平对转录本进行定量分析和/或进行单细胞分型分析的方法。
典型地，该方法包括利用基于深度学习模型DeepPASS对识别的raw peak进行筛选，选取高置信度的多聚腺苷酸化的peak。然后，将计算筛选得到的peak可用于鉴定单细胞转录组测序数据中选择性表达的转录本，并对其进行单细胞水平的定量和单细胞分型分析。
在本发明一个优选例中，所述的单细胞水平对转录本进行定量分析和/或进行单细胞分型分析的方法包括：
(W1)对于多样品的peak整合分析，整合方法为：(1)源于同一多聚腺苷酸化位点的不同样品的peak，计算该区间内每个peak末端区间(上游50bp到下游25bp)覆盖的已知多聚腺苷酸位点的数量，根据覆盖位点数量的数值对peak进行递减排序；(2)提取peak末端上游100bp到下游100bp的区间序列，输入DeepPASS模型进行预测，根据预测结果从1到0排序；(3)在同一组peak中，选取综合排序最优的peak作为多样品分析中该多聚腺苷酸化位点的唯一peak。
(W2)根据上述步骤(W1)筛选的peak，通过peak的剪接位点和区间信息将其映射到对应来源的转录本上，可鉴定单细胞转录组测序样品中由可变多聚腺苷酸化(APA)产生的转录本，实现转录本表达鉴定。
(W3)将映射转录本信息的PAS注释构建GTF格式注释文件，通过featureCounts将比对过的BAM中reads重新计算PAS注释的分配情况，得到标记PAS注释的re-assigned BAM文件。利用UMI-tools可对re-assigned BAM文件中每一个细胞条形码对应的所有PAS注释的唯一分子标识(UMI)信息进行统计，最终获得单细胞水平的转录本表达矩阵，矩阵包括每个转录本在每个细胞的唯一分子标识数目，可用于单细胞转录组的分析。
(W4)对于转录本水平的单细胞分析，将上述步骤(W3)得到基于peak注释转录本水平表达矩阵作为单细胞分析工具Seurat的输入进行单细胞分析。本发明中建立的基于peak注释转录本水平的单细胞分析体系包含如下主要部分：表达质控，特征提取，降维分析，无监督细胞聚类，转录本分子标记物表达分析。
(W5)对于输入的转录本表达矩阵进行表达质控，保留至少表达细胞数占比≥0.5％的转录本。利用Seurat的内置函数对筛选后的表达矩阵进行转录本表达的标准化和表达变异计算，将表达变异最高的2000个高变异转录本作为特征转录本用于下游降维分析。
(W6)上述步骤(W5)得到的特征转录本可通过Seurat的内置函数进行降维分析和无监督细胞聚类。利用主成分分析方法(PCA)进行降维，并保留前50个主成分结果；根据该50个主成分对总体单细胞构建共享最近邻(SNN)相似度图，以无监督方式进行特定分辨率(0.3)进行细胞集群划分，最终得到基于PAS注释转录本表达的细胞分类结果。
(W7)根据步骤(W2)和(W3)中得到peak注释信息，计算筛选给定分子标记物基因表达的转录本，利用Seurat的内置函数可对不同细胞分类的分子标记物转录本的表达量进行计算和可视化，鉴定无监督细胞分类结果中的细胞类型。
预测系统和分析系统
本发明还提供了一种基于本发明方法的多聚腺苷酸化位点的预测系统，所述的系统包括：
输入单元，所述输入单元被配置为输入数据，所述数据包括待预测的输入序列，其中，需要判断在所述待预测的输入序列中是否存在多聚腺苷酸化位点；
多聚腺苷酸化位点的预测单元，所述的预测单元被配置为执行一多聚腺苷酸化位点的预测模型，从而获得所述待预测的输入序列中是否存在多聚腺苷酸化位点的预测结果；其中，所述的预测模型是用本发明第一方面所述的方法构 建的；
输出单元，所述的输出单元被配置为输出所述多聚腺苷酸化位点的预测单元的预测结果。
本发明还提供了一种单细胞转录组测序数据分析系统(SCAPTURE系统)，所述的系统包括：
(a)输入单元，所述输入单元被配置为输入序列数据，所述的输入序列数据选自下组：(i)DeepPASS预测模型和/或待预测的输入序列，其中，需要判断在所述待预测的输入序列中是否存在多聚腺苷酸化位点；(ii)单细胞转录组测序数据；(iii)所述(i)和(ii)的组合；
(b)分析模块，所述分析模块被配置为对输入的数据进行预定的分析，从而获得分析结果，并且所述分析模块选自下组：
(M1)多聚腺苷酸化位点的预测单元；
(M2)多聚腺苷酸化测序读段信号峰(peak)的分析单元，；
(M3)单细胞水平转录本定量分析和/或单细胞分型分析单元；
(M4)多聚腺苷酸化差异变化的分析单元；
(Mx)上述(M1)、(M2)、(M3)和(M4)的组合；
(c)输出单元，所述的输出单元被配置为输出所述分析模块的分析结果。
在另一优选例中，所述的输出系统包括显示屏、pad、打印机、手机等。
本发明的主要优点包括：
(1)根据发明内容1的步骤(2)，在多聚腺苷酸化位点预测模型DeepPASS的训练数据预处理时，采用了基于正态分布随机扰动的数据偏移和数据扩增策略；数据偏移可使DeepPASS达到非位置依赖的聚腺苷酸化位点预测效果；数据扩增可避免DeepPASS训练过拟合，同时提高训练和预测效果。
(2)根据发明内容2步骤(Y2)-(Y4)，在多聚腺苷酸化测序读段peak识别中，本发明整合了给定注释中的基因和转录本信息，采用了转录本水平的测序读段信号分布进行峰值计算，可提高peak识别的灵敏度。
(3)根据发明内容2步骤(Y6)，在多聚腺苷酸化测序读段peak的筛选中，本发明采用了正态分布检验多和单峰检验筛选条件，去除peak识别过程中噪音信号。
(4)根据发明内容2步骤(Y9)，本发明在多样品测序读段peak的整合过程中，采取了已知多聚腺苷酸位点的覆盖数量和DeepPASS预测的结果排序的策略，选取置信度最高的多聚腺苷酸位点peak用于多样品中的一致性分析。
(5)根据发明内容2步骤(Y10)，本发明在单细胞转录组测序数据的多聚腺 苷酸化测序读段peak的筛选中引入了非位置依赖的深度学习模型预测，可针对性筛选出高置信度的多聚腺苷酸化的peak。
(6)根据发明内容2步骤(Y4)-(Y5)和内容3步骤(W2)，本发明利用测序读段中的读段跨外显子比对位置，结合转录本注释的剪接位点位置和基因末端位置，对识别测序读段peak所来源的转录本进行计算推断，实现选择性表达转录本的鉴定。
(7)根据发明内容3步骤(W4)-(W7)，本发明可以在单细胞水平对转录本进行定量，并将转录本水平的定量应用于单细胞分型分析系统，替代了传统的基因水平定量分析。
(8)根据发明内容3步骤(W8)，本发明可以在单细胞水平，进行不同类型不同组别的细胞间差异化比较和统计，计算筛选得到单细胞水平全转录组多聚腺苷酸化趋势以及多聚腺苷酸化差异变化的基因，从而获得单细胞水平多聚腺苷酸化差异变化的分析结果。
下面结合具体实施例，进一步阐述本发明。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围。下列实施例中未注明具体条件的实验方法，通常按照常规条件，例如Sambrook等人，分子克隆：实验室手册(New York:Cold Spring Harbor Laboratory Press,1989)中所述的条件，或按照制造厂商所建议的条件。除非另外说明，否则百分比和份数是重量百分比和重量份数。
实施例1
多聚腺苷酸化位点预测模型DeepPASS的构建
在本实施例中，包括对多聚腺苷酸化位点预测模型DeepPASS的数据收集、预处理、模型搭建和训练评估。
材料：收集已发表的多聚腺苷酸位点数据库,从威斯塔研究所网站(The Wistar Institut)的PolyA_DB3(https://exon.apps.wistar.org/PolyA_DB/v3)下载多聚腺苷酸位点注释共290051条；从美国生物技术信息中心(NCBI，National Center for Biotechonlogy Information)网站下载PolyA-seq(https://www.ncbi.nlm.nih.gov/geo)多聚腺苷酸位点注释共514262条(登录号GSE30198)；从巴塞尔大学网站(University of Basel)的Poly(A)Site2.0(https://polyasite.unibas.ch/)下载多聚腺苷酸位点注释共569005条；从GENCODE多聚腺苷酸位点注释(ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_36/gencode.v36.polyAs.gtf.gz)多聚腺苷酸位点注释共49,942条；从GENCODE(https://www.gencodegenes.org/human/)下载人基因组和基因注释文件，基因注释共228048条；从GitHub分布式版 本控制软件网站下载安装DeepPASTA预测工具(https://github.com/arefeen/Deep PASTA)和APARENT预测工具(https://github.com/johli/aparent)。
步骤01：开源跨包管理与环境管理系统Conda(
https://docs.conda.io/en/latest/)下载安装bedtools注释文件分析工具(版本号：2.28.0)，编程语言Python应用程序(版本号：3.7)，TensorFlow深度学习框架(版本号：2.0.0)，编程语言R应用程序(版本号：4.0.0)
步骤02：将来自于PolyA_DB3、PolyA-seq和Poly(A)Site2.0的多聚腺苷酸化位点数据库文件转换成BED格式文件，并在基因组坐标区间向上游和下游各拓展6bp距离
步骤03：将拓展位置后的多聚腺苷酸化位点用bedtools intersect计算不同位点之间是否有基因组区间交集，将有交集的多聚腺苷酸化位点标记为同一编号。
步骤04：将计算同一编号的多聚腺苷酸化位点的数据库数量进行统计，保留三个数据库中共有的高可信位点，再与GENCODE的多聚腺苷酸化位点合并作为高可信多聚腺苷酸化位点库，共251,072条。
步骤05：用R语言的rnorm函数(μ＝0,σ＝10)产生正态分布的数据点，将高可信多聚腺苷酸化位点的基因组坐标数值与正态分布的数据点相加，产生正态分布偏移扩增后的多聚腺苷酸化位点。
步骤06：用bedtools getfasta获取正态分布偏移扩增后的多聚腺苷酸化位点的上游100bp到下游100bp长度共200bp的基因组序列，将这些序列作为模型训练的正集。
步骤07：根据基因注释计算出基因组没有基因注释位点的基因间区，用PolyA_DB3、PolyA-seq、Poly(A)Site2.0和GENCODE的所有多聚腺苷酸化位点进行反向筛选，得到既不含基因注释也不含多聚腺苷酸化位点注释的基因间区，用bedtools getfasta获取这些基因间区内随机区域的200bp的基因组序列，将这些序列作为模型训练的负集。
步骤08：用TensorFlow搭建DeepPASS深度学习模型框架，基本结构为：(1)卷积神经网络层包括两个卷积层和最大池化层(max pooling layer)。第一个卷积层包括12-mer宽，4通道，共128filter；第二个卷积层包括6-mer，64filter可以覆盖第一个卷积层输出的128通道。卷积层输出结果通过最大池化层(Max pooling layer)进行降维压缩，选择4×4邻域内的最大值输出到下游RNN层。(2)循环神经网络层采用双向长短期记忆网络结构(BiLSTM)，设置128个单元处理卷积神经网络模块输出内容，并将结果传递到下游全连接层。全连接层设置1024个隐藏单元(hidden units)和30％的丢码率(dropout ratio)，并通过归一化指数函数 (softmax activation function)对最终二分类的分类概率进行输出。整体结构可总结如下：
步骤09：利用步骤06和07得到的正集和负集序列训练DeepPASS。DeepPASS模型训练设置数据量批次大小为5000，模型训练100期，训练结果在10期内无提升则提前终止模型训练。
步骤10：将DeepPASS与目前发表的同类多聚腺苷酸化位点预测工具(DeepPASTA和APARENT)进行比较。将步骤06和07得到的正集和负集序列合并后随机取样50000条序列作为测试集用于DeepPASS、DeepPASTA和APARENT的模型评估。
步骤11：用测试集输入DeepPASS进行预测，计算测试集序列的原始值和预测概率，计算不同概率阈值下的真阳性率和假阳性率。
步骤12：根据DeepPASTA软件要求格式化测试集文件格式，输入DeepPASTA进行预测，计算预测概率的最大值和最小值差值，根据差值区间连续递增地取不同概率阈值，计算每个阈值下的真阳性率和假阳性率。
步骤13：根据APATENT软件要求格式化测试集文件，输入APATENT进行预测，选取单条序列返回值中的最大值作为该测试序列的预测结果，计算预测分值的最大值和最小值差值，根据差值区间连续递增地取不同分值阈值，计算每个阈值下的真阳性率和假阳性率。
步骤14：用接受者操作特征曲线(ROC)和曲线下面积(AUC)评估步骤11-13中DeepPASS、DeepPASTA和APATENT的预测表现。
步骤15：用GENCODE的多聚腺苷酸化位点进行基因组坐标偏移，分别偏移上游20bp、上游15bp、上游10bp、上游5bp、0bp、下游5bp、下游10bp、下游15bp、下游20bp，构建九个不同的数据集。
步骤16：用步骤15构建的九个数据集对DeepPASS和DeepPASTA进行测试，计算在不同位置偏移条件下DeepPASS和DeepPASTA对多聚腺苷酸化位点预测的准确率和稳定性。
结果
如图1A和1B所示，综合四个多聚腺苷酸位点数据库共筛选出251072个高可信位点，经正态分布随机偏移和扩增后提取得到1907294条正集序列，从基因间区中提取数目接近的负集序列，构建出数据量约四百万条序列的训练集用于DeepPASS模型训练(图1C)。
DeepPASS模型训练结果如图1D所示，在测试集的预测准确度上达到98％。 比较DeepPASS、DeepPASTA和APARENT的总体预测表现，计算得到DeepPASS预测的曲线下面积为0.99，DeepPASTA和APARENT分别为0.92和0.72。
该结果表明，DeepPASS在多聚腺苷酸化位点预测表现优于现有的DeepPASTA和APARENT。
为了验证DeepPASS模型的非位置依赖性，本发明还基于GENCODE数据库构建偏移上游20bp、上游15bp、上游10bp、上游5bp、0bp、下游5bp、下游10bp、下游15bp、下游20bp的九个数据集测试，并与DeepPASTA进行比较测试。
结果如图1E和1F所示，DeepPASS在九个数据集上的表现准确率分别为92.6％、95.5％、96.9％、97.5％、97.6％、97.2％、96.5％、94.8％、91.8％，DeepPASTA在五个数据集上的表现准确率分别为21.3％、29.2％、43.4％、71.4％、81.1％、59.1％、29.0％、20.4％、12.2％，这表明：DeepPASS对多聚腺苷酸化位点的预测具有显著的非位置依赖性和准确率更高的优点。
以上结果表明，本发明构建的深度学习模型DeepPASS在应用于多聚腺苷酸化位点的预测上具有非位置依赖性和高准确率的优点。运用DeepPASS模型，可帮助从单细胞转录组测序数据识别的多聚腺苷酸化位点中筛选出真正可靠的位点，降低由于测序建库技术和测序数据处理引入的假阳性，推动开展单细胞水平的多聚腺苷酸化生物学研究。
实施例2
单细胞转录组数据的多聚腺苷酸化位点测序读段peak的识别和筛选
材料：从10x Genomics公司网站(https://support.10x Genomics.com/single-ce ll-gene-expression/datasets)下载人PBMC单细胞转录组测序共六个数据集，依据每个数据检测到的细胞数量命名为3k、4k、5k、6k、8k和10k。这六个数据集测序读段数分别为179147840、379462522、383941607、143988649、784064148和638901019；从GENCODE(https://www.gencodegenes.org/human/)下载人基因组和基因注释文件，基因注释共228048条；
步骤01：从10x Genomics公司网站下载安装Cell Ranger单细胞数据预处理软件(https://support.10xGenomics.com/single-cell-gene-expression/software/pipelines/latest/installation)(版本号：4.0)；从圣克鲁斯加利福尼亚大学(University of California,Santa Cruz)下载注释格式处理套件UCSC Utilities(http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/)；从圣迭戈加利福尼亚大学(University of California San Diego)下载HOMER识别peak工具(http://homer.ucsd.edu/homer /ngs/peaks.html)(版本号：4.11.1)；从GitHub下载BAM文件处理工具samtools(https://github.com/samtools/samtools)(版本号：1.9)；从GitHub下载单细胞测序工具Sierra(https://github.com/VCCRI/Sierra)(版本号：0.99.24)；从GitHub下载单细胞测序工具scAPA(https://github.com/ElkonLab/scAPA)(版本号：0.1.0)；
步骤02：用10x Genomics公司发布的单细胞转录组数据分析工具Cell Ranger对GENCODE的基因组和基因注释文件建立索引文件。
步骤03：用Cell Ranger将6个PBMC单细胞转录本测序数据比对到索引文件，获得每个样品的测序读段比对结果BAM文件和细胞×基因表达的表达量矩阵文件。
步骤04：将来自GENCODE的基因水平注释文件，根据其基因和转录本信息，转换成“基因-转录本”标记的基因转录本水平注释文件，用UCSC Utilities工具分别构建genepred、GTF和BED格式文件。
步骤05：根据基因注释文件提取每一个基因的基因组坐标范围和链方向信息，用bedtools和samtools从BAM文件中提取每一个基因对应基因组范围内的测序读段比对结果。
步骤06：根据构建的转录本水平注释文件和提取的每个基因范围内的测序读段比对结果，对于同一个基因的每条转录本，用HOMER工具中的makeTagDirectory和findPeaks识别该转录本测序读段富集的峰值区域，即测序读段peak。
步骤07：对计算得到的基于转录本水平坐标的peak区间，将转录本位置系统映射到基因组坐标位置系统，计算出对应转录本识别的peak基因组坐标系统注释文件。
步骤08：对转换坐标后的peak注释，根据其基因和转录本水平注释信息，识别转录本注释中的外显子剪接位点，计算出对应转录本的peak剪接位点信息。
步骤09：对测序样品中的每一个表达的基因执行步骤05-08，计算得到该样品全基因组范围内的测序读段peak注释。
步骤10：对获得的peak进行初步筛选，去除噪音和低丰度的peak。筛选条件包括：在R中用nortest对peak的测序读段分布曲线做正态性检验，去掉不符合正态分布的peak；在R中用diptest对peak的测序读段分布曲线做单峰检验，去掉不符合多峰分布且不符合正态分布的peak。
步骤11：对于同一个多聚腺苷酸化位点产生的测序读段分布峰，在基因注释中覆盖该加尾位点的多个转录本会计算出多个peak，对这些peak的去冗余方法为：通过计算peak区间重合比例，将重合比例大于50％的peak进行归类分组，对于同一组重合的所有peak仅保留测序读段覆盖深度最高的peak注释。
步骤12：同一基因内去重后的peaks重新进行测序读段丰度计算，去除读段数量占比小于同一基因上总读段数量1％的peak。
步骤13：对测序样品中的找到的所有peak，执行步骤10-12，得到全基因组水平初步识别和筛选的peak用于多聚腺苷酸化位点的评估。
步骤14：对于已有的工具Sierra，用相同的六个PBMC单细胞测序数据集和基因注释进行分析，根据Sierra的使用说明分析得到结果的peak注释文件。
步骤15：对于已有的工具scAPA，用相同的六个PBMC单细胞测序数据集，以及scAPA内建的基因注释进行分析，根据scAPA的使用说明分析得到结果的peak注释文件。
步骤16：统计六个PBMC单细胞测序数据集的SCAPTURE、Sierra和scAPA分析得到的peak注释，分别统计每个样品中识别的内含子区和外显子区的peak数量，并进行工具间的比较。
步骤17：将本发明构建的DeepPASS深度学习模型用于测序读段peak筛选，具体做法为：将测序读段peak注释区间的末端位置作为潜在的多聚腺苷酸化位点，取其上游100bp到下游100bp的基因组序列输入DeepPASS模型进行预测，可判断每一个测序读段peak是否为高可信多聚腺苷酸化的peak，筛选出的高可信多聚腺苷酸化peak。
结果
来源10x Genomics公司的人PBMC的六个单细胞测序数据集有良好的比对效果，基因组比对效率分别为：92.2％、89.8％、90.2％、92.0％、88.8％和89.0％。根据其基因组比对结果，利用本发明建立的SCAPTURE系统，能够利用单细胞测序数据集识别和筛选多聚腺苷酸化测序读段peak(如图2A所示)。
参见图2A，以NADK基因为例，展示了SCAPTURE计算系统在人PBMC 10k的单细胞转录组测序数据集中多聚腺苷酸化测序读段peak的识别和筛选步骤。根据本发明的计算步骤(如图2A所示)，SCAPTURE先提取NADK基因范围内的测序读段比对结果，结合构建的转录本水平的注释信息，对每一条转录本进行测序读段peak的识别，随后对NADK内的不同转录本来源的peak进行去冗余和表达质控，鉴定出NADK的多聚腺苷酸化测序读段peak信息，结果为Peak1和Peak2(参见图2A)。
如图2B所示，本发明将SCAPTURE与已有的Sierra和scAPA流程进行比较。对于六个PBMC单细胞测序样品的外显子区的测序读段peak位点，SCAPTURE能够分别鉴定出32569、45135、50970、33578、48906和54236个peak位点；Sierra鉴定结果分别为11684、18635、19956、10025、24106和24590个peak位点；scAPA 鉴定结果分别为11684、18635、19956、10025、24106和24590个peak位点。对于内含子区的外显子区的测序读段peak位点，SCAPTURE能够分别鉴定出54303、75816、92108、57772、82056和98328个peak位点；Sierra鉴定结果分别为10392、36095、49495、8659、59925和86912个peak位点；scAPA鉴定结果分别为7067、16947、19747、5631、18987和25423个peak位点。该结果表明，对于同样的单细胞测序样品，SCAPTURE能够识别出更多外显子区和内含子的测序读段peak位点，识别测序读段peak的灵敏度明显优于现有的工具。
如图2C所示，本发明利用构建的DeepPASS模型筛选高可信的测序读段peak，降低多聚腺苷酸化位点识别的假阳性。经过DeepPASS模型筛选，外显子peak和内含子peak中，经典多聚腺苷酸加尾信号占比的富集倍数分别为：外显子peak富集6.2～7.0倍，内含子peak富集2.9～7.4倍。
结果表明，本发明构建的DeepPASS模型能够用于单细胞转录组测序的多聚腺苷酸化peak置信度筛选，提高单细胞水平的多聚腺苷酸化研究的准确性。
如图2C所示，利用构建的DeepPASS模型筛选去除假阳性多聚腺苷酸化位点后，本发明构建的DeepPASS保留外显子和内含子peak数量为25247和12286，显著多于Sierra(14039和9830)和scAPA(8378和2429)。
综合以上结果，相比于现有的单细胞转录组多聚腺苷酸化研究工具，SCAPTURE具有显著更高的灵敏度和准确性。
实施例3
基于多聚腺苷酸化位点鉴定转录本表达和应用于单细胞分型分析。
材料：实施例2使用的从10x Genomics公司网站下载人PBMC单细胞转录组测序共六个数据集；实施例2中使用的从GENCODE下载的基因注释。
步骤01：从GitHub下载安装Seurat单细胞分析工具(https://github.com/satijalab/seurat)(版本号：3.2.2)；从SourceForge(https://sourceforge.net/)下载featureCounts测序读段计数工具(http://subread.sourceforge.net/)(版本号：1.6.2)；从GitHub下载分子标识计数软件UMI-tools(https://github.com/CGATOxford/UMI-tools)(版本号：1.0.1)。
步骤02：根据实施例2得到的peak注释信息和GENCODE转录本注释信息，对peak来源的特定转录本进行鉴定，鉴定方法为：转录本外显子剪接位点需与peak剪接位点信息完全匹配，且选取peak距离注释末端最近的转录本作为peak来源的转录本，最终得到单细胞转录组测序样品中基因选择性表达的转录本。
步骤03：将步骤02映射转录本信息后的peak注释进行格式转化，用UCSC Utilities工具构建GTF格式的映射转录本信息peak注释文件，用于下游的单细胞 转录本定量分析。
步骤04：用Cell Ranger将6个PBMC单细胞转录本测序数据比对到索引文件，获得每个样品的测序读段比对结果BAM文件和细胞×基因表达的表达量矩阵文件。
步骤05：将六个样品的细胞×基因表达的表达量矩阵输入Seurat进行细胞质量控制，细胞筛选标准为：线粒体基因组唯一分析标识占比≤20％、细胞唯一分析标识数量≥800以及细胞检测到的基因数≥500，保留六个样品筛选过后的细胞条形码。
步骤06：将映射转录本信息的peak注释文件，通过featureCounts将比对过后的六个样品的BAM文件中的测序读段进行重新分配，得到标记转录本信息peak注释的re-assigned BAM文件。
步骤07：利用UMI-tools可对re-assigned BAM文件中筛选过后的细胞条形码的转录本表达进行定量，计算得到细胞×转录本表达矩阵。具体做法为：将分配到每一个转录本peak注释的测序读段上唯一分子标识信息进行统计，计算每个转录本peak在每个细胞的唯一分子标识数目，该数值可作为转录本在细胞中的表达水平。
步骤08：利用步骤05和步骤08中的分别得到的基因水平表达矩阵和转录本水平表达矩阵进行单细胞转录组分析，具体包括质控、标准化、降维分析、无监督分类和细胞类型鉴定。
步骤09：对单细胞表达数据进行质控，去掉转录本表达的细胞数≤0.5％的转录本，并对基因表达量和转录本表达量进行标准化。
步骤10：对基因和转录本计算它们的细胞间表达变异度，选择排序前2000个高变异特征用于下游降维分析。
步骤11：通过Seurat的内置函数RunPCA对基因和转录本的高变异特征进行主成分分析方法降维，并保留前50个主成分结果。
步骤12：通过Seurat的内置函数FindNeighbors分别利用基因和转录本分析得到的50个主成分，构建总体单细胞构建共享最近邻相似度图。
步骤13：通过Seurat的内置函数FindClusters对基因和转录本分析分别以无监督方式进行分辨率为0.3的细胞集群划分，分别得到细胞无监督分类结果。
步骤14：寻找文献报道的人PBMC中不同细胞类型的标记基因，分别对基因和转录本表达水平鉴定细胞集群进行标记基因表达分析，通过不同细胞类型的标记基因表达情况将细胞集群注释为对应细胞类型。
结果：
如图3A和3B所示，本发明的计算流程可以基于多聚腺苷酸化测序读段peak信息鉴定基因的选择性转录本表达。以RASGRP2为例，通过Peak1和Peak2可以鉴定该基因表达转录本1(ENST00000377486.7)和转录本2(ENST00000354024.7)，且转录本1和转录本2在单细胞水平能够表征基因表达量。该结果表明，本发明能够基于多聚腺苷酸化测序读段peak进行转录本的鉴定和定量。
如图3C所示，对于筛选过后的多聚腺苷酸化测序读段peak，可以根据其剪接位点信息，鉴定其选择性表达的转录本。在六个PBMC单细胞测序样品中，传统方法得到的表达基因数为13626个，本发明构建的SCAPTURE系统得到的表达转录本数目为34316个，将单细胞分析从基因水平扩展提高到转录本水平。
如图3D、3E、3F、3G和3H所示，在单细胞的样品整合、无监督聚类分析中，基于转录本表达的分析效力与基因表达接近。参考图3G和3H，在基因水平被鉴定为无分类的细胞类群，其部分细胞能够基于转录本表达重新归类到已知的细胞类群中，如重分类的树突状细胞和重分类的自然杀伤细胞。该结果表明，本发明建立的转录本表达分析能够替代基因表达来进行常规的单细胞计算分析；此外，基于转录本表达的优势还在于，能够更准确的区分单细胞群体中的不同细胞类型。
结果表明，本发明建立的SCAPTURE计算系统能根据多聚腺苷酸化位点鉴定选择性表达的转录本，并进行单细胞水平的转录本定量和单细胞分型分析。
实施例4
不同条件下多种细胞类群的多聚腺苷酸化表达差异分析
材料：NCBI网站下载来自病毒感染研究的人PBMC单细胞转录组测序共九个公开数据集，其中五个正常人数据集，四个重症病人数据集(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi？acc＝GSE155673)。
步骤01：参考实施例2的步骤01进行软件及工具安装。
步骤02：参考实施例2的步骤02-13和17，对九个数据集进行处理，得到该样品的高可信多聚腺苷酸化peak注释文件。
步骤03：参考实施例3的步骤03-14，利用步骤02得到的注释文件对九个数据集进行单细胞水平的转录本定量和单细胞分型分析。
步骤04：利用步骤02得到的注释文件，选取含有2个以上多聚腺苷酸化位点的基因，并根据基因组注释坐标的远近，分为“近端”、“远端”和“中间”三种不同多聚腺苷酸化位点。
步骤05：建立近端位点使用率表达矩阵，用以评估每个基因的多聚腺苷酸 化位点使用倾向性。方法为：该基因的“近端”位点表达值，除以该基因其余位点表达值的加和，结果再取2为底的对数(公式如下)。对每个细胞的每个基因进行计算，即可得到近端位点使用率×细胞的数值矩阵。
步骤06：评估细胞整体的多聚腺苷酸趋势，方法为：利用近端位点使用率×细胞的数值矩阵，对每个细胞的所有基因的计算近端位点使用率的平均值，再根据正常人/重症病人进行分组，或根据免疫细胞类型进行分组，对不同组别的细胞的全基因组水平计算近端位点使用率进行累计分布统计，评估细胞整体的多聚腺苷酸化倾向。
步骤07：筛选多聚腺苷酸化调控存在差异的基因，方法为：将步骤05得到的近端位点使用率×细胞的数值矩阵构建Seurat对象，通过Seurat的内置函数FindMarkers对正常人和重症病人的细胞进行差异分析，得到不同条件下的基因的近端位点使用率差异。
结果：
如图4A所示，本发明的SCAPTURE计算流程可以基于多聚腺苷酸化位点的定量信息，来筛选不同细胞群体及不同条件下的多聚腺苷酸化差异调控。对下载的九个样品进行数据分析，首先分析得到单细胞的疾病分群以及细胞类型分群(图4B，左，中)。在细胞水平，利用SCAPTURE流程可评估每个细胞全基因组的多聚腺苷酸趋势，在单细胞水平观察到细胞间的多聚腺苷酸化选择性差异(图4B，右)。进一步对细胞类型分群比较，发现单核细胞倾向于使用远端多聚腺苷酸化位点，B淋巴细胞倾向于使用近端多聚腺苷酸化位点，揭示了不同免疫细胞类型的差异性(图4C)。更重要的是，对正常人及重症患者的比较分析发现，在感染病毒后免疫细胞整体呈现出使用近端多聚腺苷酸化位点的调控趋势(图4D)。在基因水平，利用SCAPTURE流程可鉴定不同条件不同细胞群体的多聚腺苷酸化差异基因。在单核细胞向树突细胞的分化过程中，发现PHB,HMGB1,CD47,HLA-DQB1,PSMF1,PSMD7等基因存在细胞类型相关的免疫功能差异(图4E)。在正常人和病毒感染重症患者的比较中，也发现与抗病毒生理功能直接相关的免疫球蛋白基因存在多聚腺苷酸化的显著调控(图4F)。以IGHM为例，该基因远端多聚腺苷酸化位点编码膜结合形式蛋白，而近端多聚腺苷酸化编码分泌型蛋白抗体，在病毒感染后B细胞分化出的浆细胞中，该基因的近端多聚腺苷酸化位点表达显著上调，表明浆细胞表达分泌型抗体与抗病毒免疫功能相适应(图4F)。
讨论
本发明提供了一种融合深度学习模型的多聚腺苷酸化位点(PAS)检测方法，用于检测单细胞测序数据的多聚腺苷酸化位点和转录本差异表达，及应用于单细胞分型分析和多聚腺苷酸化分析(即SCAPTURE系统)。
本发明建立的SCAPTURE计算系统其意义在于：构建的DeepPASS模型实现高准确性和非位置依赖的多聚腺苷酸化位点预测，该预测工具仅需要基因组序列作为输入，因此也可以广泛用于其它关于多聚腺苷酸化的研究体系和生物学问题；SCAPTURE可以高灵敏地从头鉴定全基因组水平的单细胞多聚腺苷酸化位点，克服了单细胞转录组测序存在的低丰度信号丢失、假阳性高等缺点；SCAPTURE还可以基于多聚腺苷酸化位点对基因选择性表达的不同转录本进行识别，将单细胞分析从基因水平提高到转录本水平，在单细胞分析中具有优势；SCAPTURE对单细胞转录组测序中多聚腺苷酸化位点的识别、筛选、定量和差异分析，还可以应用于组织特异性、细胞发育和疾病发生等体系中多聚腺苷酸化研究，使多聚腺苷酸化研究提高至单细胞水平。
在本发明中，可采用的已知多聚腺苷酸化位点收集自PolyA_DB3、PolyA-seq、Poly(A)Site2.0和GENCODE，这些多聚腺苷酸化位点也可以来自其它数据库或数据集。
在构建基因-转录本水平注释文件时，下载人的GENCODE注释信息，该注释信息也可为其它研究物种、数据库来源或手动建立的基因注释。
在本发明中，所述的3'tag-based单细胞转录组测序可以由10x Genomics公司的10x Chromium技术获得测序数据，也适用于其它同类的基于3'tag-based单细胞转录组测序原理的测序数据。
在本发明中，测序数据可以采用Cell Ranger或其它合适的比对流程进行单细胞转录组数据的测序读段基因组比对位置、细胞条形码和唯一分子标识的标记，获得结果储存于BAM文件中。
本发明中，在识别多聚腺苷酸化形成的测序读段peak时也可利用同类的基于高通量测序读段分布的peak识别工具。
在本发明中，在计算基于peak注释的转录本表达情况时，也可利用其它同类的能处理测序读序的细胞条形码和唯一分子标识的工具。
应理解，在本发明的上述实施例中是以人PBMC的单细胞测序数据进行说明。显然，本发明还可用于其它物种的3'tag-based单细胞转录组测序数据的PAS的鉴定和转录本的定量。
综上所述，通过本发明开发的SCAPTURE方法和系统，可以高效准确地识 别筛选多聚腺苷酸化位点，进一步可对转录本进行定量分析和/或进行单细胞分型分析，并且可以鉴定不同细胞群体及不同条件下的多聚腺苷酸化差异调控，从而细胞水平和基因水平揭示多聚腺苷酸化参与的免疫调控功能。
在本发明提及的所有文献都在本申请中引用作为参考，就如同每一篇文献被单独引用作为参考那样。此外应理解，在阅读了本发明的上述讲授内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。
Claims (12)
- 一种构建多聚腺苷酸化位点的预测模型的方法，其特征在于，所述方法包括步骤：(S1)提供用于模型训练的数据集，所述数据集包括正集和负集，其中，所述的正集中包含的序列为含多聚腺苷酸化位点的序列，即正序列；而所述的负集中包含的序列为不含多聚腺苷酸化位点的序列，即负序列；其中，所述的正序列包括：对于多聚腺苷酸化原始位点而言的基础正序列和基于所述基础正序列通过偏移扩增所获得的偏移正序列，并且所述的基础正序列和偏移正序列均含有所述的多聚腺苷酸化原始位点；(S2)将(S1)中所述的数据集用于模型架构进行训练；其中，所述的模型架构包括基于卷积神经网络(CNN)模块和循环神经网络(RNN)模块的模型架构；所述CNN模块包括两个卷积层和最大池化层(max pooling layer)，其中，所述卷积层的输出结果通过最大池化层(max pooling layer)进行降维压缩，并选择4×4邻域内的最大值输出到下游RNN模块；所述RNN模块设有双向长短期记忆网络结构(BiLSTM)，并将CNN模块的输出内容传递到下游全连接层；所述的全连接层设置多个隐藏单元(hidden units)和预定的丢码率(dropout ratio)，并通过归一化指数函数(softmax activation function)对最终二分类的分类结果进行输出；(S3)当训练结果达到预定终止条件，则终止模型训练，并获得所述的多聚腺苷酸化位点的预测模型(即DeepPASS模型)。
- 如权利要求1所述的方法，其特征在于，在步骤(S1)中，所述的偏移扩增包括：以多聚腺苷酸化原始基因组位点作为坐标零点，用正态分布对位点进行随机偏移，对一个多聚腺苷酸化位点进行4-20次偏移，提取偏移扩增后位点的上游100±20bp到下游100±20bp的基因组序列作为模型训练的正集中的偏移正序列。
- 一种多聚腺苷酸化位点的预测系统，其特征在于，所述的系统包括：输入单元，所述输入单元被配置为输入数据，所述数据包括待预测的输入序列，其中，需要判断在所述待预测的输入序列中是否存在多聚腺苷酸化位点；多聚腺苷酸化位点的预测单元，所述的预测单元被配置为执行一多聚腺苷酸化位点的预测模型，从而获得所述待预测的输入序列中是否存在多聚腺苷酸化位点的预测结果；其中，所述的预测模型是用权利要求1所述的方法构建的；输出单元，所述的输出单元被配置为输出所述多聚腺苷酸化位点的预测单元的预测结果。
- 一种基于单细胞转录组测序数据识别多聚腺苷酸化测序读段信号峰(peak)的分析方法，其特征在于，包括步骤：(Y1)将单细胞转录组测序数据与基因组序列进行比对，得到BAM格式存储的测序读段比对结果文件，BAM文件中应包含：BAM格式要求的测序读段比对结果信息、测序读段细胞条形码(CB)和测序读段分子标识(UMI)信息；(Y2)将基因水平注释文件，根据其基因和转录本信息，转换成“基因-转录本”标记的基因转录本水平注释文件；(Y3)根据测序读段比对结果文件和基因转录本水平注释文件，利用HOMER findPeaks在转录本水平识别测序读段分布peak，识别每个转录本对应的所有peak，并基于对应转录本的坐标信息标记peak的区间位置；(Y4)结合上述步骤(Y2)-(Y3)得到的基因转录本水平注释文件和基于转录本坐标的peak区间位置，将转录本位置系统映射到基因组坐标位置系统，计算出peak的基因组坐标系统注释文件；(Y5)结合上述步骤(Y2)-(Y4)得到的基因转录本水平注释文件和基于基因组坐标系统的peak文件，通过识别基因转录本水平注释文件中的转录本外显子剪接位点，进一步计算出对应转录本的peak剪接位点信息；(Y6)根据上述步骤(Y5)中得到的peak注释结果，对识别的peak的测序读段信号分布曲线进行正态分布检验(normality test)，去除可能源于异常扩增、多重比对或信号丰度较低的peak；并进行单峰检验(unimodality test)，保留来自邻近多聚腺苷酸化信号形成的相互重叠的peak；(Y7)上述步骤(Y2)-(Y6)计算可得到来自于同一基因的不同转录本来源的peaks，根据其基因信息可对peak进行去重；优选地，去重方法为：对于区间重合比例超过50％的peak，对其测序读段丰度进行排序，保留表达量最高的peak；(Y8)对上述步骤(Y7)去重后的peak重新计算测序读段丰度，过滤掉同一基因内测序读段数量占比≤1％的peak，得到初步筛选的peak(raw peak)；(Y9)对于多个样品中每个样品经过上述步骤(Y3)-(Y8)得到的raw peak，通过标记peak来源的样品信息，并将不同样品之间重合比例超过60％的peak标记为同一组，对同一组的peak中选取唯一值作为该多聚腺苷酸化位点的peak；(Y10)将上述步骤识别的raw peak注释区间的末端近似为多聚腺苷酸化位点，提取位点上游100bp到下游100bp序列并输入DeepPASS模型进行预测， 判断每一个raw peak是否为高可信多聚腺苷酸化的peak(PAS)，从而获得单细胞转录组测序数据的多聚腺苷酸化测序读段信号峰分析结果。
- 如权利要求4所述的分析方法，其特征在于，在步骤(Y10)中，所述的预测模型为用权利要求1所述的方法构建的预测模型，即DeepPASS模型。
- 一种在单细胞水平对转录本进行定量分析和/或进行单细胞分型分析的方法，其特征在于，所述方法包括：(Z1)提供待分析的单细胞样本的单细胞转录组测序数据；(Z2)用权利要求4所述的分析方法，对单细胞转录组测序数据进行分析，从而获得高可信peak；(Z3)基于上一步骤中获得的高可信peak，进行转录本进行定量分析和/或进行单细胞分型分析。
- 如权利要求6所述的方法，其特征在于，在步骤(Z3)中，包括：(W1)对于多样品的peak整合分析，整合方法为：将不同样品之间重合的peak标记为同一组，对同一组内的peak进行置信度排序，选取置信度最高的peak作为该多聚腺苷酸化位点的唯一peak；(W2)根据上述步骤(W1)筛选的peak，通过peak的剪接位点和基因组区间信息将其映射到对应来源的转录本上，可鉴定单细胞转录组测序样品中由可变多聚腺苷酸化(APA)产生的转录本，实现转录本表达鉴定；(W3)将映射转录本信息的PAS注释构建GTF格式注释文件，通过featureCounts将比对过的BAM中reads重新计算PAS注释的分配情况，得到标记PAS注释的re-assigned BAM文件；利用UMI-tools可对re-assigned BAM文件中每一个细胞条形码对应的所有PAS注释的唯一分子标识(UMI)信息进行统计，获得单细胞水平的转录本表达矩阵，矩阵包括每个转录本在每个细胞的唯一分子标识数目，可用于单细胞转录组的分析；(W4)对于转录本水平的单细胞分析，将上述步骤(W3)得到基于peak注释转录本水平表达矩阵作为单细胞分析工具Seurat的输入进行单细胞分析；优选地，基于peak注释转录本水平的单细胞分析体系包含如下主要部分：表达质控，特征提取，降维分析，无监督细胞聚类，转录本分子标记物表达分析；(W5)对于输入的转录本表达矩阵进行表达质控，保留至少表达细胞数占比≥0.5％的转录本；优选地，利用Seurat的内置函数对筛选后的表达矩阵进行转录本表达的标准化和表达变异计算，将表达变异最高的2000个高变异转录本作为特征转录本用于下游降维分析；(W6)上述步骤(W5)得到的特征转录本可通过Seurat的内置函数进行降维分析和无监督细胞聚类；优选地，利用主成分分析方法(PCA)进行降维，并保 留主要(如前50个)主成分结果；根据该50个主成分对总体单细胞构建共享最近邻(SNN)相似度图，以无监督方式进行特定分辨率(如0.3)进行细胞集群划分，最终得到基于PAS注释转录本表达的细胞分类结果；(W7)根据步骤(W2)和(W3)中得到peak注释信息，计算筛选给定分子标记物基因表达的转录本，利用Seurat的内置函数可对不同细胞分类的分子标记物转录本的表达量进行计算和可视化，鉴定无监督细胞分类结果中的细胞类型。
- 如权利要求7所述的方法，其特征在于，在步骤(W1)中，对于多样品的peak整合分析，整合方法为：(1)源于同一多聚腺苷酸化位点的不同样品的peak，计算该区间内每个peak末端区间(上游50bp到下游25bp)覆盖的已知多聚腺苷酸位点的数量，根据覆盖位点数量的数值对peak进行递减排序；(2)提取peak末端上游100bp到下游100bp的区间序列，输入DeepPASS模型进行预测，根据预测结果从1到0排序；(3)在同一组peak中，选取综合排序最优的peak作为多样品分析中该多聚腺苷酸化位点的唯一peak。
- 如权利要求7所述的方法，其特征在于，在步骤(W2)中，利用peak剪接位点和基因组注释信息，对符合peak来源的转录本进行鉴定；优选地，鉴定方法为：转录本外显子剪接位点需与peak剪接位点信息完全匹配，且选取末端注释距离最近的转录本作为peak来源的转录本，最终得到全基因组水平的单细胞转录组测序样品中的不同转录本表达。
- 一种单细胞转录组测序数据分析系统，其特征在于，所述的系统包括：(a)输入单元，所述输入单元被配置为输入序列数据，所述的输入序列数据选自下组：(i)DeepPASS预测模型和/或待预测的输入序列，其中，需要判断在所述待预测的输入序列中是否存在多聚腺苷酸化位点；(ii)单细胞转录组测序数据；(iii)所述(i)和(ii)的组合；(b)分析模块，所述分析模块被配置为对输入的数据进行预定的分析，从而获得分析结果，并且所述分析模块选自下组：(M1)多聚腺苷酸化位点的预测单元，所述的(M1)预测单元被配置为执行一多聚腺苷酸化位点的预测模型，从而获得所述待预测的输入序列中是否存在多聚腺苷酸化位点的预测结果；其中，所述的预测模型是用权利要求1所述的方法构建的；(M2)多聚腺苷酸化测序读段信号峰(peak)的分析单元，其中，所述的(M1)多聚腺苷酸化测序读段信号峰(peak)的分析单元被配置为执行权利要求4所述的分析方法，从而获得单细胞转录组测序数据的多聚腺苷酸化测序读段信号峰分析结果；(M3)单细胞水平转录本定量分析和/或单细胞分型分析单元，所述的 (M3)单细胞水平转录本定量分析和/或单细胞分型分析单元被配置为执行权利要求6所述的分析方法，从而获得单细胞水平转录本定量分析结果和/或单细胞分型分析结果；(Mx)上述(M1)、(M2)和(M3)的组合；(c)输出单元，所述的输出单元被配置为输出所述分析模块的分析结果。
- 如权利要求10所述的系统，其特征在于，所述的分析模块(b)还包括：(M4)多聚腺苷酸化差异变化的分析单元，所述的(M4)多聚腺苷酸化差异变化的分析单元被配置为执行以下步骤W8：根据单细胞分型和多聚腺苷酸化位点定量结果，进行不同类型不同组别的细胞间差异化比较和统计，计算筛选得到单细胞水平全转录组多聚腺苷酸化趋势以及多聚腺苷酸化差异变化的基因，从而获得单细胞水平多聚腺苷酸化差异变化的分析结果；(Mx)上述(M1)、(M2)、(M3)和(M4)的任意组合。
- 如权利要求11所述的系统，其特征在于，所述步骤W8包括以下子步骤：(L1)利用所述的注释文件，对多聚腺苷酸化位点进行分类；(L2)利用所述的注释文件和定量结果，计算位点相对使用率表达矩阵；(L3)利用(L2)得到的位点相对使用率，计算细胞整体的多聚腺苷酸趋势；和(L4)利用(L2)得到的位点相对使用率，筛选多聚腺苷酸化存在差异的基因。
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110251298.3A CN115050416A (zh) | 2021-03-08 | 2021-03-08 | 融合深度学习模型的单细胞转录组计算分析方法和系统 |
| CN202110251298.3 | 2021-03-08 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| WO2022188785A1 true WO2022188785A1 (zh) | 2022-09-15 |
Family
ID=83156113
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| PCT/CN2022/079788 WO2022188785A1 (zh) | 2021-03-08 | 2022-03-08 | 融合深度学习模型的单细胞转录组计算分析方法和系统 |
Country Status (2)
| Country | Link |
|---|---|
| CN (1) | CN115050416A (zh) |
| WO (1) | WO2022188785A1 (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115579060A (zh) * | 2022-12-08 | 2023-01-06 | 国家超级计算天津中心 | 基因位点检测方法、装置、设备及介质 |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116072217B (zh) * | 2022-11-02 | 2023-07-25 | 杭州联川基因诊断技术有限公司 | 一种单细胞转录组数据可用性处理方法、介质及设备 |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110010201A (zh) * | 2019-04-16 | 2019-07-12 | 山东农业大学 | 一种rna选择性剪接位点识别方法及系统 |
| CN110322925A (zh) * | 2019-07-18 | 2019-10-11 | 杭州纽安津生物科技有限公司 | 一种预测融合基因产生新生抗原的方法 |
| CN110910950A (zh) * | 2019-11-18 | 2020-03-24 | 广州竞远生物科技有限公司 | 一种联合分析单细胞scRNA-seq和scATAC-seq的流程方法 |
| US20200098448A1 (en) * | 2018-09-24 | 2020-03-26 | Tempus Labs, Inc. | Methods of normalizing and correcting rna expression data |
| CN111081311A (zh) * | 2019-12-26 | 2020-04-28 | 青岛科技大学 | 基于深度学习的蛋白质赖氨酸丙二酰化位点预测方法 |
| CN111192631A (zh) * | 2020-01-02 | 2020-05-22 | 中国科学院计算技术研究所 | 用于构建用于预测蛋白质-rna相互作用结合位点模型的方法和系统 |
| CN111755071A (zh) * | 2019-03-29 | 2020-10-09 | 中国科学技术大学 | 基于峰聚类的单细胞染色质可及性测序数据分析方法和系统 |
-
2021
- 2021-03-08 CN CN202110251298.3A patent/CN115050416A/zh active Pending
-
2022
- 2022-03-08 WO PCT/CN2022/079788 patent/WO2022188785A1/zh active Application Filing
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20200098448A1 (en) * | 2018-09-24 | 2020-03-26 | Tempus Labs, Inc. | Methods of normalizing and correcting rna expression data |
| CN111755071A (zh) * | 2019-03-29 | 2020-10-09 | 中国科学技术大学 | 基于峰聚类的单细胞染色质可及性测序数据分析方法和系统 |
| CN110010201A (zh) * | 2019-04-16 | 2019-07-12 | 山东农业大学 | 一种rna选择性剪接位点识别方法及系统 |
| CN110322925A (zh) * | 2019-07-18 | 2019-10-11 | 杭州纽安津生物科技有限公司 | 一种预测融合基因产生新生抗原的方法 |
| CN110910950A (zh) * | 2019-11-18 | 2020-03-24 | 广州竞远生物科技有限公司 | 一种联合分析单细胞scRNA-seq和scATAC-seq的流程方法 |
| CN111081311A (zh) * | 2019-12-26 | 2020-04-28 | 青岛科技大学 | 基于深度学习的蛋白质赖氨酸丙二酰化位点预测方法 |
| CN111192631A (zh) * | 2020-01-02 | 2020-05-22 | 中国科学院计算技术研究所 | 用于构建用于预测蛋白质-rna相互作用结合位点模型的方法和系统 |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115579060A (zh) * | 2022-12-08 | 2023-01-06 | 国家超级计算天津中心 | 基因位点检测方法、装置、设备及介质 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115050416A (zh) | 2022-09-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| WO2022188785A1 (zh) | 融合深度学习模型的单细胞转录组计算分析方法和系统 | |
| US11954614B2 (en) | Systems and methods for visualizing a pattern in a dataset | |
| Jongeneel et al. | An atlas of human gene expression from massively parallel signature sequencing (MPSS) | |
| CA3049682C (en) | Methods for non-invasive assessment of genetic alterations | |
| CA3050055C (en) | Methods and processes for assessment of genetic variations | |
| Debey et al. | A highly standardized, robust, and cost-effective method for genome-wide transcriptome analysis of peripheral blood applicable to large-scale clinical trials | |
| US20210381056A1 (en) | Systems and methods for joint interactive visualization of gene expression and dna chromatin accessibility | |
| CA3049457C (en) | Methods for non-invasive assessment of copy number alterations | |
| CN111128299A (zh) | 一种结直肠癌预后显著相关ceRNA调控网络的构建方法 | |
| CA3194557A1 (en) | Sequencing adapter manufacture and use | |
| CN112599198A (zh) | 一种用于宏基因组测序数据的微生物物种与功能组成分析方法 | |
| US20110106739A1 (en) | Method for determining the presence of disease | |
| WO2019242186A1 (zh) | 确定检测靶点的方法、装置、计算机设备和存储介质 | |
| KR102124193B1 (ko) | 기계 학습을 이용한 우울증 또는 자살 위험 예측용 마커 발굴 방법, 우울증 또는 자살 위험 예측용 마커, 및 기계 학습을 이용한 우울증 또는 자살 위험 예측 방법 | |
| US20140058682A1 (en) | Nucleic Acid Information Processing Device and Processing Method Thereof | |
| WO2006048264A2 (en) | Gene expression profiling in acute lymphoblastic leukemia (all), biphenotypic acute leukemia (bal), and acute myeloid leukemia (aml) m0 | |
| CN114974432A (zh) | 一种生物标志物的筛选方法及其相关应用 | |
| JP2022534236A (ja) | 多重オミックス分析を利用した鬱病または自殺危険の予測用マーカー発掘方法、鬱病または自殺危険の予測用マーカー、及び多重オミックス分析を利用した鬱病または自殺危険の予測方法 | |
| CN113257354B (zh) | 基于高通量实验数据挖掘进行关键rna功能挖掘的方法 | |
| Nash et al. | Pipeline for Integrated Microarray Expression Normalization Tool Kit (PIMENTo) for Tumor Microarray Profiling Experiments | |
| Sundarrajan et al. | 5 Big Data and | |
| Sundarrajan et al. | Big Data and Transcriptomics | |
| CN117219172A (zh) | 三阴性乳腺癌肿瘤细胞通路分群方法、应用及系统 | |
| WO2023161482A1 (en) | Epigenetic biomarkers for the diagnosis of thyroid cancer | |
| KR101244543B1 (ko) | 17-β 에스트라디올에 대한 노출 여부 판단용판단용 유전자 마커군, 마이크로어레이 칩 및 이를 이용한 판단 방법 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| 121 | Ep: the epo has been informed by wipo that ep was designated in this application |
Ref document number: 22766308 Country of ref document: EP Kind code of ref document: A1 |
|
| NENP | Non-entry into the national phase |
Ref country code: DE |
|
| 122 | Ep: pct application non-entry in european phase |
Ref document number: 22766308 Country of ref document: EP Kind code of ref document: A1 |