[go: up one dir, main page]

CN114875118B - 确定细胞谱系的方法、试剂盒和装置 - Google Patents

确定细胞谱系的方法、试剂盒和装置 Download PDF

Info

Publication number
CN114875118B
CN114875118B CN202210756234.3A CN202210756234A CN114875118B CN 114875118 B CN114875118 B CN 114875118B CN 202210756234 A CN202210756234 A CN 202210756234A CN 114875118 B CN114875118 B CN 114875118B
Authority
CN
China
Prior art keywords
cell
mitochondrial
transcriptome
cdna
generation sequencing
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210756234.3A
Other languages
English (en)
Other versions
CN114875118A (zh
Inventor
马玉
李丕栋
默芳
贺照人
武玉胜
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Baitu Shengke Beijing Intelligent Technology Co ltd
Original Assignee
Beijing Baitu Zhijian Technology Service Co ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing Baitu Zhijian Technology Service Co ltd filed Critical Beijing Baitu Zhijian Technology Service Co ltd
Priority to CN202210756234.3A priority Critical patent/CN114875118B/zh
Publication of CN114875118A publication Critical patent/CN114875118A/zh
Application granted granted Critical
Publication of CN114875118B publication Critical patent/CN114875118B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • C—CHEMISTRY; METALLURGY
    • C12—BIOCHEMISTRY; BEER; SPIRITS; WINE; VINEGAR; MICROBIOLOGY; ENZYMOLOGY; MUTATION OR GENETIC ENGINEERING
    • C12Q—MEASURING OR TESTING PROCESSES INVOLVING ENZYMES, NUCLEIC ACIDS OR MICROORGANISMS; COMPOSITIONS OR TEST PAPERS THEREFOR; PROCESSES OF PREPARING SUCH COMPOSITIONS; CONDITION-RESPONSIVE CONTROL IN MICROBIOLOGICAL OR ENZYMOLOGICAL PROCESSES
    • C12Q1/00—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions
    • C12Q1/68—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions involving nucleic acids
    • C12Q1/6813—Hybridisation assays
    • C12Q1/6827—Hybridisation assays for detection of mutation or polymorphism
    • C—CHEMISTRY; METALLURGY
    • C12—BIOCHEMISTRY; BEER; SPIRITS; WINE; VINEGAR; MICROBIOLOGY; ENZYMOLOGY; MUTATION OR GENETIC ENGINEERING
    • C12Q—MEASURING OR TESTING PROCESSES INVOLVING ENZYMES, NUCLEIC ACIDS OR MICROORGANISMS; COMPOSITIONS OR TEST PAPERS THEREFOR; PROCESSES OF PREPARING SUCH COMPOSITIONS; CONDITION-RESPONSIVE CONTROL IN MICROBIOLOGICAL OR ENZYMOLOGICAL PROCESSES
    • C12Q1/00—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions
    • C12Q1/68—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions involving nucleic acids
    • C12Q1/6869—Methods for sequencing
    • 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
    • G16B10/00—ICT specially adapted for evolutionary bioinformatics, e.g. phylogenetic tree construction or analysis
    • 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
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/20—Allele or variant detection, e.g. single nucleotide polymorphism [SNP] detection
    • 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
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/50—Mutagenesis
    • 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
    • G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
    • G16B25/10—Gene or protein expression profiling; Expression-ratio estimation or normalisation
    • 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
    • G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
    • G16B25/20—Polymerase chain reaction [PCR]; Primer or probe design; Probe optimisation
    • 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
    • G16B35/00—ICT specially adapted for in silico combinatorial libraries of nucleic acids, proteins or peptides
    • G16B35/10—Design of libraries
    • 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
    • G16B40/30—Unsupervised data analysis

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Chemical & Material Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Biotechnology (AREA)
  • Biophysics (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Molecular Biology (AREA)
  • Genetics & Genomics (AREA)
  • Organic Chemistry (AREA)
  • Analytical Chemistry (AREA)
  • Zoology (AREA)
  • Wood Science & Technology (AREA)
  • Biochemistry (AREA)
  • Library & Information Science (AREA)
  • Data Mining & Analysis (AREA)
  • Immunology (AREA)
  • General Engineering & Computer Science (AREA)
  • Microbiology (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physiology (AREA)
  • Artificial Intelligence (AREA)
  • Bioethics (AREA)
  • Animal Behavior & Ethology (AREA)
  • Databases & Information Systems (AREA)
  • Epidemiology (AREA)
  • Evolutionary Computation (AREA)
  • Public Health (AREA)
  • Software Systems (AREA)
  • Chemical Kinetics & Catalysis (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明涉及确定多细胞真核生物中的单细胞线粒体转录组中的变异的方法，并涉及基于细胞中线粒体变异的存在确定多细胞真核生物中的细胞谱系的方法，包括富集线粒体转录组以及进行二代测序和三代测序。本发明还涉及线粒体核酸文库和用于制备线粒体转录组文库的试剂盒、以及基于二代测序和三代测序来识别线粒体测序数据中线粒体变异的装置和确定多细胞真核生物中细胞谱系的装置。

Description

确定细胞谱系的方法、试剂盒和装置
技术领域
本发明涉及基因检测领域。具体而言，本发明涉及一种确定多细胞真核生物中的细胞谱系的方法，包括富集线粒体全长转录组、确定线粒体全长转录组的序列和其变异，以用于细胞谱系追踪。本发明还涉及用于富集线粒体全长转录组的试剂盒、以及基于二代测序和三代测序来快速地高通量确定线粒体测序数据中线粒体突变的方法和装置。
背景技术
细胞谱系(cell lineage)追踪允许人们跟踪连续的细胞后代。追踪细胞谱系对于了解控制多细胞生物发育的规则和描绘涉及具有不同谱系层次的多种细胞类型分化的复杂生物过程至关重要。
细胞谱系追踪在肿瘤和免疫系统研究中尤其具有重要意义。通过肿瘤细胞谱系追踪可以对肿瘤克隆的性质精确定义并为寻找靶点提供信息。细胞谱系追踪的一个典型应用是对比用药前、响应用药、不响应用药、用药后肿瘤进展等各阶段的单个肿瘤细胞谱系信息来判断药物对哪些细胞克隆的杀伤力更强，了解耐药性克隆、免疫逃逸克隆的演化，进行单细胞合成致死、合成拯救研究，促进肿瘤干细胞寻找、驱动突变寻找等。通过免疫细胞谱系追踪可以探索没有T细胞受体(TCR)或B细胞受体(BCR)的免疫细胞例如髓系细胞、NK细胞是如何在组织间转移的；也可以实现对T细胞进行高分辨率的追踪，从而对于具有相同TCR序列但空间分布或表型不同的T细胞，通过对时间信息的了解而获知T细胞的转移分化过程，并帮助寻找潜在靶点；也可以在免疫治疗后，研究T细胞来源(血液/淋巴结/癌旁)及T细胞如何分化或转移。
现有技术中，细胞谱系追踪的方法包括例如基因组标记方法(artificiallabeling)，其是通过病毒感染或者转基因的手段人为在目的祖先细胞基因组中插入细胞标记，比如GFP蛋白、β-半乳糖苷酶等(Kretzschmar K,Watt FM. 2012. Lineage tracing.Cell 148:33–45)。这种方法常被用在细胞系或者模式动物中，而无法应用在人体上。对人体进行实验性的谱系追踪不符合伦理要求。
在人体研究过程中大多采用检测自然变异的方式进行谱系追踪。关于自然变异，主要存在核基因组层面的变异和线粒体基因组层面的变异。
对于核基因组层面的变异，本领域已知，每个细胞的核基因组在分裂时会产生约60个突变，后代细胞会携带这些突变，因此，本领域技术人员可以利用这些突变信息来追溯细胞分裂分化的轨迹。如果对单个细胞进行基因组测序(MDA、MALBAC等方法)，虽然可以得到比较全面的基因组突变信息，但是不管是在细胞通量上还是成本上，都很难满足复杂疾病需要进行大量单细胞测序的要求。
对于线粒体基因组层面的变异，本领域已知，线粒体是细胞的动力源，除了红细胞缺乏线粒体外，几乎在身体的每个细胞中都存在100到2000个拷贝的线粒体，线粒体DNA(mitochondrial DNA，mtDNA)长约16kb，远远小于核基因组3Gb的大小，而且线粒体DNA的突变频率是核基因组的大约30倍，每个细胞有100到2000份线粒体基因组(Datlinger, P.等人, (2021). "Ultra-high-throughput single-cell RNA sequencing andperturbation screening with combinatorial fluidic indexing." Nat Methods 18(6): 635-642)，因此，利用线粒体基因组变异进行细胞谱系追踪是一个很有吸引力的选择，它将能够提供例如随着细胞增殖和老化而在细胞内产生的自然突变标记而追踪细胞谱系。
线粒体谱系研究既可以依靠线粒体基因组信息进行研究，也可以依靠线粒体转录组信息进行研究。
现有技术中直接检测线粒体基因组的方法有例如S. Ludwig 等人(LineageTracing in Humans Enabled by Mitochondrial Mutations and Single-CellGenomics, Cell, 2019, 176(6): p1325-1339)公开的方法，但是，该方法存在实验通量低，操作时间久，只能获得单一的线粒体DNA信息的缺点。
mtscATAC-seq利用单细胞ATAC-seq测序技术来捕获线粒体基因组(Lareau, C.A.等人, (2021). "Massively parallel single-cell mitochondrial DNA genotypingand chromatin profiling." Nat Biotechnol 39(4): 451-461)。在该方法中，先将细胞固定透化，以保证线粒体不会丢失，然后利用转座酶Tn5切割标记核基因组的同时会同样切割线粒体基因组。但是，该方法需要先将细胞固定透化，实验步骤较为复杂。如果在相同的测序数据量下进行分析，不同样本的线粒体基因有效信息量及测序数据深度均不同，不利于大规模的谱系研究。
现有技术中检测线粒体转录组的方法有例如MAESTER技术(Miller, T. E.等人,"Mitochondrial variant enrichment from high-throughput single-cell RNA-seqresolves clonal populations", Nat Biotechnol (2022)，https://doi.org/10.1038/s41587-022-01210-8)。该方法无需对细胞进行固定透化，但是，其实验过程复杂且所建的文库仅能通过二代测序获得较短reads(每个read约250bp)。由于线粒体基因组与核基因组有较大相似性，因此短reads测序的方式可能难以区分其来源于核基因组还是线粒体基因组。
因此，需要提供一种克服现有可用方法学局限性的、能够高通量、准确地获得细胞转录组信息和线粒体突变信息，以利于细胞谱系追踪的改进方法。
发明内容
本发明人通过深入研究，开发了一种通过确定线粒体转录组中的变异来确定多细胞真核生物中的细胞谱系的方法，其包括富集线粒体转录组以及进行短读长Bulk线粒体转录组测序、长读长单细胞线粒体全长转录组测序和任选地进行短读长单细胞转录组测序，以用于细胞谱系追踪。
因此，在第一方面，本发明提供了一种确定线粒体转录组中的变异的方法，该方法基于高通量测序，包括以下步骤：
a) 由来自受试者的样本中的单细胞RNA制备单细胞cDNA，优选地，通过全转录组扩增产生所述单细胞cDNA，其中所制备的单细胞cDNA中的每一cDNA包含通用PCR引物结合位点、用于确定转录组来源细胞的特征序列、用于确定所述来源细胞中的各转录本的特征序列，所述用于确定转录组来源细胞的特征序列优选细胞条形码(barcode)序列，所述用于确定所述来源细胞中的各转录本的特征序列优选唯一分子标识符(UMI)序列，例如，所述通用PCR引物结合位点、用于确定转录组来源细胞的特征序列、用于确定所述来源细胞中的各转录本的特征序列组合地位于cDNA的5’端或者组合地位于cDNA的3’端；
b) 当步骤a)的单细胞cDNA上的用于确定转录组来源细胞的特征序列位于5’端时，对步骤a)的单细胞cDNA使用线粒体cDNA 3’端特异性引物和能够结合步骤a)的单细胞cDNA 5’端通用PCR引物结合位点的通用PCR引物通过PCR富集线粒体cDNA，优选地，所述线粒体cDNA 3’端特异性引物位于线粒体cDNA 3’端约200bp范围内，例如，约150bp范围内、约100bp范围内、约50bp范围内、约30bp范围内；或者
当步骤a)的单细胞cDNA上的用于确定转录组来源细胞的特征序列位于3’端时，对步骤a)的单细胞cDNA使用线粒体cDNA 5’端特异性引物和能够结合步骤a)的单细胞cDNA3’端通用PCR引物结合位点的通用PCR引物通过PCR富集线粒体cDNA，优选地，所述线粒体cDNA 5’端特异性引物位于线粒体cDNA 5’端约200bp范围内，例如，约150bp范围内、约100bp范围内、约50bp范围内、约30bp范围内；
c) 对步骤b)富集的线粒体cDNA进行Bulk线粒体转录组二代测序和单细胞线粒体全长转录组三代测序，其中，所述Bulk线粒体转录组二代测序是对一个所述样本的多个单细胞的获自步骤b)的线粒体cDNA富集产物实施二代测序，以获得Bulk线粒体转录组二代测序数据；所述单细胞线粒体全长转录组三代测序是对单细胞线粒体转录组中的各转录本实施三代测序，以获得单细胞线粒体全长转录组三代测序数据；
优选地，还对步骤a)的单细胞cDNA进行单细胞转录组二代测序，例如，对靠近转录组3’端或5’端序列进行单细胞转录组二代测序，以获得用于确定转录组来源细胞的特征序列和用于确定所述来源细胞中的各转录本的特征序列的二代测序数据；
d) 用步骤c)获得的Bulk线粒体转录组二代测序数据中的用于确定转录组来源细胞的特征序列和用于确定所述来源细胞中的各转录本的特征序列的二代测序数据，对步骤c)获得的单细胞线粒体全长转录组三代测序数据中的用于确定转录组来源细胞的特征序列和用于确定所述来源细胞中的各转录本的特征序列的三代测序数据进行校正，然后将校正后的单细胞线粒体全长转录组三代测序数据与步骤c)获得的Bulk线粒体转录组二代测序数据中的变异位点数据(例如，添加、缺失、取代和／或删除，例如，SNP)进行整合，获得线粒体转录组中的变异信息；
优选地，用步骤c)获得的单细胞转录组二代测序数据中的用于确定转录组来源细胞的特征序列和用于确定所述来源细胞中的各转录本的特征序列的二代测序数据，对步骤c)获得的单细胞线粒体全长转录组三代测序数据中的用于确定转录组来源细胞的特征序列和用于确定所述来源细胞中的各转录本的特征序列的三代测序数据进行校正，然后将校正后的单细胞线粒体全长转录组三代测序数据与步骤c)获得的Bulk线粒体转录组二代测序数据中的变异位点数据(例如，添加、缺失、取代和／或删除，例如，SNP)进行整合，获得线粒体转录组中的变异信息；
所述线粒体转录组中的变异信息是例如，添加、缺失、取代和／或删除，例如，SNP。
在一些实施方案中，所述步骤b)的线粒体cDNA 3’端特异性引物或线粒体cDNA 5’端特异性引物是一条或多条引物，例如，一条或多条选自表1的引物，以及在1个或多个PCR管中实施步骤b)，由此，通过PCR富集线粒体cDNA。在一个实施方案中，在1个PCR管中实施步骤b)，优选地，根据线粒体中各转录本的表达水平比例来混合相应量的特异性引物，由此，通过PCR富集线粒体cDNA。
在一些实施方案中，所述步骤c)的三代测序是将额外的多个样本(例如，额外的1个样本-10个样本，例如，额外的2个、3个、4个、5个、6个、7个、8个、9个样本)分别经步骤a)和步骤b)后富集的线粒体cDNA进行混合，然后进行三代测序实施的，优选地，在所述三代测序时实施的建库还包括添加至少一个样本索引的核苷酸序列，例如，所述添加至少一个样本索引的核苷酸序列选自至少一对SEQ ID NO: 16和17、SEQ ID NO: 18和19、SEQ ID NO: 20和21、SEQ ID NO: 22和23、SEQ ID NO: 24和25、SEQ ID NO: 26和27样本索引引物。
在一些实施方案中，所述步骤d)包括：
(1) 从步骤c)生成的单细胞线粒体全长转录组三代测序数据，获得单细胞线粒体全长转录组的三代测序结果；
(2) 从步骤c)生成的二代Bulk线粒体转录组测序数据，获得线粒体转录组的高准确性二代测序序列，例如，其中包括barcode/UMI序列和SNP/InDel位点；
(3) 用所述(2)中获得的线粒体转录组的高准确性二代测序序列对所述(1)中获得的单细胞线粒体全长转录组的三代测序序列进行校正，获得具有高准确性单细胞变异位点(例如，SNP/InDel位点)的单细胞线粒体全长转录组序列。
在一些优选的实施方案中，所述步骤d)包括：
(1) 从步骤c)生成的三代单细胞线粒体全长转录组测序数据，获得单细胞线粒体全长转录组的变异位点，例如，单核苷酸多态性位点以及插入和／或删除(也简称为SNP/InDel)位点；
(2) 从步骤c)生成的二代Bulk线粒体转录组测序数据，获得线粒体转录组的短读长中的高准确性单细胞变异位点，例如，SNP/InDel位点；
(3) 从步骤c)生成的二代单细胞转录组测序数据，获得高准确性细胞barcode/UMI序列；
(4) 将所述(1)和(3)中分别获得的结果进行barcode/UMI映射，获得单细胞线粒体全长转录组序列，其具有高准确性细胞barcode/UMI序列和较高准确性的变异位点，例如，SNP/InDel位点；
(5) 对(4)得到的barcode/UMI校正后的单细胞线粒体全长转录组序列与(2)得到的高准确性单细胞变异位点(例如，SNP/InDel位点)进行整合，获得具有高准确性单细胞变异位点(例如，SNP/InDel位点)的单细胞线粒体全长转录组序列。
在一些实施方案中，所述步骤d)获得的单细胞变异位点在Bulk线粒体转录组二代测序数据中的至少5个读长中被检测到。
在第二方面，本发明提供了一种确定多细胞真核生物中细胞谱系的方法，其包括上述第一方面的确定线粒体转录组中的变异的方法所述的步骤a)- d)，并进一步包括：
步骤e)基于细胞中线粒体变异的存在对细胞进行聚类，并推断细胞的谱系；优选地，通过计算单细胞变异位点(例如，SNP/InDel位点)的变异频率，构建细胞-变异(例如，SNP/InDel变异) 矩阵，由此推断细胞的进化谱系。
在第三方面，本发明提供了一种线粒体核酸文库，其为通过上述第一方面的确定线粒体转录组中的变异的方法中所述的步骤a)- c)制备的二代Bulk线粒体转录组文库和三代单细胞线粒体全长转录组文库的组合文库。
在第四方面，本发明提供了用于制备线粒体转录组文库的试剂盒，其包含含有通用PCR引物结合位点、用于确定转录组来源细胞的细胞条形码(barcode)序列和用于确定所述来源细胞中的各转录本的唯一分子标识符(UMI)序列的核苷酸序列，用于添加至线粒体RNA反转录获得的cDNA上，以及与线粒体RNA反转录获得的cDNA退火的一个或多个引物组，且每个引物组包含一条正向引物和一条反向引物，用于一条线粒体核酸链的靶向扩增，由此扩增一条或多条线粒体核酸链。
在一些实施方案中，所述每个引物组中的一条引物是针对线粒体转录组的3’端特异性引物或5’端特异性引物，优选地，所述线粒体转录组的3’端特异性引物和/或5’端特异性引物位于线粒体转录组的3’端或5’端约200bp范围内，例如，约150bp范围内、约100bp范围内、约50bp范围内、约30bp范围内，例如，表1所示的针对线粒体转录组的一条或多条3’端特异性引物；所述每个引物组中的另一条引物是通用PCR引物，用于结合线粒体RNA反转录获得的cDNA上的通用PCR引物结合位点。
在一些实施方案中，本发明的用于制备线粒体转录组文库的试剂盒还包含向扩增的一条或多条线粒体核酸链添加至少一个样本索引的核苷酸序列(例如，所述添加至少一个样本索引的核苷酸序列选自至少一对SEQ ID NO: 16和17、SEQ ID NO: 18和19、SEQ IDNO: 20和21、SEQ ID NO: 22和23、SEQ ID NO: 24和25、SEQ ID NO: 26和27样本索引引物)，用于混合样本的高通量三代测序。
在第五方面，本发明提供了一种确定线粒体测序数据中线粒体变异的装置，其包含一个或多个以下模块：
模块1. 接收单细胞线粒体全长三代测序数据，并映射到靶线粒体DNA样本；
模块2. 接收Bulk线粒体转录组二代测序数据，获得高准确性的基因水平的变异和barcode/UMI序列，例如，获得高准确性单细胞变异位点(例如，SNP/InDel位点)和barcode/UMI序列；
模块3. 用模块2中获得的高准确性barcode/UMI序列校正模块1的barcode/UMI序列， 再将经barcode/UMI序列校正的模块1序列与模块2得到的高准确性的变异位点(例如，SNP/InDel位点)进行整合，获得高准确性的线粒体变异位点(例如，SNP/InDel位点)。
在一些实施方案中，本发明的识别线粒体测序数据中线粒体变异的装置包含一个或多个以下模块：
模块1. 接收单细胞线粒体全长三代测序数据，并映射到靶线粒体DNA样本；
模块2. 接收单细胞转录组二代测序数据并分析，获得高准确性barcode/UMI序列；
模块3. 接收Bulk线粒体转录组二代测序数据并分析，获得高准确性的基因水平的变异，从而获得高准确性单细胞变异位点(例如，SNP/InDel位点)；
模块4. 将所述模块1和模块2中分别获得的获得的单细胞线粒体全长三代测序数据、单细胞转录组二代测序数据barcode/UMI映射，获得高准确性barcode/UMI序列及较高准确性单细胞变异位点(例如，SNP/InDel位点)；
模块5. 对模块4得到的barcode/UMI校正后的单细胞全长转录组序列与模块3得到的高准确性的变异位点(例如，SNP/InDel位点)进行整合，获得高准确性的线粒体变异位点(例如，SNP/InDel位点)。
在第六方面，本发明提供了一种确定多细胞真核生物中细胞谱系的装置，其包含本发明在第五方面所述的装置中的模块；以及一个额外的模块：基于第五方面的模块得到的高准确性的线粒体变异位点(例如，SNP/InDel位点)的存在对细胞进行聚类，并推断细胞的谱系；优选地，通过计算单细胞变异位点(例如，SNP/InDel位点)的变异频率，构建细胞-变异(例如，SNP/InDel变异) 矩阵，由此获得细胞的进化谱系。
发明的效果
本发明涉及通过确定线粒体转录组中的变异来确定多细胞真核生物中的细胞谱系的方法，所述方法是通过富集线粒体转录组以及进行二代Bulk线粒体转录组测序、三代单细胞线粒体全长转录组测序和任选地二代单细胞转录组测序而实施的。由本发明方法获得的线粒体转录组变异数据更准确，在分析细胞谱系方面有更广泛的应用。本发明的该方法不仅与多个单细胞组学平台兼容，而且操作简单，不存在细胞固定透化等损伤核酸的实验步骤，更加贴合本发明的目的。
本发明还提供了用于富集线粒体转录组的试剂盒，用于防止扩增线粒体cDNA期间的PCR偏差(PCR bias)，从而使得Bulk线粒体转录组二代测序和长读长单细胞线粒体全长转录组三代测序获得的数据量稳定可控，极大地提高了细胞谱系分析中有效信息的占比。
进一步地，使用本发明的对短读长Bulk线粒体转录组二代测序和长读长单细胞线粒体全长转录组三代测序，以及任选地短读长单细胞转录组二代测序，并使用本发明的对获得的数据进行识别和处理的生物信息分析装置，能够准确地进行单细胞中线粒体突变的检出和多细胞真核生物中细胞谱系的确定。
附图说明
结合以下附图一起阅读时，将更好地理解以下详细描述的本发明的优选实施方案。出于说明本发明的目的，图中显示了目前优选的实施方案。然而，应当理解本发明不限于图中所示实施方案的精确安排和手段。
图1显示了确定多细胞真核生物中的细胞谱系方法的一个示例性流程图。
图2例示了通过10x Genomics 5’试剂盒(10x Genomics，Chromium Next GEMSingle Cell 5ʹ Kit v2，1000263和10x Genomics，Chromium Next GEM Chip K SingleCell Kit，1000287)获得的单细胞转录组的cDNA产物的结构。图中，“read1”表示“read 1测序引物的互补序列，是read 1测序引物的结合位点，read1测序引物是一种通用引物”，“10x细胞条形码(也称为10xbarcode)”表示使用10x Genomics 5’试剂盒进行反转录获得的cDNA产物上的细胞条形码序列，“UMI”表示“唯一分子标识符(unique molecularidentifier)”，“TSO”表示模板转换寡核苷酸(template switching oligo)。
图3例示了对线粒体转录组cDNA的富集、以及对线粒体转录组cDNA富集产物进行长读长的单细胞线粒体全长转录组三代测序和短读长的Bulk线粒体转录本二代测序中测序文库构建的流程图。
图4显示了使用单细胞转录组的总全长cDNA进行短读长的单细胞转录组二代测序文库构建的流程图。
图5例示了生物信息分析框架图。
图6显示了对实施例1中的细胞样本利用高通量的单细胞线粒体转录组三代测序数据和Bulk线粒体转录组二代测序数据进行生物信息学分析得到的细胞-变异频率层次聚类图。
图7显示了对实施例1中的细胞样本利用高通量的单细胞线粒体转录组三代测序数据和Bulk线粒体转录组二代测序数据进行生物信息学分析推断的细胞谱系图。
图8显示了对实施例1中的细胞样本利用高通量的单细胞线粒体转录组三代测序数据、Bulk线粒体转录组二代测序数据和单细胞转录组二代测序数据进行生物信息学分析得到的细胞-变异频率层次聚类图。
图9显示了对实施例1中的细胞样本利用高通量的单细胞线粒体转录组三代测序数据、Bulk线粒体转录组二代测序数据和单细胞转录组二代测序数据进行生物信息学分析推断的细胞谱系图。
图10显示了对实施例7中的细胞样本利用高通量的单细胞线粒体转录组三代测序数据、Bulk线粒体转录组二代测序数据和单细胞转录组二代测序数据进行生物信息学分析得到的细胞-变异频率层次聚类图。
图11显示了对实施例7中的细胞样本利用高通量的单细胞线粒体转录组三代测序数据、Bulk线粒体转录组二代测序数据和单细胞转录组二代测序数据进行生物信息学分析推断的细胞谱系图。
图12显示了对实施例8中的细胞样本利用高通量的单细胞线粒体转录组三代测序数据、Bulk线粒体转录组二代测序数据和单细胞转录组二代测序数据进行生物信息学分析得到的细胞-变异频率层次聚类图。
图13显示了对实施例8中的细胞样本的短读长二代测序单细胞转录组表达数据进行聚类分析后，获得的细胞亚群。
图14显示了对实施例8中的细胞样本基于谱系聚类的细胞亚群分类结果。
具体实施方式
除非另外限定，否则本文中所用的全部技术与科学术语具有如本发明所属领域的普通技术人员通常理解的相同含义。本文所提及的全部出版物、专利申请、专利和其他参考文献通过引用的方式完整地并入。此外，本文中所述的材料、方法和例子仅是说明性的并且不意在是限制性的。本发明的其他特征、目的和优点将从本说明书及附图并且从后附的权利要求书中显而易见。
定义
为了解释本说明书，将使用以下定义，并且只要适当，以单数形式使用的术语也可以包括复数，并且反之亦然。要理解，本文所用的术语仅是为了描述具体的实施方案，并且不意欲是限制性的。
术语“约”在与数字数值联合使用时意为涵盖具有比指定数字数值小5%的下限和比指定数字数值大5%的上限的范围内的数字数值。
如本文所用，术语“和/或”意指可选项中的任一项或可选项的两项或多项。
在本文中，当使用术语“包含”或“包括”时，除非另有指明，否则也涵盖由所述及的要素、整数或步骤组成的情形。例如，当提及“包含”某个具体序列的引物时，也旨在涵盖由该具体序列组成的引物。
除非另有说明，否则“一个”、“一种”、“该”和“至少一个”可互换使用并且表示一个或多于一个。
“核酸序列”指核苷酸(即，例 如，腺嘌呤(A)，胸腺嘧啶(T)，胞嘧啶(C)，鸟嘌呤(G)，和/或尿嘧啶(U))的任何聚合物。将包含脱氧核糖核苷的核酸序列称为脱氧核糖核酸(DNA)。将包含核糖核苷的核酸序 列被称为核糖核酸(RNA)。RNA可进一步表征为以下类型，(1)蛋白质编码RNA，也称为信使RNA (mRNA)；(2)非编码RNA(non-codingRNA，ncRNA)，ncRNA的核苷酸一般不长，共同特点是转录自基因组，但是不翻译成蛋白质，在RNA水平上行使各自的生物学功能，可以从长度上来划分为小于50nt的ncRNA，例如微量RNA(miRNA)，小干扰RNA(siRNA)等；50nt到500nt的ncRNA，例如，转运RNA(tRNA)、核糖体RNA(rRNA)、小核(snRNA)和小核仁RNA(snoRNA)等；大于500nt的ncRNA，包括长的mRNA-like的非编码RNA、长的不带polyA尾巴的非编码RNA等。
如本文所用，术语“测序”通常是指用于确定一种或多种多核苷酸中的核苷酸碱基的序列的方法和技术。多核苷酸可以是核酸分子例如脱氧核糖核酸(DNA)或核糖核酸(RNA)，包括它们的变体或衍生物(例如，单链DNA)。可以通过目前可用的各种装置执行测序，例如但不限于通过Illumina®、Pacific Biosciences(PacBio®)、Oxford Nanopore®或Life Technologies(IonTorrent®)的测序装置。可以使用核酸扩增、聚合酶链式反应(PCR)或等温扩增来执行测序。在一些实施方案中，此类装置提供测序读长(本文中也称为“read”)。读长(read)包括一串核酸碱基，其对应于已经测序的核酸分子的序列。
“二代测序”与“下一代测序(Next-generation sequencing，NGS)”可互换地使用，是一种基于PCR和基因芯片发展而来的测序技术。与一代测序通过合成终止测序不同的是：二代测序开创性的引入了可逆终止末端，从而实现边合成边测序(Sequencing bySynthesis)。二代测序在DNA复制过程中通过捕捉新添加的碱基所携带的特殊标记(一般为荧光分子标记)来确定核酸的序列。二代测序的出现将分子生物学研究推向了一个高通量发展的时代，利用NGS产生了大量基因组和转录组数据。但是，NGS测序平台的读长基本不超过500个碱基，是短读长测序，主要适用于检测单个细胞中的单核苷酸变异(SNV)、小的插入缺失(<50bp，Indel)、以及拷贝数变异(CNV)等。通过NGS测序获得的序列准确性高，但覆盖度低。
三代测序(third-generation sequencing, TGS)是指单分子实时测序技术，也叫“从头测序技术”。与一代测序技术和二代测序技术相比，三代测序技术最大的特点就是单分子测序，测序过程中无需进行PCR扩增，实现了对每一条DNA分子的单独测序，克服了NGS测序平台读长短的缺点。目前TGS技术主要分为两类，一类是由美国太平洋生物科学公司(PacBio)推出的SMRT测序(Single Molecule Real-Time Sequencing, SMRT)，即单分子实时测序，其中，核酸分子在被合成过程中，通过读取荧光信号识别碱基。另一类是由英国Oxford Nanopore公司推出的基于核苷酸通过纳米孔时合成不同碱基所产生的电信号差异以实现单分子测序得方法。三代测序是一种长读长(long-read)测序技术，现有技术中存在利用例如Nanopore MinION测序仪或PacBio平台实施三代测序。
术语“单细胞测序”是指对单个细胞的基因组进行测序。单个细胞是生物体的基本结构和功能单位。在相同的外部刺激或生理条件下，源自同一类型细胞的细胞可能会表现出细胞间差异，即，异质性(heteroplasmy)。构成例如大脑、血液系统、免疫系统的细胞之间存在异质性。如果能够从组织系统中挑选多个不同的单细胞进行研究，则可以提供该组织中多个不同单细胞的有价值的信息，甚至能够重建出所述系统。
术语“测序深度(Sequencing Depth)”是指测序得到的总碱基数与待测基因组或转录组或染色体区段大小的比值。假设一个基因组大小为2M，测序深度为10X，那么获得的总数据量为20M。测序深度能够用测序得到的碱基总量(bp)与待测基因组或转录组或染色体区段大小的比值来表示。
术语“覆盖度(coverage)”是指，基因组或转录组或染色体区段上获知序列信息的序列部分占整个组或区段的比例。在一些实施方案中，覆盖度是指，(例如通过序列检测手段，如测序)测到序列信息的碱基数占所检测区域的总碱基数的比例。例如，在测序检测全基因组序列时，由于存在大片段拼接的缺口(gap)、测序读长有限、重复序列等问题，测序后得到的基因组序列通常无法完全覆盖基因组的所有区域，此时，覆盖度就是最终得到的测序碱基数占整个基因组碱基数的比例。例如，对人的基因组测序获得的覆盖度为98.5%，这表明该基因组还有1.5%的区域的序列没有得到。在另一些实施方案中，覆盖度是指，就所检测的区域而言，例如通过测序分析测到序列信息的基因位点(例如SNP位点或基因变异位点)的数目，占该区域中所检测的总基因位点数的比例。所检测区域可以是全基因组、特定染色体、或特定染色体区段、或转录物组、或特定转录区域。
“文库构建”是指在DNA的一端或两端加上可以用于测序以及文库拆分的接头(adaptor)结构。在一个实施方案中，所构建的文库由P5/Truseq Read/TSO/barcode/UMI和/或DNA插入片段(DNA Insert)组成，其中，P5与流动池(flow cell)上的接头引物序列互补和相同，进行簇(cluster)生成；Truseq Read与全长Illumina TruSeq Read 1 测序引物互补和相同。
单细胞测序的基础流程是：获得单细胞，单细胞裂解，对单细胞内的核酸进行扩增，以达到测序的最小上样量的要求，然后上机，读取数据并分析。
可以通过流式细胞术(FACS)、激光捕获显微切割技术(Laser capturemicrodissection，LCM)、微流控技术(Microfluidics)获得单细胞，其中微流控技术是一种用于精确控制微量液体的技术，微流控芯片是实施该技术的平台。微流控芯片通过细微的管道对液体实施操控，微流控对液体的操控程度刚好适合于单细胞样品的处理操作。
细胞裂解方法可以是物理方法、化学方法或酶法。存在三种主要的物理细胞裂解形式：机械裂解、热裂解和电裂解。化学裂解应用裂解缓冲液并诱导高效裂解以破坏细胞。酶促细胞裂解是减少DNA断裂的最温和的方法。可以使用蛋白酶，如胃蛋白酶和胰蛋白酶。
可以通过全基因组扩增(Whole Genome Amplification，WGA)和/或全转录组扩增(Whole Transcriptome Amplification，WTA)对细胞内的核酸进行扩增。全基因组扩增方法包括例如简并寡核苷酸引发的PCR(Degenerate Oligonucleotide Primed PCR，DOC-PCR)、多重置换扩增(Multiple Displacement Amplification，MAD)、基于多重退火和成环的扩增循环(MALBAC)。全转录组扩增包括对单细胞中提取的RNA逆转录出cDNA，然后对逆转录生成的cDNA进行扩增。
单细胞基因组测序，即scDNA-seq，用来获得细胞基因组的突变和结构变化。
单细胞转录组测序，即 scRNA-seq，可以以单细胞分辨率来测量转录组中的基因表达，识别细胞簇中的生物学相关差异。单细胞RNA测序方法包括例如Smart-seq2、CEL-seq2、sci-RNA-seq、10x Chromium、Drop-seq、Seq-Well、inDrops，其中Smart-seq2和CEL-seq2是基于微孔板的低通量的方法，其他5种方法是高通量方法。
如本文所用，术语“条形码(barcode)”是指能够提供分析物的信息的标记或标识符。条形码可以是分析物的一部分。条形码也可以独立于分析物。除了分析物的内源特征(例如，分析物的大小或一个或多个末端序列)之外，条形码可以是连接至分析物(例如，核酸分子)的一个标签或者是多个标签的组合。条形码是独特的。条形码可以有多种不同的形式。例如，条形码可包括多核苷酸条形码；随机核酸和/或氨基酸序列；和合成核酸和/或氨基酸序列。条形码可以以可逆或不可逆的方式连接于分析物。在样品测序之前、期间和/或之后，可以将条形码添加到例如脱氧核糖核酸(DNA)或核糖核酸(RNA)样品的片段中。条形码可以允许确定和/或定量各个测序读长。在一些实施方案中，当条形码是一段核苷酸序列时，任意两个条形码之间至少存在两个或者两个以上碱基的不同，因为测序存在对碱基的误读，这样的设计可以避免将两个条形码混淆。对于一段长16bp的核苷酸序列，大约存在350万种条形码(barcode)。在一个实施方案中，barcode用来区别GEMs，是对一个单细胞分配的标记。
术语“固定”是指固定细胞的任何方法或过程。在本文中与“交联”可互换地使用。许多化学品能够用来固定细胞，包括但不限于，甲醛、福尔马林、或戊二醛。
术语“富集”指对线粒体核酸的任何分离或目的线粒体核酸区段的比例相对于核酸组合物中的其他核酸区段增加。例如，富集或分离目的线粒体核酸核酸或核酸区段可以包括正向方法，如“分出”目的线粒体核酸核酸或核酸区段，或可以包括负向方法，如排除不作为目的线粒体核酸的核酸或排除不包含线粒体核酸区段的核酸区段。备选地，富集包括选择性或定向扩增目的线粒体核酸或核酸区段。这类选择性或定向扩增目的线粒体核酸将增加线粒体核酸在核酸组合物中的比例(即富集所述线粒体核酸)。
“单核苷酸多态性”或“SNP”是基因组内的单核苷酸变化(即A，C，G或T)，其在生物物种的成员之间或成对染色体之间不同。
本文提及的“纯化”，可以指核酸已经历处理(即，例如分级)以除去各种其它组分。当使用术语“基本上纯化的”时，表示目的核酸构成组合物中的主要组分，诸如构成组合物的约50％、约60％、约70％、 约80％、约90％、约95％或更多(即，例如，重量/重量和/或重量/体积)。
多重 PCR (Multiplex polymerase chain reaction，MPCR) 是将不同的引物对和模板分散于相对独立的空间中进行PCR扩增。例如，在同一个反应体系中加入不同的引物对，针对不同的模板或同一模板的不同区域进行特异性的PCR扩增，从而得到多个目的扩增产物，结合一定的检测手段能够实现同时对多个靶标进行检测。MPCR 具有高效率、高通量、低成本的特性。
术语“通用引物结合位点”是指存在于(通常，通过人工添加至)不同的许多目标核酸中的、与通用引物结合的位点。“通用引物”通过结合“通用引物结合位点”指导引物延伸。测序时所使用的通用引物也称为通用测序引物。PCR时所使用的通用引物也称为通用PCR引物。通用PCR引物和通用测序引物可以相同或不同。
术语受试者的“样本”是指例如来自所述受试者个体的血液、唾液、口腔拭子、尿液、指甲、毛囊、皮屑、组织、体液样本。
术语“混合样本(pooling)”是指在进行高通量测序时，对多个样本进行混合，以最大化利用测序仪实施高通量测序。在本发明的方法中可以混合任意数目的样本，不局限于混合几个样本。
术语“样本索引(sample index)”是在三代测序文库构建中使用的样本标记物，用于区分不同样本的三代测序数据。样本索引(sample index)通常是一段碱基组成的DNA序列。在三代测序文库构建过程中，每个样本的文库都要分别带上不同的sample index序列，像样本的“身份证号”一样，能将流动池(flow cell)里所测得的海量序列正确地拆分。选择合适的样本索引集合(sample index sets)以确保在多重测序(multiplexed sequencing)运行中多个样本索引之间没有重叠。
术语“流动池(flow cell)”是指固体表面的腔室，一种或多种流体物质可以流过该固体表面。流动池以及相关的流体系统和检测平台的实例例如在US 7,329,492；US 7,211,414中描述。
术语线粒体“核酸文库”是本发明创建的线粒体核酸的集合，可以通过生物合成来制备。
术语“Cell Ranger”是10X genomics公司提供的一套针对单细胞 RNA 测序输出结果进行比对、定量、聚类及基因表达分析的分析流程，它包含有与单细胞基因表达分析相关的四个pipelines，分别是：
cellranger mkfastq 流程：其功能为将 Illumina 测序仪产生的raw base call(BCL) 文件解析成 FASTQ 文件；
cellranger count 流程：其功能为将 cellranger mkfastq 产生的或其他来源的FASTQ 文件进行比对、过滤、barcode 计数以及 UMI 计数，并可以生成 feature-barcode 定量矩阵，随后确定细胞群并进行基因表达分析；
cellranger aggr 流程：其功能为将多个cellranger count 产生的数据进行整合、标准化，并可以对整合后的数据进行分析；
cellranger reanalyze 流程：其功能为使用 cellranger count 或 cellrangeraggr 产生的表达矩阵重新进行降维、聚类等后续分析。
以上四个管线(pipeline)均将转录组常用比对软件STAR封装其中，所述STAR用来判断reads是否比对到了基因组上。
术语“GATK”是 Genome Analysis ToolKit的缩写，是一款从高通量测序数据中分析变异信息的软件(下载链接：https://software.broadinstitute.org/gatk/download/)，可以用于寻找SNP和Indel，一个标准的分析流程称为GATK Best Practise，主要包括以下几个步骤：
数据预处理：对从测序仪下机后的数据进行质控，去除低质量的reads，将过滤后的reads比对到参考基因组上，产生BAM格式的比对文件。
寻找变异：进行变异召回(variant calling)，寻找例如SNP和Indel，将比对数据存储在VCF格式的文件中。
使用寻找出的变异位点进行后续的分析。
术语“Fastq”是序列数据存储的标准格式之一，每4行为一条read的信息，包含测序read名、序列、正反链标示、序列质量值。
本发明的确定线粒体转录组中的变异和确定多细胞真核生物中的单细胞谱系的方法
线粒体在细胞功能中起着重要的作用，不仅提供细胞所需90%以上的能量，还是细胞凋亡启动和执行的主要细胞器之一。人类线粒体基因组包括16569个碱基对，其中包含2个转录为核糖体RNA的基因、13个转录为线粒体信使RNA并编码蛋白质的基因、以及22个tRNA基因。人类线粒体基因组中的突变会导致多种人类疾病，例如，癌症、心脏病、糖尿病等。
本发明提供了使用线粒体基因组中的体细胞突变来推断细胞谱系的方法，所述线粒体基因组中的体细胞突变也可以用作细胞的遗传条形码(genetic barcode)。在本发明的方法中，线粒体基因组中的体细胞突变是通过对单细胞转录组富集线粒体cDNA并进行二代测序和三代测序确定的，基于细胞中突变的存在对细胞进行聚类，由此推断细胞的谱系。本发明的方法可以确定线粒体DNA中的缺失(deletion)、插入(insertion)、重复(duplication)和倒位(inversion)等突变。
本发明提供了确定线粒体转录组中的变异和确定多细胞真核生物中的细胞谱系的方法，所述方法的一个示例性流程图如图1所示，其涉及以下方面。
1) 利用酶法或者物理方法制备多细胞真核生物的单细胞悬液。
2) 捕获单细胞，然后裂解细胞，使细胞mRNA自细胞中游离。
在一个实施方案中，使用10x Genomics公司的仪器捕获细胞；将获得的细胞按照10x Genomics 操作说明书进行上机，可捕获的细胞数为500个细胞-10000个细胞。该方案也同样适用于其他单细胞WTA(Whole transcriptome analysis)捕获平台如Drop-seq、InDrop-seq、BD Rhapsody及smart-seq等。
3) 逆转录细胞mRNA，获得3’端或5’端细胞条形码标记的全长cDNA。
在一些实施方案中，获得如图2所示的3’/5’细胞条形码标记的全长cDNA。
在一些实施方案中，使用了来自10x Genomics公司的产品制备单细胞转录组，利用微流控技术将单细胞与反应试剂包裹形成GEMs(Gel bead in emulsion)，每一个GEM里面除了包含一个细胞和试剂外，还具备带有双层标签的Gel beads oligo，第一层标签是10x barcode，用来标记各个细胞，即，每个GEM里面的barcode是唯一的；第二层标签是用来进行表达定量的UMI，用来标记逆转录的cDNA分子。例如，使用10x Chromium Next GEMSingle Cell 5' Kit v2试剂盒(10x Genomics公司，目录号：PN-1000263)，在10xChromium Controller上机结束后会形成GEMs，按照相应的说明书操作流程将GEMs进行RT反应程序孵育及cDNA扩增，即可获得带有单细胞barcode的全长cDNA。具体地，在GEM生成后，将凝胶珠(Gel Bead)溶解并使任何共同分配在该凝胶珠中的细胞裂解。释放出包含(i)Illumina R1序列(结合read 1测序引物)、(ii) 16 nt 10x Barcode、(iii) 10 nt UMI和(iv) 13 nt TSO的寡核苷酸，将其与细胞裂解物和含有逆转录 (RT)试剂和poly(dT) RT引物的Master Mix混合。孵育GEM，获得反转录自聚腺苷酸化mRNA的包含10x Barcoded的全长cDNA。
通过所述步骤1)至步骤3)可以获得细胞cDNA，其包含核基因转录组和线粒体基因转录组的cDNA。在一个实施方案中，获得的细胞全长转录组cDNA如图2所示，其包含条形码(barcode)标记。该部分适配于市面上大多数单细胞测序技术，包括10x genomics 3’RNAkit、10x genomics 5’RNA kit、10x genomics multiome kit、smart-seq、smart-seq2、smart-seq3、drop-seq、indrop、microwell-seq等可产生单细胞经标记的全长cDNA的产品或者技术。
4) 设计线粒体特异性引物，通过PCR富集线粒体cDNA。
在一些实施方案中，当步骤3)中的细胞条形码序列标记在cDNA的5’端时，对于每个线粒体基因设计一条或者多条位于cDNA 3’端的特异性引物，利用3’端特异性引物与5’端通用PCR引物进行该线粒体基因的PCR反应。在特异性引物结合位置方面要求尽量靠近基因的3’端，同时所述引物具有较高的特异性，能够区分线粒体基因组和核基因组序列。
优选地，设计并合成针对线粒体DNA(mtDNA)上各线粒体基因的引物。线粒体DNA包含如下线粒体基因：13个编码基因，2个rRNA基因，以及22个tRNA基因。本发明的方法由于仅需要对每个线粒体基因合成1条下游反向引物，由此降低了成本。由于所述引物是特异性靶向线粒体转录组的引物，因此具有富集线粒体基因的作用。
又在一些实施方案中，当步骤3)中的细胞条形码序列标记在cDNA的3’端时，对于每个基因设计一条或者多条位于cDNA 5’端的特异性引物，利用5’端特异性引物与3’端通用PCR引物进行该基因的PCR反应。在特异性引物结合位置方面要求尽量靠近基因的5’端，同时所述引物具有较高的特异性，能够区分线粒体基因组和核基因组序列。
优选地，设计并合成针对线粒体DNA(mtDNA)上各线粒体基因的引物。线粒体DNA包含如下线粒体基因：13个编码基因，2个rRNA基因，以及22个tRNA基因。由于仅需要对每个线粒体基因合成1条上游反向引物，由此降低了成本。由于所述引物是特异性靶向线粒体转录组的引物，因此具有富集线粒体基因的作用。
在一些实施方案中，当步骤3)中的细胞条形码序列标记在cDNA的5’端且获得的全长转录组cDNA如图2所示时，为获得带有单细胞条形码(Single cell barcode)标记的线粒体全长cDNA，线粒体基因的靶向富集用上游引物使用的是通用PCR引物。在一个实施方案中，所述通用PCR引物为Read 1引物序列(5’-CTACACGACGCTCTTCCGATCT-3’(SEQ ID NO：28))。因此仅需要设计每个线粒体基因的下游引物以特异性富集线粒体基因的全长cDNA。在一个实施方案中，设计的下游反向引物的位置尽量靠近mtDNA基因的3’端200bp范围内，例如，尽量靠近mtDNA基因的3’端约200bp范围内，例如，约150bp范围内、约100bp范围内、约50bp范围内、约30bp范围内，避开SNP位点，且能够与核基因组进行区分，以获得mtDNA基因的全长转录组。
在一些实施方案中，本发明设计了如表1所示的具有线粒体基因特异性的下游反向引物。对于每种mtRNA反转录获得的cDNA，仅需设计一条针对其的线粒体特异性反向引物(reverse primer)即可特异性扩增该cDNA。鉴于线粒体基因与细胞核基因具有较高的同源性，因此，需要设计仅对特定的线粒体基因序列具有更强的特异性，从而能够与核基因组进行区分的引物。
现有技术中没有公开对线粒体基因的全长cDNA进行特异性富集，然后用于测序的技术方案。例如，MEASTER技术(参见 Miller, T. E.等人, "Mitochondrial variantenrichment from high-throughput single-cell RNA-seq resolves clonalpopulations", Nat Biotechnol (2022). https://doi.org/10.1038/s41587-022-01210-8)需要设计65对引物(参见 Miller, T. E.等人，第14页图1)以满足短序列测序时覆盖全长线粒体cDNA的要求。与MEASTER技术针对每个线粒体基因设计多个引物相比，本发明的引物能够特异性扩增出全长线粒体cDNA，而不是如MEASTER技术那样仅扩增约250 bp长的线粒体cDNA片段。
进一步地，在富集线粒体cDNA时，本发明人为了解决构建的测序文库中各线粒体转录本信息可能存在覆盖度不均一的问题，对步骤4)中用于PCR扩增各线粒体cDNA的富集用特异性引物的比例进行调整，旨在得到数量均等的线粒体各转录本来源信息，以防止引入PCR偏差(PCR bias)，从而有利于后续分析的准确性。
在一些实施方案中，将用于PCR扩增各线粒体cDNA的富集用特异性引物的比例参照HPA数据库https://www.proteinatlas.org/中线粒体各基因的表达水平进行调整。
在一些实施方案中，当通过PCR扩增特异性富集各线粒体cDNA，例如，富集选自线粒体13个编码基因、2个rRNA基因的cDNA时，将各特异性引物依照各基因表达水平来调整所述特异性引物的使用浓度，例如通过HPA数据库https://www.proteinatlas.org/调取线粒体各基因的表达水平，由此设定针对不同基因的特异性引物的添加量，具体如下表A所示。
表A. 线粒体基因表达量与扩增用引物的终浓度范围
待扩增的线粒体cDNA所对应的基因名称 线粒体基因表达量(覆盖深度)占比(%) 扩增用引物的终浓度范围uM
MT-RNR1 0.07~0.55 0.07～1.02
MT-RNR2 0.25~2.62 0.05～1.7
MT-ND1 0.19~2.12 0.3～2.68
MT-ND2 0.26~1.81 0.27～0.48
MT-CO1 0.83~4.86 0.05～0.12
MT-CO2 0.44~3.09 0.06～0.47
MT-ATP8 0.42~4.15 0.16～0.5
MT-ATP6 0.5~3.72 0.08～0.34
MT-CO3 0.57~3.93 0.07～0.2
MT-ND3 0.12~1.84 0.09～0.4
MT-ND4L 0.23~2.32 0.17～0.57
MT-ND4 0.39~3.07 0.3～0.76
MT-ND5 0.12~1.23 0.3～1.55
MT-ND6 0.13~1.63 0.5～2.82
MT-CYB 0.24~2.11 0.2～0.6
使用表A所示的针对各线粒体基因的扩增用引物的添加量对样本进行线粒体cDNA富集。经过测序分析不同线粒体基因的覆盖深度，由此对各特异性引物的添加量优化而获得各特异性引物的经调整的添加量。对其他样本同样地获得各特异性引物的经调整的添加量，也能实现稳定且相对均一的线粒体基因富集。
按所获得的针对线粒体cDNA的各特异性引物经调整的添加量将多条特异性引物(例如，针对线粒体13个编码基因、2个rRNA基因的cDNA的15条特异性下游反向引物)混合在1个PCR管中，补加无核酸酶的ddH2O至总体积为100ul，获得下游反向引物Mix。在一个实施方案中，如下进行线粒体转录组扩增PCR1，为1管反应/样本，PCR1扩增程序为：
预变性：95℃ 3min
6个循环(98℃ 30s；65℃ 30s；72℃ 3min)
终延伸：72℃ 5min
保存温度：4℃。
线粒体转录组扩增PCR1反应结束之后，对PCR1反应产物进行纯化，纯化产物作为后续反应的模板。在一些实施方案中，使用1.0XSPRI磁珠(Solid Phase ReversibleImmobilization，Beckman，B23317)或1.0X Ampure XP磁珠(Beckman，A63880)等进行纯化，得到靶向富集的线粒体cDNA。
与本发明可以便捷、省时地得到引物Mix相比，在使用MAESTER技术(参见Miller,T. E.等人，见上文)时，将65条引物溶解后，按固定比例配制成12管不同的下游反向引物Mix(参见Miller, T. E.等人，第14页，图1)，通过补加不同体积的缓冲液，保证每条引物在扩增用混合物中的浓度为1μM。随后，每个样本需要配制12管扩增用混合物进行线粒体转录组扩增的PCR1反应。MAESTER技术的PCR1反应为：
预变性：95℃ 3min
6个循环(98℃ 20s；65℃ 15s；72℃ 3min)
终延伸：72℃ 5min
保存温度：4℃。
MAESTER技术在PCR1反应结束后，基于获自测序的覆盖率而凭经验确定需要从12管PCR产物分别取不同体积的产物量混合为1管(参见Miller, T. E.等人，第14页，图1)进行下一步纯化，纯化产物作为PCR2的反应模板，进行PCR2扩增反应。MAESTER技术的PCR2反应为：
预变性：95℃ 3min
8个循环(98℃ 20s；60℃ 30s；72℃ 3min)
终延伸：72℃ 5min
保存温度：4℃。
因此，与采用MAESTER技术相比较，本发明的步骤4)具有耗材用量少，节省成本，配制时间短的优点；以及在样本数量较多时，操作更快，且更节省时间的优点。
5) 对经纯化的步骤4)的线粒体转录组富集产物实施下述i)和ii)，进行二代Bulk线粒体转录组片段测序和三代单细胞线粒体全长转录组测序。
i) 取一部分步骤4)的线粒体转录组富集产物进行三代单细胞线粒体全长转录组测序。
在一个实施方案中，如下实施步骤5)的i)：取一部分步骤4)的线粒体转录组富集产物进行单一样本长读长三代测序建库并实施三代测序。例如，按照ONT连接测序试剂盒(SQK-LSK110)操作说明书执行末端修复和测序接头连接，实施单一样本三代测序文库构建，然后上机测序。
在一个实施方案中，如下实施步骤5)的i)：取一部分步骤4)的线粒体转录组富集产物作为来自目的样本的线粒体转录组富集产物，与作为来自其他非目的样本的线粒体转录组富集产物混合，进行多个样本(即，包括目的样本和非目的样本)的线粒体转录组长读长三代测序建库并实施三代测序。具体地，通过PCR扩增的方式向线粒体转录组富集产物的混合物(即，此轮PCR的模板)添加样本索引(Sample index)序列，由此，获得带样本索引序列的多个样本单细胞线粒体全长转录组扩增产物，然后，进行长读长混合样本三代测序建库并实施混合样本三代测序。在一些实施方案中，混合多个样本，例如，混合1、2、3…、N个样本的线粒体转录组富集产物，其中N是正整数。例如，按照ONT连接测序试剂盒(SQK-LSK110)操作说明书执行末端修复和测序接头连接，实施文库构建，然后上机测序。
ii) 取另一部分步骤4)的线粒体转录组富集产物进行Bulk线粒体转录组片段构建，例如，使用NEBNext® Ultra™ II FS DNA Library Prep Kit for Illumina(NEB，E7805)按照产品说明书分别执行片段化、末端修复和加A尾、接头(adaptor)连接、及PCR扩增，构建二代测序文库，然后使用Illumina测序平台对文库进行测序，以获得Bulk线粒体转录组短读长二代测序数据。
具体而言，将步骤4)获得的线粒体cDNA富集产物进行片段化，筛选满足二代测序读长的片段范围进行片段选择，然后通过接头连接及PCR扩增的方式得到所有的Bulk线粒体转录组二代文库，经二代测序获得的数据中是单个样本的一群细胞的线粒体转录组短读长的数据。在本发明中，将对单个样本中多个细胞的具有细胞条形码(barcode)序列和UMI序列的线粒体cDNA富集产物进行二代建库和测序也简称为Bulk线粒体转录组二代建库和测序，所获得的测序数据也称为“Bulk线粒体转录组二代测序数据”，与之相关的二代测序也称为“短读长Bulk线粒体转录组二代测序”、“二代Bulk线粒体转录组测序” 、“Bulk线粒体转录组短读长测序”或‘“Bulk线粒体转录组短读长测序”。
在一个实施方案中，步骤5)的i)采用的是ONT(Oxford Nanopore Technology公司)三代长读长的测序技术，可以用于分析单细胞全长线粒体cDNA，测序数据准确度也随着测序技术的发展和分析方法的优化得到改善，结合步骤5)的ii)Illumina Novaseq6000平台的二代测序数据共同分析，能够得到更多更准确的数据。
现有技术例如MAESTER采用的是二代测序技术分析线粒体基因组突变，虽然二代测序数据具有较高的测序数据准确度，成本低，但同时也存在读长短，信息丢失的弊端，因此，本发明的方法克服了现有技术中的不足。
6) 任选地，以步骤3)中制备的单细胞全长cDNA转录组为模板，进行二代短读长单细胞转录组测序。
在一个实施方案中，按照10x Genomics CG000331 Rev B中的Gene Expression(GEX) Library Construction制备单细胞转录组文库，其流程如图4所示，将单细胞全长转录组cDNA进行片段化，然后经过末端修复、加A碱基及接头连接，即可在后续PCR过程中富集含有cell barcode及UMI信息的基因表达片段。也可以使用其他能够制备单细胞转录组二代测序文库的试剂。单细胞转录组文库随后进行二代短读长测序，与步骤5)的ii) 短读长Bulk线粒体转录组二代测序数据及步骤5)的i)三代长读长测序数据共同分析。
所述步骤6)是任选实施的步骤，用于获得靠近转录组3’端或5’端序列信息 ，以进一步校正三代测序reads中的barcode和UMI的信息。
另外，所述步骤步骤5)的i)和步骤5)的ii)可以以任何顺序实施；所述步骤6)也可以在步骤4)之前实施，并且所述步骤4)至步骤6)可以以不同的组合二代测序和三代测序的方式进行实施。
在一些实施方案中，如图3和图4所示实施所述步骤4)至步骤6)。在一个具体实施方案中，图3中的Bulk线粒体短读长序列(Short reads sequence)二代测序文库是以线粒体转录组富集产物为模板，使用NEBNext® Ultra™ II FS DNA Library Prep Kit forIllumina(NEB，E7805)试剂盒构建的，因此，所获得的Bulk线粒体转录组二代测序数据中仅含有线粒体基因的数据，并且通过本发明的线粒体转录组富集步骤，可以确保样本的线粒体转录组文库中各片段的二代测序深度一致，更利于生物信息的分析。在一个具体实施方案中，图4中的单细胞转录组短读长序列(Short reads sequence)二代测序文库是以单细胞全长cDNA转录组为模板，使用10x Genomics CG000331 Rev B中的Gene Expression(GEX) Library Construction Kit, PN-1000190构建的。由于单细胞转录组二代测序数据中主要包含的是核基因表达的信息，其中也包含部分线粒体数据，但通常是否包含所有线粒体基因数据是未知的，且样本的线粒体各基因测序深度也难以一致。在常规单细胞转录组分析流程中，在并不需要线粒体数据时，线粒体数据往往作为背景噪音被去除。因此，如果仅采用单细胞转录组二代测序数据进行线粒体信息分析，可能会存在有效的线粒体转录组二代测序数据极少且不全面等技术问题。将本发明的所述步骤4)至步骤6)组合实施克服了该现有技术中存在的技术问题。
7) 对步骤5)和任选地步骤6)生成的测序数据进行生物信息分析。
本发明利用长读长三代测序获得的读长更长、且得到全长转录组信息的优点，对单细胞线粒体全长转录组进行了三代测序(所获得的reads中也包含barcode/UMI信息)。另一方面，由于长读长三代测序，特别是Oxford Nanopore测序的单碱基准确率低，在barcode、UMI及全长转录组都会出现测序错误，因此，在一个实施方案中，使用步骤5)的ii)中得到的Bulk线粒体转录组二代测序数据，对步骤5)的i)的长读长三代测序单细胞线粒体转录组的barcode/UMI和变异信息(SNP/InDel)进行校正，以得到携带有高准确度的barcode/UMI、SNP/InDel的全长单细胞线粒体转录组序列。
又在一个实施方案中，使用步骤5)的ii)中得到的Bulk线粒体转录组二代测序数据和进一步地使用步骤6)中得到的单细胞转录组短读长二代测序数据，对步骤5)的i)的长读长三代测序单细胞线粒体转录组的barcode/UMI和变异信息(SNP/InDel)进行校正，以得到携带有高准确度的barcode/UMI、SNP/InDel的全长单细胞线粒体转录组序列。图5例示了该技术方案中的一种生物信息分析框架图，具体如下所述。
a模块. 长读长测序单细胞线粒体转录组数据分析，以期获得单细胞线粒体全长转录组水平的变异信息(Nanopore 测序因为单碱基准确性低，所以会携带低准确性barcode/UMI序列、低准确性 SNP/InDel位点)：
(a) 采用Smith-Waterman局部比对算法(https://github.com/mengyao/Complete-Striped-Smith-Waterman-Library) 将样本index序列与read两端 300bp的序列进行比对，允许4 个碱基错配，将 read分配到具有最佳比对得分的样本；
(b) 单个样本采用 minimap2软件(https://github.com/lh3/minimap2) 将reads比对到参考基因组(例如，GRCh38)；
(c) 提取read 5’端未比对到参考基因组的序列(定义为softclip5, 包含例如P5/Truseq Read/TSO/barcode/UMI等通过10X Genomics技术添加的序列);
(d) 在比对完成的bam文件中提取转录本比对到参考基因组的起始位置plrs。
b模块. 短读长测序单细胞转录组数据提供高准确性barcode/UMI序列：
(a) 使用 fastp软件(https://github.com/OpenGene/fastp)选择默认参数对单细胞转录组短读长测序下机的原始数据进行质控，去除低质量 reads；
(b) 质控后的数据使用 Cell Ranger进行处理，所述Cell Ranger 是 10Xgenomics公司提供的一套针对单细胞 RNA 测序输出结果进行比对、定量、聚类及基因表达分析的分析流程，由此获得高质量的barcode/UMI序列；
(c) 质控后的序列使用 STAR软件(https://github.com/alexdobin/STAR)比对到参考基因组(例如，GRCh38)，提取转录本比对到参考基因组的起始位置psrs。
c模块. 短读长测序Bulk 线粒体转录组数据提供高准确性的barcode/UMI序列及基因水平的变异：
(a) 使用 fastp软件，选择默认参数对短读长测序Bulk 线粒体转录组的原始数据进行质控，去除低质量 reads；
(b) 采用STAR软件将序列比对到参考基因组(例如，GRCh38)，随后采用GATK4Mitochondria Pipeline (https://github.com/gatk-workflows/gatk4-mitochondria-pipeline)进行变异检测，得到高准确性的SNP/InDel;
(c) 质控后的数据使用 Cell Ranger进行处理，由此获得高质量的barcode/UMI序列。
d模块. 长读长测序单细胞线粒体转录组数据与短读长测序单细胞转录组数据整合，用来校正线粒体全长转录组的barcode/UMI序列：
(a) 根据基因名称将单细胞短读长测序的barcode/UMI 及长读长测序softclip5进行分组。同一基因的barcode/UMI与softclip5具有相同的转录本来源，进行比对时，可以减少比对时间，提高比对效率；
(b) 将长读长测序单细胞线粒体转录本数据的softclip5序列分为P5/TruseqRead/TSO等模板序列(有效 reads应携带这些序列，且所有 reads的这些序列相同)与barcode/UMI序列。将模板序列作为锚点序列，锚点序列允许更多的错配(15%的错配率)，可以获得更多有效 reads，同时定位到 barcode/UMI序列在softclip5的位置。根据此位置，获得softclip5的barcode/UMI序列，进而与短读长测序单细胞转录本的相同基因来源的barcode/UMI序列进行两两比对；
(c) 整合长读长、短读长barcode/ UMI 比对的结果(mismatches_gaps)与转录本在参考基因组的起始位置(pos_diff =|psrs-plrs|)， 比较不同mismatches_gaps与pos_diff下的长读长测序reads数目，最终选择(pos_diff<30 & mismatches_gaps<=3) |(pos_diff<20 & mismatches_gaps<=4) |(pos_diff<10 & mismatches_gaps<=5)作为长读长测序与短读长测序 barcode/UMI 映射的标准。
e模块. 对d模块得到的barcode/UMI校正后的单细胞全长转录本序列与c模块得到的高准确性的SNP/InDel进行整合，预期得到高准确性的SNP/InDel的全长转录本：
(a) 以短读长测序 Bulk 线粒体转录本数据得到的SNP/InDel为依据，确定 SNP/InDel 在长读长测序线粒体转录本单条reads 中是否存在；
(b) 携带相同 barcode/UMI的不同 reads 来自于相同细胞的相同转录本，对这样的reads进行SNP/InDel的相互校正,得到高准确性的变异结果。
f模块. 根据上述各模块得到的携带高准确度 barcode/UMI、高准确度 SNP/InDel的全长转录本计算单个细胞(barcode)水平 SNP/InDel的变异频率，构建细胞-变异(SNP/InDel) 矩阵，并根据细胞间的欧式距离 (Euclidean distance)，进行层次聚类分析，据此推断细胞的进化谱系。
本发明中采用的线粒体测序技术是二代测序与三代测序的结合。在线粒体转录组三代测序数据中需要单细胞条形码信息，以便进行三代单细胞线粒体全长转录组的序列分析。在线粒体三代测序数据中，单细胞条形码(Single cell barcode)能够将三代测序数据归类到单细胞的水平。另外，在Bulk线粒体二代测序数据中能够提供高准确性单细胞变异位点(例如，SNP/InDel位点)和单细胞条形码(Single cell barcode)信息，因此，本发明的单细胞线粒体转录组的三代测序数据可以被Bulk线粒体转录组二代测序数据校正。任选地，对单细胞转录组靠近3’端或5’端的序列进行二代测序，也可以进一步校正三代测序reads中的barcode和UMI的信息。
本发明中将步骤5)的ii)的Bulk线粒体转录组二代短读长测序与将步骤5)的i)单细胞线粒体全长转录组的长读长三代测序(例如，ONT纳米孔测序技术、Pacbio单分子测序技术或者其他长读长测序技术)组合，任选地，还与步骤6)的单细胞转录组靠近3’端或5’端的序列的短读长二代测序组合，能够将短读长的二代测序技术与长读长的三代测序技术的各自优势结合起来。
对于同一样本同时执行Bulk线粒体转录组二代测序和单细胞线粒体全长转录组三代测序，以及任选地，还执行单细胞转录组靠近3’端或5’端的序列的短读长二代测序，具有如下优势。其一，能够准确地获得线粒体全长转录组信息；其二，利用同一样本中多个单细胞的Bulk线粒体转录组进行的二代测序结果作为线粒体变异的参照，可显著地降低三代测序变异位点的假阳性，能够降低三代测序所需测序数据量和测序成本。
现有的MAESTER技术仅利用了二代测序平台。二代测序准确度高，但读长短，往往会丢失全长信息，MAESTER技术可能会造成全长线粒体转录组覆盖不完全、测序深度不均一的问题，进而漏检部分真实变异。
本发明的方法将二代测序平台和三代测序平台的测序结果结合起来进行分析，二代测序数据可以用于纠正三代测序数据，降低线粒体变异位点的假阳性；二代测序数据结合三代测序数据的全长特点，可以获得更多线粒体转录组的信息。
因此，本发明的二代测序平台和三代测序平台的组合能获得更加准确的线粒体全长转录组序列信息。
本发明的试剂盒
本发明提供了试剂盒。试剂盒可以包括一种或多种实施本发明方法的物品和/或试剂。
在一个实施方案中，试剂盒用于制备线粒体转录组测序文库，其包含含有通用PCR引物结合位点、用于确定转录组来源细胞的细胞条形码(barcode)序列和用于确定所述来源细胞中的各转录本的唯一分子标识符(UMI)序列的核苷酸序列，用于添加至线粒体RNA反转录获得的cDNA上，以及与线粒体RNA反转录获得的cDNA退火的一个或多个引物组，且每个引物组用于一条线粒体核酸链的靶向扩增，由此扩增一条或多条线粒体核酸链。在一个实施方案中，所述每个引物组中的一条引物是针对线粒体转录本的3’端特异性引物或5’端特异性引物，例如，表1所示的针对线粒体转录本的一条或多条3’端特异性引物；所述每个引物组中的另一条引物是通用PCR引物，用于结合线粒体RNA反转录获得的cDNA上的通用PCR引物结合位点。
在进一步的实施方案中，试剂盒还包含向扩增的一条或多条线粒体核酸添加至少一个样本索引的核苷酸序列(例如，所述添加至少一个样本索引的核苷酸序列选自至少一对SEQ ID NO: 16和17、SEQ ID NO: 18和19、SEQ ID NO: 20和21、SEQ ID NO: 22和23、SEQID NO: 24和25、SEQ ID NO: 26和27样本索引引物)，用于多样本三代测序文库的构建。
试剂盒还可以包含用于产生线粒体转录组测序文库所需的一种或多种其他组分。例如，引物延伸或扩增的酶、缓冲剂等。
试剂盒的组分通常存在于合适的包装材料中。
术语“包装材料”是指用于容纳试剂盒内容物的一种或多种物理结构。包装材料通过常规方法构造，一般提供无菌的、无污染的环境。例如，保护内容物免于外部环境的合适容器如小瓶。
包装材料可具有标签，其表明可以使用产生线粒体转录组测序文库的组分。
另外，包装材料包含说明如何使用试剂盒中的材料的说明。“使用说明”通常包括描述试剂浓度或至少一种测定方法参数的有形表达，如待混合的试剂和样品的相对量、试剂/样品混合物的维持时间、温度、缓冲条件等等。
识别线粒体测序数据中线粒体变异和确定多细胞真核生物中细胞谱系的装置
本发明提供了基于二代测序和三代测序来识别单细胞线粒体测序数据中线粒体变异和确定多细胞真核生物中细胞谱系的装置，所述装置被配置成能够实现本发明的确定单细胞线粒体转录组中的变异的方法和/或确定多细胞真核生物中细胞谱系的方法，例如，所述装置被配置成：
－接收输入，所述输入包括二代测序数据和三代测序数据，例如，单细胞线粒体全长三代测序数据和Bulk线粒体转录组二代测序数据，以及任选地，单细胞转录组二代测序数据，
－对输入的各测序数据进行质控和过滤，并进行生物信息学分析。
在一些实施方案中，所述装置包含一个或多个以下模块：
模块1. 接收单细胞线粒体全长三代测序数据，并映射到靶线粒体DNA样本；
模块2. 接收Bulk线粒体转录组二代测序数据，获得高准确性的基因水平的变异和barcode/UMI序列，例如，获得高准确性单细胞变异位点(例如，SNP/InDel位点)和barcode/UMI序列；
模块3. 用模块2中获得的高准确性barcode/UMI序列校正模块1的barcode/UMI序列， 再与模块2得到的高准确性的变异位点(例如，SNP/InDel位点)进行整合，获得高准确性的线粒体变异位点(例如，SNP/InDel位点)；以及
模块4. 基于模块3得到的高准确性的线粒体变异位点(例如，SNP/InDel位点)的存在对细胞进行聚类，并推断细胞的谱系；优选地，通过计算单细胞变异位点(例如，SNP/InDel位点)的变异频率，构建细胞-变异(例如，SNP/InDel变异) 矩阵，由此获得细胞的进化谱系。
任一所述模块1-4可以通过计算机执行。因此，本发明还提供经编程可以执行任一所述模块1-4的计算机。计算机典型地包括：与计算机通讯界面连接的CPU，系统记忆(RAM),非暂时性存储器(ROM)，和一个或多个其他存储装置如硬板、软盘、CD ROM drive。计算机还可以包括展示装置，例如打印机、CRT监测器或LCD展示器，以及输入装置，例如键盘、鼠标、笔、触摸屏或语音激活系统。输入装置可以接收数据，例如通过界面直接从测序仪接收。
在一些实施方案中，所述装置包含一个或多个以下模块：
模块1. 接收单细胞线粒体全长三代测序数据，并映射到靶线粒体DNA样本；
模块2. 接收单细胞转录组二代测序数据并分析，获得高准确性barcode/UMI序列；
模块3. 接收Bulk线粒体转录组二代测序数据并分析，获得高准确性的基因水平的变异，从而获得高准确性单细胞变异位点(例如，SNP/InDel位点)；
模块4. 将所述模块1和模块2中分别获得的长、短读长单细胞数据barcode/UMI映射，获得高准确性barcode/UMI序列及较高准确性单细胞变异位点(例如，SNP/InDel位点)；
模块5. 对模块4得到的barcode/UMI校正后的单细胞全长转录组序列与模块3得到的高准确性的变异位点(例如，SNP/InDel位点)进行整合，获得高准确性的线粒体变异位点(例如，SNP/InDel位点)；以及
模块6. 基于模块5得到的高准确性的线粒体变异位点(例如，SNP/InDel位点)的存在对细胞进行聚类，并推断细胞的谱系；优选地，通过计算单细胞变异位点(例如，SNP/InDel位点)的变异频率，构建细胞-变异(例如，SNP/InDel变异) 矩阵，由此获得细胞的进化谱系。
任一所述模块1-6可以通过计算机执行。因此，本发明还提供经编程可以执行任一所述模块1-6的计算机。计算机典型地包括：与计算机通讯界面连接的CPU，系统记忆(RAM),非暂时性存储器(ROM)，和一个或多个其他存储装置如硬板、软盘、CD ROM drive。计算机还可以包括展示装置，例如打印机、CRT监测器或LCD展示器，以及输入装置，例如键盘、鼠标、笔、触摸屏或语音激活系统。输入装置可以接收数据，例如通过界面直接从测序仪接收。
本发明的方法、试剂盒和装置的应用
本发明的方法、试剂盒和装置能够用于单细胞线粒体相关突变的高通量检测。
常规的线粒体相关突变的检测包括：PCR-RFLP和DNA一代测序、芯片方法等，但这些方法通常集中在对少数常见的线粒体基因位点进行突变和缺失筛查，如对线粒体相关性糖尿病常见突变A3243G的筛查，但阳性检出率低，这是因为线粒体疾病具有量效现象，即小量的线粒体DNA突变往往不出现临床症状，但随着突变线粒体比例增高，才会出现临床表现，且临床严重程度可能和突变比例成正相关，而现有的常规检测方法对于单细胞的线粒体突变都不能有效地检测。本发明的方法、试剂盒和装置能够克服现有技术中存在的对线粒体相关突变的检测结果假阳性高、无法高通量检测的技术问题，故本发明对单细胞中线粒体转录组的准确检测具有非常重要的临床指导意义。
进一步地，本发明的方法、试剂盒和装置能够用于追踪细胞谱系。细胞谱系追踪在肿瘤和免疫系统研究中尤其具有重要意义。通过肿瘤细胞谱系追踪可以对肿瘤克隆的性质精确定义并为寻找靶点提供信息。
实施例
描述以下实施例以辅助对本发明的理解。不意在且不应当以任何方式将实施例解释成限制本发明的保护范围。
实施例1. 细胞条形码标记的全长cDNA的制备
将来自健康个体的5mL全血收集在含肝素的采血管中。通过聚蔗糖密度梯度离心分离外周血单个核细胞(PBMC)(供应商：迈顺生物，货号：P121111505C)。计数细胞，获得约10000个单细胞的细胞稀释液。
使用10x Chromium Next GEM Single Cell 5' Kit v2试剂盒(10x Genomics公司，目录号：PN-1000263)，芯片为Chromium Next GEM Chip K(10x Genomics公司，目录号：PN-1000287)，按照生产商的操作说明书所述，将获得的细胞稀释液加载到Chromium NextGEM Chip K芯片进样孔内。
将加载有细胞稀释液的芯片装载到10× Chromium Controller仪器(购自10×Genomics公司)上，运行程序，形成油包水乳液微滴(GEMs，Gel Beads-in-emulsion)。在所产生的GEMs中大部分(约90–99%)不含有细胞，剩余的GEMs中大部分含有仅一个细胞，GEMs中的凝胶微珠上的寡核苷酸包含(i) 一个Illumina R1序列(结合read 1测序引物)，(ii)10x条形码(10x Barcode)，具有16 nt长度，一个微珠对应于一种条形码；(iii) 一个10 nt唯一分子标识符(UMI)，所述UMI是一段随机序列，10个碱基长的UMI有100多万种序列的变化(410 = 1,048,576)；和(iv) 13 nt模板转换寡核苷酸(TSO)。
在GEMs产生后，通过10x Chromium Next GEM Single Cell 5' Kit v2试剂盒中的含有逆转录(RT)试剂和多聚(dT)逆转录引物的Master Mix逆转录细胞mRNA，并通过模板转换和转录物延伸获得如图2所示的10x 细胞条形码标记的全长cDNA。
对GEMs进行破油和纯化，然后对经纯化的10x 细胞条形码标记的全长cDNA使用10x Chromium Next GEM Single Cell 5' Kit v2试剂盒中的引物进行全转录组扩增PCR(whole transcriptome amplification (WTA) PCR)，其中所用引物结合单细胞转录组逆转录期间在各全长cDNA两端添加的共有5'和3'末端引物结合位点，正向引物为：5’-CTACACGACGCTCTTCCGATCT-3’(SEQ ID NO：28)，Poly-dT反向引物为：5’- Poly-dTAAGCAGTGGTATCAACGCAGAG-3’(SEQ ID NO：29)，由此获得各单细胞的全转录组扩增产物。
实施例2. 线粒体cDNA的富集
对一部分实施例1中制备的细胞条形码标记的全转录组全长cDNA进行了线粒体cDNA的富集。
对线粒体如下基因的cDNA进行了富集：MT-RNR1、MT-RNR2、MT-ND1、MT-ND2、MT-CO1、MT-CO2、MT-ATP8、MT-ATP6、MT-CO3、MT-ND3、MT-ND4L、MT-ND4、MT-ND5、MT-ND6、MT-CYTB。具体而言，设计并合成了下表1所示的15条下游反向引物。将15条下游反向引物各溶解至100μM，然后按表1所示引物终浓度条件混合在1管PCR管中，补加无核酸酶的ddH2O至总体积为100ul，即得到下游反向引物Mix；向该PCR管中添加上游引物，为10x Chromium NextGEM Single Cell 5' Kit v2试剂盒中的Read1引物序列(5’-CTACACGACGCTCTTCCGATCT-3’(SEQ ID NO：28))，并添加实施例1获得的WTA产物，实施线粒体转录组的PCR扩增。
表1. 用于特异性扩增线粒体cDNA的下游反向引物序列信息：
SEQ ID NO: 核苷酸序列(5’- 3’) 扩增的转录本的来源基因 引物终浓度uM
1 CCATGTTACGACTTGTCTCCTCTATAT MT-RNR1 0.1
2 GTATAATACTAAGTTGAGATGATATCATTTACGG MT-RNR2 0.05
3 GATTGTAATGGGTATGGAGACATATCATATAAG MT-ND1 0.3
4 CGTGGTAAGGGCGATGAGT MT-ND2 0.3
5 GTTCTTCGAATGTGTGGTAGGGTG MT-CO1 0.1
6 AGATTTTTAGGGGAATTAATTCTAGGACGA MT-CO2 0.2
7 TGAAGCGAACAGATTTTCGTTCA MT-ATP8 0.5
8 CTAGAAGTGTGAAAACGTAGGCTTG MT-ATP6 0.2
9 AGACATACAGAAATAGTCAAACCACATCTA MT-CO3 0.2
10 CTCATAGGCCAGACTTAGGGCTAG MT-ND3 0.3
11 GTCTAGGCCATATGTGTTGGAGAT MT-ND4L 0.4
12 GGATAGGAGGAGAATGGGGGATAG MT-ND4 0.3
13 GAGTAGGGTTAGGATGAGTGGG MT-ND5 0.3
14 TGTGGGTTTAGTAATGGGGTTTGT MT-ND6 0.5
15 TAGGGAGATAGTTGGTATTAGGATTAGGAT MT-CYTB 0.3
扩增线粒体转录组的PCR反应以1个PCR管/1个全血个体样本的线粒体cDNA实施，具体PCR扩增程序为：
预变性：95℃ 3min
6个循环(98℃ 30s；65℃ 30s；72℃ 3min)
终延伸：72℃ 5min
保存温度：4℃。
使用1.0X SPRI磁珠(Solid Phase Reversible Immobilization，Beckman，B23317)纯化PCR扩增产物，获得富集的线粒体转录组cDNA产物。
实施例3. 单细胞线粒体转录组的三代测序
将一部分实施例2获得的线粒体转录组富集产物进行长读长的高通量三代测序。
具体地，采用ONT(Oxford Nanopore Technology公司)的长读长三代测序技术，对实施例2获得的单细胞线粒体转录组构建Nanopore文库，使用Nanopore MinION测序仪(Oxford Nanopore Technology公司)对构建的Nanopore文库进行测序，获得如下结果。
长读长三代测序线粒体全长转录组的原始下机数据reads数目为95,123,378，采用 minimap2软件(https://github.com/lh3/minimap2) 将reads比对到参考基因组GRCh38，结果表明比对到参考基因组的reads数目为94,930,502，比对率为99.08%，目标区域reads数目为50,892,494，占比53.61%；在bam文件中提取softclip5，并提取plrs。
实施例4. 线粒体转录组的Bulk二代测序
将一部分实施例2获得的线粒体转录组富集产物进行短读长的高通量Bulk二代测序。
具体地，按照NEBNext® Ultra™ II FS DNA Library Prep Kit for Illumina(NEB，目录号：E7805)产品说明书对实施例2的单细胞线粒体转录组富集产物执行快速可靠的线粒体cDNA片段化、末端修复、接头连接及Bulk二代测序文库构建。
Bulk二代测序文库构建后，使用Illumina Novaseq6000平台上机测序，获得Bulk线粒体转录组短读长的二代测序数据。
Bulk线粒体转录组二代测序的原始下机reads数目为83,164,940，使用 fastp软件，选择默认参数对该下机reads进行质控, 经质控过滤后的剩余reads数目为78,310,624, 占原始下机reads数目的94.16%。
采用STAR软件将质控后的reads序列比对到参考基因组GRCh38，随后使用GATK4Mitochondria Pipeline进行线粒体变异检测，得到高准确性的SNP/InDel，结果如表2所示。
表2. Bulk线粒体转录组短读长二代测序后的变异信息统计
项目 对项目的描述 数值
raw_reads 原始reads数目 83,164,940
clean_reads 质控后reads数目 78,310,624
primary mapped 比对到参考基因组的reads数目 65,894,825
mean depth of variant 目标基因变异的平均覆盖深度 3847
total_variant_count 目标基因的变异数目 335
snp_count 目标基因SNP数目 196
indel_count 目标基因InDel数目 96
mean depth of variant[PASS] 目标基因变异的平均覆盖深度[高质量] 4110.24
total_variant_count[PASS] 目标基因的变异数目[高质量] 125
snp_count[PASS] 目标基因SNP数目[高质量] 78
indel_count[PASS] 目标基因InDel数目[高质量] 28
实施例5. 单细胞全转录组的短读长二代测序
对一部分实施例1中制备的细胞条形码标记的单细胞全转录组全长cDNA进行了短读长高通量二代测序。
使用10x Chromium Next GEM Single Cell 5' Kit v2试剂盒(10x Genomics公司，CG000331 Rev B)中的Gene Expression (GEX) Library Construction如图4所示制备单细胞转录组的短读长测序文库。
具体而言，将实施例1中制备的单细胞全长转录组cDNA进行片段化，然后经过末端修复、加A碱基及接头连接，即可在PCR过程中富集含有细胞条形码和UMI信息的基因表达片段。将制备的单细胞转录组二代测序文库随后使用Illumina Novaseq6000平台进行短读长二代测序，获得单细胞转录组的短读长二代测序数据。
单细胞转录组的短读长二代测序原始下机reads数目为672,118,966，使用 fastp软件，选择默认参数对该下机reads进行质控, 经质控过滤后剩余reads数目为672,118,966, 占原始下机reads数目的100%。
将该质控后的数据使用 Cell Ranger进行处理，统计结果如下表3所示。
表3. Cell Ranger 统计结果
项目 对项目的描述 数值
Estimated Number of Cells 捕获的细胞数目 10,929
Mean Reads per Cell 每个细胞的平均reads数目 30,749
Median Genes per Cell 每个细胞中位基因数目 1,614
Number of Reads reads数目 336,059,483
Valid Barcodes 有效 barcode 占比 89.70%
Reads Mapped to Genome 比对到参考基因组的比例 91.20%
Fraction Reads in Cells 细胞内 reads的比例 90.30%
Total Genes Detected 检测到的基因总数 22,447
Median UMI Counts per Cell 每个细胞的中位 UMI数目 4,107
将该质控后的reads序列使用 STAR软件比对到参考基因组GRCh38，提取转录组reads序列比对到参考基因组的起始位置psrs，由此，提供了高准确性的单细胞条形码/UMI序列和其后接的基因表达片段序列。
实施例6. 对测序数据的生物信息学分析
本实施例对通过二代测序和三代测序获得的测序数据进行了生物信息学分析。
实施例6.1. 利用线粒体转录组Bulk二代测序数据和单细胞线粒体全长转录组三代测序数据获得高准确度 barcode/UMI、高准确度 SNP/InDel的线粒体全长转录组序列
将实施例3获得的长读长三代测序单细胞线粒体转录组数据与实施例4获得的Bulk线粒体转录组短读长二代测序数据进行整合，旨在进一步校正所述长读长三代测序单细胞线粒体转录组的序列。
具体而言，根据线粒体基因名称将实施例4中的Bulk线粒体转录组短读长二代测序数据通过Cellranger软件获得的barcode/UMI和实施例3中的长读长三代测序单细胞线粒体转录组数据的softclip5进行分组，获得下表4。
表4. Bulk线粒体转录组短读长二代测序和单细胞线粒体转录组长读长三代测序的reads数目统计
目标线粒体基因 Bulk线粒体转录组短读长二代测序barcode/UMI数目 单细胞线粒体转录组长读长三代测序softclip5数目
MT-ATP6 20,434 7,722,652
MT-ATP8 95,735 3,079,504
MT-CO1 46,962 5,280,040
MT-CO2 149,669 9,610,192
MT-CO3 96,974 5,546,072
MT-CYB 82,807 4,978,790
MT-ND1 12,291 1,467,608
MT-ND2 31,765 2,756,783
MT-ND3 31,103 2,179,502
MT-ND4 26,954 2,075,683
MT-ND4L 58,752 2,113,493
MT-ND5 14,624 1,130,849
MT-ND6 7,983 790,002
MT-RNR1 191,186 4,702,531
MT-RNR2 408,461 3,811,227
采用blast软件，将长读长三代测序的单细胞线粒体转录本softclip5比对到短读长二代测序的Bulk线粒体转录本barcode/UMI，选择最佳比对的barcode/UMI。在进行该最佳比对时，不考虑错配(mismatches)和空位(gaps)。
整合短读长二代测序的Bulk线粒体转录本、长读长三代测序的单细胞线粒体转录本barcode/ UMI比对的结果(mismatches_gaps)与转录本在参考基因组的起始位置(pos_diff =|psrs-plrs|)，比较不同mismatches_gaps与pos_diff下的长读长测序reads数目(见表5所示)，最终选择(pos_diff<30 & mismatches_gaps<=3) |(pos_diff<20 &mismatches_gaps<=4) |(pos_diff<10 & mismatches_gaps<=5)作为长读长三代测序单细胞线粒体转录组与短读长二代测序Bulk线粒体转录组barcode/UMI映射的标准。
表5. 不同比对结果与位置差异条件下的短读长Bulk线粒体转录本、长读长测序单细胞线粒体转录本barcode/UMI匹配情况
过滤条件 长读长测序单细胞线粒体转录本softclip5数目 匹配到短读长Bulk线粒体转录本barcode/UMI的长读长测序单细胞线粒体转录本reads数目 占比
mismatches_gaps<=3 47,751,920 19,563,044 0.4097
pos_diff<10 & mismatches_gaps<=3 47,751,920 18,351,810 0.3843
pos_diff<20 & mismatches_gaps<=3 47,751,920 18,967,567 0.3972
pos_diff<30 & mismatches_gaps<=3 47,751,920 19,119,653 0.4004
mismatches_gaps<=4 47,751,920 21,538,332 0.4510
pos_diff<10 & mismatches_gaps<=4 47,751,920 19,528,885 0.4090
pos_diff<20 & mismatches_gaps<=4 47,751,920 20,226,372 0.4236
pos_diff<30 & mismatches_gaps<=4 47,751,920 20,403,109 0.4273
mismatches_gaps<=5 47,751,920 22,671,477 0.4748
pos_diff<10 & mismatches_gaps<=5 47,751,920 20,228,268 0.4236
pos_diff<20 & mismatches_gaps<=5 47,751,920 20,977,926 0.4393
pos_diff<30 & mismatches_gaps<=5 47,751,920 21,175,084 0.4434
(pos_diff<30 & mismatches_gaps<=3)| (pos_diff<20 &mismatches_gaps<=4) | (pos_diff<10 & mismatches_gaps<=5) 47,751,920 30,016,535 0.4414
将匹配到短读长测序Bulk线粒体转录本barcode/UMI的长读长测序单细胞线粒体转录本reads序列(即，经短读长测序Bulk线粒体barcode/UMI校正后的单细胞线粒体全长转录本序列)与实施例4获得的短读长Bulk测序线粒体转录本reads序列进行整合，并以实施例4获得的短读长Bulk测序线粒体转录本reads序列中的SNP/InDel为依据，获得包含高准确性SNP/InDel的长读长测序单细胞线粒体转录本序列。
由于携带相同barcode/UMI的不同reads来自于相同细胞的相同转录本，通过对这样的reads进行SNP/InDel的相互校正，最终获得28个高质量线粒体变异。
根据以上生物信息学分析获得的携带高准确度barcode/UMI、高准确度SNP/InDel的线粒体全长转录本序列计算单个细胞(即，barcode标记的单细胞)水平的变异频率，构建了10795个细胞－28个变异矩阵，并根据细胞间欧式距离进行层次聚类分析(图6)，据此推断出细胞的进化谱系(图7，不同颜色表示不同的细胞群分类)。
实施例6.2. 利用线粒体转录组Bulk二代测序数据、单细胞转录组二代测序数据和线粒体转录组三代测序数据获得高准确度barcode/UMI、高准确度SNP/InDel的线粒体全长转录组序列
将实施例3获得的长读长测序单细胞线粒体转录组数据与实施例5获得的短读长测序单细胞转录组数据进行整合，旨在进一步校正长读长测序单细胞线粒体转录组的barcode/UMI序列。
根据线粒体基因名称将实施例5中的短读长测序单细胞转录组数据中的barcode/UMI和实施例3中的长读长测序单细胞线粒体转录组数据的softclip5进行分组，获得下表6。
表6. 短读长测序单细胞转录本与长读长测序单细胞线粒体转录本reads数目统计
目标线粒体基因 短读长测序单细胞转录本barcode/UMI数目 长读长测序单细胞线粒体转录本 softclip5数目
MT-ATP6 325,189 7,722,652
MT-ATP8 164,661 3,079,504
MT-CO1 558,951 5,280,040
MT-CO2 246,860 9,610,192
MT-CO3 188,977 5,546,072
MT-CYB 158,675 4,978,790
MT-ND1 86,467 1,467,608
MT-ND2 72,650 2,756,783
MT-ND3 62,762 2,179,502
MT-ND4 85,837 2,075,683
MT-ND4L 78,658 2,113,493
MT-ND5 88,469 1,130,849
MT-ND6 37,381 790,002
MT-RNR1 574,902 4,702,531
MT-RNR2 2,069,266 3,811,227
采用 blast软件，将长读长测序的单细胞线粒体转录本softclip5比对到短读长测序的单细胞转录本barcode/UMI，选择最佳比对的barcode/UMI。在进行该最佳比对时，不考虑错配(mismatches)和空位(gaps)。
整合短读长测序单细胞转录本、长读长测序单细胞线粒体转录本barcode/ UMI比对的结果(mismatches_gaps)与转录本在参考基因组的起始位置(pos_diff =|psrs-plrs|)，比较不同mismatches_gaps与pos_diff下的长读长测序reads数目(见表7所示)，最终选择(pos_diff<30 & mismatches_gaps<=3) |(pos_diff<20 & mismatches_gaps<=4) |(pos_diff<10 & mismatches_gaps<=5)作为长读长测序单细胞线粒体转录本与短读长测序单细胞转录本barcode/UMI映射的标准。
表7. 不同比对结果与位置差异条件下的短读长测序单细胞转录本、长读长测序单细胞线粒体转录本barcode/UMI匹配情况
过滤条件 长读长测序单细胞线粒体转录本softclip5数目 匹配到短读长测序单细胞转录本barcode/UMI的长读长测序单细胞线粒体转录本reads数目 占比
mismatches_gaps<=3 47,751,920 28,008,938 0.5866
pos_diff<10 & mismatches_gaps<=3 47,751,920 26,047,223 0.5455
pos_diff<20 & mismatches_gaps<=3 47,751,920 26,991,453 0.5652
pos_diff<30 & mismatches_gaps<=3 47,751,920 27,225,384 0.5701
mismatches_gaps<=4 47,751,920 30,725,553 0.6434
pos_diff<10 & mismatches_gaps<=4 47,751,920 27,695,975 0.5800
pos_diff<20 & mismatches_gaps<=4 47,751,920 28,755,991 0.6022
pos_diff<30 & mismatches_gaps<=4 47,751,920 29,026,108 0.6079
mismatches_gaps<=5 47,751,920 32,475,194 0.6801
pos_diff<10 & mismatches_gaps<=5 47,751,920 28,722,588 0.6015
pos_diff<20 & mismatches_gaps<=5 47,751,920 29,860,299 0.6253
pos_diff<30 & mismatches_gaps<=5 47,751,920 30,160,527 0.6316
(pos_diff<30 & mismatches_gaps<=3) | (pos_diff<20 &mismatches_gaps<=4) | (pos_diff<10 & mismatches_gaps<=5) 47,751,920 30,016,535 0.6286
将匹配到短读长测序单细胞转录本barcode/UMI的长读长测序单细胞线粒体转录本reads序列(即，经短读长测序单细胞转录本barcode/UMI校正后的单细胞线粒体全长转录本序列)与实施例4获得的短读长Bulk测序线粒体转录本reads序列进行整合，并以实施例4获得的短读长Bulk测序线粒体转录本reads序列中的SNP/InDel为依据，获得包含高准确性 SNP/InDel的长读长测序单细胞线粒体转录本序列。
由于携带相同barcode/UMI的不同reads来自于相同细胞的相同转录本，通过对这样的reads进行SNP/InDel的相互校正，最终获得63个高质量线粒体变异。
根据以上生物信息学分析获得的携带高准确度 barcode/UMI、高准确度 SNP/InDel的线粒体全长转录本序列计算单个细胞(即，barcode标记的单细胞)水平的变异频率，构建了10871个细胞－63个变异矩阵，并根据细胞间欧式距离，进行层次聚类分析(图8)，据此推断出细胞的进化谱系(图9，不同颜色表示不同的细胞群分类)。
实施例7. 人源肝细胞肝癌的细胞谱系的确定
将来自人源肝细胞肝癌的细胞(中美冠科 Liver-Cancer，批号LI0050)制备单细胞悬液，计数细胞。除了所使用的细胞样本不同之外，与实施例1同样地制备细胞条形码标记的全长cDNA；与实施例2同样地进行线粒体cDNA的富集。
实施例7.1. 混合细胞样本长读长测序线粒体全长转录组
采用ONT(Oxford Nanopore Technology公司)三代长读长的测序技术，使用Nanopore MinION测序仪(Oxford Nanopore Technology公司)对本实施例的肝细胞肝癌的经富集的单细胞线粒体全长转录组进行测序。
将本实施例上述富集的单细胞线粒体全长转录组的肝癌样本(即，目标肝细胞肝癌样本)与自两例人源PBMC样本(迈顺生物，货号：P121111505C(注：该公司此货号仅表示提供的是人源PBMC样本，所述人源PBMC样本可以为自不同的人类个体(donor ID)获得的))和三例人源肝细胞肝癌样本(中美冠科，批号：LI0050)(注：该公司此货号仅表示提供的是人源肝细胞肝癌样本，所述人源肝细胞肝癌样本可以为自不同的人类个体(donor ID)获得的)分别富集的单细胞线粒体全长转录组血液和肝癌样本(共6个样本)混合在1个PCR管中，添加下述样本索引(sample index)引物，进行进一步的PCR扩增并建库，由此实施混合样本三代测序。样本索引的设计原理：首先，要满足引物设计的基本要求，碱基平衡，无发卡结构，无串联重复序列等，其次，在序列上需要满足与线粒体基因无互补序列，避免造成非特异性扩增。使用了下述样本索引(sample index)正向（F）和反向（R）引物：
SIF-1：ATAGTGATACTGACCACCGAGATCTACACATATACGCAGTCGACAACTTTCTTGCGCCATGGGACTACACGACGCTC （SEQ ID NO: 16）；
SIF-2：ATAGTGATACTGACCACCGAGATCTACACCGCACGACGTACAAACGGAATCGAGCGCCATGGGACTACACGACGCTC （SEQ ID NO: 17）；
SIF-3：ATAGTGATACTGACCACCGAGATCTACACGGGACAGAAGGGGACACAAGACTCGCGCCATGGGACTACACGACGCTC （SEQ ID NO: 18）；
SIF-4：ATAGTGATACTGACCACCGAGATCTACACTTAACAAGTATGATAGAATCCGAAGCGCCATGGGACTACACGACGCTC （SEQ ID NO: 19）；
SIF-5：ATAGTGATACTGACCACCGAGATCTACACACTCCTGGGTAAACCCTGGACAAGGCGCCATGGGACTACACGACGCTC （SEQ ID NO: 20）；
SIF-6：ATAGTGATACTGACCACCGAGATCTACACCAGAAAACCCCGCCTTTGCGAGAAGCGCCATGGGACTACACGACGCTC （SEQ ID NO:21）；
SIR-1：CATCAGCGAACCGGCATACGAGATGCACGAGCAGTCCCACGGTAACACTGTTGAACTCCAGTTCAGACGTGT （SEQ ID NO: 22）；
SIR-2：CATCAGCGAACCGGCATACGAGATTATTATGCTTCTTTGTTCCCTGAATGTTGAACTCCAGTTCAGACGTGT （SEQ ID NO: 23）；
SIR-3：CATCAGCGAACCGGCATACGAGATAGCGACTGTATGCTGTGCCTAGTTTGTTGAACTCCAGTTCAGACGTGT （SEQ ID NO: 24）；
SIR-4：CATCAGCGAACCGGCATACGAGATCCTTGCGACAGGGTTTCAACGCTTTGTTGAACTCCAGTTCAGACGTGT （SEQ ID NO: 25）；
SIR-5：CATCAGCGAACCGGCATACGAGATGTTAGATTGCACGATAGATGAAACTGTTGAACTCCAGTTCAGACGTGT （SEQ ID NO: 26）；
SIR-6：CATCAGCGAACCGGCATACGAGATTAAATACCATCTTCTTTCTACCTGTGTTGAACTCCAGTTCAGACGTGT （SEQ ID NO: 27）；
PCR反应体系
试剂 体积 终浓度
2X KAPA HiFi HotStart ReadyMix 25 ul 1X
SIF各引物(分别为10uM) 5 ul 1uM
SIR各引物(分别为10uM) 5 ul 1uM
DNA模板 (为线粒体富集产物) 15 µL
总体积 50 µL
PCR反应程序：
预变性：98℃ 3min
8个循环(98℃ 20s；54℃ 30s；72℃ 3min)
终延伸：72℃ 5min
保存温度：4℃。
将PCR产物使用0.8X SPRI磁珠进行纯化。使用Nanopore MinION测序仪(OxfordNanopore Technology公司)进行测序。获得如下结果：长读长测序混合样本线粒体转录组的原始下机数据reads数目为100,532,000。采用Smith-Waterman局部比对算法将样本index序列与read两端 300bp的序列进行比对，允许 4 个碱基错配，将read分配到具有最佳比对得分的样本，从而将本实施例中的目标肝细胞肝癌样本从混合的6个样本中拆分得到reads数目为16,375,154(占混合样本reads数目比为16.29%)。
采用minimap2软件(https://github.com/lh3/minimap2) 将拆分完成后的reads比对到参考基因组GRCh38，结果表明：比对到参考基因组的reads数目为16,317,110，比对率为99.64%，目标区域reads数目为16,212,816 ，占比99.36%。
提取bam文件的softclip5序列，并提取plrs。
实施例7.2. 短读长Bulk线粒体转录组二代测序
按照NEBNext® Ultra™ II FS DNA Library Prep Kit for Illumina(NEB，目录号：E7805)产品说明书对本实施例的肝细胞肝癌的经富集的单细胞线粒体全长转录组执行快速可靠的线粒体cDNA片段化、末端修复、接头连接及Bulk二代测序文库构建。
Bulk二代测序文库构建后，使用Illumina Novaseq6000平台上机测序，获得Bulk线粒体转录组短读长的二代测序数据。
Bulk线粒体转录组二代测序原始下机reads数目为66,946,990, 使用 fastp软件，选择默认参数对该下机reads进行质控, 经质控过滤后的剩余reads数目为59,857,212, 占原始下机reads数目的89.41%。
采用STAR软件将质控后的reads序列比对到参考基因组GRCh38，随后使用GATK4Mitochondria Pipeline进行线粒体变异检测，得到高准确性的SNP/InDel，结果如表8所示。
表8. Bulk线粒体转录组短读长二代测序后的变异信息统计
项目 对项目的描述 数值
raw_reads 原始reads数目 66,946,990
clean_reads 质控后reads数目 59,857,212
primary mapped 比对到参考基因组的reads数目 53,204,720
mean depth of variant 目标基因变异的平均覆盖深度 4,958
total_variant_count 目标基因的变异数目 321
snp_count 目标基因SNP数目 175
indel_count 目标基因InDel数目 115
mean depth of variant[PASS] 目标基因变异的平均覆盖深度[高质量] 4832.144
total_variant_count[PASS] 目标基因的变异数目[高质量] 104
snp_count[PASS] 目标基因SNP数目[高质量] 43
indel_count[PASS] 目标基因InDel数目[高质量] 40
实施例7.3. 单细胞全转录组的短读长二代测序
对本实施例制备的细胞条形码标记的单细胞全转录组全长cDNA进行了短读长高通量二代测序。
使用10x Chromium Next GEM Single Cell 5' Kit v2试剂盒(10x Genomics公司，CG000331 Rev B)中的Gene Expression (GEX) Library Construction如图4所示制备单细胞转录组文库。
具体而言，将本实施例制备的单细胞全长转录组cDNA进行片段化，然后经过末端修复、加A碱基及接头连接，即可在PCR过程中富集含有细胞条形码和UMI信息的基因表达片段。将制备的单细胞转录组二代测序文库随后使用Illumina Novaseq6000平台进行短读长二代测序，获得单细胞转录组的短读长二代测序数据。
单细胞转录组短读长二代测序的原始下机reads数目为423,788,290，使用 fastp软件，选择默认参数对该下机reads进行质控, 经质控过滤后剩余reads数目为423,788,290, 占原始下机reads数目的100%。
将该质控后的数据使用 Cell Ranger进行处理，统计结果如下表9所示。
表9. Cell Ranger 统计结果
项目 对项目的描述 数值
Estimated Number of Cells 捕获的细胞数目 18,187
Mean Reads per Cell 每个细胞的平均reads数目 11,651
Median Genes per Cell 每个细胞中位基因数目 929
Number of Reads reads数目 211,894,145
Valid Barcodes 有效 barcode 占比 91.70%
Reads Mapped to Genome 比对到参考基因组的比例 78.6%
Fraction Reads in Cells 细胞内 reads的比例 73.80%
Total Genes Detected 检测到的基因总数 24,037
Median UMI Counts per Cell 每个细胞的中位 UMI数目 2,001
将该质控后的序列使用 STAR软件(https://github.com/alexdobin/STAR)比对到参考基因组GRCh38，提取转录本比对到参考基因组的起始位置psrs，由此，提供了高准确性的单细胞条形码/UMI序列和其后接的基因表达片段序列。
实施例7.4. 对测序数据的生物信息学分析
将实施例7.1获得的长读长三代测序单细胞线粒体转录组数据与实施例7.3获得的短读长二代测序单细胞转录组数据进行整合，旨在进一步校正长读长测序单细胞线粒体转录组的barcode/UMI序列。
根据线粒体基因名称将实施例7.3中的短读长二代测序单细胞转录组数据中的barcode/UMI和实施例7.1中的长读长三代测序单细胞线粒体转录组数据的softclip5进行分组，获得下表10。
表10. 短读长测序单细胞转录本与长读长测序单细胞线粒体转录本reads数目统计
目标线粒体基因 短读长测序单细胞转录本barcode/UMI数目 长读长测序单细胞线粒体转录本 softclip5数目
MT-ATP6 325,189 1,862,460
MT-ATP8 164,661 1,623,516
MT-CO1 558,951 1,628,099
MT-CO2 246,860 824,904
MT-CO3 188,977 1,540,157
MT-CYB 158,675 742,738
MT-ND1 86,467 195,213
MT-ND2 72,650 1,135,344
MT-ND3 62,762 3,106,842
MT-ND4 85,837 1,329,685
MT-ND4L 78,658 1,687,007
MT-ND5 88,469 121,260
MT-ND6 37,381 168,240
MT-RNR1 574,902 22,065
MT-RNR2 2,069,266 44,434
采用 blast软件，将长读长测序单细胞线粒体转录本softclip5比对到短读长测序单细胞转录本barcode/UMI，选择最佳比对的barcode/UMI。在进行该最佳比对时，不考虑错配(mismatches)和空位(gaps)。
整合短读长测序单细胞转录本、长读长测序单细胞线粒体转录本barcode/ UMI比对的结果(mismatches_gaps)与转录本在参考基因组的起始位置(pos_diff =|psrs-plrs|)， 比较不同mismatches_gaps与pos_diff下的长读长测序单细胞线粒体转录本reads数目(表11所示)，最终选择(pos_diff<30 & mismatches_gaps<=3) |(pos_diff<20 &mismatches_gaps<=4) |(pos_diff<10 & mismatches_gaps<=5)作为长读长测序单细胞线粒体转录本与短读长测序单细胞转录本barcode/UMI 映射的标准。
表11. 不同比对结果与位置差异条件下的短读长测序单细胞转录本、长读长测序单细胞线粒体转录本barcode/UMI 匹配情况
过滤条件 长读长测序单细胞线粒体转录本 softclip5数目 匹配到短读长测序单细胞转录本barcode/UMI的长读长测序单细胞线粒体转录本reads数目 占比
mismatches_gaps<=3 15,184,700 2,819,913 0.1857
pos_diff<10 & mismatches_gaps<=3 15,184,700 2,240,071 0.1475
pos_diff<20 & mismatches_gaps<=3 15,184,700 2,320,856 0.1528
pos_diff<30 & mismatches_gaps<=3 15,184,700 2,347,948 0.1546
mismatches_gaps<=4 15,184,700 3,620,135 0.2384
pos_diff<10 & mismatches_gaps<=4 15,184,700 2,641,501 0.1740
pos_diff<20 & mismatches_gaps<=4 15,184,700 2,741,979 0.1806
pos_diff<30 & mismatches_gaps<=4 15,184,700 2,782,715 0.1833
mismatches_gaps<=5 15,184,700 4,713,299 0.3104
pos_diff<10 & mismatches_gaps<=5 15,184,700 3,149,672 0.2074
pos_diff<20 & mismatches_gaps<=5 15,184,700 3,277,754 0.2159
pos_diff<30 & mismatches_gaps<=5 15,184,700 3,339,542 0.2199
(pos_diff<30 & mismatches_gaps<=3) | (pos_diff<20 & mismatches_gaps<=4) | (pos_diff<10 &mismatches_gaps<=5) 15,184,700 3,277,242 0.2158
将匹配到短读长测序单细胞转录本barcode/UMI的长读长测序单细胞线粒体转录本reads序列(即，经短读长测序单细胞转录本barcode/UMI校正后的单细胞线粒体全长转录本序列)与实施例7.2获得的短读长Bulk测序线粒体转录本reads序列进行整合，并以实施例7.2获得的短读长Bulk测序线粒体转录本reads序列中的SNP/InDel为依据，获得包含高准确性 SNP/InDel的长读长测序单细胞线粒体转录本序列。
由于携带相同barcode/UMI的不同reads来自于相同细胞的相同转录本，通过对这样的reads进行SNP/InDel的相互校正，最终获得37个高质量线粒体变异。
根据以上生物信息学分析获得的携带高准确度 barcode/UMI、高准确度 SNP/InDel的线粒体全长转录本序列计算单个细胞(即，barcode标记的单细胞)水平的变异频率，构建了17700个细胞－37个变异的矩阵，并根据细胞间欧式距离，进行层次聚类分析(图10)，据此推断出细胞的进化谱系(图11，不同颜色表示不同的细胞群分类)。
实施例8. 使用商业化细胞系检验本发明的方法
在本实施例中，使用3种具有已知线粒体基因突变的细胞系验证了本发明的方法可以有效识别线粒体基因突变。与实施例1－6类似地，对3种细胞系(HepG2(供应商：ATCC，货号：HB-8065)、THP-1(供应商：ATCC，货号：TIB-202)、AGS(供应商：ATCC，货号：CRL-1739))实施了单细胞线粒体全长转录组三代测序、单细胞全转录组二代测序和Bulk线粒体转录组二代测序，并对获得的测序数据进行了分析。
实施例8.1. 单细胞线粒体全长转录组三代测序
与实施例3类似地，对3种细胞系(HepG2、THP-1、AGS)实施了单细胞线粒体全长转录组三代测序。该长读长三代测序线粒体转录组的原始下机数据reads数目为63,334,135。采用 minimap2软件(https://github.com/lh3/minimap2) 将reads比对到参考基因组GRCh38，结果表明：比对到参考基因组的reads数目为56,414,636，比对率为89.07%，目标区域reads数目为38,774,532，占比68.73%。
提取bam 文件的softclip5序列，并提取plrs。
实施例8.2. 短读长Bulk线粒体转录组二代测序
与实施例4类似地，按照NEBNext® Ultra™ II FS DNA Library Prep Kit forIllumina(NEB，目录号：E7805)产品说明书对3种细胞系(HepG2、THP-1、AGS)的经富集的单细胞线粒体全长转录组执行快速可靠的线粒体cDNA片段化、末端修复、接头连接及Bulk二代测序文库构建，并测序。
Bulk线粒体转录组二代测序原始下机reads数目为184,685,770, 使用 fastp软件，选择默认参数对该下机reads进行质控, 过滤后剩余reads数目为184,589,764, 占原始下机reads数目的99.94%。
采用STAR软件将质控后的reads序列比对到参考基因组GRCh38，随后使用GATK4Mitochondria Pipeline进行线粒体变异检测，得到11个高准确性的SNP(表12)，用于后续分析。
表12. Bulk线粒体转录组短读长二代测序后的变异信息统计
细胞系名称 疾病名称 基因 基因组改变 覆盖深度 变异频率
AGS 胃癌 MT-ATP6 g.chrM:9055G>A 5125 0.233
AGS 胃癌 MT-ND4L g.chrM:10550A>G 2261 0.238
AGS 胃癌 MT-CO3 g.chrM:9698T>C 3019 0.298
AGS 胃癌 MT-ND2 g.chrM:4561T>C 2803 0.232
HepG2 肝癌 MT-CO3 g.chrM:9950T>C 110 0.607
HepG2 肝癌 MT-ND3 g.chrM:10373G>A 1116 0.114
HepG2 肝癌 MT-CYB g.chrM:14757T>C 1015 0.602
THP-1 白血病 MT-ND1 g.chrM:3970C>T 5455 0.199
THP-1 白血病 MT-ND2 g.chrM:4732A>G 4062 0.087
THP-1 白血病 MT-ND5 g.chrM:13928G>C 1160 0.176
THP-1 白血病 MT-CO1 g.chrM:6962G>A 2565 0.208
实施例8.3. 单细胞全转录组的短读长二代测序
与实施例5类似地，对3种细胞系(HepG2、THP-1、AGS)分别制备了细胞条形码标记的单细胞全转录组全长cDNA，并进行了短读长高通量二代测序。
单细胞转录组短读长二代测序的原始下机reads数目为483,079,508，使用 fastp软件，选择默认参数对该下机reads进行质控, 经质控过滤后剩余reads数目为483,079,508, 占原始下机reads数目的100%。
将该质控后的数据使用 Cell Ranger进行处理，统计结果如下表13所示。
表13. Cell Ranger统计结果
项目 对项目的描述 数值
Estimated Number of Cells 捕获的细胞数目 16,670
Mean Reads per Cell 每个细胞的平均reads数目 14,489
Median Genes per Cell 每个细胞中位基因数目 1,164
Number of Reads reads数目 241,539,754
Valid Barcodes 有效 barcode 占比 92.30%
Reads Mapped to Genome 比对到参考基因组的比例 93.1%
Fraction Reads in Cells 细胞内 reads的比例 97.20%
Total Genes Detected 检测到的基因总数 25,401
Median UMI Counts per Cell 每个细胞的中位 UMI数目 1,833
将该质控后的序列使用 STAR软件(https://github.com/alexdobin/STAR)比对到参考基因组GRCh38，提取转录本比对到参考基因组的起始位置psrs，由此，提供了高准确性的单细胞条形码/UMI序列和其后接的基因表达片段序列。
实施例8.4. 对测序数据的生物信息学分析
将实施例8.1获得的长读长三代测序单细胞线粒体转录组数据与实施例8.3获得的短读长二代测序单细胞转录组数据进行整合，旨在进一步校正长读长测序单细胞线粒体转录组的barcode/UMI序列。
根据线粒体基因名称将实施例8.3中的短读长二代测序单细胞转录组数据中的barcode/UMI和实施例8.1中的长读长三代测序单细胞线粒体转录组数据的softclip5进行分组，获得下表14。
表14. 短读长测序单细胞转录本与长读长测序单细胞线粒体转录本reads数目统计
目标线粒体基因 短读长测序单细胞转录本barcode/UMI数目 长读长测序单细胞线粒体转录本 softclip5数目
MT-ATP6 1,151,299 1,427,969
MT-ATP8 505,344 2,133,519
MT-CO1 1,404,862 1,689,369
MT-CO2 1,127,308 4,319,140
MT-CO3 723,077 3,714,381
MT-CYB 388,633 2,068,604
MT-ND1 333,244 2,238,233
MT-ND2 196,785 1,367,499
MT-ND3 129,169 1,210,243
MT-ND4 310,451 1,682,299
MT-ND4L 263,600 1,466,838
MT-ND5 230,940 683,674
MT-ND6 124,823 511,105
MT-RNR1 2,089,213 5,099,172
MT-RNR2 5,968,084 5,504,621
采用 blast软件，将长读长测序单细胞线粒体转录本softclip5比对到短读长测序单细胞转录本barcode/UMI，选择最佳比对的barcode/UMI。在进行该最佳比对时，不考虑错配(mismatches)和空位(gaps)。
整合短读长测序单细胞转录本、长读长测序单细胞线粒体转录本barcode/ UMI比对的结果(mismatches_gaps)与转录本在参考基因组的起始位置(pos_diff =|psrs-plrs|)， 比较不同mismatches_gaps与pos_diff下的长读长测序单细胞线粒体转录本reads数目，最终选择(pos_diff<30 & mismatches_gaps<=3) |(pos_diff<20 & mismatches_gaps<=4) |(pos_diff<10 & mismatches_gaps<=5)作为长读长测序单细胞线粒体转录本与短读长测序单细胞转录本barcode/UMI 映射的标准，获得了8,219,127 reads, 占比21.19%。
由于携带相同barcode/UMI的不同reads来自于相同细胞的相同转录本，通过对这样的reads进行SNP/InDel的相互校正，最终获得7个高质量线粒体变异。
根据以上生物信息学分析获得的携带高准确度 barcode/UMI、高准确度 SNP/InDel的线粒体全长转录本序列计算单个细胞(即，barcode标记的单细胞)水平的变异频率，构建了16,602个细胞－7个变异的矩阵，并根据细胞间欧式距离，进行层次聚类分析，获得如图12所示的细胞-变异频率层次聚类图。
将实施例 8.3短读长二代测序单细胞转录组表达数据进行聚类分析(分析方法参考https://scanpy-tutorials.readthedocs.io/en/latest/pbmc3k.html)，并将HepG2细胞系特异性表达的TF基因、THP-1细胞系特异性表达的AZU1基因、AGS细胞系特异性表达的PRSS23基因作为标志基因，对这三种细胞系进行细胞群标记，结果如图13所示，图中所示“Mix”亚群表示因为细胞周期等原因造成的一部分细胞无法通过基因表达水平进行准确分类。将谱系结果映射到短读长二代测序单细胞转录本表达数据，得到基于谱系聚类的细胞亚群分类结果(图14)。
将基于短读长二代测序单细胞转录组表达的细胞亚群结果作为标准，基于谱系聚类的HepG2、THP-1、AGS细胞系的召回率(recall)分别为99.48%、99.37%、98.62%，由此，本发明的方法可以有效识别各细胞中的线粒体基因突变。
以上描述了本发明的示例性实施方案，本领域技术人员应当理解的是，这些公开内容仅是示例性的，在本发明的范围内可以进行各种其它替换、适应和修改。因此，本发明不限于文中列举的具体实施方案。
序列表
<110> 北京百图智检科技服务有限公司
<120> 确定细胞谱系的方法、试剂盒和装置
<130>
<160> 28
<170> PatentIn版本3.3
<210> 1
<211> 27
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 1
ccatgttacg acttgtctcc tctatat 27
<210> 2
<211> 34
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 2
gtataatact aagttgagat gatatcattt acgg 34
<210> 3
<211> 33
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 3
gattgtaatg ggtatggaga catatcatat aag 33
<210> 4
<211> 19
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 4
cgtggtaagg gcgatgagt 19
<210> 5
<211> 24
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 5
gttcttcgaa tgtgtggtag ggtg 24
<210> 6
<211> 30
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 6
agatttttag gggaattaat tctaggacga 30
<210> 7
<211> 23
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 7
tgaagcgaac agattttcgt tca 23
<210> 8
<211> 25
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 8
ctagaagtgt gaaaacgtag gcttg 25
<210> 9
<211> 30
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 9
agacatacag aaatagtcaa accacatcta 30
<210> 10
<211> 24
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 10
ctcataggcc agacttaggg ctag 24
<210> 11
<211> 24
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 11
gtctaggcca tatgtgttgg agat 24
<210> 12
<211> 24
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 12
ggataggagg agaatggggg atag 24
<210> 13
<211> 22
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 13
gagtagggtt aggatgagtg gg 22
<210> 14
<211> 24
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 14
tgtgggttta gtaatggggt ttgt 24
<210> 15
<211> 30
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 15
tagggagata gttggtatta ggattaggat 30
<210> 16
<211> 77
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 16
atagtgatac tgaccaccga gatctacaca tatacgcagt cgacaacttt cttgcgccat 60
gggactacac gacgctc 77
<210> 17
<211> 77
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 17
atagtgatac tgaccaccga gatctacacc gcacgacgta caaacggaat cgagcgccat 60
gggactacac gacgctc 77
<210> 18
<211> 77
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 18
atagtgatac tgaccaccga gatctacacg ggacagaagg ggacacaaga ctcgcgccat 60
gggactacac gacgctc 77
<210> 19
<211> 77
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 19
atagtgatac tgaccaccga gatctacact taacaagtat gatagaatcc gaagcgccat 60
gggactacac gacgctc 77
<210> 20
<211> 77
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 20
atagtgatac tgaccaccga gatctacaca ctcctgggta aaccctggac aaggcgccat 60
gggactacac gacgctc 77
<210> 21
<211> 77
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 21
atagtgatac tgaccaccga gatctacacc agaaaacccc gcctttgcga gaagcgccat 60
gggactacac gacgctc 77
<210> 22
<211> 72
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 22
catcagcgaa ccggcatacg agatgcacga gcagtcccac ggtaacactg ttgaactcca 60
gttcagacgt gt 72
<210> 23
<211> 72
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 23
catcagcgaa ccggcatacg agattattat gcttctttgt tccctgaatg ttgaactcca 60
gttcagacgt gt 72
<210> 24
<211> 72
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 24
catcagcgaa ccggcatacg agatagcgac tgtatgctgt gcctagtttg ttgaactcca 60
gttcagacgt gt 72
<210> 25
<211> 72
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 25
catcagcgaa ccggcatacg agatccttgc gacagggttt caacgctttg ttgaactcca 60
gttcagacgt gt 72
<210> 26
<211> 72
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 26
catcagcgaa ccggcatacg agatgttaga ttgcacgata gatgaaactg ttgaactcca 60
gttcagacgt gt 72
<210> 27
<211> 72
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 27
catcagcgaa ccggcatacg agatgttaga ttgcacgata gatgaaactg ttgaactcca 60
gttcagacgt gt 72
<210> 28
<211> 22
<212> DNA
<213> 人工序列
<220>
<223> 引物
<400> 28
ctacacgacg ctcttccgat ct 22

Claims (32)

1.一种确定多细胞真核生物中的单细胞线粒体转录组中的变异的方法，其特征在于，该方法基于高通量测序，由以下步骤组成：
a)由来自受试者的样本中的单细胞RNA制备单细胞cDNA，其中所制备的单细胞cDNA中的每一cDNA包含通用PCR引物结合位点、用于确定转录组来源细胞的特征序列、用于确定所述来源细胞中的各转录本的特征序列；
b)当步骤a)的单细胞cDNA上的用于确定转录组来源细胞的特征序列位于5’端时，对步骤a)的单细胞cDNA使用线粒体cDNA 3’端特异性引物和能够结合步骤a)的单细胞cDNA 5’端通用PCR引物结合位点的通用PCR引物通过PCR富集线粒体cDNA；或者
当步骤a)的单细胞cDNA上的用于确定转录组来源细胞的特征序列位于3’端时，对步骤a)的单细胞cDNA使用线粒体cDNA 5’端特异性引物和能够结合步骤a)的单细胞cDNA 3’端通用PCR引物结合位点的通用PCR引物通过PCR富集线粒体cDNA；
c)对步骤b)富集的线粒体cDNA进行Bulk线粒体转录组二代测序和单细胞线粒体全长转录组三代测序，其中，所述Bulk线粒体转录组二代测序是对一个所述样本的多个单细胞的获自步骤b)的线粒体cDNA富集产物实施二代测序，以获得Bulk线粒体转录组二代测序数据；所述单细胞线粒体全长转录组三代测序是对单细胞线粒体转录组中的各转录本实施三代测序，以获得单细胞线粒体全长转录组三代测序数据；
d) 用步骤c)获得的Bulk线粒体转录组二代测序数据中的用于确定转录组来源细胞的特征序列和用于确定所述来源细胞中的各转录本的特征序列的测序数据，对步骤c)获得的单细胞线粒体全长转录组三代测序数据中的用于确定转录组来源细胞的特征序列和用于确定所述来源细胞中的各转录本的特征序列的测序数据进行校正，然后将校正后的单细胞线粒体全长转录组三代测序数据与步骤c)获得的Bulk线粒体转录组二代测序数据中的变异位点数据进行整合，获得线粒体转录组中的变异信息。
2.一种确定多细胞真核生物中的单细胞线粒体转录组中的变异的方法，其特征在于，该方法基于高通量测序，由以下步骤组成：
a)由来自受试者的样本中的单细胞RNA制备单细胞cDNA，其中所制备的单细胞cDNA中的每一cDNA包含通用PCR引物结合位点、用于确定转录组来源细胞的特征序列、用于确定所述来源细胞中的各转录本的特征序列；
b)当步骤a)的单细胞cDNA上的用于确定转录组来源细胞的特征序列位于5’端时，对步骤a)的单细胞cDNA使用线粒体cDNA 3’端特异性引物和能够结合步骤a)的单细胞cDNA 5’端通用PCR引物结合位点的通用PCR引物通过PCR富集线粒体cDNA；或者
当步骤a)的单细胞cDNA上的用于确定转录组来源细胞的特征序列位于3’端时，对步骤a)的单细胞cDNA使用线粒体cDNA 5’端特异性引物和能够结合步骤a)的单细胞cDNA 3’端通用PCR引物结合位点的通用PCR引物通过PCR富集线粒体cDNA；
c)对步骤b)富集的线粒体cDNA进行Bulk线粒体转录组二代测序和单细胞线粒体全长转录组三代测序，其中，所述Bulk线粒体转录组二代测序是对一个所述样本的多个单细胞的获自步骤b)的线粒体cDNA富集产物实施二代测序，以获得Bulk线粒体转录组二代测序数据；所述单细胞线粒体全长转录组三代测序是对单细胞线粒体转录组中的各转录本实施三代测序，以获得单细胞线粒体全长转录组三代测序数据；以及对步骤a)的单细胞cDNA进行单细胞转录组二代测序，以获得单细胞转录组二代测序数据；
d) 用步骤c)获得的单细胞转录组二代测序数据中的用于确定转录组来源细胞的特征序列和用于确定所述来源细胞中的各转录本的特征序列的测序数据，对步骤c)获得的单细胞线粒体全长转录组三代测序数据中的用于确定转录组来源细胞的特征序列和用于确定所述来源细胞中的各转录本的特征序列的测序数据进行校正，然后将校正后的单细胞线粒体全长转录组三代测序数据与步骤c)获得的Bulk线粒体转录组二代测序数据中的变异位点数据进行整合，获得线粒体转录组中的变异信息。
3.根据权利要求1或权利要求2所述的方法，其中，步骤a)的单细胞cDNA是通过全转录组扩增产生的。
4.根据权利要求1或权利要求2所述的方法，其中，步骤a)的所述用于确定转录组来源细胞的特征序列是细胞条形码序列，所述用于确定来源细胞中的各转录本的特征序列是唯一分子标识符UMI序列。
5.根据权利要求1或权利要求2所述的方法，其中，步骤b)中所述线粒体cDNA 3’端特异性引物位于线粒体cDNA 3’端约200bp范围内；所述线粒体cDNA 5’端特异性引物位于线粒体cDNA 5’端约200bp范围内。
6.根据权利要求5所述的方法，其中，步骤b)中所述线粒体cDNA 3’端特异性引物位于线粒体cDNA 3’端约150bp范围内；所述线粒体cDNA 5’端特异性引物位于线粒体cDNA 5’端约150bp范围内。
7.根据权利要求6所述的方法，其中，步骤b)中所述线粒体cDNA 3’端特异性引物位于线粒体cDNA 3’端约100bp范围内；所述线粒体cDNA 5’端特异性引物位于线粒体cDNA 5’端约100bp范围内。
8.根据权利要求7所述的方法，其中，步骤b)中所述线粒体cDNA 3’端特异性引物位于线粒体cDNA 3’端约50bp范围内；所述线粒体cDNA 5’端特异性引物位于线粒体cDNA 5’端约50bp范围内。
9.根据权利要求8所述的方法，其中，步骤b)中所述线粒体cDNA 3’端特异性引物位于线粒体cDNA 3’端约30bp范围内；所述线粒体cDNA 5’端特异性引物位于线粒体cDNA 5’端约30bp范围内。
10.根据权利要求1或权利要求2所述的方法，其中，步骤d)的所述线粒体转录组中的变异信息是添加、缺失、取代和／或删除。
11.根据权利要求10所述的方法，其中，步骤d)的所述线粒体转录组中的变异信息是SNP。
12.根据权利要求1或权利要求2所述的方法，其中，步骤b)的线粒体cDNA 3’端特异性引物或线粒体cDNA 5’端特异性引物是一条或多条引物，以及在1个或多个PCR管中实施步骤b)，由此，通过PCR富集线粒体cDNA。
13.根据权利要求12所述的方法，其中，步骤b)的线粒体cDNA 3’端特异性引物是SEQID NO: 1－15中的一条或多条引物，以及在1个PCR管中实施步骤b)，由此，通过PCR富集线粒体cDNA。
14.根据权利要求13所述的方法，其中，步骤b)的线粒体cDNA 3’端特异性引物是SEQID NO: 1－15中的一条或多条引物，以及在1个PCR管中根据线粒体中各转录本的表达水平比例来混合相应量的所述特异性引物实施步骤b)，由此，通过PCR富集线粒体cDNA。
15.根据权利要求1或权利要求2所述的方法，其中，步骤c)的三代测序是将额外的多个样本分别经步骤a)和步骤b)后富集的线粒体cDNA在一管中进行混合，然后进行三代测序实施的。
16.根据权利要求15所述的方法，其中，所述额外的多个样本是额外的1个、2个、3个、4个、5个、6个、7个、8个、9个样本或10个样本。
17.根据权利要求16所述的方法，其中，在步骤c)的所述三代测序时实施的建库还包括添加至少一个样本索引的核苷酸序列。
18.根据权利要求17所述的方法，其中，所述添加至少一个样本索引的核苷酸序列选自至少一对SEQ ID NO: 16和17、SEQ ID NO: 18和19、SEQ ID NO: 20和21、SEQ ID NO: 22和23、SEQ ID NO: 24和25、SEQ ID NO: 26和27样本索引引物。
19.根据权利要求1所述的方法，其中，步骤d)包括：
(1) 从步骤c)生成的单细胞线粒体全长转录组三代测序数据，获得单细胞线粒体全长转录组的三代测序序列；
(2) 从步骤c)生成的Bulk线粒体转录组二代测序数据，获得Bulk线粒体转录组的高准确性二代测序序列，其中包括细胞条形码/唯一分子标识符UMI序列和单细胞变异位点；
(3)用所述(2)中获得的Bulk线粒体转录组的高准确性二代测序短读长序列对所述(1)中获得的单细胞线粒体全长转录组的三代测序序列进行校正，获得具有高准确性单细胞变异位点的单细胞线粒体全长转录组序列。
20.根据权利要求2所述的方法，其中，步骤d)包括：
(1) 从步骤c)生成的三代单细胞线粒体全长转录组测序数据，获得单细胞线粒体全长转录组的变异位点；
(2) 从步骤c)生成的二代Bulk线粒体转录组测序数据，获得线粒体转录组的短读长中的高准确性单细胞变异位点；
(3) 从步骤c)生成的二代短读长单细胞转录组测序数据，获得高准确性细胞条形码/唯一分子标识符UMI序列；
(4) 将所述(1)和(3)中分别获得的结果进行细胞条形码/唯一分子标识符UMI映射，获得单细胞线粒体全长转录组序列，其具有高准确性细胞条形码/唯一分子标识符UMI序列和较高准确性的变异位点；
(5) 对(4)得到的细胞条形码/唯一分子标识符UMI校正后的单细胞线粒体全长转录组序列与(2)得到的高准确性单细胞变异位点进行整合，获得具有高准确性单细胞变异位点的单细胞线粒体全长转录组序列。
21.根据权利要求19或权利要求20所述的方法，其中，所获得的单细胞变异位点在二代Bulk线粒体转录组测序数据中的至少5个读长中被检测到。
22.一种确定多细胞真核生物中细胞谱系的方法，其特征在于，包括：权利要求1-21中任一项所述的步骤a)- d)，并进一步包括：
步骤e)基于细胞中线粒体变异的存在对细胞进行聚类，并推断细胞的谱系。
23.根据权利要求22所述的方法，其中，通过计算单细胞变异位点的变异频率，构建细胞-变异矩阵，由此推断细胞的进化谱系。
24.根据权利要求23所述的方法，其中，所述变异是SNP/InDel变异。
25.一种识别线粒体测序数据中线粒体变异的装置，其特征在于，由以下模块组成：
模块1. 接收单细胞线粒体全长三代测序数据，并映射到靶线粒体DNA样本；
模块2. 接收Bulk线粒体转录组二代测序数据，获得高准确性的基因水平的变异和细胞条形码/UMI序列；
模块3. 用模块2中获得的高准确性细胞条形码/UMI序列校正模块1的细胞条形码/UMI序列，再与模块2得到的高准确性的变异位点进行整合，获得高准确性的线粒体变异位点。
26.根据权利要求25所述的装置，其中，所述变异是SNP/InDel变异。
27.一种确定多细胞真核生物中细胞谱系的装置，其特征在于，由以下模块组成：
权利要求25或权利要求26所述的装置中的模块；以及
模块4. 基于模块3得到的高准确性的线粒体变异位点的存在对细胞进行聚类，并推断细胞的谱系。
28.根据权利要求27所述的装置，其中，通过计算单细胞变异位点的变异频率，构建细胞-变异矩阵，由此推断细胞的进化谱系。
29.一种识别线粒体测序数据中线粒体变异的装置，其特征在于，由以下模块组成：
模块1. 接收单细胞线粒体全长三代测序数据，并映射到靶线粒体DNA样本；
模块2. 接收单细胞转录组二代测序数据并分析，获得高准确性细胞条形码/UMI序列数据；
模块3. 接收Bulk线粒体转录组二代测序数据并分析，获得高准确性的基因水平的变异，从而获得高准确性单细胞变异位点；
模块4. 将所述模块1和模块2中分别获得的长、短读长单细胞数据细胞条形码/UMI映射，获得高准确性细胞条形码/UMI序列及较高准确性单细胞变异位点；
模块5. 对模块4得到的细胞条形码/UMI校正后的单细胞全长转录组序列与模块3得到的高准确性的变异位点进行整合，获得高准确性的线粒体变异位点。
30.根据权利要求29所述的装置，其中，所述变异是SNP/InDel变异。
31.一种确定多细胞真核生物中细胞谱系的装置，其特征在于，由以下模块组成：
权利要求29或权利要求30所述的装置中的模块；以及
模块6. 基于模块5得到的高准确性的线粒体变异位点的存在对细胞进行聚类，并推断细胞的谱系。
32.根据权利要求31所述的装置，其中，通过计算单细胞变异位点的变异频率，构建细胞-变异矩阵，由此推断细胞的进化谱系。
CN202210756234.3A 2022-06-30 2022-06-30 确定细胞谱系的方法、试剂盒和装置 Active CN114875118B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210756234.3A CN114875118B (zh) 2022-06-30 2022-06-30 确定细胞谱系的方法、试剂盒和装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210756234.3A CN114875118B (zh) 2022-06-30 2022-06-30 确定细胞谱系的方法、试剂盒和装置

