10X单细胞 & 10XATAC 联合分析表征细胞调控网络（MIRA）

严格比较同一单细胞中的基因表达和染色质可及性，可以阐明这些机制的耦合或解耦如何调节fate commitment的逻辑 （ 看来这里分析的是单细胞多组学的RNA & ATAC ）。在这里，开发了 MIRA： 用于综合调节分析的概率多模态模型，这是一种综合方法，可以系统地对比转录和可及性，以推断沿着发育轨迹驱动细胞的调节网络 。 MIRA 利用细胞状态的joint topic modeling和单个基因位点的调节潜力建模 。 MIRA 从而在一个有效且可解释的潜在空间中代表细胞状态，推断高保真谱系树，确定分支点命运决定的关键调节因子，并揭示局部可及性对不同位点转录的可变影响 。 MIRA 应用于来自两个不同多模式平台的表皮维持分化和胚胎大脑发育，揭示早期发育基因受到局部染色质景观的严格调控，而终末命运基因的滴定无需大量染色质重塑。

对同一单细胞中的表达和染色质可及性进行分析为了解驱动细胞发育连续体的转录和表观遗传机制的相互作用提供了前所未有的机会 。 虽然许多计算方法分别分析表达和可及性， 但最近的几种算法采用了联合分析，其中基于两种数据模式将细胞投影到共享的潜在空间上，从而更好地捕获数据的生物结构 。 然而，该领域缺乏超越可视化和聚类的工具来严格对比每个细胞中的转录和可及性，以阐明驱动发育命运决定的复杂调控网络。

跨发育轨迹的全局转录和可及性状态的综合分析将能够发现控制谱系分支点命运决定的关键调节因子 。 在基因水平上， 检查基因位点附近的转录动态与染色质可及性可能会揭示这些机制如何相互作用以调节不同的基因模块 。 某些基因可能受顺式调控元件的调控，这些调控元件在它们变得可及时同时被激活，而其他基因可能受可及性和激活分离的元件的调控 。 确定哪些基因受这些不同机制调控的逻辑可以深入了解需要严格时空调控与信号响应的通路模式。

在这里，开发了 MIRA： 用于综合调节分析的概率多模态模型，这是一种综合方法，系统地对比转录和可及性，以确定沿着发育连续体驱动细胞的调节网络 。 MIRA 利用细胞状态的joint topic modeling和单个基因位点的调节潜力 (RP) 建模 。 MIRA 从而在一个有效且可解释的潜在空间中代表细胞状态，推断高保真谱系树，确定分支点命运决定的关键调节因子，并揭示局部可及性对不同位点转录的可变影响 。将 MIRA 应用于表皮维持分化和大脑发育系统，通过多模式单细胞 RNA 测序（scRNA-seq）和来自两个不同平台（SHARE-seq 和 10x ）。在每个系统中，MIRA 构建了一个高保真发展轨迹，并确定了在谱系分支点驱动关键命运决定的调控因素。此外， MIRA 区分了受局部染色质landscope严格时空调节的早期发育基因与允许在对局部染色质影响最小的因素滴定时保持可接近的终末命运基因，揭示了可变调节电路如何协调命运承诺和终末身份 。

MIRA 利用单细胞中表达和染色质可及性的joint topic modeling和 RP 建模来确定驱动发育连续体中关键命运决定的调节机制 。 Probabilistic topic modeling has been employed in natural language understanding to elucidate the abstract topics that shape the meaning of a given collection of text 。 最近，opic modeling已分别应用于 scRNA-seq 和 scATAC-seq，分别将转录或表观遗传细胞状态描述为共调控基因或顺式调控元件的“themat”组。

MIRA 的 topic model使用变分自编码器方法， 将深度学习与概率图形模型相结合 ，以学习定义每个细胞身份的表达和accessibility topics。

MIRA 通过对过度分散的 scRNA-seq 计数和稀疏的 scATAC-seq 数据使用不同的生成分布来解释每种模式的不同统计特性 。 使用稀疏约束来确保细胞的topics组成是连贯的和可解释的。 MIRA 的 贝叶斯超参数调整方案 找到了全面但非冗余地描述每个数据集所需的适当数量的topics。

MIRA 接下来将表达式和可及性topics组合成一个联合表示，用于计算 k 最近邻 (KNN) 图 。 然后利用开发的一种新方法， 利用 KNN 图构建高保真谱系树，以定义谱系之间的分支点，其中分化为一个终端状态的概率与另一个不同 。 然后，MIRA 对比了映射在此谱系树上的表达和可及性topics的出现，以 阐明在推断的谱系分支点处驱动命运决定的关键调节因素 。

接下来，MIRA 利用 RP 建模以单个基因位点的分辨率整合转录组和可及性数据，以确定每个基因周围的调控元件如何影响其表达 。 虽然染色质可及性和表达之间的相关性被归因于细胞状态的全基因组协调变化所混淆，但基因组可及性表明顺式调控元件和转录之间存在机械调控关系 。因此，顺式调控元件的感知影响被建模为随着转录起始位点 (TSS) 上游或下游的基因组距离以 MIRA 从多模式数据中学习到的独立速率呈指数衰减。每个基因的 RP 被评分为单个调控元件贡献的总和。 MIRA 通过检查预测因子中的转录因子基序富集或占有率（如果提供染色质免疫沉淀测序（ChIP-seq）数据）预测每个基因座的关键调节因子，这些元素通过概率性计算机内缺失（pISD）预测会高度影响该基因座的转录 。

此外，MIRA 通过将局部 RP 模型与第二个扩展模型进行比较，量化了局部染色质可及性对基因表达的调节影响，该模型增强了 MIRA 的可及性topics编码的全基因组可及性状态的内容 。 其转录被 RP 模型充分预测的基因仅基于局部可及性（来自 TSS 的±600 千碱基）被定义为局部染色质可及性影响转录表达 (LITE) 基因。 其表达被具有全基因组范围的模型明显更好地描述的基因被定义为非局部染色质可及性影响的 转录表达 (NITE) 基因。

虽然 LITE 基因似乎受到局部染色质可及性的严格调控，但 NITE 基因的转录似乎不需要大量的局部染色质重塑就可以确认 。 MIRA 将 LITE 模型高估或低估每个细胞中表达的程度定义为“染色质差异”， 突出显示转录与局部染色质可及性变化脱钩的细胞 。 MIRA 检查整个发育连续体的染色质差异，以揭示可变网络如何调节ate commitment和终端身份。

应用于通过 SHARE-seq3 检测的毛囊维持分化，MIRA 的joint topic表示构建了一个谱系图，其潜在结构模仿了毛囊的真实空间布局。IRA 推断的谱系树重建了毛囊谱系的祖先层次结构，外根鞘细胞导致早期基质祖细胞，随后分支为后代内根鞘 (IRS)，然后是髓质和皮质谱系。 准确的谱系树是确定在谱系分支点指导细胞命运决定的因素的关键先决条件 。

MIRA 使用Stream graphs对比了推断谱系树中的表达流和可及性topics，以揭示沿不同发育路径驱动细胞命运的调节模块 。 Stream graphs支持沿连续体的高维、多模态比较 。 表达topic e2 捕获了控制祖基质细胞的转录状态，包括细胞增殖和 Eda 和 Shh 信号。 此后，表达topic e6 描述了对应于 Notch 相关因子激活的皮层规范。 相反，表达topic e4 表征了髓质规范，包含与髓质特异性可及性topic a5 中 Smad5/Smad2/3 motif的富集一致的 Bmp/Tgf-β 相关因子。 comparison with cortex-specific accessibility topic a6 showed both lineages were enriched for motifs bound by canonical hair shaft regulators Lef1 and Hoxc, with expression implicating the influence of Hoxc13

