单细胞去除聚类的离群点

最近收到学员提问，关于单细胞转录组聚类分群后可视化发现有极个别细胞是离群点，如下，想去除掉这几个少数分子，其实我指点了应该是提取坐标即可，本质上仍然是对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)

就是根据坐标，拿到细胞的ID，然后从对象里面剔除掉这些细胞即可。

是不是非常简单啊！

其实我们两年前就在单细胞天地发布的全网第一个单细胞转录组课程，精炼了常规单细胞转录组数据分析主线，就是5大R包， scater,monocle,Seurat,scran,M3Drop，然后10个步骤：

  • step1: 创建对象

  • step2: 质量控制

  • step3: 表达量的标准化和归一化

  • step4: 去除干扰因素(多个样本整合)

  • step5: 判断重要的基因

  • step6: 多种降维算法

  • step7: 可视化降维结果

  • step8: 多种聚类算法

  • step9: 聚类后找每个细胞亚群的标志基因进行亚群命名

  • step10: 继续分类

如果你真的认真学了 ，这样的问题轻而易举就可以回答，可惜的是，知道这个课程的人寥寥无几，或者购买了也很少有人静下心来花费十几个小时听我絮叨。

“
赞 (0)

相关推荐