Publications (2)

Publication Number Publication Date
CN114875118A CN114875118A (zh) 2022-08-09
CN114875118B true CN114875118B (zh) 2022-10-11

Family

ID=82683087

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210756234.3A Active CN114875118B (zh) 2022-06-30 2022-06-30 确定细胞谱系的方法、试剂盒和装置

Country Status (1)

Country Link
CN (1) CN114875118B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2622371A (en) * 2022-09-13 2024-03-20 Agecurve Ltd Cell tree rings: Method and cell lineage tree based aging timer for calculating biological age of biological sample
CN118486378B (zh) * 2024-07-16 2024-11-29 墨卓生物科技(浙江)有限公司 用于单细胞免疫组库测序数据的降噪方法、设备和存储介质

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108460246A (zh) * 2018-03-08 2018-08-28 北京希望组生物科技有限公司 一种基于三代测序平台的hla基因分型方法
CN113035269A (zh) * 2021-04-16 2021-06-25 北京计算科学研究中心 基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法
CN113343736A (zh) * 2021-06-21 2021-09-03 天津大学合肥创新发展研究院 一种dna测序用条形码识别算法的硬件加速装置
WO2021207576A1 (en) * 2020-04-09 2021-10-14 Takeda Vaccines, Inc. Qualitative and quantitative determination of single virus haplotypes in complex samples
CN113724788A (zh) * 2021-07-29 2021-11-30 哈尔滨医科大学 一种鉴定肿瘤细胞的染色体外环状dna组成基因的方法
CN114292912A (zh) * 2021-12-24 2022-04-08 广州燃石医学检验所有限公司 一种变体核酸的检测方法
CN114540473A (zh) * 2021-08-27 2022-05-27 四川大学华西第二医院 一种新型核酸测序系统

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108460246A (zh) * 2018-03-08 2018-08-28 北京希望组生物科技有限公司 一种基于三代测序平台的hla基因分型方法
WO2021207576A1 (en) * 2020-04-09 2021-10-14 Takeda Vaccines, Inc. Qualitative and quantitative determination of single virus haplotypes in complex samples
CN113035269A (zh) * 2021-04-16 2021-06-25 北京计算科学研究中心 基于高通量测序技术的基因组代谢模型构建、优化及可视化的方法
CN113343736A (zh) * 2021-06-21 2021-09-03 天津大学合肥创新发展研究院 一种dna测序用条形码识别算法的硬件加速装置
CN113724788A (zh) * 2021-07-29 2021-11-30 哈尔滨医科大学 一种鉴定肿瘤细胞的染色体外环状dna组成基因的方法
CN114540473A (zh) * 2021-08-27 2022-05-27 四川大学华西第二医院 一种新型核酸测序系统
CN114540472A (zh) * 2021-08-27 2022-05-27 四川大学华西第二医院 一种新型三代测序方法
CN114292912A (zh) * 2021-12-24 2022-04-08 广州燃石医学检验所有限公司 一种变体核酸的检测方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Jabba: hybrid error correction for long sequencing reads;Giles Miclotte，等;《Algorithms for Molecular Biology》;20160503;第11卷;第1-12页 *
Lineage Tracing in Humans Enabled byMitochondrial Mutations and Single-Cell Genomics;Leif S. Ludwig等;《Cell》;20190307;第176卷(第6期);第1325-1339页，第1326页右列第3段，第1331页右列第2段，图1，图5 *
基因组二代测序数据与三代测序数据的混合校正和组装;马东娜，等;《基因组学与应用生物学》;20180425;第37卷(第4期);第1147-1155页 *

