单细胞去除聚类的离群点

Original 生信技能树单细胞天地 2022-06-07

最近收到学员提问，关于单细胞转录组聚类分群后可视化发现有极个别细胞是离群点，如下，想去除掉这几个少数分子，其实我指点了应该是提取坐标即可，本质上仍然是对seurat包的熟练程度罢了。

跟我们前面回答的问题类似，不过那些问题是针对于monocle包，使用monocle做拟时序分析（单细胞谱系发育）学员的问题是: 拟时序分析的热图提取基因问题 , 本质上，都对R包返回对象的了解程度罢了。

为了解释如何提取坐标，我需要使用大家都理解的数据集，然后创造出上面的聚类图。

rm(list=ls())
options(stringsAsFactors = F)

library(scRNAseq)
## ----- Load Example Data -----
data(fluidigm)
# Set assay to RSEM estimated counts
assay(fluidigm) <- assays(fluidigm)$rsem_counts
ct <- floor(assays(fluidigm)$rsem_counts)
ct[1:4,1:4]
meta <- as.data.frame(colData(fluidigm))
counts <- ct   
# 下面的代码仅仅是对 Seurat(V3) 有效
# Create Seurat object
library(Seurat)
sce_test <- CreateSeuratObject(counts = ct, 
                             meta.data = meta,
                             min.cells = 1, 
                             min.features = 0, 
                             project = 'c1_xy') # already normalized
sce_test  

# Add meta.data (nUMI and timePoints)
sce_test <- AddMetaData(object = sce_test, 
                      metadata = apply(ct, 2, sum), 
                      col.name = 'nUMI_raw')

# 步骤 ScaleData 的耗时取决于电脑系统配置（保守估计大于一分钟）
sce_test <- ScaleData(object = sce_test, 
                    vars.to.regress = c('nUMI_raw'), 
                    model.use = 'linear', 
                    use.umi = FALSE)
sce_test <- FindVariableFeatures(object = sce_test, 
                               mean.function = ExpMean, 
                               dispersion.function = LogVMR, 
                               x.low.cutoff = 0.0125, 
                               x.high.cutoff = 4, 
                               y.cutoff = 0.5)
length(VariableFeatures(sce_test))
sce_test <- RunPCA(object = sce_test, pc.genes = VariableFeatures(sce_test))
# 下面只是展现不同降维算法而已，并不要求都使用一次。
sce_test <- RunICA(sce_test )
sce_test <- RunTSNE(sce_test )
#sce_test <- RunUMAP(sce_test,dims = 1:10)
#VizPCA( sce_test, pcs.use = 1:2)
DimPlot(object = sce_test, reduction = "pca") 
DimPlot(object = sce_test, reduction = "ica")
DimPlot(object = sce_test, reduction = "tsne")
#DimPlot(object = sce_test, reduction = "umap")

# 针对PCA降维后的表达矩阵进行聚类 FindNeighbors+FindClusters 两个步骤。
sce_test <- FindNeighbors(object = sce_test, dims = 1:20, verbose = FALSE) 
sce_test <- FindClusters(object = sce_test, verbose = FALSE)
# 继续tSNE可视化
DimPlot(object = sce_test, reduction = "tsne")

出图如下：

可以看到是很明显的2个细胞亚群，但是有少数几个细胞，走错了地方，这个时候，学员突发奇想要删掉它，我这里不想评价这样做对不对，先给出解决方案吧。

其实一行代码就明白了：

> head(sce_test@reductions$tsne@cell.embeddings)
               tSNE_1    tSNE_2
SRR1275356  -7.326137 -2.427610
SRR1274090   7.266558  2.613246
SRR1275251 -11.391271 -3.700327
SRR1275287  -5.283922 -2.496858
SRR1275364  -4.570582 -1.116134
SRR1275269  -4.215317 -5.497514
plot(sce_test@reductions$tsne@cell.embeddings)