Contrasting modalities，Wnt 驱动的可及性 topics a4 描述了在髓质和皮质谱系之间的分支点处的短暂可及性状态，没有任何相应的表达topic 。 因此，祖细胞基质细胞中的细胞水平染色质重塑先于指定每个下游谱系的转录改变。

虽然毛囊中的大多数基因表现出 LITE 调节，局部可及性与转录同步增加，但表达与 Krt23 等基因的 LITE 模型预测的不同 。 尽管局部染色质可及性对 Krt23 表达的预测很差，但其转录是谱系特异性的，并且与可及性topics a5 的激活密切相关，编码了髓质全基因组的可及性模式。 一致地，NITE 模型更接近地预测了 Krt23 表达，该模型包括这些全基因组可及性状态作为特征。 LITE 基因因此受到局部染色质重塑的严格调控，而 NITE 基因的滴定不需要大量的局部染色质重塑，将转录与局部可及性分离。

MIRA 的“染色质差异”描绘了局部可及性与整个发育轨迹中的转录分离的程度 。 尽管 Krt23 局部可及性在髓质和皮质谱系之间的分支点增加并且在两者中都保持升高，但它最终仅在髓质中高度表达，导致高染色质差异，高估了其在皮质中的表达。 尽管在两个谱系中均可访问， Krt23’s lineage-specific expression despite accessibility in both lineages suggests its activation requires addition of a factor that does not primarily impact transcription via remodeling local accessibility 。

在细胞水平上，终末分化的髓质和皮质细胞中的基因表达比毛囊分化早期的基因表达表现出更多的 NITE 调节（p<005，Wilcoxon） 。 通常，终末表达基因的可及性在fate commitment之前增加并在随后的两个谱系中保持，但仅在髓质和皮质之间的分支点之后以谱系特异性方式激活表达。 在分支点使用染色质差异来识别具有这些“分支引发”动态的基因。 虽然启动表明表达的必然性，但这些基因表明启动基因座的后续表达可能是有条件的，这种模式被检测为强 NITE 调控 。

细胞级topics建模还支持fate commitment之前的启动可及性模式 。 例如，最终在前一个分支点启动可及性的皮层中表达的基因的动态由皮层特异性表达topics e6 和跨分支可及性topics a4 描述。 如前所述， 可及性topics a4 描述了染色质状态的全细胞变化，这与表达topics影响的同步变化不符 。

随后在髓质或皮质中条件性表达的分支引发基因似乎对髓质或皮质fate commitment的调节作出反应。 MIRA pISD 暗示 Notch 效应器 Rbpj 作为分支引发的皮层基因的顶部调节因子，Bmp/Tgf-β-诱导的 Smad5/Smad2/3 作为分支引发的髓质基因的调节因子，与与这些因子诱导相关的基因表达一致。 因此，MIRA 确定分支点的细胞具有允许多种命运的染色质状态，由短暂的可及性topics a4 描述，最终通过随后添加 fate-defining signal，即 Bmp/Tgf-β 或 Notch。

接下来将 MIRA 应用于同一数据集 3 中的一个单独系统，即滤泡间表皮 (IFE)。 分化的两个空间轴指定 IFE，一个控制基底干细胞分化为越来越浅的表皮层（表皮分层轴），另一个控制基底细胞内陷和滤泡形成（滤泡轴）。 MIRA joint topic representation的潜在结构再次模仿了这个分化系统的空间布局，重构了两个分化轴 。

此外，与先前报告的对该数据集 3 没有联合建模表达和可及性的分析不同，MIRA 确定了两个不同的基底-棘突-颗粒轨迹 。 一种标记为“中间”的轨迹在转录和表观遗传上与上毛囊结构更加相似，表明这些细胞在空间上靠近毛囊并受到更多的促毛囊调节。 这些“中间”基底细胞显示出 Egr2 表达和基序的激活，之前与表皮增殖和伤口愈合有关。 相比之下，远离毛囊的基底细胞表现出更强的 Thbs1 表达，这与先前的工作 一致，该工作确定了两个不同的基底细胞群，其中 Thbs1 标记远离毛囊的基底细胞。 这两个不同的基底细胞niche中的每一个都产生了自己的表皮层列，由 MIRA joint topics modeling捕获 。

从joint topic representation中推断出的谱系树揭示了共同的和谱系特定的branch points，它们塑造了两条基底棘颗粒轨迹的空间程序 。 通过可及性topics可视化状态变化确定了 Hes 对基底细胞的共同调节影响，其次是 Pou2f3 对棘细胞的影响，并在两个轨迹中以 Grhl 和 C/ebp 终止于颗粒细胞。 相比之下，谱系特定的可及性topics区分了 Klf4 基序在“中间”棘状细胞和颗粒细胞中的影响，而不是 Gata3 对颗粒细胞中的影响，这些颗粒细胞由距毛囊更远的 Thbs1+ 基底细胞产生。

观察到终端populations中的表达显着富集 NITE 调控，尤其是在谱系之间差异表达的终端基因中（p<005，Wilcoxon）。 同样， 最终命运染色质可及性似乎指定了可用的细胞状态，而转录最终取决于额外的空间或信号队列 。 总体而言， MIRA 阐明了沿两条平行轨迹的共享和谱系特异性分化机制，在 IFE 内具有不同的空间调节 。

接下来将 MIRA 应用于胚胎大脑数据集（在不同平台上进行检测，10x Multiome）以确定驱动抑制性和兴奋性神经元命运之间决定的关键因素，其平衡对正常大脑发育至关重要。 与对表达和可及性的单独分析相比， MIRA joint topic modeling 构建了一个高保真谱系树，展示了从共同祖细胞到星形胶质细胞或神经元祖细胞的连续流动，随后分支为抑制性或兴奋性神经元 。

MIRA 的高保真谱系树对于精确对比整个发展连续体中的表达动态和可及性至关重要 。 定义对神经元祖细胞而非星形胶质细胞的承诺的主要表达topics e3 富含细胞周期基因，这可能反映了在承诺抑制或兴奋性命运之前神经元群的扩张。 随后细胞周期topics e3 的失活与可及性topics a2 的激活一致，后者富含 Ascl1 的motif，Ascl1 是已知促进细胞周期退出和分化的神经祖细胞中的开创性转录因子。 Ascl1 基序也在早期抑制性可及性topic a4 中得到丰富，与 Ascl1 在抑制性神经元分化中的关键作用一致。

在通往兴奋性命运的替代轨迹中，早期兴奋性可及性topic a6 表明 Ascl1 motif的消耗与 Neurod1 motif的增加相协调。 Ascl1 和 Neurod1 属于不同的基本螺旋-环-螺旋转录因子亚组； Neurod1 促进分化为兴奋性神经元，而 Ascl1 指定抑制性神经元。

MIRA topics对比了由抑制性驱动 Ascl1 或兴奋性驱动 Neurod1 启动的规范的时间进程。 抑制轨迹激活了抑制成熟驱动 Bdnf 信号，最终激活了定义终端抑制命运的 GABA 突触成分（topic e13）。一致地，对齐的终端抑制可及性topic a10 富含 Egr1 motif，下游 Bdnf 效应器直接 激活 GABA 能神经传递基因。

发散的兴奋性分支首先激活对支持神经元代谢需求很重要的线粒体成分（topic e14），然后是谷氨酸能突触机械的末端激活，包括唯一区分兴奋性神经元的谷氨酸转运蛋白（topic e20）。 鉴于 Mef2c 在兴奋性分支中的表达，对齐的终端兴奋性可及性topic a13 被丰富了可归因于 Mef2c 的 Mef2 motif，这与其已知的通过促进兴奋性分化维持兴奋性/抑制性平衡的作用一致 。

