【细胞通讯】CellPhone

今天分享一篇关于分析细胞通讯的，发表于2020年nature protocols的相关文章。

CellPhoneDB的配体-受体数据库集成于UniProt、Ensembl、PDB、IUPHAR等，共存储978种蛋白质：其中501种是分泌蛋白，而585种是膜蛋白，这些蛋白质参与1,396种相互作用。CellPhoneDB数据库还考虑了配体和受体的亚基结构，准确地表示了异聚体复合物，有466个是异聚体，对于准确研究多亚基复合物介导的细胞通讯很关键。CellPhoneDB有474种相互作用涉及分泌蛋白，而490种相互作用涉及膜蛋白，总共有250个涉及整合素的相互作用。相比于早期的方法，改方法的优点就是考虑到多聚体受体的情况，而非简单地把受体限定为单个分子。

当然目前只支持人的数据。

依据构建的配体-受体库，对于用户输入的单细胞数据，基于一种细胞类型的受体和另一种细胞类型的配体的表达，得出两种细胞群之间丰富的受体-配体相互作用。对于细胞群所表达的基因，计算表达该基因的细胞百分比和基因表达平均值。若该基因只在该群中10%及以下的细胞中表达（默认值为0.1），则被移除。

将所有细胞的簇标签随机排列1000次（默认值），并确定簇中受体平均表达水平和相互作用簇中配体平均表达水平的平均值。对于两种细胞类型之间每对比较中的每一个受体-配体对，这产生了一个零分布（null distribution，亦称伯努利分布、两点分布，指的是对于随机变量X有, 参数为p(0<p<1)，如果它分别以概率p和1-p取1和0为值。EX= p,DX=p(1-p)。伯努利试验成功的次数服从伯努利分布,参数p是试验成功的概率）。

通过计算等于或高于实际平均值的平均值的比例，可以得到了一个P值，表示给定受体-配体复合物细胞类型特异性的可能性。换句话说，如果观察到的平均值在前5%，则相互作用的P值为0.05。

根据两种细胞类型中富集到的显著的受体-配体对的数量，对细胞类型之间高度特异性的相互作用进行排序，以便手动筛选出生物学相关的相互作用关系。

===CellPhoneDB的使用测试===

安装：

官网地址：https://github.com/Teichlab/cellphonedb

下面是官网的推荐安装方式：

1. Create python=>3.6 environment

Using conda: conda create -n cpdb python=3.7
Using virtualenv: python -m venv cpdb

Activate environment

Using conda: source activate cpdb
Using virtualenv: source cpdb/bin/activate

Install CellPhoneDB pip install cellphonedb

我是用的conda安装的：

conda create -n cpdb python=3.7

source activate cpdb

pip install cellphonedb

因为后面的画图是基于R的，所以还需要安装R以及需要的几个画图相关的包：ggplot2，pheatmap。

测试数据：

用的也是官网的测试数据：

https://raw.githubusercontent.com/Teichlab/cellphonedb/master/in/example_data/test_counts.txt

https://raw.githubusercontent.com/Teichlab/cellphonedb/master/in/example_data/test_meta.txt

counts里面是表达谱文件：

meta是细胞类型和分组信息

CellPhoneDB主要分成method、plot、query和database4个模块。

我们主要进行分析的是method（分析模块）和plot模块（画图模块）。

query是进行数据查询的模块，查询基因有哪些互作结果（get_interaction_gene结果为数据库涉及到的基因，find_interactions_by_element可以找到特定基因的受体配体作用数据对）。

database是输入数据库，一般默认可以不输入，直接使用即可。

运行程序：

cellphonedb method analysis test_meta.txt test_counts.txt

当然还有一些其它可供选择的参数：例如threads设置线程的，output-path设置输出文件夹的等。

运行完之后会有一个叫做默认out的输出文件夹，里面会有相应的pvalues.txt means.txt significant_means.txt等输出文件。

绘制气泡图：

cellphonedb plot dot_plot

当然可以通过参数：--pvalues-path --output-path --output-name来设置输出和输入的文件路径和名字。

绘制热图：

cellphonedb plot heatmap_plot test_meta.txt

当然如果觉得官方画的不好看，也可以自己提取输出文件的信息来美化图片。

library(tidyverse)

library(DT)

create_dt <- function(x)

{

DT::datatable(x,extensions='Buttons',options=list(dom='Blfrtip',buttons=c('csv','excel','pdf'),lengthMenu=list(c(10,25,50,-1),c(10,25,50,"All"))))

}

//读取文件

mypvals <- read.delim("pvalues.txt",check.names=FALSE)

mymeans <- read.delim("means.txt",check.names=FALSE)

mysigmeans <- read.delim("significant_means.txt",check.names=FALSE)

dim(mysigmeans)

mymeans[1:4,1:4]

//调整配体-受体顺序，配体在前，受体在后

library(dplyr)

order_sequence <- function(df)

