背景

受配体介导的细胞之间的通讯分析对于协调发育、分化和炎症等多种生物学过程至关重要。CellPhoneDB是一个受配体及其相互作用的数据库，有别于其他数据库的是CellPhoneDB考虑了受配体的亚单位结构，准确的表达了异聚体复合物。除此之外，该软件允许客户引入自己的相互作用分子，已达到研究目的。

CellPhoneDB

安装

这里我们首先介绍一下该软件的安装，该软件依赖于python => 3.6 ，我们可以借助conda环境来进行安装：

#创建名为cellphonedb的虚拟环境
conda create -n cellphonedb python=3.7

#激活虚拟环境
conda activate cellphonedb

#在虚拟环境中下载软件
pip install cellphonedb
conda install -n cellphonedb r #由于后面我们需要使用R语言绘制网络图，因此在这里先下载R

#获取帮助信息
cellphonedb --help #这里只能使用--help不能使用-h，不明白为什么这样干QAQ
Usage: cellphonedb [OPTIONS] COMMAND [ARGS]...

Options:
  --help  Show this message and exit.

Commands:
  database
  method
  plot
  query

其中4个模块分别是用来进行数据库更新下载、统计分析、画图以及提取信息。这里将重点介绍method模块。

cellphonedb method --help
Usage: cellphonedb method [OPTIONS] COMMAND [ARGS]...

Options:
  --help  Show this message and exit.

Commands:
  analysis
  statistical_analysis

下面我们以进行了统计分析的statistical_analysis算法为例进行介绍，首先运行程序：

#下载测试数据
curl https://raw.githubusercontent.com/Teichlab/cellphonedb/master/in/example_data/test_counts.txt --output test_counts.txt
curl https://raw.githubusercontent.com/Teichlab/cellphonedb/master/in/example_data/test_meta.txt --output test_meta.txt
cat test_counts.txt #表达谱
1 Gene    d-pos_AAACCTGAGCAGGTCA  d-pos_AAACCTGGTACCGAGA  d-pos_AAACCTGTCGCCATAA  d-pos_AAACGGGTCAGTTGAC  d-pos_AAAGATGCATTGAGCT  d-pos_AAAGATGTCCAAAGTC  d-pos_AAAGCAAAGAGGACGG  d-pos_AAAGCAACACA
2 ENSG00000238009 0       0       0       0       0       0       0       0       0       0
3 ENSG00000279457 0       0       0       0       0.723769155614119       0       1.1269975757326 1.81828622356148        0       0
4 ENSG00000228463 0       0       0       0       0       0       0       0       0.737864655764131       1.40825228390187
5 ENSG00000237094 0       0       0       0       0       0       0       0       0       0
6 ENSG00000230021 0       0       0       0       0       0       0       0       0       0
7 ENSG00000237491 0       0       0       0       0       0       0       0       0       0
8 ENSG00000177757 0       0       0       0       0       0       0       0       0       0

cat test_meta.txt  #细胞类型或者分组信息
Cell    cell_type
d-pos_AAACCTGAGCAGGTCA  NKcells_1
d-pos_AAACCTGGTACCGAGA  NKcells_0
d-pos_AAACCTGTCGCCATAA  NKcells_1
d-pos_AAACGGGTCAGTTGAC  Tcells
d-pos_AAAGATGCATTGAGCT  NKcells_0
d-pos_AAAGATGTCCAAAGTC  NKcells_0
d-pos_AAAGCAAAGAGGACGG  Myeloid
d-pos_AAAGCAACACATTCGA  NKcells_1
d-pos_AAAGTAGAGAGCCCAA  NKcells_0
d-pos_AAAGTAGCAAGCTGAG  NKcells_0