为了确定胚胎大脑中的 LITE 和 NITE 基因， 为定义每个表达和可及性topics的基因训练了 MIRA RP 模型 。 值得注意的 LITE 基因包括具有严格时空调控的命运驱动转录因子，例如常见的祖基因 Pax6 和促进兴奋的 Mef2c。 相反，NITE 基因丰富了细胞周期机制以及由神经递质和离子通道基因组成的神经元分化基因组。 先前已报道局部染色质景观对细胞周期基因的激活贡献有限，与 NITE 调节一致。 这可能反映了对控制每个细胞周期阶段的基因表达的要求，这与重塑局部染色质landscope所需的时间不相容。 同样，突触维持和可塑性可能需要神经递质和离子通道基因的快速反应调节，反映为 NITE 调节 。

Analogously to the hair follicle and IFE, expression topics describing the neuronal

progenitors were significantly enriched for LITE regulation, whereas after commitment to the excitatory or inhibitory fate, topics were significantly enriched for NITE regulation (p<005, Wilcoxon) Neuronal progenitor and early inhibitory regulator Ascl1 is known to be a pioneering transcription factor that remodels the chromatin landscape to regulate its targets By contrast, terminal inhibitory regulator Egr1 was previously reported to have non-pioneer-like properties40 Notably, targets predicted by MIRA pISD to be downstream of Ascl1 demonstrated significantly stronger LITE regulation than predicted Egr1 targets, potentially reflective of local chromatin remodeling by pioneering Ascl1 driving their expression (p<005, Wilcoxon)

总之，MIRA 利用ell-level topic modeling和基因级 RP 建模来严格对比单细胞转录的时空动态与染色质可及性，以揭示这些机制如何相互作用以协调发育轨迹中的关键命运决定 。 在我们分析的每个系统中，MIRA 展示了表达和可及性数据joint topics modeling的能力，以推断准确捕获复杂时空分化轴的高保真谱系树。 将表达和可及性topics映射到联合谱系树上，阐明了在关键谱系分支点驱动命运决定的关键branch points。

MIRA 对比了转录动力学和局部染色质可及性，以定义每个基因位点的染色质差异，揭示主要受 LITE 或 NITE 机制调节的离散基因模块 。 有趣的是，在我们从皮肤和大脑数据集测试的所有三个系统中，早期表达的基因都被丰富用于 LITE 调节。 早期表达基因的 LITE 调控可能反映了严格调控的重要性，需要对其表达进行广泛的染色质重塑，然后在其异常表达将产生破坏性后果的命运中强行沉默。 相反，对于维持终末细胞功能很重要的基因电池不太依赖局部染色质重塑来进行调节，这表明细胞信号传导等机制的影响更大，这些机制允许转录表达以满足波动的细胞需求。

在 NITE 调节的基因中，我们还注意到在谱系分支点具有启动可及性的基因，这些基因显示出随后的谱系特异性激活，以响应诸如信号的命运定义力，大概是通过结合或激活对局部可及性影响最小的因子 在这些情况下，可及性似乎反映了编码细胞可用转录状态的塑料细胞身份， 其中最终转录响应细胞的空间或信号生态位 。 未来的工作需要进一步确定细胞何时使用 LITE 与 NITE 机制来调节不同细胞过程的逻辑。

MIRA（用于综合调控分析的概率多模态模型）是一种综合方法，它系统地对比了单细胞转录和可及性，以推断沿着发育轨迹驱动细胞的调控网络。

MIRA leverages joint topic modeling of cell states and regulatory potential modeling at individual gene loci to:

MIRA harnesses a variational autoencoder approach to model both transcription and chromatin accessibility topics defining each cell’s identity while accounting for their distinct statistical properties and employing a sparsity constraint to ensure topics are coherent and interpretable MIRA’s hyperparameter tuning scheme learns the appropriate number of topics needed to comprehensively yet non-redundantly describe each dataset MIRA next combines the expression and accessibility topics into a joint representation used to calculate a k-nearest neighbors (KNN) graph This output can then be leveraged for visualization and clustering, construction of high fidelity lineage trajectories, and rigorous topic analysis to determine regulators driving key fate decisions at lineage branch points

MIRA’s regulatory potential (RP) model integrates transcriptional and chromatin accessibility data at each gene locus to determine how regulatory elements surrounding each gene influence its expression Regulatory influence of enhancers is modeled to decay exponentially with genomic distance at a rate learned by the MIRA RP model from the joint multimodal data MIRA learns independent upstream and downstream decay rates and includes parameters to weigh upstream, downstream, and promoter effects The RP of each gene is scored as the sum of the contribution of individual regulatory elements MIRA predicts key regulators at each locus by examining transcription factor motif enrichment or occupancy (if provided chromatin immunoprecipitation (ChIP-seq) data) within elements predicted to highly influence transcription at that locus using probabilistic in silico deletion (ISD)

MIRA quantifies the regulatory influence of local chromatin accessibility by comparing the local RP model with a second, expanded model that augments the local RP model with genome-wide accessibility states encoded by MIRA’s chromatin accessibility topics Genes whose expression is significantly better described by this expanded model are defined as non-local chromatin accessibility-influenced transcriptional expression (NITE) genes Genes whose transcription is sufficiently predicted by the RP model based on local accessibility alone are defined as local chromatin accessibility-influenced transcriptional expression (LITE) genes While LITE genes appear tightly regulated by local chromatin accessibility, the transcription of NITE genes appears to be titrated without requiring extensive local chromatin remodeling MIRA defines the extent to which the LITE model over- or under-estimates expression in each cell as “chromatin differential”, highlighting cells where transcription is decoupled from shifts in local chromatin accessibility MIRA examines chromatin differential across the developmental continuum to reveal how variable circuitry regulates fate commitment and terminal identity

生活很好，有你更好

机器学习算法与Python实践这个系列主要是参考《机器学习实战》这本书。因为自己想学习Python，然后也想对一些机器学习算法加深下了解，所以就想通过Python来实现几个比较常用的机器学习算法。恰好遇见这本同样定位的书籍，所以就参考这本书的过程来学习了。

一、kNN算法分析

K最近邻（k-Nearest Neighbor，KNN）分类算法可以说是最简单的机器学习算法了。它采用测量不同特征值之间的距离方法进行分类。它的思想很简单：如果一个样本在特征空间中的k个最相似（即特征空间中最邻近）的样本中的大多数属于某一个类别，则该样本也属于这个类别。

比如上面这个图，我们有两类数据，分别是蓝色方块和红色三角形，他们分布在一个上图的二维中间中。那么假如我们有一个绿色圆圈这个数据，需要判断这个数据是属于蓝色方块这一类，还是与红色三角形同类。怎么做呢？我们先把离这个绿色圆圈最近的几个点找到，因为我们觉得离绿色圆圈最近的才对它的类别有判断的帮助。那到底要用多少个来判断呢？这个个数就是k了。如果k=3，就表示我们选择离绿色圆圈最近的3个点来判断，由于红色三角形所占比例为2/3，所以我们认为绿色圆是和红色三角形同类。如果k=5，由于蓝色四方形比例为3/5，因此绿色圆被赋予蓝色四方形类。从这里可以看到，k的值还是很重要的。

KNN算法中，所选择的邻居都是已经正确分类的对象。该方法在定类决策上只依据最邻近的一个或者几个样本的类别来决定待分样本所属的类别。由于KNN方法主要靠周围有限的邻近的样本，而不是靠判别类域的方法来确定所属类别的，因此对于类域的交叉或重叠较多的待分样本集来说，KNN方法较其他方法更为适合。