{

da<-data.frame()

for(i in 1:length(df$gene_a))

{

sub_data <- df[i,]

if(sub_data$receptor_b=="False")

{

if(sub_data$receptor_a=="True")

{

old_names<- colnames(sub_data)

my_list<-strsplit(old_names[-c(1:11)],split="\\|")

my_character <- paste(sapply(my_list,'[[',2L),sapply(my_list,'[[',1L),sep='|')

new_names<-c(names(sub_data)[1:4],"gene_b","gene_a","secreted","receptor_b","receptor_a","annotation_strategy","is_integrin",my_character)

sub_data=dplyr::select(sub_data,new_names)

names(sub_data)<-old_names

da=rbind(da,sub_data)

}

else

{

da=rbind(da,sub_data)

}

return(da)

}

df<-subset(mymeans,receptor_a=="True"&receptor_b=="False"|receptor_a=="False"&receptor_b=="True")

df <- df %>% dplyr::mutate(na_count=rowSums(is.na(df)|df=="")) %>% subset(na_count==0) %>% dplyr::select(-na_count)

dim(df)

means_order <- order_sequence(df) %>% tidyr::unite(Pairs,gene_a,gene_b)

pvals_order <- order_sequence(mypvals) %>% tidyr::unite(Pairs,gene_a,gene_b)

//合并表达量文件和pvalue文件

means_sub <- means_order[,c('Pairs',colnames(mymeans)[-c(1:11)])]

pvals_sub <- pvals_order[,c('Pairs',colnames(mymeans)[-c(1:11)])]

means_gather <- tidyr::gather(means_sub,celltype,mean_exprssion,names(means_sub)[-1])

pvals_gather <- tidyr::gather(pvals_sub,celltype,pval,names(pvals_sub)[-1])

mean_pval <- dplyr::left_join(means_gather,pvals_gather,by=c('Pairs','celltype'))

create_dt(mean_pval)

//提取显著性表达的受体配体对

a <- mean_pval %>% dplyr::select(Pairs,celltype,pval) %>% tidyr::spread(key=celltype,value=pval)

sig_pairs <- a[which(rowSums(a<=0.05)!=0),]

dim(sig_pairs)

mean_pval_sub <- subset(mean_pval,Pairs %in% sig_pairs$Pairs)

//可视化显著性表达的受体配体对

library(RColorBrewer)

library(scales)

library(ggplot2)

library(cowplot)

p1 <- mean_pval_sub %>% ggplot(aes(x=mean_exprssion)) +geom_density(alpha=0.2,color='red')

p2 <- mean_pval_sub %>% ggplot(aes(x=log2(mean_exprssion))) +geom_density(alpha=0.2,color='red')

mean_distribution <- plot_grid(p1,p2,labels="AUTO",nrow=1)

p1 <- mean_pval_sub %>% ggplot(aes(x=pval)) +geom_density(alpha=0.2,color='red')

p2 <- mean_pval_sub %>% ggplot(aes(x=-log10(pval+1*10^-3))) +geom_density(alpha=0.2,color='red')

pval_distribution <- plot_grid(p1,p2,labels="AUTO",nrow=1)

plot_grid(mean_distribution,pval_distribution,labels="AUTO",ncol=1)

p <- mean_pval_sub %>% dplyr::arrange(Pairs) %>% head(16*10) %>% ggplot(aes(x=Pairs,y=celltype)) +

geom_point(aes(color=log2(mean_exprssion),size=-log10(pval+1*10^-3))) +

guides(colour=guide_colourbar(order=1),size=guide_legend(order=2)) +

labs(x='',y='') +

scale_color_gradientn(name='Expression level',colours=terrain.colors(100)) +

theme(axis.text.x=element_text(angle=45,hjust=1)) +

theme(panel.border=element_rect(color='black',fill=NA),panel.grid.major.x=element_blank(),panel.grid.major.y=element_blank(),panel.grid.minor.x=element_blank(),panel.grid.minor.y=element_blank(),panel.background=element_blank(),axis.title.x=element_blank(),axis.title.y=element_blank(),axis.ticks=element_blank())+

theme(legend.key.size=unit(0.4,'cm'),legend.title=element_text(size=9),legend.text=element_text(size=8))

save_plot("ligands_Receptor_dotplot.png",p,base_height=4,base_aspect_ratio=3,base_width=NULL,dpi=600)

save_plot("ligands_Receptor_dotplot.pdf",p,base_height=4,base_aspect_ratio=3,base_width=NULL)

====真实例子测试====

还用最常用的pbmc3k，上次做轨迹分析的例子。因为我们已经保存了，所以直接load进来。

pbmc3k <- readRDS("pbmc.rds")

//准备cellphone的输入文件

write.table(as.matrix(pbmc3k@assays$RNA@data), 'cellphonedb_count.txt', sep='\t', quote=F)

meta_data <- cbind(rownames(pbmc3k@meta.data), pbmc3k@meta.data[,'cell_type', drop=F])

meta_data <- as.matrix(meta_data)

meta_data[is.na(meta_data)] = "Unkown"

write.table(meta_data, 'cellphonedb_meta.txt', sep='\t', quote=F, row.names=F)

//运行cellphone

cellphonedb method statistical_analysis cellphonedb_meta.txt cellphonedb_count.txt --counts-data=gene_name --threads 20

cellphonedb plot dot_plot

cellphonedb plot heatmap_plot cellphonedb_meta.txt

//读入网络数据，构建网络，并绘制最简单的一张网络图：

pbmc='out/'

library(psych)

library(qgraph)

library(igraph)

library(tidyverse)

netf<- "count_network.txt"

mynet <- read.delim(paste0(pbmc,"count_network.txt"), check.names = FALSE)

table(mynet$count)

mynet %>% filter(count>0) -> mynet

head(mynet)

net<- graph_from_data_frame(mynet)

plot(net)

//为了给这个网络图的边点mapping上不同的属性我们引入一串颜色。

allcolour=c("#DC143C","#0000FF","#20B2AA","#FFA500","#9370DB",

"#98FB98","#F08080","#1E90FF","#7CFC00","#FFFF00",

"#808000","#FF00FF","#FA8072","#7B68EE","#9400D3",

"#800080","#A0522D","#D2B48C","#D2691E","#87CEEB",

"#40E0D0","#5F9EA0","#FF1493",

"#FFE4B5","#8A2BE2","#228B22","#E9967A","#4682B4",

"#32CD32","#F0E68C","#FFFFE0","#EE82EE","#FF6347",

"#6A5ACD","#9932CC","#8B008B","#8B4513","#DEB887")

karate_groups <- cluster_optimal(net)

coords <- layout_in_circle(net, order =order(membership(karate_groups))) #设置网络布局

E(net)$width <- E(net)$count/10 # 边点权重（粗细）

plot(net, edge.arrow.size=.1,

edge.curved=0,

vertex.color=allcolour,

vertex.frame.color="#555555",

vertex.label.color="black",

layout = coords,

vertex.label.cex=.7)

经过以上一番修饰后，我们得到的网络图如下：

但是边的颜色和点的颜色还是对应不上，我们来修改一番边的属性。

net2 <- net # 复制一份备用

for (i in 1: length(unique(mynet$SOURCE)) ){

E(net)[map(unique(mynet$SOURCE),function(x) {

get.edge.ids(net,vp = c(unique(mynet$SOURCE)[i],x))

})%>% unlist()]$color <- allcolour[i]

}

plot(net, edge.arrow.size=.1,

edge.curved=0,

vertex.color=allcolour,

vertex.frame.color="#555555",

vertex.label.color="black",

layout = coords,

vertex.label.cex=.7)

//把线条换成曲线

plot(net, edge.arrow.size=.1,

edge.curved=0.2, # 只是调了这个参数

vertex.color=allcolour,

vertex.frame.color="#555555",

vertex.label.color="black",

layout = coords,

vertex.label.cex=.7)

//接下来绘制，另一种贝壳一样的图，相当于查看每个单独的source的网络图

length(unique(mynet$SOURCE)) # 查看需要绘制多少张图，以方便布局

par(mfrow=c(2,5), mar=c(.3,.3,.3,.3))

for (i in 1: length(unique(mynet$SOURCE)) ){

net1<-net2

E(net1)[map(unique(mynet$SOURCE),function(x) {

get.edge.ids(net,vp = c(unique(mynet$SOURCE)[i],x))

})%>% unlist()]$color <- allcolour[i]

plot(net1, edge.arrow.size=.1,

edge.curved=0.4,

vertex.color=allcolour,

vertex.frame.color="#555555",

vertex.label.color="black",

layout = coords,

vertex.label.cex=1)

}

//可以加上频数

length(unique(mynet$SOURCE))

par(mfrow=c(2,5), mar=c(.3,.3,.3,.3))

for (i in 1: length(unique(mynet$SOURCE)) ){

net1<-net2

E(net1)$count <- ""

E(net1)[map(unique(mynet$SOURCE),function(x) {

get.edge.ids(net,vp = c(unique(mynet$SOURCE)[i],x))

})%>% unlist()]$count <- E(net2)[map(unique(mynet$SOURCE),function(x) {

get.edge.ids(net,vp = c(unique(mynet$SOURCE)[i],x))

})%>% unlist()]$count # 故技重施

E(net1)[map(unique(mynet$SOURCE),function(x) {

get.edge.ids(net,vp = c(unique(mynet$SOURCE)[i],x))

})%>% unlist()]$color <- allcolour[i]

plot(net1, edge.arrow.size=.1,

edge.curved=0.4,

edge.label = E(net1)$count, # 绘制边的权重

vertex.color=allcolour,

vertex.frame.color="#555555",

vertex.label.color="black",

layout = coords,

vertex.label.cex=1

)

}

//也可以只挑自己感兴趣的配体受体等进行点图的展示

chemokines <- grep("^CXC|CCL|CCR|CX3|XCL|XCR", mymeans$interacting_pair,value = T)

th1 <- grep("IL2|IL12|IL18|IL27|IFNG|IL10|TNF$|TNF |LTA|LTB|STAT1|CCR5|CXCR3|IL12RB1|IFNGR1|TBX21|STAT4",