在UMAP图上标识我们感兴趣的基因所在的类群(单细胞数据)

参考链接: https://www.jianshu.com/p/37d2e8d68c91
https://satijalab.org/seurat/articles/pbmc3k_tutorial.html
https://satijalab.org/seurat/articles/visualization_vignette.html
我们这个代码要解决的需求,就是将我们从GEO数据库中下载的表达矩阵(.csv文件)使用seurat这个包进行处理。期望的目标是绘制出UMAP图,将我们感兴趣的基因标记在上面。
(1)目前,我对于seurat这个包的认识几乎为零。相当于是从头开始学。
(2)另一方面,我们所使用的这个处理矩阵的处理方式是相对比较模糊的。
这是我遇到的两个难题,我相信我每次都具有化险为夷的能力,我相信自己可以接下去克服难关。我希望我的分析结果,足够可靠,能够给他们接下来的实验提供指导。

其实,我们没必要使用seurat这个封装好的包,虽然每一步的处理都是傻瓜式的,用起来很爽。但是这样处理得到的结果的可信度我是需要打一个问号的。
我们就单纯的搜索一下,直接使用矩阵画UMAP图可不可以?
昨晚,我们使用umap直接处理了原始的count矩阵,结果计算了一整晚都没有计算完成(这样的方法肯定是有错误的)。后来,查看官网的pipeline,发现在计算主成分之前,是先利用高变基因来寻找的。所以,现在大致的思路应该是按照官网的流程,重新来走一遍。

1。首先第一步,将我们的.csv数据转换为seurat的对象。

参考链接:https://www.jianshu.com/p/41d7fdae0484 #这个作者的代码是有问题的,不用她的
https://www.jianshu.com/p/37d2e8d68c91

data<-read.csv("process.csv",header=T)
dim(data)
#36920   454
#首先将数据转换为secrut对象,我们这个矩阵就是count矩阵
#将数据转换为稀疏矩阵
dataan <- as(as.matrix(data),"dgCMatrix")

出错显示:

Error in as(as.matrix(data), “dgCMatrix”) :
没有可以用来強制转换“matrix”成“dgCMatrix”的方法或默认函数

但是,昨天好像是成功了,这是为什么呢?
找到了主要的原因是,缺少依赖的包。

library(dplyr)
library(Seurat)
library(patchwork)

加载完成之后,继续运行这一步就没有出错。
但是又warning的出现:

Warning message:
In storage.mode(from) <- “double” : NAs introduced by coercion
#翻译一下,就是说*引入了NA的值。

我们来看一下,是否产生了NA的值。
找到了主要的原因,是我的矩阵有一列是字符,也是我在读入数据的时候,没有将这一列(下图中的X列,也即基因名)设为行名。
在UMAP图上标识我们感兴趣的基因所在的类群(单细胞数据)
我重新读入矩阵试试看。

library(dplyr)
library(Seurat)
library(patchwork)
data<-read.csv("process.csv",header=T)
row.names(data)<-data[,1]
data<-data[,-1]
dim(data)
#36920   453
#首先将数据转换为secrut对象,我们这个矩阵就是count矩阵
#将数据转换为稀疏矩阵
data_1<- as.matrix(data)
data_2 <- as(as.matrix(data),"dgCMatrix") #主要原因是没有加载包。

重新读入矩阵之后,又出现了新的错误。

Error in asMethod(object) :
‘Realloc’ could not re-allocate memory (132335648 bytes)

以关键词“ ‘Realloc’ could not re-allocate memory”,在搜索引擎上进行检索。
翻译一下,就是无法重新分配内存(每次遇到这种类型的问题,我就立马傻掉)。
查了一下,是需要清除R的运行内存,可能我前前后后有处理很多很大的矩阵,所以把内存占满了吧。
参考链接:https://www.zhihu.com/question/390053502
解决方式:

rm(list=ls())
gc()
           used  (Mb) gc trigger   (Mb)  max used
Ncells  3632349 194.0    6540962  349.4   6540962
Vcells 40283137 307.4  139136965 1061.6 173921129
         (Mb)
Ncells  349.4
Vcells 1327.0

rstudioapi::restartSession() #重启R

Restarting R session...