该算法在分类时有个主要的不足是，当样本不平衡时，如一个类的样本容量很大，而其他类样本容量很小时，有可能导致当输入一个新样本时，该样本的K个邻居中大容量类的样本占多数。因此可以采用权值的方法（和该样本距离小的邻居权值大）来改进。该方法的另一个不足之处是计算量较大，因为对每一个待分类的文本都要计算它到全体已知样本的距离，才能求得它的K个最近邻点。目前常用的解决方法是事先对已知样本点进行剪辑，事先去除对分类作用不大的样本。该算法比较适用于样本容量比较大的类域的自动分类，而那些样本容量较小的类域采用这种算法比较容易产生误分[参考机器学习十大算法]。

总的来说就是我们已经存在了一个带标签的数据库，然后输入没有标签的新数据后，将新数据的每个特征与样本集中数据对应的特征进行比较，然后算法提取样本集中特征最相似（最近邻）的分类标签。一般来说，只选择样本数据库中前k个最相似的数据。最后，选择k个最相似数据中出现次数最多的分类。其算法描述如下：

1）计算已知类别数据集中的点与当前点之间的距离；

2）按照距离递增次序排序；

3）选取与当前点距离最小的k个点；

4）确定前k个点所在类别的出现频率；

5）返回前k个点出现频率最高的类别作为当前点的预测分类。

二、Python实现

对于机器学习而已，Python需要额外安装三件宝，分别是Numpy，scipy和Matplotlib。前两者用于数值计算，后者用于画图。安装很简单，直接到各自的官网下载回来安装即可。安装程序会自动搜索我们的python版本和目录，然后安装到python支持的搜索路径下。反正就python和这三个插件都默认安装就没问题了。

另外，如果我们需要添加我们的脚本目录进Python的目录（这样Python的命令行就可以直接import），可以在系统环境变量中添加：PYTHONPATH环境变量，值为我们的路径，例如：E:\Python\Machine Learning in Action

21、kNN基础实践

一般实现一个算法后，我们需要先用一个很小的数据库来测试它的正确性，否则一下子给个大数据给它，它也很难消化，而且还不利于我们分析代码的有效性。

首先，我们新建一个kNNpy脚本文件，文件里面包含两个函数，一个用来生成小数据库，一个实现kNN分类算法。代码如下：

[python] view plain copy

#########################################

# kNN: k Nearest Neighbors

# Input: newInput: vector to compare to existing dataset (1xN)

# dataSet: size m data set of known vectors (NxM)

# labels: data set labels (1xM vector)

# k: number of neighbors to use for comparison

# Output: the most popular class label

#########################################

from numpy import

import operator

# create a dataset which contains 4 samples with 2 classes

def createDataSet():

# create a matrix: each row as a sample

group = array([[10, 09], [10, 10], [01, 02], [00, 01]])

labels = ['A', 'A', 'B', 'B'] # four samples and two classes

return group, labels

# classify using kNN

def kNNClassify(newInput, dataSet, labels, k):

numSamples = dataSetshape[0] # shape[0] stands for the num of row

## step 1: calculate Euclidean distance

# tile(A, reps): Construct an array by repeating A reps times

# the following copy numSamples rows for dataSet

diff = tile(newInput, (numSamples, 1)) - dataSet # Subtract element-wise

squaredDiff = diff 2 # squared for the subtract

squaredDist = sum(squaredDiff, axis = 1) # sum is performed by row

distance = squaredDist 05

## step 2: sort the distance

# argsort() returns the indices that would sort an array in a ascending order

sortedDistIndices = argsort(distance)

classCount = {} # define a dictionary (can be append element)

for i in xrange(k):

## step 3: choose the min k distance

voteLabel = labels[sortedDistIndices[i]]

## step 4: count the times labels occur

# when the key voteLabel is not in dictionary classCount, get()

# will return 0

classCount[voteLabel] = classCountget(voteLabel, 0) + 1

## step 5: the max voted class will return

maxCount = 0

for key, value in classCountitems():

if value > maxCount:

maxCount = value

maxIndex = key

return maxIndex

然后我们在命令行中这样测试即可：

[python] view plain copy

import kNN

from numpy import

dataSet, labels = kNNcreateDataSet()

testX = array([12, 10])

k = 3

outputLabel = kNNkNNClassify(testX, dataSet, labels, 3)

print "Your input is:", testX, "and classified to class: ", outputLabel

testX = array([01, 03])

outputLabel = kNNkNNClassify(testX, dataSet, labels, 3)

print "Your input is:", testX, "and classified to class: ", outputLabel

这时候会输出：

[python] view plain copy

Your input is: [ 12 10] and classified to class: A

Your input is: [ 01 03] and classified to class: B

22、kNN进阶

这里我们用kNN来分类一个大点的数据库，包括数据维度比较大和样本数比较多的数据库。这里我们用到一个手写数字的数据库，可以到这里下载。这个数据库包括数字0-9的手写体。每个数字大约有200个样本。每个样本保持在一个txt文件中。手写体图像本身的大小是32x32的二值图，转换到txt文件保存后，内容也是32x32个数字，0或者1，如下：

数据库解压后有两个目录：目录trainingDigits存放的是大约2000个训练数据，testDigits存放大约900个测试数据。

这里我们还是新建一个kNNpy脚本文件，文件里面包含四个函数，一个用来生成将每个样本的txt文件转换为对应的一个向量，一个用来加载整个数据库，一个实现kNN分类算法。最后就是实现这个加载，测试的函数。

[python] view plain copy

#########################################

# kNN: k Nearest Neighbors

# Input: inX: vector to compare to existing dataset (1xN)

# dataSet: size m data set of known vectors (NxM)

# labels: data set labels (1xM vector)

# k: number of neighbors to use for comparison

# Output: the most popular class label

#########################################

from numpy import

import operator

import os

# classify using kNN

def kNNClassify(newInput, dataSet, labels, k):

numSamples = dataSetshape[0] # shape[0] stands for the num of row

## step 1: calculate Euclidean distance

# tile(A, reps): Construct an array by repeating A reps times

# the following copy numSamples rows for dataSet

diff = tile(newInput, (numSamples, 1)) - dataSet # Subtract element-wise

squaredDiff = diff 2 # squared for the subtract

squaredDist = sum(squaredDiff, axis = 1) # sum is performed by row

distance = squaredDist 05

## step 2: sort the distance

# argsort() returns the indices that would sort an array in a ascending order

sortedDistIndices = argsort(distance)

classCount = {} # define a dictionary (can be append element)

for i in xrange(k):

## step 3: choose the min k distance

voteLabel = labels[sortedDistIndices[i]]

## step 4: count the times labels occur

# when the key voteLabel is not in dictionary classCount, get()

# will return 0

classCount[voteLabel] = classCountget(voteLabel, 0) + 1

## step 5: the max voted class will return

maxCount = 0

for key, value in classCountitems():

if value > maxCount:

maxCount = value

maxIndex = key

return maxIndex

# convert image to vector

def img2vector(filename):

rows = 32

cols = 32

imgVector = zeros((1, rows cols))

fileIn = open(filename)

for row in xrange(rows):

lineStr = fileInreadline()

for col in xrange(cols):

imgVector[0, row 32 + col] = int(lineStr[col])

return imgVector

# load dataSet

def loadDataSet():

## step 1: Getting training set

print "---Getting training set"

dataSetDir = 'E:/Python/Machine Learning in Action/'

trainingFileList = oslistdir(dataSetDir + 'trainingDigits') # load the training set

numSamples = len(trainingFileList)

train_x = zeros((numSamples, 1024))

train_y = []

for i in xrange(numSamples):

filename = trainingFileList[i]

# get train_x

train_x[i, :] = img2vector(dataSetDir + 'trainingDigits/%s' % filename)

# get label from file name such as "1_18txt"

label = int(filenamesplit('_')[0]) # return 1

train_yappend(label)

## step 2: Getting testing set