#运行程序
cellphonedb method statistical_analysis --counts-data ensembl --output-path test_output test_meta.txt test_counts.txt
#气泡图
cellphonedb plot dot_plot --pvalues-path test/pvalues.txt --output-path test_output  --output-name test.dotplot.pdf
#热图
cellphonedb plot heatmap_plot --pvalues-path test/pvalues.txt --output-path test_output --pvalue 0.05 --count-name test.heatmap_count.pdf --log-name test.heatmap_log_count.pdf --count-network-name test.count_network.txt --interaction-count-name test.interaction_count.txt test_meta.txt

ll test_output
-rw-r--r-- 1 yangyupeng research  18258 Jun 24 15:45 deconvoluted.txt
-rw-r--r-- 1 yangyupeng research  22772 Jun 24 15:45 means.txt
-rw-r--r-- 1 yangyupeng research  22651 Jun 24 15:45 pvalues.txt
-rw-r--r-- 1 yangyupeng research  17616 Jun 24 15:45 significant_means.txt
-rw-r--r-- 1 yangyupeng research    344 Jun 24 15:46 test.count_network.txt
-rw-r--r-- 1 yangyupeng research  76816 Jun 24 15:46 test.dotplot.pdf
-rw-r--r-- 1 yangyupeng research  10533 Jun 24 15:46 test.heatmap_count.pdf
-rw-r--r-- 1 yangyupeng research  10530 Jun 24 15:46 test.heatmap_log_count.pdf
-rw-r--r-- 1 yangyupeng research     56 Jun 24 15:46 test.interaction_count.txt

气泡图：

气泡图

热图：

热图

网络图

不同于气泡图和热图，cellphonedb软件只计算了受配体的网络关系，如上述结果中的test.count_network.txt：

SOURCE  TARGET  count
NKcells_1       NKcells_1       5
NKcells_1       NKcells_0       14
NKcells_1       Tcells  8
NKcells_1       Myeloid 12
NKcells_0       NKcells_1       14
NKcells_0       NKcells_0       15
NKcells_0       Tcells  11
NKcells_0       Myeloid 32
Tcells  NKcells_1       8
Tcells  NKcells_0       11
Tcells  Tcells  1
Tcells  Myeloid 17
Myeloid NKcells_1       12
Myeloid NKcells_0       32
Myeloid Tcells  17
Myeloid Myeloid 20

因此网络图需要我们自己绘制，用于绘制网络图的常用软件软件有2个：Cytoscape和igraph，这里我们选择igraph作为我们的绘图软甲。igraph软件有R、python、Mathematica和C/C++4个版本，这里我们采用R语言的版本。
rigraph的安装有许多坑需要踩，这里我就帮助各位踩了。首先需要安装一个叫做GLPK的软件作为igraph的依赖。

#使用conda安装，方便快捷无污染
conda install glpk
conda install r-psych r-tidyverse r-ggsci

然后在R端安装igraph，这里有一个坑需要注意一下，你可以直接使用conda安装igraph，但是在安装完之后的使用有没有问题没法确定，没有问题最好。如果有问题，作者在GitHub上的issue里也提到，由于CRPN不允许作者把GLPK嵌套进模块中，所以使用者需要额外安装上述提到的GLPK，而这一个过程也导致了后续可能会发生的一些错误，所以作者建议使用源代码编译（虽然这里有坑，当然这里以及后面提到的一些坑更多的是针对centos linux的conda环境，如果你使用的是window那么恭喜你，这没啥坑可踩），所以这里我们使用的是github上的源码安装。（就比如悲催的小编，鬼知道我找了多少种方法，在这里卡了好几天找不到解决的办法，最后误打误撞在用源码安装GLPK的时候在安装步骤里借鉴到了一个解决思路，顺利解决，撒花）

if (!required(remotes)){
  install.packages("remotes")
}
if (!required(qgraph)){
  install.packages("qgraph")
}