最后解决。
我们看一下这个稀疏矩阵的结构。
在UMAP图上标识我们感兴趣的基因所在的类群(单细胞数据)
这个矩阵应该就类似于我比较熟悉的表达矩阵。
然后,将数据转换为seurat对象。

data_3<- CreateSeuratObject(counts = data_2, project = "tang", min.cells = 3, min.features = 200)

在转换的时候,又出现了新的问题。

Error in h(simpleError(msg, call)) :
在为’colSums’函数选择方法时评估’x’参数出了错: Cholmod error ‘out of memory’ at file …/Core/cholmod_memory.c, line 146

还是说out of memory,就很奇怪。

这个问题的解决很简单,是因为我前面处理,把空间给占满了,几个超级大的matrix,不满才怪。那这个时候呢,最简单的办法就是,重启Rstudio。

这个问题解决之后,我们的seurat对象就已经创建完成了。那么接下来,我们要面对的事情就是如何处理这个对象。

2。处理seurat对象。

(忘记记录了,接下来的步骤是seurat包封装好的,只需要傻瓜式的运行即可)

data_3<-FindVariableFeatures(data_3,selection.method = "vst",nfeatures = 5633) #寻找高变基因
all.genes<-rownames(data_3)
data_3<-ScaleData(data_3,features = all.genes) #归一化
data_3<-RunPCA(data_3,features = VariableFeatures(object = data_3)) #PCA
data_3<-FindNeighbors(data_3,dims = 1:40) #聚类
data_3<-FindClusters(data_3,resolution = 0.5)
data_3<-RunUMAP(data_3,dims = 1:40) #UMAP

3。单独标记我们感兴趣的细胞群

DimPlot(data_3,reduction = "umap")
#cell.highlight
#cols.highlight

g1_treat <- WhichCells(data_3, idents = c("1","2","3"))
DimPlot(data_3,label=T,cells.highlight= list(g1_treat), cols.highlight = c("pink"),cols = "grey")

g2_treat <- WhichCells(data_3, idents = c("5"))
DimPlot(data_3,label=T,cells.highlight= list(g2_treat), cols.highlight = c("pink"),cols = "grey")

g3_treat <- WhichCells(data_3, idents = c("5","0"))
DimPlot(data_3,label=T,cells.highlight= list(g3_treat), cols.highlight = c("pink"),cols = "grey")

g4_treat <- WhichCells(data_3, idents = c("0"))
DimPlot(data_3,label=T,cells.highlight= list(g4_treat), cols.highlight = c("pink"),cols = "grey")

得到的图如下所示。

在UMAP图上标识我们感兴趣的基因所在的类群(单细胞数据)
在UMAP图上标识我们感兴趣的基因所在的类群(单细胞数据)

4。标记我们感兴趣的基因表达的细胞

FeaturePlot(data_3,features = c("Gpr65"))
FeaturePlot(data_3,features = c("Vip"))
FeaturePlot(data_3,features = c("Glp1r"))
FeaturePlot(data_3,features = c("Oxtr"))

FeaturePlot(data_3,features = c("Gpr65","Vip","Glp1r","Oxtr"))

在UMAP图上标识我们感兴趣的基因所在的类群(单细胞数据)

当然,还有其他的细节可以继续的探究。到这里基本上完成了我们的需求。


但是,对于我们的项目而言,要明白的是,作者绘制的是tsne图。
在UMAP图上标识我们感兴趣的基因所在的类群(单细胞数据)
在tsne图中,作者将细胞分成了27类,而我们只有6类,不知道是否会有什么区别。
如果只是想看看,那没有什么关系。

data_3<-RunTSNE(data_3,dims = 1:40)
DimPlot(data_3,reduction = "tsne")

在UMAP图上标识我们感兴趣的基因所在的类群(单细胞数据)
文章中又说,有一部分的细胞是由于操作的过程中细胞损伤,要被移除的,我们在图中将其标出来。
(这里上传不了图片了)
大致的意思,这些损伤的细胞,不在我们感兴趣的基因所在的类。故而我们可以将其移除。或者这部分损伤的细胞不会对我们的实验产生影响。

上一篇:Opencv-GPU 编译错误 CUDA_nppicom_LIBRARY (ADVANCED) 未定义


下一篇:出现Head https://registry-1.docker.io/v2/library/node/manifests/14-alpine的解决方法