print "---Getting testing set"

testingFileList = oslistdir(dataSetDir + 'testDigits') # load the testing set

numSamples = len(testingFileList)

test_x = zeros((numSamples, 1024))

test_y = []

for i in xrange(numSamples):

filename = testingFileList[i]

# get train_x

test_x[i, :] = img2vector(dataSetDir + 'testDigits/%s' % filename)

# get label from file name such as "1_18txt"

label = int(filenamesplit('_')[0]) # return 1

test_yappend(label)

return train_x, train_y, test_x, test_y

# test hand writing class

def testHandWritingClass():

## step 1: load data

print "step 1: load data"

train_x, train_y, test_x, test_y = loadDataSet()

## step 2: training

print "step 2: training"

pass

## step 3: testing

print "step 3: testing"

numTestSamples = test_xshape[0]

matchCount = 0

for i in xrange(numTestSamples):

predict = kNNClassify(test_x[i], train_x, train_y, 3)

if predict == test_y[i]:

matchCount += 1

accuracy = float(matchCount) / numTestSamples

## step 4: show the result

print "step 4: show the result"

print 'The classify accuracy is: %2f%%' % (accuracy 100)

测试非常简单，只需要在命令行中输入：

[python] view plain copy

import kNN

kNNtestHandWritingClass()

输出结果如下：

[python] view plain copy

step 1: load data

---Getting training set

---Getting testing set

step 2: training

step 3: testing

step 4: show the result

The classify accuracy is: 9884%

给定测试实例，基于某种距离度量找出训练集中与其最靠近的k个实例点，然后基于这k个最近邻的信息来进行预测。

通常，在分类任务中可使用“投票法”，即选择这k个实例中出现最多的标记类别作为预测结果；在回归任务中可使用“平均法”，即将这k个实例的实值输出标记的平均值作为预测结果；还可基于距离远近进行加权平均或加权投票，距离越近的实例权重越大。

k近邻法不具有显式的学习过程，事实上，它是懒惰学习（lazy learning）的著名代表，此类学习技术在训练阶段仅仅是把样本保存起来，训练时间开销为零，待收到测试样本后再进行处理。

KNN一般采用欧氏距离，也可采用其他距离度量，一般的Lp距离：

KNN中的K值选取对K近邻算法的结果会产生重大影响。如果选择较小的K值，就相当于用较小的领域中的训练实例进行预测，“学习”近似误差（近似误差：可以理解为对现有训练集的训练误差）会减小，只有与输入实例较近或相似的训练实例才会对预测结果起作用，与此同时带来的问题是“学习”的估计误差会增大，换句话说，K值的减小就意味着整体模型变得复杂，容易发生过拟合；

如果选择较大的K值，就相当于用较大领域中的训练实例进行预测，其优点是可以减少学习的估计误差，但缺点是学习的近似误差会增大。这时候，与输入实例较远（不相似的）训练实例也会对预测器作用，使预测发生错误，且K值的增大就意味着整体的模型变得简单。

在实际应用中，K值一般取一个比较小的数值，例如采用交叉验证法来选择最优的K值。经验规则：k一般低于训练样本数的平方根

1、计算测试对象到训练集中每个对象的距离

2、按照距离的远近排序

3、选取与当前测试对象最近的k的训练对象，作为该测试对象的邻居

4、统计这k个邻居的类别频率

5、k个邻居里频率最高的类别，即为测试对象的类别

输入X可以采用BallTree或KDTree两种数据结构，优化计算效率，可以在实例化KNeighborsClassifier的时候指定。

KDTree

基本思想是，若A点距离B点非常远，B点距离C点非常近， 可知A点与C点很遥远，不需要明确计算它们的距离。 通过这样的方式，近邻搜索的计算成本可以降低为O[DNlog(N)]或更低。 这是对于暴力搜索在大样本数N中表现的显著改善。KD 树的构造非常快，对于低维度 (D<20) 近邻搜索也非常快, 当D增长到很大时，效率变低： 这就是所谓的 “维度灾难” 的一种体现。

KD 树是一个二叉树结构，它沿着数据轴递归地划分参数空间，将其划分为嵌入数据点的嵌套的各向异性区域。 KD 树的构造非常快：因为只需沿数据轴执行分区, 无需计算D-dimensional 距离。 一旦构建完成, 查询点的最近邻距离计算复杂度仅为O[log(N)]。 虽然 KD 树的方法对于低维度 (D<20) 近邻搜索非常快, 当D增长到很大时, 效率变低。

KD树的特性适合使用欧氏距离。

BallTree

BallTree解决了KDTree在高维上效率低下的问题，这种方法构建的树要比 KD 树消耗更多的时间，但是这种数据结构对于高结构化的数据是非常有效的， 即使在高维度上也是一样。

KD树是依次对K维坐标轴，以中值切分构造的树；ball tree 是以质心C和半径r分割样本空间，每一个节点是一个超球体。换句简单的话来说，对于目标空间(q, r)，所有被该超球体截断的子超球体内的所有子空间都将被遍历搜索。

BallTree通过使用三角不等式减少近邻搜索的候选点数:|x+y|<=|x|+|y|通过这种设置, 测试点和质心之间的单一距离计算足以确定距节点内所有点的距离的下限和上限 由于 ball 树节点的球形几何, 它在高维度上的性能超出 KD-tree, 尽管实际的性能高度依赖于训练数据的结构。

BallTree适用于更一般的距离。

1、优点

非常简单的分类算法没有之一，人性化，易于理解，易于实现

适合处理多分类问题，比如推荐用户

可用于数值型数据和离散型数据，既可以用来做分类也可以用来做回归

对异常值不敏感

2、缺点

属于懒惰算法，时间复杂度较高，因为需要计算未知样本到所有已知样本的距离

样本平衡度依赖高，当出现极端情况样本不平衡时，分类绝对会出现偏差，可以调整样本权值改善

可解释性差，无法给出类似决策树那样的规则

向量的维度越高，欧式距离的区分能力就越弱

样本空间太大不适合，因为计算量太大，预测缓慢

文本分类

用户推荐

回归问题

1）所有的观测实例中随机抽取出k个观测点，作为聚类中心点，然后遍历其余的观测点找到距离各自最近的聚类中心点，将其加入到该聚类中。这样，我们就有了一个初始的聚类结果，这是一次迭代的过程。

2）我们每个聚类中心都至少有一个观测实例，这样，我们可以求出每个聚类的中心点（means），作为新的聚类中心，然后再遍历所有的观测点，找到距离其最近的中心点，加入到该聚类中。然后继续运行2）。

3）如此往复2），直到前后两次迭代得到的聚类中心点一模一样。

本算法的时间复杂度：O(tkmn)，其中，t为迭代次数，k为簇的数目，m为记录数，n为维数；

空间复杂度：O((m+k)n)，其中，k为簇的数目，m为记录数，n为维数。

适用范围：

K-menas算法试图找到使平凡误差准则函数最小的簇。当潜在的簇形状是凸面的，簇与簇之间区别较明显，且簇大小相近时，其聚类结果较理想。前面提到，该算法时间复杂度为O(tkmn)，与样本数量线性相关，所以，对于处理大数据集合，该算法非常高效，且伸缩性较好。但该算法除了要事先确定簇数K和对初始聚类中心敏感外，经常以局部最优结束，同时对“噪声”和孤立点敏感，并且该方法不适于发现非凸面形状的簇或大小差别很大的簇。

1）首先，算法只能找到局部最优的聚类，而不是全局最优的聚类。而且算法的结果非常依赖于初始随机选择的聚类中心的位置。我们通过多次运行算法，使用不同的随机生成的聚类中心点运行算法，然后对各自结果C通过evaluate(C)函数进行评估，选择多次结果中evaluate(C)值最小的那一个。k-means++算法选择初始seeds的基本思想就是：初始的聚类中心之间的相互距离要尽可能的远

