CN111292807B - 一种单细胞转录组数据中分析双细胞的方法 - Google Patents
一种单细胞转录组数据中分析双细胞的方法 Download PDFInfo
- Publication number
- CN111292807B CN111292807B CN201811487856.0A CN201811487856A CN111292807B CN 111292807 B CN111292807 B CN 111292807B CN 201811487856 A CN201811487856 A CN 201811487856A CN 111292807 B CN111292807 B CN 111292807B
- Authority
- CN
- China
- Prior art keywords
- cells
- cell
- heterogeneous
- synthetic
- double
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 17
- 230000014509 gene expression Effects 0.000 claims abstract description 41
- 238000005094 computer simulation Methods 0.000 claims abstract description 7
- 230000002194 synthesizing effect Effects 0.000 claims abstract description 7
- 238000004422 calculation algorithm Methods 0.000 claims description 24
- 239000000592 Artificial Cell Substances 0.000 claims description 16
- 238000001228 spectrum Methods 0.000 claims description 10
- 238000007476 Maximum Likelihood Methods 0.000 claims description 4
- 238000005070 sampling Methods 0.000 claims description 4
- 230000002093 peripheral effect Effects 0.000 claims description 3
- 238000012163 sequencing technique Methods 0.000 abstract description 8
- 210000004027 cell Anatomy 0.000 description 185
- 238000005516 engineering process Methods 0.000 description 6
- 238000003064 k means clustering Methods 0.000 description 6
- 239000011159 matrix material Substances 0.000 description 6
- 238000004458 analytical method Methods 0.000 description 5
- 238000004364 calculation method Methods 0.000 description 3
- 210000004460 N cell Anatomy 0.000 description 2
- 239000013598 vector Substances 0.000 description 2
- 238000009825 accumulation Methods 0.000 description 1
- 101150010487 are gene Proteins 0.000 description 1
- 238000003491 array Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000024245 cell differentiation Effects 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 239000013604 expression vector Substances 0.000 description 1
- 239000003550 marker Substances 0.000 description 1
- 108020004999 messenger RNA Proteins 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012174 single-cell RNA sequencing Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Images
Landscapes
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本发明公开了一种单细胞转录组数据中分析双细胞的方法，属于单细胞转录组测序技术领域。本发明通过计算机模拟，基于实验数据，合成细胞表达谱，通过鉴别合成异类双细胞以鉴别实验数据中的双细胞。本发明既不依赖于对细胞类型的现有知识，也不依赖于对细胞总表达量的假设，而是在计算机模拟的基础上就可实现全自动化双细胞的鉴别。
Description
技术领域
本发明属于单细胞转录组测序技术领域，具体地说本发明涉及一种单细胞转录组数据中分析双细胞的方法。
背景技术
单细胞转录组测序技术(scRNA-Seq)在近年来发展迅速，相比于常规转录组测序而言，该技术可以将基因表达的分辨率从群体水平的平均值精确到单个细胞水平（具体可参见文献1、2）。尤其是从2015年以来刚刚兴起的基于微流控技术的各类海量单细胞测序技术(Massively Parallel Single Cell Sequencing)（具体可参见文献3-5），可以将大量细胞分离成单细胞并给每个单细胞加上分子标签后联合扩增建库测序，避免了人工分离单细胞的繁琐实验流程。
具体而言，海量单细胞测序技术通过液滴或微孔阵列的微流控系统分离单细胞。根据泊松分布的原理，当细胞数远少于液滴数或微孔数时，大部分液滴或微孔为空，一部分液滴或微孔含有一个细胞，只有极少比例的的液滴或微孔含有2个或2个以上的细胞。由于总液滴数或微孔数一定，当细胞数增多时，双细胞（doublet）比例也随之提高。当分析10000个细胞时，双细胞比例通常情况下约为5%。该技术可以在双细胞率较低的基础上实现同时分析成千上万个单细胞。
然而，即使如上所述双细胞率相对较低，但在通量较高的实验中，双细胞个数也会达到几百个。这些双细胞基因表达谱在后续的数据分析中也可能会影响细胞亚型聚类、细胞分化状态分析、细胞亚型标志物分析等的结果，造成对生物学意义理解上的偏差或误导。因此，在数据中鉴别出双细胞可以进一步提高单细胞转录组实验的精度，降低实验噪音，减少后续分析错误的可能性。
在现有的实践之中，通常会采用两种方法鉴定并去除单细胞表达谱中的双细胞：i). 去除总表达分子数偏高的细胞；ii). 通过已知细胞类型的标志物，去除表达两种或两种以上细胞亚群标志物的细胞（如文献6、7）。但是第一种方法假设任意两个细胞表达的分子总数大于任意单个细胞的表达分子数。而事实上，单细胞表达谱的动态范围很大，且不同细胞类型的mRNA总表达水平可相差几倍。因此，采用这种方法去除双细胞所基于的假设并不成立。而第二种方法依赖于已有的生物学知识积累，通过手动注释的方法，去除已知标志物的细胞类型混合而成的双细胞，但有可能错误的去除了细胞类型转化的中间状态。因此，第二种方法需要极高的知识背景，步骤繁琐，而且假阳性和假阴性都很高。
上述文献1-7具体为：
1. Kolodziejczyk, A.A., et al., The technology and biology of single-cell RNA sequencing. Mol Cell, 2015. 58(4): p. 610-20.
2. Shapiro, E., T. Biezuner, and S. Linnarsson, Single-cellsequencing-based technologies will revolutionize whole-organism science. NatRev Genet, 2013. 14(9): p. 618-30.
3. Fan, H.C., G.K. Fu, and S.P. Fodor, Expression profiling.Combinatorial labeling of single cells for gene expression cytometry.Science, 2015. 347(6222): p. 1258367.
4. Klein, A.M., et al., Droplet barcoding for single-celltranscriptomics applied to embryonic stem cells. Cell, 2015. 161(5): p. 1187-1201.
5. Macosko, E.Z., et al., Highly Parallel Genome-wide ExpressionProfiling of Individual Cells Using Nanoliter Droplets. Cell, 2015. 161(5):p. 1202-1214.
6. Griffiths JA, Scialdone A, Marioni JC: Using single-cell genomicsto understand developmental processes and cell fate decisions. Mol Syst Biol2018, 14:e8046.
7. Stoeckius, M. et al. Cell "hashing" with barcoded antibodiesenables multiplexing and doublet detection for single cell genomics. bioRxiv(2017)。
发明内容
本发明目的是：为了克服现有技术鉴定并去除单细胞表达谱中的双细胞存在的问题，提供一种自动化鉴定双细胞并拆分双细胞中细胞表达谱的方法。
具体地说，本发明是采用以下技术方案实现的，包括以下步骤：
合成细胞的步骤：通过计算机模拟，基于实验数据合成细胞表达谱，得到合成细胞，所述实验数据为检测到的细胞数据；
鉴别合成异类双细胞的步骤：将合成细胞与实验数据合并，采用k近邻算法计算合并的数据集中每个合成细胞的近邻及近邻中合成双细胞的比例，并基于该比例采用聚类算法将合成细胞分成两类，其中高于分类阈值的那类为合成异类双细胞；
鉴别实验数据中的双细胞的步骤：将合成异类双细胞和实验数据合并，采用k近邻算法计算合并的数据集中实验数据中每个检测到的细胞周围近邻及近邻中合成异类双细胞的比例，并基于该比例采用聚类算法将检测到的细胞分为两类，其中高于分类阈值的那类为实验数据中的异类双细胞。
进一步而言，上述方案还包括拆分鉴别出的异类双细胞为两个单细胞表达谱的步骤：对每个鉴别出的实验数据中的异类双细胞，选取该异类双细胞所有近邻中的合成异类双细胞及各合成异类双细胞的来源单细胞，使用聚类算法把来源单细胞分为两类，得到两个单细胞亚型表达谱的分布，通过在两个单细胞亚群表达谱分布中多轮次随机抽样，采用最大似然估计寻找两个叠加之后和该异类双细胞表达谱距离最小的单细胞，作为该异类双细胞的两个替换对应单细胞。
进一步而言，所述鉴别合成异类双细胞的步骤或者鉴别实验数据中的双细胞的步骤中，所述聚类算法为k-means、DBSCAN、层次聚类或其他聚类算法。
本发明的有益效果如下：本发明既不依赖于对细胞类型的现有知识，也不依赖于对细胞总表达量的假设，而是在计算机模拟的基础上就可实现全自动化双细胞的鉴别。在此基础上，在双细胞鉴别出之后，通过拆分双细胞表达谱，可以将检测到的双细胞表达谱拆分成为2个单细胞表达谱。因此，采用本发明，可以无需去除鉴别出的双细胞，在最大程度上利用数据，提高实验中检测到细胞数的效率。
附图说明
图1：本发明流程示意图。
图2：本发明数据示例的同类双细胞与个异类双细胞的分类结果图。
图3：本发明数据示例的待检测细胞的近邻中合成双细胞的比例计算结果图。
具体实施方式
下面结合实施例并参照附图对本发明作进一步详细描述。
实施例1：
本实施例公开了一种自动化鉴定双细胞并拆分双细胞中细胞表达谱的方法，其原理是基于计算机模拟进行双细胞的鉴定分析，主要包括以下四个过程：
合成细胞：通过计算机模拟，基于实验数据，合成细胞表达谱。实验数据为检测到的细胞数据，由于实验数据中双细胞的比例很低，绝大多数为单细胞，所以绝大多数合成细胞将由检测到的单细胞数据中的两个随机来源单细胞叠加而成，即形成合成双细胞。
双细胞分为两种，同类双细胞由2个属于同种细胞类型的细胞混合而成，而异类双细胞由2个不同种类的细胞类型混合而成。故形成的合成双细胞也分为合成同类双细胞和合成异类双细胞两类。
而同类双细胞的表达谱和其来源细胞类型的单细胞表达谱高度相似，并不会对下游分析解读造成很大的影响；而异类双细胞具有两种细胞类型的表达特征，因此是需要专注的对象。
鉴别合成异类双细胞：合并合成细胞和实验数据，采用k近邻算法计算每个合成细胞的近邻及近邻中合成双细胞的比例，并以此给每个双细胞打分。而同类双细胞的近邻中将既有合成的同类双细胞，也有检测到的单细胞。而异类双细胞的近邻中将有合成异类双细胞和实验数据中待鉴定的异类双细胞。由于待鉴定的双细胞比例很小，异类双细胞的近邻将绝大多数为合成异类双细胞。可以基于这一原理，基于近邻中合成双细胞的比例采用聚类算法将合成细胞分成两类，其中高于分类阈值的那类为合成异类双细胞。
鉴别实验数据中的双细胞：通过合并合成异类双细胞和实验数据，采用k近邻算法计算实验数据中每个检测到的细胞周围近邻及近邻中合成异类双细胞的比例，实验数据中每个检测到的细胞即为每个待鉴定的细胞，并以此给每个待鉴定的细胞打分。和上面同样的原理，真实的单细胞的近邻大部分为其他单细胞，而检测到的异类双细胞的近邻中应含有较多的合成异类双细胞。再次采用聚类算法将检测到的细胞分为两类，其中高于分类阈值的那类为实验数据中的异类双细胞。
拆分鉴别出的异类双细胞为两个单细胞表达谱：对每个鉴别出的异类双细胞，选取该鉴别出的异类双细胞所有近邻中的合成异类双细胞，并找到每个近邻中的合成异类双细胞的两个来源单细胞。采用聚类算法将所有近邻的来源单细胞分为两类，并得到两个单细胞亚型表达谱的分布。这里的假设是表达谱相似的合成双细胞的来源细胞是来自于相同的两个细胞亚型。通过在两个细胞亚群表达谱分布中多轮次随机抽样，采用最大似然估计寻找两个叠加之后和该鉴别出的异类双细胞表达谱距离最小的单细胞，作为该鉴别出的异类双细胞的两个替换对应单细胞。这两个单细胞可在原检测到的实验数据中替换对应的异类双细胞。
上述的聚类算法可以为k-means、DBSCAN或层次聚类等聚类算法，其中对于k-means算法，采用k=2的k-means聚类。
以下给出一个具体的计算实例：
1）设检测到的细胞数据可以表述成一个N*M的表达矩阵A，其中N为总细胞数，M为总基因数。矩阵A的元素为基因表达量，为整数值。从A中随机按行抽取2个1*M的向量（有放回），将2个抽取的向量按列求和，得到一个合成细胞表达向量。重复S次后得到一个S*M的合成细胞矩阵B，即得到S个合成细胞。
2）将矩阵A和B按列合并，并计算合并后的矩阵中所有S个合成细胞的k近邻。这里采用pearson相关系数或其他类似指标作为衡量两点之间距离远近的指标。k为用户可设参数。
3）计算S个合成细胞的近邻中合成双细胞所占比例{f1, f2, …, fs}。
4）通过k=2的k-means聚类的方法计算最佳的分类阈值c将{f1, f2, …, fs}分为两部分，从而将合成细胞表达矩阵B分成两个矩阵B1、B2，其中B1包含近邻中合成双细胞比例低于阈值c的合成细胞，B2包含近邻中合成双细胞比例高于阈值c的双细胞（即合成异类双细胞）。
5）将矩阵A和B2按列合并，并计算合并后的矩阵中所有N个细胞的k近邻，计算采用的参数与上述步骤2）中参数相同。
6）计算N个细胞的近邻中合成异类双细胞所占比例{g1, g2, …, gn}。
7）通过k=2的k-means聚类计算最佳的分类阈值c’将{g1, g2, …, gn}分为两部分。其中近邻中合成异类双细胞比例高于阈值c’的细胞注释为异类双细胞。
8）对于每一个表达矩阵A中的异类双细胞x，提取其所有近邻合成异类双细胞{n1,n2, …, nk}和各近邻合成异类双细胞的来源单细胞对{(p1, q1), (p2, q2), …,(pk,qk)}。使用k=2的k-means聚类算法将单细胞{p1, p2, …, pk, q1, q2, …, qk}分成2类D1、D2。从D1和D2两个分布中随机抽样合成双细胞y，并计算其与x的距离。采用最大似然估计，在1000次抽样中选择yi使得 yi-x = min{y1-x, y2-x, …, y1000-x}。此时可以认为yi的两个来源细胞pyi, qyi的表达谱可近似为双细胞x拆分后的单细胞表达谱。
数据示例如下：
实验数据采用10x genomics平台公开数据：50%:50% Jurkat:293T CellMixture，其中总细胞数为3388，总基因数为32738。各步骤的相关情况如下：
1、通过1000次抽样，每次随机抽取2个细胞并求和，得到1000个合成双细胞。
2、合并实验检测到的细胞和合成双细胞后得到总共4388个细胞。
3、通过k=10的k近邻算法计算1000个合成双细胞的近邻。
4、计算1000个合成双细胞的近邻中，其他合成双细胞的比例。通过k-means聚类将1000个合成双细胞分成2类，得到220个同类双细胞，780个异类双细胞。如图2所示，图中最左侧两个立柱的上部（浅灰色部分）为3388个待检测细胞的近邻中的合成双细胞比例，其余部分为1000个合成双细胞的近邻中的合成双细胞比例。近邻中合成双细胞比例的阈值为0.4，左侧为检测到的单细胞和合成同类双细胞，右侧为合成异类双细胞。
5、合并实验检测到的3388个细胞，和步骤4中鉴定出的780个异类双细胞，得到共4168个细胞。通过k=10的k近邻算法计算3388个待检测细胞的近邻。计算3388个待检测细胞的近邻中合成双细胞的比例。如图3所示，为共4168个细胞的近邻中的合成双细胞比例。通过k-means聚类将待检测细胞分为2类，得到3349个单细胞，39个双细胞。近邻中合成双细胞比例的阈值为0.5，左侧为鉴定出的真实单细胞，右侧为鉴定出的双细胞和合成异类双细胞。
经过单细胞拆分步骤，将39个已鉴定出的的双细胞拆分为78个单细胞。将原始实验数据的3388个待检测细胞替换为3349个真实单细胞和78个拆分单细胞的集合。
虽然本发明已以较佳实施例公开如上，但实施例并不是用来限定本发明的。在不脱离本发明之精神和范围内，所做的任何等效变化或润饰，同样属于本发明之保护范围。因此本发明的保护范围应当以本申请的权利要求所界定的内容为标准。
Claims (3)
1.一种单细胞转录组数据中分析双细胞的方法，其特征在于，包括以下步骤：
合成细胞的步骤：通过计算机模拟，基于实验数据合成细胞表达谱，得到合成细胞，所述实验数据为检测到的细胞数据；
鉴别合成异类双细胞的步骤：将合成细胞与实验数据合并，采用k近邻算法计算合并的数据集中每个合成细胞的近邻及近邻中合成双细胞的比例，并基于该比例采用聚类算法将合成细胞分成两类，其中高于分类阈值的那类为合成异类双细胞；
鉴别实验数据中的双细胞的步骤：将合成异类双细胞和实验数据合并，采用k近邻算法计算合并的数据集中实验数据中每个检测到的细胞周围近邻及近邻中合成异类双细胞的比例，并基于该比例采用聚类算法将检测到的细胞分为两类，其中高于分类阈值的那类为实验数据中的异类双细胞；
拆分鉴别出的异类双细胞为两个单细胞表达谱的步骤：对每个鉴别出的实验数据中的异类双细胞，选取该异类双细胞所有近邻中的合成异类双细胞及各合成异类双细胞的来源单细胞，使用聚类算法把来源单细胞分为两类，得到两个单细胞亚型表达谱的分布，通过在两个单细胞亚群表达谱分布中多轮次随机抽样，采用最大似然估计寻找两个叠加之后和该异类双细胞表达谱距离最小的单细胞，作为该异类双细胞的两个替换对应单细胞。
2.根据权利要求1所述的单细胞转录组数据中分析双细胞的方法，其特征在于：所述鉴别合成异类双细胞的步骤中，所述聚类算法为k-means、DBSCAN或层次聚类。
3.根据权利要求1所述的单细胞转录组数据中分析双细胞的方法，其特征在于：所述鉴别实验数据中的双细胞的步骤中，所述聚类算法为k-means、DBSCAN或层次聚类。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811487856.0A CN111292807B (zh) | 2018-12-06 | 2018-12-06 | 一种单细胞转录组数据中分析双细胞的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811487856.0A CN111292807B (zh) | 2018-12-06 | 2018-12-06 | 一种单细胞转录组数据中分析双细胞的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111292807A CN111292807A (zh) | 2020-06-16 |
| CN111292807B true CN111292807B (zh) | 2021-10-08 |
Family
ID=71022794
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811487856.0A Active CN111292807B (zh) | 2018-12-06 | 2018-12-06 | 一种单细胞转录组数据中分析双细胞的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111292807B (zh) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113257365B (zh) * | 2021-05-26 | 2022-07-12 | 南开大学 | 面向非标准化单细胞转录组测序数据的聚类方法及系统 |
| CN113611368B (zh) * | 2021-07-26 | 2022-04-01 | 哈尔滨工业大学（深圳） | 基于2d嵌入的半监督单细胞聚类方法、装置、计算机设备 |
| CN116486916A (zh) * | 2022-11-03 | 2023-07-25 | 杭州联川生物技术股份有限公司 | 一种单细胞转录组濒死细胞和多细胞过滤方法、介质和设备 |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106701995A (zh) * | 2017-02-20 | 2017-05-24 | 元码基因科技（北京）有限公司 | 通过单细胞转录组测序进行细胞质量控制的方法 |
-
2018
- 2018-12-06 CN CN201811487856.0A patent/CN111292807B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106701995A (zh) * | 2017-02-20 | 2017-05-24 | 元码基因科技（北京）有限公司 | 通过单细胞转录组测序进行细胞质量控制的方法 |
Non-Patent Citations (2)
| Title |
|---|
| DoubletFinder: Doublet detection in single-cell RNA sequencing data using artificial nearest neighbors;Christopher S. McGinnis等;《https://www.biorxiv.org/content/10.1101/352484v3.full.pdf》;20180719;第2-10、13页以及图1-2 * |
| Palladium-based mass tag cell barcoding with a doublet-filtering scheme and single-cell deconvolution algorithm;Eli R Zunder等;《nature protocols》;20150122;全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111292807A (zh) | 2020-06-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Mereu et al. | Benchmarking single-cell RNA-sequencing protocols for cell atlas projects | |
| Lähnemann et al. | Eleven grand challenges in single-cell data science | |
| US10347365B2 (en) | Systems and methods for visualizing a pattern in a dataset | |
| Seyfferth et al. | Advances and opportunities in single-cell transcriptomics for plant research | |
| Xu et al. | Identification of cell types from single-cell transcriptomes using a novel clustering method | |
| Zhao et al. | Bartender: a fast and accurate clustering algorithm to count barcode reads | |
| Ames et al. | Scalable metagenomic taxonomy classification using a reference genome database | |
| Liu et al. | Unsupervised embedding of single-cell Hi-C data | |
| Abdelaal et al. | Predicting cell populations in single cell mass cytometry data | |
| CN111292807B (zh) | 一种单细胞转录组数据中分析双细胞的方法 | |
| Patruno et al. | A review of computational strategies for denoising and imputation of single-cell transcriptomic data | |
| US20190332963A1 (en) | Systems and methods for visualizing a pattern in a dataset | |
| US20210381056A1 (en) | Systems and methods for joint interactive visualization of gene expression and dna chromatin accessibility | |
| CN112466404B (zh) | 一种宏基因组重叠群无监督聚类方法及系统 | |
| Gambardella et al. | The impact of microRNAs on transcriptional heterogeneity and gene co-expression across single embryonic stem cells | |
| Todorov et al. | Computational approaches for high‐throughput single‐cell data analysis | |
| Jeong et al. | PRIME: a probabilistic imputation method to reduce dropout effects in single-cell RNA sequencing | |
| CN113470743A (zh) | 一种基于bd单细胞转录组和蛋白组测序数据的差异基因分析方法 | |
| Massey et al. | High-throughput analysis of single human cells reveals the complex nature of DNA replication timing control | |
| Chen et al. | STEEL enables high-resolution delineation of spatiotemporal transcriptomic data | |
| Trapnell et al. | Monocle: Cell counting, differential expression, and trajectory analysis for single-cell RNA-Seq experiments | |
| Der et al. | Single-cell RNA sequencing for the study of lupus nephritis | |
| Chen et al. | DoRC: Discovery of rare cells from ultra-large scRNA-seq data | |
| CN107609348B (zh) | 高通量转录组数据样本分类数目估计方法 | |
| Kawaguchi et al. | Learning single-cell chromatin accessibility profiles using meta-analytic marker genes |
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 |