hi，各位道友，上次我们介绍了R包DoubletFinder用于去除多细胞 那么python是否也有类似的模块去除多细胞呢，答案是有的。这次我们就来使用一下python模块去除多细胞
Single-Cell Remover of Doublets
Python code for identifying doublets in single-cell RNA-seq data
给定一个原始的（未归一化的）UMI，以细胞为行，基因为列的矩阵counts_matrix计数，计算每个单元的多细胞得分。
scr.scrub_doublets（）从观察到的数据模拟双峰，并使用k最近邻分类器为每个转录组计算一个连续doublet_score（介于0和1之间）。 分数将自动设置为阈值以生成predicted_doublets，这是一个布尔数组，对于预测的doublets为True，否则为False。
最佳做法：
一、处理来自多个样本的数据时，请分别对每个样本运行Scrublet。 因为Scrublet旨在检测由两个细胞的随机共封装形成的多细胞捕获，所以它在多个样本的合并数据集上可能表现不佳（原因大家都懂的）。
二、检查doublet分数阈值是否合理（在理想情况下，如本例所示，将双峰模拟doublet分数直方图的两个峰分开），并在必要时进行手动调整。例子在本文的后面展示。
三、可视化二维嵌入中的多细胞预测（例如UMAP或t-SNE）。 预测的双峰应该大体上共定位（可能在多个群集中）。 如果不是，则可能需要调整doublet得分阈值，或更改预处理参数以更好地解析数据中存在的单元格状态。
接下来我们看一下如何使用
第一步，导入必要的模块
第二步：读入矩阵，要求如上述所讲,计算多细胞比率
这一步包括
Initialize Scrublet object
相关参数是：
expected_doublet_rate:预期多细胞的比率，通常为0.05-0.1。 结果对该参数不是特别敏感。
sim_doublet_ratio：相对于观察到的转录组数量，要模拟的双峰数量。 此值应该足够高，以使所有的doublet状态都能通过模拟doublet很好地表示。 设置得太高在计算上是耗时的。 默认值是2，尽管低至0.5的值会为已测试的数据集提供非常相似的结果。
n_neighbors：用于构造观察到的转录组和模拟多细胞的KNN分类器的邻居数。 通常，round（0.5 * sqrt（n_cells））的默认值效果很好。
运行默认pipeline，其中包括：
双重模拟
标准化，基因过滤，重新缩放，PCA
多细胞计算
多细胞得分阈值检测和双峰调用
绘制观察到的转录组和模拟多细胞的多细胞得分直方图
模拟的多细胞直方图通常是双峰的。左模式对应于由具有相似基因表达的两个细胞产生的“嵌入”多细胞。 右边的的模式对应于“新型”多细胞，其由具有不同基因表达的细胞产生。 Scrublet只能检测”新型“双峰，这一点和doubleFinder的R包一样。
要比较单细胞与多细胞，我们必须设置一个阈值多细胞得分，理想情况下，应在模拟的双峰直方图的两种模式之间设置最小值。 scrub_doublets（）尝试自动识别这一点，并且在本示例中做得很好。 但是，如果自动阈值检测效果不佳，则可以使用call_doublets（）函数调整阈值。 例如：
scrub.call_doublets(threshold=0.25)
接下来我们画一下这个多细胞分布的直方图：
获取二维嵌入以可视化结果 (Tsne同理)
将空间位置信息和转录组分析相结合，对于癌症、免疫、肿瘤免疫相互作用，组织微环境，神经和发育等领域，有着令人期待的应用前景。
单细胞的一切分析，加前缀Spatial 都是一个新的分析点，因此Scanpy 扩展后也可用于空间转录组数据分析。 https://scanpy-tutorials.readthedocs.io/en/latest/spatial/integration-scanorama.html
python 包Scanpy很多函数是借鉴了R包Seurat，所以这两种方法分析结果差异不大，大家可以对照Seurat分析，上面网址也提供了Seurat包处理单细胞空间转录分析过程。
和分析单细胞转录组数据一样，单细胞空间转录组主要包括了：质控（QC），标准化（Normalization），降维聚类（Dimensional reduction and clustering），Cluster marker genes， Spatially variable genes
为了和Seurat结果比较，我们使用了相同的一套数据集： https://support.10xgenomics.com/spatial-gene-expression/datasets/1.1.0/V1_Mouse_Brain_Sagittal_Anterior ， 新鲜的冷冻小鼠脑组织， 前牙矢状切面，可以参考前面讲述的ABA大脑图谱： https://www.jianshu.com/p/5d087fffeb35
导入相关包
读取数据
预处理
我们根据总counts和表达的genes对spots进行一些基本的过滤:
这儿和Seurat得到的QC小提琴图一样，只是形式不同。
标准化，HVG
PCA，聚类，可视化
每一簇marker genes
空间特异性genes
保存数据
我们可以发现与Seurat相比，分类结果还是有差异，不过大的区域识别两种方法都没有什么问题。
单细胞转录组数据分析在阐述多细胞生物发育与疾病进程方面已经开发了多种新的方法，如比较有名的轨迹推断(TI，trajectory inference)。但是，我们知道，各种轨迹推断方法只是一种利用表达量的排序手段而已，而且严重依赖先验的知识，如根节点的选择。有没有一种技术可以真正的在RNA转录的时候为转录的RNA打上时间的标签呢？
今天，就为大家介绍一种基于RNA代谢标记的单细胞分析方法：dynamo。我们发现当我们真的可以标记RNA新旧之后，其实基于数学上我们对时间的定义，可以开发出许多新的单细胞动力学模型。
单细胞RNA测序提供了整个转录组的快照，但掩盖了RNA生成动态。Qi Qiu等提出了一种单细胞代谢标记的新RNA标记测序(scNT-seq，基于微流控微孔板单细胞RNA标记技术为新格元商业化的DynaSCOPE技术)，一种用于大规模平行分析来自同一细胞的新转录和已存在mRNA的方法。该方法基于液滴微流体的方法能够在条形码珠上进行高通量化学转换，有效地用T-to-C替换标记新转录的mRNA。
Qi Qiu等在文章中的数据揭示了时间序列的转录因子活性和在单细胞水平上响应神经元激活的细胞状态轨迹。作者进一步确定了RNA生物发生和衰变的速率，以揭示在多能和罕见的全能双细胞胚胎(2C)干细胞状态之间逐步转换的RNA调控策略。最后，将RNA代谢标记与遗传扰动相结合，确定DNA甲基胞嘧啶双加氧酶作为进入like-2C细胞状态的表观遗传屏障。时间分辨率下的单细胞转录组分析，打开了关于细胞类型特异性RNA调节机制的新线索。
我们知道，RNA水平的动态变化是由RNA转录、加工和降解的相互作用调控的。因此，要了解转录组在不同功能细胞类型中调控的机制，需要对基因表达的时间动态进行细胞类型特异性测量。单细胞RNA测序(scRNA-Seq)的发展使人们对细胞类型和状态的异质性有了更全面的了解。然而，标准的scRNA-Seq方法捕获到的是新转录(“新”)和已存在(“旧”)RNA的混合物，而不能获得RNA的时间动态。
单细胞转录组测序，连同RNA速度(velocity )和代谢标记(metabolic labeling)方法，以前所未有的分辨率揭示细胞状态及其之间的转换关系。然而，充分利用这些数据需要开发出能够预测细胞命运并揭示调控机制的动态模型。在这里，python 分析工具dynamo是一个分析框架，整合固有的剪接与否和标记与否，以估计绝对的RNA速度，重建速度矢量场(vector fields )，以预测未来的细胞命运，并最终使用微分几何分析阐明潜在的调控网络。
在文章 Massively parallel and time-resolved RNA sequencing in single cells with scNT-seq 中，作者将dynamo应用不同的生物学过程，包括预测分化造血干细胞谱系的未来状态，从细胞周期进程中糖皮质激素反应的反卷积、驱动斑马鱼色素沉着的调控网络的表征、确定对SARS-CoV-2感染的可能耐药性途径。
因此，基于RNA代谢标记的工作代表了从定性的、隐式的分化概念到定量和预测理论的重要一步探索。新模型的提出也为之前定性的和描述性质的单细胞数据分析朝着定量和预测迈出一大步。
开发dynamo的初心
在相同的转录本群体中区分新的和旧的RNA，常用的方法依赖于RNA代谢标记，利用外源性核苷类似物，如4-硫嘌呤(4sU)和标记RNA的生化富集。尽管这些方法对RNA动力学调控有重要的见解，但它们需要充足的起始材料，并对富集足够的细胞信号提出了挑战。最近发展了几种方法将4sU化学转化为胞嘧啶类似物，产生尿嘧啶-胞嘧啶复合物，在逆转录后标记新转录的RNA。这些化学方法允许直接测量细胞RNA的时间信息，而不需要生化富集。最近的研究表明，通过将Smart-Seq/plate-based scRNA-Seq与其中一种化学方法(如用于RNA代谢测序的巯基(SH)链烷基化)结合，可以在单细胞水平上联合分析新旧转录组。然而，这些基于Smart-Seq/plate的方法存在一些局限性。首先，它们成本高、耗时长，很难对高度异质性的细胞群进行大规模分析。其次，这些方法缺乏独特的分子标识符(UMIs)，无法准确定量新的转录水平。
为了克服这些限制，Qi Qiu等开发了scNT-Seq，一种基于UMI的高通量scRNA-Seq方法，结合了代谢RNA标记、液滴微流体和化学诱导的s4U到胞嘧啶模拟物的编码，同时测量来自同一细胞的新和旧转录组。作者证明，scNT-Seq能够在单细胞水平上对细胞RNA动力学、基因调控网络(GRN)活性和细胞状态轨迹进行时间分辨分析，同时它极大地提高了可扩展性并降低了成本。
Dynamo是一个python库，目前提供了一个完整的解决方案来分析传统scRNA-seq或基于scRNA-seq的时间代谢标记的动力学分析。它渴望成为持续集成机器学习、系统生物学、信息论、随机物理等最激动人心的发展的领先工具，以建模、理解和解释各种尖端单细胞基因组学技术生成的数据集(dynamo 2/3的开发正在进行中)。我们希望这些模型、理解和解释不仅能促进你的研究，而且最终可能导致新的生物学发现。
可以在单细胞分辨率下揭示细胞类型特异性TF调节活性的时间动态。
转录因子(TF)的调控活性可以通过将顺式调控序列与单细胞基因表达联系起来，以单细胞分辨率进行量化。通过scNT-Seq联合分析新的和旧的转录组，可以并行分析由外部刺激诱导的动态调控和与细胞特性相关的稳定调控。通过将单细胞调控网络推断和聚类(SCENIC)应用于成对的单细胞新/旧转录组。作者鉴定出79个协同调节的TF调控，且至少在一种细胞类型中具有显著的顺式调节基序富集。
最近的研究表明，基因表达的时间导数被称为“RNA velocity”，可以通过在scRNA-Seq数据集中区分未剪接(内含子reads)和剪接(外显子reads)的mrna来估计，并用于告知单个细胞的转录状态如何随时间(以小时为尺度)变化。我们首先检验了RNA velocity分析是否能够预测单个细胞在短暂(分钟)和持续(小时)神经元激活时的转录状态轨迹。因为代谢标记可以捕捉RNA水平的快速变化，而通过3 '标记的UMIs检测新的RNA在很大程度上独立于基因结构，我们推断，从scNT-Seq中对新RNA和总RNA的单细胞配对测量可以用来计算按标记时间(单位时间的分子)缩放的基于代谢标记的RNA velocity。作者根据富集靶基因的新转录RNA水平计算这些活性调节转录因子在每个细胞中的调控活性。通过将这些转录因子的调控活性投射到RNA velocity 流上，我们构建了一个针对不同类别转录因子的单细胞分辨率、时间分辨率的调控活性图谱。
通过将TimeLapse化学与高通量液滴微流体平台相结合，scNT-Seq能够共同分析同一细胞的新合成和已存在的转录组，在单细胞水平捕获mRNA的时间信息。传统的RNA速度分析使用内源性RNA剪接动力学来告知细胞的未来轨迹。因此，它受到RNA剪接时间不受控制和许多基因内含子reads稀少性的限制。由于代谢标记周期的时间和长度可以通过实验控制，在scNT-Seq中通过3 '标记的UMIs直接计数新的和旧的转录本，提供了一种无偏性的方法来计算所有可检测基因的RNA动力学参数。
使用计算模型明确地结合了基于代谢标记的单细胞测序，我们可以计算高度动态过程(数分钟或者数小时间)的时间分辨RNA速度。此外，在外部刺激或细胞状态转变期间，测量与TF相连接的靶基因的新RNA水平可以在单细胞水平上刻画TF调节活性。最后，通过脉冲追踪实验，scNT-Seq可以更准确地估计RNA动力学参数，揭示RNA调控策略在罕见细胞群体。
https://jef.works/blog/2020/01/14/rna_velocity_analysis_tutorial_tips/
https://www.kallistobus.tools/velocity_tutorial
https://velocyto.org/velocyto.py/index.html
https://pypi.org/project/scvelo/
https://dynamo-release.readthedocs.io/en/latest/index.html
https://github.com/wulabupenn/scNT-seq/tree/master/TC_calling
https://zhuanlan.zhihu.com/p/359100974
https://www.global-sci.org/intro/article_detail/csiam-am/18653.html
https://www.biorxiv.org/content/10.1101/2020.09.19.304584v1.full
https://zhuanlan.zhihu.com/p/351104418
RNA velocity, methods
La Manno et al (2018, Nature) [ steady-state model ] Implemented in velocyto
Bergen et al (2019, bioRxiv) [ dynamical model ] Implemented in scVelo
Soneson et al (2020 bioRxiv) finds that “preprocessing choices affect RNA velocity”
Velociraptor from Rue / Lun / Soneson uses basilisk to embed scVelo in R/Bioc:
https://github.com/kevinrue/velociraptor Works great !