2）关于初始k值选择的问题。首先的想法是，从一个起始值开始，到一个最大值，每一个值运行k-means算法聚类，通过一个评价函数计算出最好的一次聚类结果，这个k就是最优的k。我们首先想到了上面用到的evaluate(C)。然而，k越大，聚类中心越多，显然每个观测点距离其中心的距离的平方和会越小，这在实践中也得到了验证。第四节中的实验结果分析中将详细讨论这个问题。

3）关于性能问题。原始的算法，每一次迭代都要计算每一个观测点与所有聚类中心的距离。有没有方法能够提高效率呢？是有的，可以使用k-d tree或者ball tree这种数据结构来提高算法的效率。特定条件下，对于一定区域内的观测点，无需遍历每一个观测点，就可以把这个区域内所有的点放到距离最近的一个聚类中去。这将在第三节中详细地介绍。

相似点：都包含这样的过程，给定一个点，在数据集中找离它最近的点。即二者都用到了NN(Nears Neighbor)算法，一般用KD树来实现NN。

k-d tree 与 ball tree

1）k-d tree[5]

把n维特征的观测实例放到n维空间中，k-d tree每次通过某种算法选择一个特征(坐标轴)，以它的某一个值作为分界做超平面，把当前所有观测点分为两部分，然后对每一个部分使用同样的方法，直到达到某个条件为止。

上面的表述中，有几个地方下面将会详细说明：（1）选择特征（坐标轴）的方法 （2）以该特征的哪一个为界 （3）达到什么条件算法结束。

(1)选择特征的方法

计算当前观测点集合中每个特征的方差，选择方差最大的一个特征，然后画一个垂直于这个特征的超平面将所有观测点分为两个集合。

（2)以该特征的哪一个值为界 即垂直选择坐标轴的超平面的具体位置。

第一种是以各个点的方差的中值（median）为界。这样会使建好的树非常地平衡，会均匀地分开一个集合。这样做的问题是，如果点的分布非常不好地偏斜的，选择中值会造成连续相同方向的分割，形成细长的超矩形(hyperrectangles)。

替代的方法是计算这些点该坐标轴的平均值，选择距离这个平均值最近的点作为超平面与这个坐标轴的交点。这样这个树不会完美地平衡，但区域会倾向于正方地被划分，连续的分割更有可能在不同方向上发生。

（3）达到什么条件算法结束

实际中，不用指导叶子结点只包含两个点时才结束算法。你可以设定一个预先设定的最小值，当这个最小值达到时结束算法。

图6中，星号标注的是目标点，我们在k-d tree中找到这个点所处的区域后，依次计算此区域包含的点的距离，找出最近的一个点（黑色点），如果在其他region中还包含更近的点则一定在以这两个点为半径的圆中。假设这个圆如图中所示包含其他区域。先看这个区域兄弟结点对应区域，与圆不重叠；再看其双亲结点的兄弟结点对应区域。从它的子结点对应区域中寻找（图中确实与这个双亲结点的兄弟结点的子结点对应区域重叠了）。在其中找是否有更近的结点。

k-d tree的优势是可以递增更新。新的观测点可以不断地加入进来。找到新观测点应该在的区域，如果它是空的，就把它添加进去，否则，沿着最长的边分割这个区域来保持接近正方形的性质。这样会破坏树的平衡性，同时让区域不利于找最近邻。我们可以当树的深度到达一定值时重建这棵树。

然而，k-d tree也有问题。矩形并不是用到这里最好的方式。偏斜的数据集会造成我们想要保持树的平衡与保持区域的正方形特性的冲突。另外，矩形甚至是正方形并不是用在这里最完美的形状，由于它的角。如果图6中的圆再大一些，即黑点距离目标点点再远一些，圆就会与左上角的矩形相交，需要多检查一个区域的点，而且那个区域是当前区域双亲结点的兄弟结点的子结点。

为了解决上面的问题，我们引入了ball tree。

2）ball tree[4]

解决上面问题的方案就是使用超球面而不是超矩形划分区域。使用球面可能会造成球面间的重叠，但却没有关系。ball tree就是一个k维超球面来覆盖这些观测点，把它们放到树里面。图7（a)显示了一个2维平面包含16个观测实例的图,图7（b）是其对应的ball tree，其中结点中的数字表示包含的观测点数。

不同层次的圆被用不同的风格画出。树中的每个结点对应一个圆，结点的数字表示该区域保含的观测点数，但不一定就是图中该区域囊括的点数，因为有重叠的情况，并且一个观测点只能属于一个区域。实际的ball tree的结点保存圆心和半径。叶子结点保存它包含的观测点。

使用ball tree时，先自上而下找到包含target的叶子结点，从此结点中找到离它最近的观测点。这个距离就是最近邻的距离的上界。检查它的兄弟结点中是否包含比这个上界更小的观测点。方法是：如果目标点距离兄弟结点的圆心的距离大于这个圆的圆心加上前面的上界的值，则这个兄弟结点不可能包含所要的观测点。（如图8）否则，检查这个兄弟结点是否包含符合条件的观测点。

那么，ball tree的分割算法是什么呢？

选择一个距离当前圆心最远的观测点i1，和距离i1最远的观测点 i2，将圆中所有离这两个点最近的观测点都赋给这两个簇的中心，然后计算每一个簇的中心点和包含所有其所属观测点的最小半径。对包含n个观测点的超圆进行分割，只需要线性的时间。

与k-d tree一样，如果结点包含的观测点到达了预先设定的最小值，这个顶点就可以不再分割了。

起因： 最近有个问题样本，跑完cellranger，样本的cellranger结果如下，细胞数目极高(3W+)。在后续数据质控分析中，线粒体基因占比和双细胞率均很高，用scDblFinder进行双细胞预测，双细胞占率竟然高达34%。我很好奇，双细胞率为什么会这么高，如何审视这个结果？决定看看scDblFinder的细节。

问题1：如何理解Doublets？

在scRNA-seq的细胞捕获步骤中，两个或多个细胞聚集成单个液滴（双联体/多联体）会导致混合的转录组，也就是两个或多个细胞共用一个barcode，称为doublets或multiplets（后面统称为doublets）。它是基于液滴的单细胞测序的技术副产品。双细胞会造成每个“细胞”的高UMI计数，改变cluster的细胞类型鉴定干扰到下游分析。这会导致对稀有细胞类型、中间细胞状态和疾病相关转录组学特征的人为错误发现。双细胞率已被证明与捕获的细胞数量成正比（Bloom 2018; Kang et al 2018）。

双细胞可以分为同型（相同细胞类型）或异型（不同细胞类型）。

问题2： Cell Ranger可以自动剔除doublets和multiplets吗

答：目前没有方法可以识别与双细胞中单个细胞相关的转录本信息。

10X官网对双细胞率的 相关答复 ，我们目前没有一种方法通过算法识别单个物种的单细胞基因表达数据，barcode是否包含多个细胞。

目前，Cell Ranger 软件仅在barnyard实验或多物种实验中估计双细胞率。

对此，10X也给出三条参考意见：

他的意思是：1) 通过已知细胞类型的marker基因来鉴别双细胞，比如T/B细胞的marker基因，在同一个barcode细胞中同时高表达就可判定为双细胞；2) seurat标准分析流程，质控环节通过UMI和gene指标过滤；3）运用scDblFinder双细胞预测软件。

问题3：10X Genomics 单细胞实验中估计的双细胞率是多少？

假设不存在细胞结团，可使用下表（取自 10X基因组学用户指南 ）来估计单细胞实验中估计的双细胞率。

