关于什么是真正的科研精神

文摘   2024-11-03 15:51   中国澳门  

hello呀 各位bro and sis,我是不是好久没有分享错误了,还记得上一次我们debug是在6月份了,后面就开始做一些有的没的东西。这次我们把这cellchat继续更新完吧。

起初我在分析数据的时候并没有遇到很多错误,毕竟我的代码都是在这个包的作者那搞来的。

#倒库
library(dplyr)
library(Seurat)
library(patchwork)
library(harmony)
library(ggsci)
library(ggplot2)
library(stringr)
library(DDRTree)
library(pheatmap)
library(monocle)
library(CellChat)####细胞互作核心R包
load("../merge.relabel2.Rdata")

到这里应该都没啥问题,慢慢的加载数据,然后查看一下自己merge后有多少个亚群,然后把每个亚群都提取出来,进行交互分析。

table(merge@meta.data[["celltype"]])
###
Tcell <- subset(merge,celltype == "Tcell")
Epithelial <- subset(merge,celltype == "Epithelial")
Plasma <- subset(merge,celltype == "Plasma")
Macrophage <- subset(merge,celltype == "Macrophage")
Fibroblast <- subset(merge,celltype == "Fibroblast")
Endothelial <- subset(merge,celltype == "Endothelial")
NK_cell<- subset(merge,celltype == "NK_cell")
B_cell <- subset(merge,celltype == "B_cell")
Pericyte <- subset(merge,celltype == "Pericyte")
Mast <- subset(merge,celltype == "Mast")
Pit <- subset(merge,celltype == "Pit")
####合并###
cellchat.use <- merge(x=Tcell,y=c(Epithelial,Plasma,Macrophage,Fibroblast,Endothelial,NK_cell,B_cell,Pericyte,Mast,Pit))
cellchat.use <- JoinLayers(cellchat.use)
cellchat.use<-NormalizeData(cellchat.use)

labels <- Idents(cellchat.use)
#查看一下数据状态
table(cellchat.use@meta.data[["celltype"]])
table(cellchat.use@meta.data[["group"]])
######载入cellchat对象#######
###导入分组数据###
data<-subset(cellchat.use,group=="Normal")###根据组别要修改!!!
table(data$group)
#创建CellChat专用文件
#注意CellChat需要用normalize之后的数据
data.input <- GetAssayData(data[["RNA"]],layer = "data")
#此处注意,CellChat分析时取出的是Normalize之后的数据
#分组信息
data@meta.data$labels <- data$celltype
table(data@meta.data$labels)
meta <- data@meta.data
cellchat <- createCellChat(object = data.input,
                           meta = meta,
                           group.by = "labels",
                           assay = "RNA")
# 添加meta.data信息
cellchat <- addMeta(cellchat, meta = meta)
## 设置默认的labels
cellchat <- setIdent(cellchat,ident.use = "labels")
### # number of cells in each cell group
groupSize <- as.numeric(table(cellchat@idents))
groupSize

基本上到这步就是人尽皆知的内容了,查看R包作者给的教程也可以看到这些内容。然后我们需要倒入cellchat的数据库,并且将这个数据做一系列跟seurat一样的标准化处理。

#####载入CellChatDB数据集#####
CellChatDB <- CellChatDB.mouse
CellChatDB <- CellChatDB.human
#展示一下CellChatDB的数据库结构
showDatabaseCategory(CellChatDB)
dplyr::glimpse(CellChatDB$interaction)
CellChatDB.use <- subsetDB(CellChatDB,search = "Secreted Signaling",key = "annotation"# use Secreted Signaling
CellChatDB.use <- subsetDB(CellChatDB,search = "Cell-Cell Contact")
CellChatDB.use<-CellChatDB
CellChatDB.use <- subsetDB(CellChatDB)
CellChatDB.use <- CellChatDB
#将数据库加入到cellchat文件当中
cellchat@DB <- CellChatDB.use
######基于CellChat的DB细胞通讯数据分析#######
#数据预处理,用于细胞间通讯分析
#这一步一定要有,不然会报错的
#如果不想筛选就用NULL
cellchat <- subsetData(cellchat)
#下面2步类似于Seurat里面的FindVariableFeatures
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)
#投射数据到PPI,模拟转录过程
cellchat <- projectData(cellchat,PPI.human)#这里需要注意一下是人还是小鼠的,我之前用小鼠的那个做了一下结果就失败了,因为数据集不匹配。
###
#这一步特别慢
cellchat <- computeCommunProb(cellchat, type = "triMean")
#删除含有细胞量<10的通路
cellchat <- filterCommunication(cellchat,min.cells = 10)
#重复一下上面的操作
#这一步和上一步不同,上一步是在配体水平分析
#这一步是在通路水平上分析,这一步你不做的话,netp里面没信息
cellchat <- computeCommunProbPathway(cellchat)
cellchat <- aggregateNet(cellchat)#汇总之前的分析
cellchat <- aggregateNet(cellchat)
###保存###
cellchat.normal<-cellchat
save(cellchat.normal,file = "cellchat.normal.Rdata")

到这里呢基本上normal组的数据就分析完了,同样的步骤再提取一次tumor进行一次分析就好了。然后根据内容要求把两种得到了的内容合并一下,然后继续分析,作图。