Also Published As

Publication number Publication date
CN114875118A (zh) 2022-08-09

Similar Documents

Publication Publication Date Title
US11814678B2 (en) Universal short adapters for indexing of polynucleotide samples
US11788139B2 (en) Optimal index sequences for multiplex massively parallel sequencing
US11161087B2 (en) Methods and compositions for tagging and analyzing samples
CA2983935C (en) Error suppression in sequenced dna fragments using redundant reads with unique molecular indices (umis)
CN111808854B (zh) 带有分子条码的平衡接头及快速构建转录组文库的方法
CN113463202B (zh) 一种新的rna高通量测序的方法、引物组和试剂盒及其应用
CN114875118B (zh) 确定细胞谱系的方法、试剂盒和装置
CN112703253A (zh) 液滴单细胞表观基因组谱分析用于患者分层的用途
CN112996924A (zh) 液滴单细胞表观基因组谱分析用于患者分层的用途
US20210324454A1 (en) Systems and methods for correcting sample preparation artifacts in droplet-based sequencing
El Amrani Computational methods for the identification and characterization of tissues and cells
Sos The Single Cell Transposome Hypersensitive Sites Sequencing (scTHS-seq) assay for Chromatin Accessibility and Assessment of Epigenetic States in the Human Adult Brain
CN109097481A (zh) 基于高通量测序技术鉴定鱼卵仔稚鱼的方法
Coden Computational detection of doublets in single-cell RNA sequencing

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20240320

Address after: 101, 2nd Floor, Building 10, Yard 9, Yongteng North Road, Haidian District, Beijing, 100089

Patentee after: Baitu Shengke (Beijing) Intelligent Technology Co.,Ltd.

Country or region after: Zhong Guo

Address before: Room 614, room 909, floor 9, block B, No. 18, Zhongguancun Street, Haidian District, Beijing 100086

Patentee before: Beijing Baitu Zhijian Technology Service Co.,Ltd.

Country or region before: Zhong Guo