双细胞率（08%/1000cells），如果细胞数为1W，双细胞率为76%，约8%。

单细胞实验双细胞率为10-20%，这个数值明显高于上面10X给出的双细胞率(~8%)。这个怎么理解呢？

我想到的几个原因：

1）10X给出的是理想情况的双细胞率(细胞不结团)，用标准样本做基准比较；

2）双细胞率跟实验环节中的样本处理和细胞上样量都有关系；

3）还取决如何计算双细胞，双细胞率=双细胞数目/总细胞数；因计算的细胞总体不同而不同；

我们一般会进行QC细胞质控，用UMI/genes指标过滤掉低质量细胞和异常值细胞；

如果QC细胞质控后计算双细胞率，和QC细胞质控前计算双细胞率，预测的双细胞率会不一致；

最近，在网上找到一个10x Genomics 提供的估计双细胞率是：

比如1W个细胞，双细胞率为：0008(10000/1000)=008=8%。

下面我们看看scDblFinder软件是如何具体执行的。

双细胞在单细胞测序数据中很普遍，可能会导致人为错误的发现。

目前，实验层面还是无法检测同一样本的细胞形成的双细胞，包括异型双细胞。

算法层，已经开发了许多计算方法来根据转录谱识别双细胞。大多数这些方法依赖于通过对真实细胞求和或求平均来生成人造双细胞，并对它们与真实细胞之间的相似性进行评分。 例如，DoubletFinder在真实细胞和人工双细胞的合集上生成k最近邻近图(kNN) ，并估计每个细胞附近人造双细胞的密度（McGinnis, Murrow, and Gartner 2019）。 以类似的方式，Bais和Kostka (2020) 提出的bcds算法和共表达评分cxds算法。

Xi 和 Li (2021a) 最近发表的文章中对双细胞检测方法进行基准测试，使用模拟数据和包含双细胞实验标记的真实数据数据集，发现DoubletFinder的算法最为准确。 但是，基准测试也发现，没有一种方法在所有数据集上都是系统性最优，强调在各种数据集上测试和基准测试方法的必要性，并表明某些算法在不同情况下可能具有优势和劣势。

没有一种算法是完美不缺的，特别是预测模型算法。

下图比较了scDblFinder 这个包中一些方法（以粗体显示）与其他方法：

因此，我们在单细胞转录组数据质控过滤时，会考虑到双细胞的因素，通过相关软件进行预测双细胞。常用的软件有 scDblFinder (R语言)和 Scrublet (python)。这里仅讨论scDblFinder。

安装scDblFinder需要满足R >= 40 和 Bioconductor >= 312

scDblFinder的输入数据是 SingleCellExperiment 对象（空的drops已经移出），至少要包含counts矩阵(assay ‘counts’)。即sce对象都不应该包含空滴，但不应该经过非常严格的过滤（这会影响双细胞率的估计）。

如果还包含归一化矩阵 (assay ‘logcounts’) 和PCA (reducedDim ‘PCA’)，可以使用scDblFinder的cluster模式（不常见）。

对于 10x 数据，通常将dbr留空是安全的，它会自动估计。 scDblFinder的输出会在sce的colData中添加一些以‘scDblFinder’为前缀的列，其中最重要的是：

如果你有多个样本（理解为不同的细胞捕获），那么最好为每个样本分别进行双细胞识别（对于cell hashes实验中的多重样本，那意味着每个批次）。 可以通过简单地向scDblFinder的samples参数提供样本 id来完成，或者，将样本信息存储在colData列中，提供列名即可。另外，还可以考虑使用BPPARAM参数对其进行多线程处理（假设有足够的RAM）。 例如：

我们用之前案例中的数据测试下scDblFinder函数。

单个样本

该样本共10194个细胞，其中968个细胞被预测为双细胞，双细胞率为95%

那么该如何审视这个结果，选择怎样的参数？

我想到的是，对于一个预测模型来说，调参的意义不大，我们对结果没有预判，修改参数，都会出现不同的预测值。另外，对于一个常规的10X单细胞转录组数据，我们对双细胞率是有一定预判的，10X的实验步骤大致固定，cellranger的细胞数大致1W，双细胞率大致10%，我们知道预测的边界。我们有粗略的“标尺”。

但是现在，cellranger给出的细胞数是3W+，我们其实是不清楚双细胞率的边界，细胞数“超纲”了，只知道细胞数越多，双细胞也会越多。这类样本太少，我们没有横向可参考的实例。如果出现这种结果，最应该审视的是实验端出了什么问题，线粒体基因占比也非常高。

scDblFinder有两种生成人造双细胞的主要模式： 随机模式 （scDblFinderrandom，clusters=FALSE, 默认方式）和 基于cluster的模式 （scDblFinderclusters，clusters=TRUE 或提供你自定义的cluster - 以前版本的方法）。在实际中，我们观察到两种方法都表现良好（比其他方法要好）。当数据集被分成清晰的cluster时，教程建议使用基于cluster的方法（例如发展轨迹），否则使用随机模式。

双细胞分为 同型双细胞 （'Homotypic' doublets）和 异性双细胞 （ 'Heterotypic' doublets）。同型双细胞由相同类型的细胞（即相似的转录状态）组成，仅根据它们的转录组信息很难辨识。而且，它们对于大多数分析来说也相对无害，因为它们看起来与单细胞高度相似。 相反，异型双细胞（由具有不同转录状态的细胞形成）表现为一种人为的新型细胞类型，会影响下游分析。

scDblFinder只关注异性双细胞。

step1： 将数据集缩减到仅高表达的基因（默认为 1000）； 如果使用基于cluster的方法，则会选择每个cluster的表达靠前的基因。另外使用基于cluster的方法（而不是人为指定cluster），将会执行快速聚类（请参阅fastcluster）。

step2： 通过组合不同cluster的细胞来创建人工双细胞，创建的人工双细胞数与cluster的数目成比例。 我们主要关注不同cluster间的双细胞，我们不会试图识别同型双细胞，无论如何，它们实际上无法识别且对下游分析相对无害。因此，我们减少了人工双细胞的必要数量， 也防止分类器被训练以识别与单细胞无法区分的细胞（因此将单细胞称为双细胞）。 scDblFinder另一种策略是生成完全随机的人工双细胞，并使用迭代程序从训练中排除无法识别的人工双细胞。 在实践中，这两种方法具有相当的性能，它们也可以结合使用。

step3： 然后对真实细胞和人工双细胞的合集进行降维，并生成最近邻网络。 接着使用网络来估计每个细胞的许多特征，特别是最近邻居中人工双细胞的比例。 该比率不是选择特定的邻域大小，该比率是在不同的 k 值下计算的，通过使用多个预测变量创建分类器。预测变量还包括距离加权比，进一步添加了的细胞层面上的预测变量：对主成分的预测；文库大小； 和共表达分数（基于Bais 和 Kostka2020 的变体）。 然后 scDblFinder训练梯度提升树（ gradient boosted trees ），以根据这些特征区分来自真实细胞的人工双细胞。 最后，阈值程序通过同时最小化错误分类率和预期的双细胞率来决定调用细胞的分数（参见Thresholding）。

step4： 基于分类器方法的一个关键问题是一些真实细胞被错误标记，从某种意义上说，它们实际上是双细胞，但被标记为单细胞。这些会误导分类器。出于这个原因，分类和阈值处理以迭代方式执行：在每一轮中，从下一轮的训练集中删除识别为双细胞的真实细胞。

双细胞只能出现在给定的样本或某次捕获中，因此需要为每个样本单独进行双细胞判别，这也加快了分析速度。如果给定samples参数，scDblFinder将利用该参数将细胞拆分为单个样本/捕获，并在给出BPPARAM参数的情况下并行分析。分类器将在全局范围内进行训练，但阈值将在每个样本的基础上进行优化。如果你的样品是多标签，即不同的样品混合在不同的批次中，那么需要提供批次信息。

