【细胞通讯】CellCall
今天的分享是做细胞通讯的另外一种方法---cellcall，文章发表在Nucleic Acids Research杂志。
CellCall 是一个通过整合细胞内和细胞间信号来推断细胞间通讯网络和内部调节信号的工具包。
======数据集收集======
CellCall 基于 KEGG 通路收集配体-受体-转录因子 (L-R-TF) 轴数据集。和可能也是这个文章最大的亮点了。
（1）人类的L-R对是从NATMI, Cellinker, CellTalkDB, CellChat and STRING v11 databases (只提取有文献支持的验证的)。如果有配体-受体complex的也收集。
（2）从KEGG中收集L-R下游互作的TF。只有在同一个branch上面的才收集。最后收集到19 144 个L-R-TF互作的对。
（3）人类的TF-TG是从TRANSFAC, JASPAR, RegNetwork,Pathway Commons, TRRUST, Ontogenet, ReMap, EVEX, HTRIDB, CHEA, ENCODE and MOTIFMAP中收集的。最后收集到了587 248对有实验支持的TF-TG对。
（4）通过同源map，在小鼠中获得了小鼠中12 069 L-R-TF对，554 207对TF-TG对。
======CellCall测试======
安装：
library(devtools)
devtools::install_github("ShellyCoder/cellcall")
数据格式：
简单地说，cellcall需要的分析数据是一个行名为基因，列名为细胞类型的表达矩阵。列名不能包含逗号、句号、破折号等标点符号。建议使用下划线连接条形码和单元格类型。由于这个包的方法依赖于细胞类型信息，因此正确获取细胞类型信息很重要。
当然也可以直接输入单细胞的对象。这个包的作者提供了一个函数，用来直接将seurat对象转化为cellcall需要输入的数据格式，所以直接用函数即可。准备好一个定义好细胞群的单细胞seurat对象。
构建cellcall object：
test <- CreateObject_fromSeurat(Seurat.object=Seurat_single,#seurat对象
slot="counts",
cell_type="celltype", #细胞类型
data_source="UMI",
scale.factor =10^6,
Org = "Mus musculus") #物种信息
利用包里自带的数据集进行测试：
f.tmp <- system.file("extdata", "example_Data.Rdata", package="cellcall")
load(f.tmp)
## gene expression stored in the variable in.content
dim(in.content)
in.content[1:4, 1:4]
table(str_split(colnames(in.content), "_", simplify = T)[,2])
mt <- CreateNichConObject(data=in.content, min.feature = 3,
names.field = 2,
names.delim = "_",
source = "TPM",
scale.factor = 10^6,
Org = "Homo sapiens",
project = "Microenvironment")
接下来推断细胞-细胞通讯得分。参数我们都使用默认值，具体意义可以查看函数，根据需要修改！
mt <- TransCommuProfile(object = mt,
pValueCor = 0.05,
CorValue = 0.1,
topTargetCor=1,
p.adjust = 0.05,
use.type="median",
probs = 0.9,
method="weighted",
IS_core = TRUE,
Org = 'Homo sapiens')
CellCall 嵌入了一种通路活性分析方法，以帮助探索某些细胞之间通信所涉及的主要通路。这个整合在其他工具中不常见，所以做一下看看。
n <- mt@data$expr_l_r_log2_scale
pathway.hyper.list <- lapply(colnames(n), function(i){
print(i)
tmp <- getHyperPathway(data = n, object = mt, cella_cellb = i, Org="Homo sapiens")
return(tmp)
})
myPub.df <- getForBubble(pathway.hyper.list, cella_cellb=colnames(n))
plotBubble(myPub.df)
还是和别的工具一样，来一个圈图，展示所有细胞间的互作联系。
cell_color <- data.frame(color=c('#e31a1c','#1f78b4','#e78ac3','#ff7f00'), stringsAsFactors = FALSE)
rownames(cell_color) <- c("SSC", "SPGing", "SPGed", "ST")
ViewInterCircos(object = mt, font = 2, cellColor = cell_color,
lrColor = c("#F16B6F", "#84B1ED"),
arr.type = "big.arrow",arr.length = 0.04,
trackhight1 = 0.05, slot="expr_l_r_log2_scale",
linkcolor.from.sender = TRUE,
linkcolor = NULL, gap.degree = 2,
order.vector=c('ST', "SSC", "SPGing", "SPGed"),
trackhight2 = 0.032, track.margin2 = c(0.01,0.12), DIY = FALSE)
这里作者提出两种object画图，所以也可以根据另外一种object画图。
ViewInterCircos(object = mt@data$expr_l_r_log2_scale, font = 2,
cellColor = cell_color,
lrColor = c("#F16B6F", "#84B1ED"),
arr.type = "big.arrow",arr.length = 0.04,
trackhight1 = 0.05, slot="expr_l_r_log2_scale",
linkcolor.from.sender = TRUE,
linkcolor = NULL, gap.degree = 5,
order.vector=c('ST', "SSC", "SPGing", "SPGed"),
trackhight2 = 0.032, track.margin2 = c(0.05,0.12), DIY = T)
还可以采用热图来呈现不同细胞类型之间 L-R 相互作用的详细通信分数。
viewPheatmap(object = mt, slot="expr_l_r_log2_scale", show_rownames = T,
show_colnames = T,treeheight_row=0, treeheight_col=10,
cluster_rows = T,cluster_cols = F,fontsize = 12,angle_col = "45",
main="score")
也可以用桑基图展示：
mt <- LR2TF(object = mt, sender_cell="ST", recevier_cell="SSC",
slot="expr_l_r_log2_scale", org="Homo sapiens")
head(mt@reductions$sankey)
if(!require(networkD3)){
BiocManager::install("networkD3")
}
sank <- LRT.Dimplot(mt, fontSize = 8, nodeWidth = 30, height = NULL, width = 1200,
sinksRight=FALSE, DIY.color = FALSE)
networkD3::saveNetwork(sank, "~/ST-SSC_full.html")
其中：第一列是ligand，第二列是receptor，最后一列是TF。左边和右边的颜色是和ligand和TF一致的。
library(magrittr)
library(dplyr)
tmp <- mt@reductions$sankey
tmp1 <- dplyr::filter(tmp, weight1>0) ## filter triple relation with weight1 (LR score)
tmp.df <- trans2tripleScore(tmp1) ## transform weight1 and weight2 to one value (weight)
head(tmp.df)
## set the color of node in sankey graph
mycol.vector = c('#5d62b5','#29c3be','#f2726f','#62b58f','#bc95df', '#67cdf2', '#ffc533', '#5d62b5', '#29c3be')
elments.num <- tmp.df %>% unlist %>% unique %>% length()
mycol.vector.list <- rep(mycol.vector, times=ceiling(elments.num/length(mycol.vector)))
sankey_graph(df = tmp.df, axes=1:3, mycol = mycol.vector.list[1:elments.num], nudge_x = NULL, font.size = 4, boder.col="white", isGrandSon = F)
其中：左边和右边的颜色是和ligand和receptor一致的。
tmp <- mt@reductions$sankey
tmp1 <- dplyr::filter(tmp, weight1>0) ## filter triple relation with weight1 (LR score)
tmp.df <- trans2tripleScore(tmp1) ## transform weight1 and weight2 to one value (weight)
## set the color of node in sankey graph
mycol.vector = c('#9e0142','#d53e4f','#f46d43','#fdae61','#fee08b','#e6f598','#abdda4','#66c2a5','#3288bd','#5e4fa2')
elments.num <- length(unique(tmp.df$Ligand))
mycol.vector.list <- rep(mycol.vector, times=ceiling(elments.num/length(mycol.vector)))
sankey_graph(df = tmp.df, axes=1:3, mycol = mycol.vector.list[1:elments.num],
isGrandSon = TRUE, nudge_x = nudge_x, font.size = 2, boder.col="white",
set_alpha = 0.8)
其中：左边和右边的颜色都是和ligand一致的。
采用转录因子富集图来呈现受体细胞中的转录因子活性，这个很好，其他工具很少有整转录因子的，而且转录因子的分析费时费力。
str(mt@data$gsea.list$SSC@geneSets)
ssc.tf <- names(mt@data$gsea.list$SSC@geneSets)
ssc.tf
getGSEAplot(gsea.list=mt@data$gsea.list, geneSetID=c("MYB", "TCF7", "FOXO3"),
myCelltype="SSC", fc.list=mt@data$fc.list,
selectedGeneID = mt@data$gsea.list$SSC@geneSets$CREBBP[1:10],
mycol = NULL)
也可以山脊图展示变化倍数
## gsea object
egmt <- mt@data$gsea.list$SSC
## filter TF
egmt.df <- data.frame(egmt)
head(egmt.df[,1:6])
flag.index <- which(egmt.df$p.adjust < 0.05)
ridgeplot.DIY(x=egmt, fill="p.adjust", showCategory=flag.index, core_enrichment = T,orderBy = "NES", decreasing = FALSE)
本文使用 文章同步助手 同步
