背景简介：
环境游离RNA污染是单细胞测序中可能存在的情况，他对细胞测序质量的影响较大，因此，有效地计算和预测游离RNA污染，去除污染严重的低质量细胞对单细胞测序分析具有重要意义。在生信技能树中已经有关于相关的软件包soupX介绍。详见：
这篇笔记我们介绍另一个R工具decontX。笔者认为，他的主要优点是使用方便，步骤简便，且结果容易使用。
原文链接：Decontamination of ambient RNA in single-cell RNA-seq with DecontX
原文摘要：基于液滴的微流控装置（主要是10X Genomics）已被广泛用于单细胞 RNA 测序(scRNA-seq)。然而，存在于细胞悬浮液中的环境 RNA 可以与细胞的天然 mRNA 一起被异常计数，并导致不同细胞群之间转录物的交叉污染。DecontX 是一种新的贝叶斯方法，用于估计和去除单个细胞中的污染。DecontX 准确地预测了鼠-人混合数据集中的污染水平，并去除了 PBMC 数据集中标记基因的异常表达。我们还比较了四种不同的 scRNA-seq 方案的污染水平。总的来说，DecontX 可以整合到 scRNA-seq 工作流中，以改进下游分析。
心路历程：
笔者以前并没有对游离RNA污染这个事非常在意，而包括Seurat/Scanpy等知名分析流程中，也没有明确其重要性，生信技能树关于soupX的推文虽然看了，但是由于步骤看起来比较繁琐，因此一直没有尝试，直到我最近尝试整合多个样本数据时，发现一个难以解决的问题：为什么我的UMAP可视化时候总是出现毛毛躁躁的边缘和大量散在细胞，还有很多细胞亚群之间有连续的细胞（maybe可能时真是存在的过渡态细胞）。
更让我不解的是，为什么顶刊大佬们的图就这么圆滑、这么清晰？比如：
【Naftali Kaminski, et al. Circulation. 2021】
HLCA_UMAP【而我的图】可以看到有很多散在细胞，边缘不干净，细胞出现亚群之间界限不明确，小亚群特别多
my_umap_before_decontX这时候我做了两个事情：
- 提高质控的严格程度：由于我是整合数据，细胞比较多，于是采取了更加严格的质控条件。要求每个细胞检测的基因数量 > 1000，检测到的 UMI > 1000（之前是基因数量>500）。这一操作使得我的细胞数量从14.8万直接下降到8.8万（这其中很重要的原因是某个数据集的测序质量很差）
- 评估环境游离RNA：这实际上只是一个尝试，我并不抱有很大希望。我使用了
decontX预测了计算RNA污染程度（contamination）。decontX官方并没有推荐筛选的阈值（maybe是我没看到），我参考了生信技能树 soupX 推文中提到的，soupX 官方推荐污染程度控制在20%以下能够去除99%的污染。
经过上面两个步骤后，我分别重新跑了scVI和harmony，umap的结果就变成这样了：【scVI】
scvi_umap_after_decontX【Harmony】
harmony_umap_after_decontXWow ! 看起来非常不错~
这里我需要解释一下为什么scVI看起来仍然有批次效应（source这个图），我问了生物信息学的老师，他的解释是：scVI不会强行拟合，他会更好地保留生物学差异。后来我查了数据集对应的文献，并结合我的图，好像确实是这样。
因此，未来我将会把这两个步骤纳入我的常规分析流程中。
代码实操：
讲了这么多，我们现在开始实操吧：
1. 下载decontX包
# install.packages("devtools")
devtools::install_github("campbio/decontX")
2. decontX评估
decontX接收两类数据：SingleCellExperiment对象， 或counts矩阵。
在这里我推荐后者，因为counts我们可以直接从seurat对象中获取，不必进行对象转换，也减轻了我们运行内存的负担，提高运行速度。
library(Seurat)
library(decontX)
# 如果你已经有一个处理好的seurat对象，但不想转化数据类型，可以直接使用这个方法
# 读取seurat对象
scobj <- readRDS('test.rds')
# 提取矩阵(genes x cells)
# 需要注意的是，你应该查看你的矩阵是否已经注释了行名和列明，以及如果你的矩阵是scanpy对象中提取的矩阵应该行列转置
counts <- scobj@assays$RNA@counts
# 计算
library(decontX)
# 需要把结果储存在新变量
decontX_results <- decontX(counts)
没错，就是这么简单就搞定了。
接下来我们提取结果，需要注意的是，decontX还会提供一个矫正后的矩阵，但我通常不会使用，有兴趣的话你可以提出来探索一下是否好用。
# 你可以使用str()查看结果的形式
# RNA污染的计算结果储存在：decontX_results$contamination
# 他是一个数值型向量，长度与细胞数量一致，顺序也与矩阵的colnames一致
# 我们可以直接把他写入metadata
scobj$Contamination =decontX_results$contamination
head(scobj@meta.data) # 可以查看一下结果，这里就不展示了
# 我们对他进行可视化
library(ggplot2)
FeaturePlot(scobj,
features = 'Contamination',
raster=FALSE # 细胞过多时候需要加这个参数
) +
scale_color_viridis_c()+
theme_bw()+
theme(panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank())+
xlab('scVI_UMAP_1')+
ylab('scVI_UMAP_2')
# 在这里看到，对比前文最开始的状态，我们经过了严格的质控之后，散在的细胞已经少了很多，但是仍然有亚群之间不清晰的界限
# 我们可视化Contamination，惊奇地发现这些边缘的毛躁就是Contamination较高地细胞
contamination现在我们根据Contamination进行筛选，decontX官方并没有给予明确的建议，这里参考的了soupX的官方文档，并基于自己摸索，根据这个数据的情况，选择 Contamination < 0.2是比较好的选择。
当然，你需要结合你的数据情况去评估！
low_con_scobj = scobj[,scobj$Contamination < 0.2] #保留小于等于0.2的细胞
# 这里大概删除了1700+细胞
删除之后，如开始所说的，我分别重新进行scVI和Harmony整合，并进行了umap降维，这些流程都是官方文档，没有什么特别的。效果如开头所见，这里只展示了scVI结果。
scvi_umap_after_decontX那么本次分享就到这里了。此外，本次画图除了decontX结果的可视化，主要用到的工具为omicverse（python工具包），开发作者上个月在技能树中也分享了一个超详细的系列教程，非常推荐大家阅读和使用，至于效果大家也可以看到了。
最后，感谢健明老师的支持，笔者非生信专业，如果错误之处，欢迎各位朋友提出建议和批评。
补充
如果你对SingleCellExperiment对象如果做感兴趣，我这里也贴了官方给的代码作为参考
# 如果你只有seurat对象，我推荐你用counts方法
# 如果你非要用sce对象，可以自己查一下如何转换
# 我们以10X的pbmc4k为例
# BiocManager::install("celda")
# BiocManager::install("TENxPBMCData")
# Load PBMC data
library(TENxPBMCData)
pbmc4k <- TENxPBMCData("pbmc4k") # 首次下载可能会比较慢
colnames(pbmc4k) <- paste(pbmc4k$Sample, pbmc4k$Barcode, sep = "_")
rownames(pbmc4k) <- rowData(pbmc4k)$Symbol_TENx
library(decontX)
# 这里与counts不同，可以直接储存在sce对象
pbmc4k <- decontX(x = pbmc4k) # 一行代码即可
# --------------------------------------------------
# Starting DecontX
# --------------------------------------------------
# Wed Aug 30 17:14:28 2023 .. Analyzing all cells
# Wed Aug 30 17:14:28 2023 .... Converting to sparse matrix
# Wed Aug 30 17:14:33 2023 .... Generating UMAP and estimating cell types
# Wed Aug 30 17:14:47 2023 .... Estimating contamination
# Wed Aug 30 17:14:48 2023 ...... Completed iteration: 10 | converge: 0.02716
# Wed Aug 30 17:14:49 2023 ...... Completed iteration: 20 | converge: 0.01045
# Wed Aug 30 17:14:49 2023 ...... Completed iteration: 30 | converge: 0.005854
# Wed Aug 30 17:14:50 2023 ...... Completed iteration: 40 | converge: 0.00389
# Wed Aug 30 17:14:51 2023 ...... Completed iteration: 50 | converge: 0.002784
# Wed Aug 30 17:14:52 2023 ...... Completed iteration: 60 | converge: 0.002249
# Wed Aug 30 17:14:52 2023 ...... Completed iteration: 70 | converge: 0.001869
# Wed Aug 30 17:14:53 2023 ...... Completed iteration: 80 | converge: 0.001577
# Wed Aug 30 17:14:54 2023 ...... Completed iteration: 90 | converge: 0.001344
# Wed Aug 30 17:14:55 2023 ...... Completed iteration: 100 | converge: 0.001125
# Wed Aug 30 17:14:55 2023 ...... Completed iteration: 109 | converge: 0.0009944
# Wed Aug 30 17:14:55 2023 .. Calculating final decontaminated matrix
# Wed Aug 30 17:14:56 2023 .. Converting decontaminated matrix to DelayedMatrix
# --------------------------------------------------
# Completed DecontX. Total time: 56.66165 secs
# --------------------------------------------------
# 查看数据
pbmc4k
# class: SingleCellExperiment
# dim: 33694 4340
# metadata(1): decontX
# assays(2): counts decontXcounts
# rownames(33694): RP11-34P13.3 FAM138A ... AC213203.1 FAM231B
# rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
# colnames(4340): pbmc4k_AAACCTGAGAAGGCCT-1
# pbmc4k_AAACCTGAGACAGACC-1 ... pbmc4k_TTTGTCAGTTAAGACA-1
# pbmc4k_TTTGTCATCCCAAGAT-1
# colData names(13): Sample Barcode ... decontX_contamination
# decontX_clusters # 结果添加在colData
# reducedDimNames(1): decontX_UMAP # 添加了UMAP
# mainExpName: NULL
# altExpNames(0):
# 现在我们把计算结果提取出来看看
# decontX_contamination为计算的游离RNA污染率(contamination rate)
coldata = colData(pbmc4k)
head(coldata)
# 内容很多，只列出来有关的两列
# decontX_contamination decontX_clusters
# <numeric> <factor>
# pbmc4k_AAACCTGAGAAGGCCT-1 0.0646881 1
# pbmc4k_AAACCTGAGACAGACC-1 0.1464950 1
# pbmc4k_AAACCTGAGATAGTCA-1 0.0141301 1
# pbmc4k_AAACCTGAGCGCCTCA-1 0.0802633 2
# pbmc4k_AAACCTGAGGCATGGT-1 0.0237364 2
# pbmc4k_AAACCTGCAAGGTTCT-1 0.0351587 2
# 最后，如果你需要把他写入seurat对象，可以提取出来之后直接添加在metadata中
# 需要注意的是，你的seurat对象的细胞顺序需要与这里一致