通过将数据集减少到仅高表达的基因（由nfeatures参数控制），可以大大加快分析速度，即使会稍微影响到准确度。然后，根据cluster参数，将执行最终的PCA和聚类（使用内部 fastcluster函数）。基于cluster方法的基本原理是同型双细胞几乎不可能根据它们的转录组进行区分，因此创建这种双细胞是一种计算资源的浪费，而且还会误导分类器标记为单细胞。

然而，另一种方法是随机生成双细胞（将clusters设置为 FALSE 或 NULL），并使用迭代方法从训练中排除无法识别的人工双细胞。

根据cluster和propRandom参数，将通过合并随机细胞和/或不相同的cluster对的细胞合并，形成人工双细胞（这可以使用getArtificialDoublets函数手动执行）。 一部分双细胞将简单地使用组成细胞的counts总和，而其余的将进行调整文库大小和进行泊松重采样，数据校正。

对真实细胞和人工细胞的组合执行新的PCA，从中生成 kNN网络。 使用这个 kNN，为每个细胞收集了许多参数，例如 KNN中双细胞的比例、到最近双细胞和最近非双细胞的距离之比等。在输出中报告了一些带有“scDblFinder”前缀的功能，例如：

scDblFinder有相当多的参数来控制预处理、双细胞的生成、分类等（参见scDblFinder）。 我们仅对重要参数进行说明。

双细胞的期望检出率对邻域中人造双细胞的密度没有影响，但会影响分类器的分数，特别是分类临界值。是通过dbr和dbrsd参数指定（dbrsd指定dbr周围的 +/- 范围，在该范围内与dbr的偏差将被视为空）。

对于10x数据，捕获的细胞越多，产生双细胞的概率越大，Chromium文档表明每1000个细胞捕获的双细胞率大约为 1%（因此对于 5000 个细胞，(001 5) 5000 = 250 个双细胞） ，scDblFinder默认的预期双细胞率将设置为01（默认标准偏差为 0015）。但是请注意，不同的实验方案可能会产生更多的双细胞率，因此需要相应地更新。如果不确定双细胞率，您可能会考虑增加 dbrsd，以便大多数/纯粹从错误分类错误中估计它。

那么你很可能有错误的双细胞率。 如果你没有提供dbr参数，双细胞率将使用10X Genomics预期双细胞率自动计算，这意味着捕获的细胞越多，双细胞率就越高。 如果你认为不适用于你的数据，可手动设置dbr。

如果出现意外高的双细胞率最常见原因是，1）你有一个多样本数据集并且没有按样本进行拆分。 scDblFinder会认为数据是具有大量细胞的单次捕获，因此具有非常高的双细胞率。 按样本拆分应该可以解决问题。

阈值根据预期双细胞数量和错误分类（即人造双细胞）试图最小化方差，这意味着有效（即最终）双细胞率将与给定的不同。 scDblFinder还认为假阳性比假阴性对后续的分析问题要小些。你可以通过设置 [dbrsd=0]在一定程度上减少与输入双细胞率的偏差。

虽然这两种方法在基准测试中的表现非常相似，但随机生成方法在复杂数据集中通常略胜一筹。 如果你的数据被非常清晰地划分为cluster，或者你对双细胞的起源感兴趣，则基于cluster的方法更可取。 这也将能够更准确地计算同型双细胞率，因此略好于阈值法(thresholding)。 否则，特别是如果你的数据没有非常清楚地划分为cluster，则随机方法（例如clusters= FALSE）更可取。

如果你在多样本数据集上运行scDblFinder但是未提供cluster标签，而是基于特定样本的标签（意味着一个样本中的标签“1”可能与另一个样本中的标签“1”无关），并且 在 tSNE上绘制它们看起来没有意义。 出于这个原因，当运行多个样本时，建议首先将所有样本聚集在一起（例如使用 sce$cluster <- fastcluster(sce)），然后将cluster信息提供给 scDblFinder。

如果某些细胞的读数为零（或非常接近于零），则会出现‘Size factors should be positive’此错误。 过滤掉这些细胞后，错误应该消失了。 但是请注意，我们建议在运行 scDblFinder之前不要进行过于严格的过滤。

由于它依赖于人工双细胞的部分随机生成，因此对同一数据多次运行scDblFinder会产生略有不同的结果。你可以使用 setseed() 确保可重复性，但是在多线程时使用 setseed() 是不行的。请使用以下程序：

如果输入的sce对象已经包含归一化矩阵( logcounts)或名为“PCA”的reducedDim数据，scDblFinder将使用它们进行聚类分析。 此外，可以使用 scDblFinder() 函数的 cluster参数手动指定。 通过这种方式，seurat聚类可以例如用于创建人造双细胞（参见 Seurat::asSingleCellExperimentSeurat for conversion to SCE）。

在人造双细胞生成之后，真实和人造双细胞的计数必须一起重新处理（即归一化和 PCA），在内部使用scater执行的。 如果您希望以不同的方式执行此步骤，您可以提供自定义函数来执行此 *** 作（请参阅 scDblFinder 中的处理参数）。 然而，我们注意到，这一步的变化对双细胞检测的影响不大。 事实上，例如，根本不执行任何归一化会降低双峰识别的准确性，但也是一点点。

可以，专门处理峰值数据。由于单细胞ATAC-seq数据的稀疏性比转录组要大得多，而且scDblFinder需要处理一系列基因，因此使用默认的标准参数，运行的性能较差(执行慢)。因此，我们推荐使用aggregateFeatures=TRUE，这将在正常的scDblFinder 过程之前聚合相似的基因（而不是选择基因），会产生不错的结果。如果基因足够少，我们推荐直接基于距离计算而不是通过SVD步骤获得，如下所示：

cDblFinder的输入数据不应包含空液滴，并且可能需要移除覆盖率非常低的细胞以避免错误（例如 <200 reads）。

进一步的细胞质控应该在双细胞识别的下游进行，有两个理由：

1默认的预期双细胞率是根据给定的细胞计算的，如果你排除了很多质量低的细胞，scDblFinder可能会认为双细胞率应该低于实际值。

2剔除所有低质量细胞可能会妨碍我们检测由高质量细胞和低质量细胞组合形成的双细胞的能力。

话虽如此，这些主要是理论依据，除非你的QC过滤非常严格（而且不应该如此！），否则结果不太可能有很大的不同。

scDblFinder运行之前要做一些初次过滤，但不要太严格（例如 <200 UMI）

质控粗略过滤->运行scDblFinder->较严格过滤

后记：

之前在群里看到有人问，双细胞质控后，拿质控后的数据重新跑scDblFinder，还是会有大量双细胞被检出？

这个是必然的，由scDblFinder的算法决定，它是由输入数据建立起预测分类模型，不像singleR有另外一套参考数据集，拿输入数据去map到参考数据集。

只要你喂给scDblFinder数据，它都会给输入的细胞一个score值，然后设置阈值进行分类，值高的为双细胞。

scDblFinder强烈依赖输入数据，有它使用的范畴，所以反复进行双细胞scDblFinder质控，用法是不对的。

拿到单细胞转录组数据，先质控粗略过滤->运行scDblFinder进行双细胞质控->较严格质控过滤。

其实最好有一个相互验证的过程。

以上就是关于10X单细胞 & 10XATAC 联合分析表征细胞调控网络（MIRA）全部的内容，包括:10X单细胞 & 10XATAC 联合分析表征细胞调控网络（MIRA）、如何在实践中学习Python、KNN算法常见问题总结等相关内容解答，如果想了解更多相关内容，可以关注我们，你们的支持是我们更新的动力！