seurat标识我们感兴趣的基因所在的类群（单细胞数据）

参考链接： 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”)
1
2
3
4
5
6
出错显示：

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

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

library(dplyr)
library(Seurat)
library(patchwork)
1
2
3
加载完成之后，继续运行这一步就没有出错。
但是又warning的出现：

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

我们来看一下，是否产生了NA的值。
找到了主要的原因，是我的矩阵有一列是字符，也是我在读入数据的时候，没有将这一列（下图中的X列，也即基因名）设为行名。

我重新读入矩阵试试看。

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”) #主要原因是没有加载包。
1
2
3
4
5
6
7
8
9
10
11
12
重新读入矩阵之后，又出现了新的错误。

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…
1
2
3
4
5
6
7
8
9
10
11
12
最后解决。
我们看一下这个稀疏矩阵的结构。

这个矩阵应该就类似于我比较熟悉的表达矩阵。
然后，将数据转换为seurat对象。

data_3<- CreateSeuratObject(counts = data_2, project = “tang”, min.cells = 3, min.features = 200)
1
在转换的时候，又出现了新的问题。

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
1
2
3
4
5
6
7
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”)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
得到的图如下所示。

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”))
1
2
3
4
5
6

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

但是，对于我们的项目而言，要明白的是，作者绘制的是tsne图。

在tsne图中，作者将细胞分成了27类，而我们只有6类，不知道是否会有什么区别。
如果只是想看看，那没有什么关系。

data_3<-RunTSNE(data_3,dims = 1:40)
DimPlot(data_3,reduction = “tsne”)
1
2

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

另外，经过偶然的尝试，我换用了作者文章中的指标，resolution=3，用36个维度来划分，结果真的划分成了12个类群。我终于弄明白这里面的主要的机制了，非常的开心。

接下来，我应该知道如何更加精确的来复现这张图了。
先完成这个结果，然后看怎样更加的完美。
我们现在先调整受损的细胞的数目，主要考虑的方式是去掉，高度表达与受损细胞相关的基因sprr1a,Ecel1的细胞。

现在把原先的模式关掉，重新退出。

调整了受损的细胞，但是不尽如人意。
————————————————
版权声明：本文为CSDN博主「今天也是个妖精头子呀」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/weixin_40640700/article/details/118944675