###载入数据###
list.files()
load("cellchat.normal.Rdata")
load("cellchat.Tumor.Rdata")
load("cellchat.Rdata")
cc.list=list(Normal=cellchat.normal,Tumor=cellchat.Tumor)
###第一张图####
pdf("1.pdf",width =5,height = 5)
par(mfrow=c(1,2),xpd=T)#构图1行2列
##所有细胞群总体观:通讯数量与强度对比
gg1<-compareInteractions(cellchat,show.legend = F,group = c(1,2),measure = "count")
gg2<-compareInteractions(cellchat,show.legend = F,group = c(1,2),measure = "weight")
gg1+gg2
dev.off()

得到这样的图,应该没啥问题。

par(mfrow=c(1,2),xpd=T)#构图1行2列
netVisual_diffInteraction(cellchat,weight.scale = T,measure = "count")
netVisual_diffInteraction(cellchat,weight.scale = T,measure = "weight")
##红色是case相对于control上调的,蓝色是下调的,两组数据相减
names(cc.list)
dev.off()

但是到这里发现报错了??

Error in i_set_edge_attr(x, attr(value, "name"), index = value, value = attr(value,  : 
   Length of new attribute value must be 1 or 111, the number of target edges, not 102

我查看了所有的内容,并且用我们之前的一个步骤排查了错误,

obj1_raw <- cellchat@net[["Normal"]][["count"]]
> obj2_raw <- cellchat@net[["Tumor"]][["count"]]
> shared_label = intersect(rownames(obj1_raw), rownames(obj2_raw))
> obj1 = obj1_raw[shared_label, shared_label, drop = FALSE]
> obj2 = obj2_raw[shared_label, shared_label, drop = FALSE]
> net.diff <- obj2 - obj1
> View(net.diff)

其实我们可以发现这个是没有任何问题的,但我也不问题出现在什么地方。于是我询问一下Ai,

针对你的问题,netVisual_diffInteraction 函数在设置边属性时出现了长度不匹配的错误。以下是精确的修改方案,以确保 cellchat@net$count 数据正确并能够被该函数使用:

解决方案步骤:
检查网络矩阵的一致性: 确保 obj1 和 obj2 是对称矩阵,并且具有相同的维度和行列名称。如果 net.diff 的边数量与 cellchat 对象所需的不匹配,就会导致错误。

确保边数匹配: 使用以下代码来检查 net.diff 中的非零边数,并与 netVisual_diffInteraction 期望的边数匹配:
`
length(which(net.diff != 0))  # 查看非零边的数量
`
确认这个数量与 netVisual_diffInteraction 函数中引用的边数一致。

创建合适的数据格式: 如果需要,你可以重新构建网络矩阵,以确保数据格式与 cellchat@net$count 插槽中的预期格式相符。可以通过以下方式处理:
`
# 使用邻接矩阵或数据框格式,确保网络中所有边的数量一致
net.diff[is.na(net.diff)] <- 0  # 处理可能的 NA 值
`
赋值到 cellchat@net$count 插槽: 确保赋值的网络矩阵格式正确,并使用相同的 rownames 和 colnames。
`
cellchat@net$count <- net.diff
`
检查并运行 netVisual_diffInteraction: 尝试再次运行可视化函数,并确认是否解决了问题。
`
netVisual_diffInteraction(cellchat, weight.scale = TRUE, measure = "count")
`

经过这个Ai的指导后我运行了一次,结果还是同样的报错,我甚至觉得可能是我在前面构建提取的时候出现了问题,是不是有的细胞很少,所以他们之间没有交流。

查看完了之后发现很正常,每个注释的细胞都很多。我实在是百撕不得骑姐(杰哥告诉我的词汇),我查看了文章,并且有勇气的给作者发了邮件,询问这个事情。

事实上作者很快就给了我回复,第二天的时候,但是我在没看到邮件之前,就在想是不是r包不兼容的问题,比如上次我们看到一个包的作者删除了一些函数,下面很多人再说国粹的那个内容,于是我尝试制定安装一下r包的版本。

install.packages("igraph",version="1.3.5",type = "binary")

但是不知道是网络的原因还是啥,一直安装不上,可是我们可以去CRAN上下载这个包呀,下载完本地安装就OK了。安装好了,我们就可以看到结果是正常的了。

我们可以正常得到结果了,一个是数量上一个是强度上。到这里基本上我们的问题就解决完了。最后把R包的作者给我的回复贴出来,非常感谢作者大大给的回复,金教授目前也在国内,武汉大学。相比于教授给的回复,这可比国内某些通讯作者好太多了,想跟他要点东西,从来都不会回你,这可能就是差距吧,科研是共同进步的,而不是闭门造车的。

最后给我的好大哥做个广告,他在帮助大家共同进步这一块真的是very nice的。各种群,只要你能想到,从细胞层面到动物层面,从生物信息学到机制层面,甚至还有找对象层面,你都可以找他,他都能帮你解决。至少我加了的就有好几个了。甚至这些群里面真的有教授再给你解答问题的。这才是科研人该有的素质,而且找他可以免费领抗体,嘿嘿嘿。扫码加关注,解决一切科研问题。

好了,今天的小文章就到这里吧,拜拜了。


生信夜班侠
一个不想搞医疗,想改行搞科研的儿科Dr