标签：seurat https 矩阵 library UMAP 图上 单细胞 data 我们
参考链接： 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列，也即基因名）设为行名。
我重新读入矩阵试试看。
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...
最后解决。
我们看一下这个稀疏矩阵的结构。
这个矩阵应该就类似于我比较熟悉的表达矩阵。
然后，将数据转换为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")
得到的图如下所示。
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"))
当然，还有其他的细节可以继续的探究。到这里基本上完成了我们的需求。
但是，对于我们的项目而言，要明白的是，作者绘制的是tsne图。
在tsne图中，作者将细胞分成了27类，而我们只有6类，不知道是否会有什么区别。
如果只是想看看，那没有什么关系。
data_3<-RunTSNE(data_3,dims = 1:40)
DimPlot(data_3,reduction = "tsne")
文章中又说，有一部分的细胞是由于操作的过程中细胞损伤，要被移除的，我们在图中将其标出来。
（这里上传不了图片了）
大致的意思，这些损伤的细胞，不在我们感兴趣的基因所在的类。故而我们可以将其移除。或者这部分损伤的细胞不会对我们的实验产生影响。
标签：seurat,https,矩阵,library,UMAP,图上,单细胞,data,我们 来源： https://blog.csdn.net/weixin_40640700/article/details/118944675
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。