#这里先不要安装igraph因为在这里有一个天坑，
#igraph在安装的时候会读取你环境的glpk，但是
#不知道为什么我在安装的时候R死活读取不到
#glpk，并且在安装过程中也不会报错，对你没看错
#读取不到这个glpk他也不会报错，如果你安装的时候
#不注意的话，后面如果报了找不到GLPK的error，你会
#一脸懵逼的不知道那里出了问题，并且找不到解决的
#办法，你以为添加了环境变量就好了，其实完全没用
#这是在安装的时候就已经埋好的坑QAQ。

#先配置GLPK的环境，告诉R GLPK要到哪里去找
Sys.setenv(CPPFLAGS="-DNDEBUG -D_FORTIFY_SOURCE=2 -O2 -isystem ~/miniconda2/envs/cellphonedb/include")
Sys.setenv(LDFLAGS="-Wl,-O2 -Wl,--sort-common -Wl,--as-needed -Wl,-z,relro -Wl,-z,now -Wl,--disable-new-dtags -Wl,--gc-sections -Wl,-rpath,~/miniconda2/envs/cellphonedb/lib -Wl,-rpath-link,~/miniconda2/envs/cellphonedb/lib -L~/miniconda2/envs/cellphonedb/lib")
remotes::install_github("igraph/rigraph@master")


library(psych)
library(igraph)
library(qgraph)
library(tidyverse)
library(ggsci)

#输入文件，并确定数据中所有行的count值均大于0
network_df <- read.table("/path/to/test.count_network.txt", sep = "\t", quote = "", header = TRUE)
network_df <- network_df %>% filter(count > 0)
#建立网络
net <- graph_from_data_frame(network_df)
#提取顶点位置节点的坐标
edge.start <- igraph::ends(net, es=igraph::E(net), names=FALSE)
#  选取ggsci保重的d3-category20c调色板
col <- pal_d3("category20c")(length(unique(network_df$SOURCE)))
names(col) <- unique(network_df$cluster)
#计算社区结构，并按照社区结构对位点进行排列
group <-  cluster_optimal(net)
coords <- layout_in_circle(net, order = order(membership(group)))
#计算edge的权重
E(net)$width <- E(net)$count / 10
vextex.size <- 20#当然这里可以取权重比如vextex.size <- max(vertex.weight) / vertex.weight * V(net)$count + 5
#调节线条和label的配色
E(net)$color <- col[E(net)$cluster]
V(net)$label.color <- "black"

#调整节点位置的线条角度，这里一开始也挺挠头的，因为默认定点的线是一致向右的，感谢cellchat，在这个软件的源代码中找到了如何调整，很巧妙的一个方法
loop.angle<-ifelse(coords[igraph::V(net),1]>0,-atan(coords[igraph::V(net),2]/coords[igraph::V(net),1]),pi-atan(coords[igraph::V(net),2]/coords[igraph::V(net),1]))
igraph::E(net)$loop.angle[which(edge.start[,2]==edge.start[,1])] <- loop.angle[edge.start[which(edge.start[,2]==edge.start[,1]),1]]
plot(net,
     edge.arrow.size = 1,
     edge.curved = 0.2,
     vertex.color = col,
     vertex.frame.color = F,
     layout = coords,
     vertex.label.cex = 1,
     vertex.size = vextex.size)

network

当然我们也可以将每个配体分开进行展示：

lapply(names(col), function(x){
    E(net1)$color[E(net1)$cluster != x] <- NA
    E(net1)$count[E(net1)$cluster != x] <- NA
    par(mfrow=c(nrow,ncol), mar = c(.5,.5,.5,.5))
    p <- plot(net1,
              edge.arrow.size = 1,
              edge.curved = 0.2,
              edge.label = E(net1)$count,
              edge.label.color = "black",
              vertex.color = col,
              vertex.frame.color = F,
              layout = coords,
              vertex.label.cex = 1)
    return(p)
  })

network

好了，通常我们文章里这两幅图用的最多，基本上也够用了。通过这次的经历发现igraph这个包还挺有意思的，后面有机会专门做一期介绍igraph，再见！！！