我们在数据分析过程中，有了一群重点关注的基因（比如差异表达基因、预后相关基因、特定通路或signature成员基因），接下来就可以分析他们之间的关联。gene-gene相互作用网络除了可以用常见的string工具画PPI网络图，还可以纳入更多的信息，从而找到重点基因或突出差异。
比如下面这两幅图：
这种图有其他差异分析结果展示图，比如热图、火山图、箱线图等无法带来的好处：我们不仅可以看出哪些基因高表达或低表达了以及每个基因差异表达水平的高低，还以知道差异表达基因之间的关联--诸如他们是否处在同一通路上、表达模式上的相关性及其程度，甚至还能帮助我们找到其中的关键基因。那么我们就可以针对这些关键基因或通路进行后续分析。
类似的还可以画下面这种图：
上图表示的是细胞类型之间的关联，可以一眼看出在肿瘤中最强烈相关的几种细胞类型。
这种网络图的重点不外乎vertices(点)的大小和颜色，edges(边)的粗细和颜色。抓住这点，就可以很容易的根据自己的数据复现出来类似的图。关键是要明白自己的数据里应该应用什么值来进行对应的设置。上面图1和图2两个图中都是用的差异表达倍数（log FoldChange）来设置每个基因节点的颜色（和大小）。那edges(边)是什么呢？可以用基因功能相似性系数也可以用表达值相关性系数。
画网络图的工具除了大名鼎鼎的Cytoscape外，还有很多其他工具，其中就包括R的igraph包。接下来我们就分别用这两个来复现这个图。
###安装包
BiocManager::install("BioCor")
###
library("BioCor")
library("org.Hs.eg.db")
#library("reactome.db")
library(msigdbr)
library(clusterProfiler)
library(enrichplot)
library(reshape2)
library(corrplot)
library(plyr)
library(igraph)
###############################数据准备
###参考文献中用的GO数据库中BP分类下的通路，当然也可使使用其他通路
hsGO <- msigdbr(species = "Homo sapiens", category = "C5",subcategory='BP')# %>% dplyr::select(gs_exact_source,gs_name, gene_symbol)
###差异分析结果
tumor_ec_markers0=read.xlsx('41467_2020_16164_MOESM8_ESM.xlsx',sheet = 3)
tumor_ec_markers0[1:5,]
Genes log2FC p-value Bonferroni PCT1 PCT2
INSR 2.147045 7.910142e-57 1.328666e-52 0.585 0.126
HSPG2 1.818224 4.165828e-120 6.997341e-116 0.864 0.321
VWA1 1.716222 2.147301e-71 3.606822e-67 0.623 0.127
PLVAP 1.695246 1.237309e-81 2.078309e-77 0.798 0.294
IGHG3 1.694869 6.457962e-20 1.084744e-15 0.302 0.060
###因为GO中很多通路是高度相似的，直接用会导致很多基因功能相似性过高，导致网络图中边的数量过多，所以先做富集得到重点通路
bp <- enrichGO(tumor_ec_markers0$Genes,keyType='SYMBOL', ont="BP", OrgDb = 'org.Hs.eg.db')
bp <- pairwise_termsim(bp)
bp2 <- simplify(bp, cutoff=0.7, by="p.adjust", select_fun=min)
###
GO_terms=bp2@result$ID#[bp2@result$pvalue<=0.05]
hsGO_tmp=hsGO[hsGO$gs_exact_source%in%GO_terms,]
hsGO_list=split(hsGO_tmp$gs_name,hsGO_tmp$gene_symbol)
###用Jaccard index系数计算基因之间的功能相似性
goSemSim <- mgeneSim(names(hsGO_list), info= hsGO_list,method = 'avg') ##avg , reciprocal
jaccard_index=D2J(goSemSim)
a=reshape2::melt(lower.tri(jaccard_index))
jaccard_index_df= reshape2::melt(jaccard_index)
jaccard_index_df=jaccard_index_df[a$value,]
#只保留值相似性较高的
jaccard_index_df0=jaccard_index_df[abs(jaccard_index_df$value)>0.15,]
############### 设置网络图的节点和边
nodes=tumor_ec_markers0
edges=jaccard_index_df0
###设置节点的特征
#用pvalue控制节点的大小
nodes$size <- c(5,6,8)[cut(abs(log10(nodes$Bonferroni)),3)]
#用FoldChange控制节点的颜色和深浅
nodes$cut=cut(nodes$log2FC,breaks = c(seq(min(nodes$log2FC)-0.1,-0.5,lengt=5),seq(0.5,max(nodes$log2FC)+0.1,lengt=5)))
# nodes$cut=cut(nodes$log2FC,10)
cols_nodes=colorRampPalette(colors = c("#143CFF","#F0F0F0","#FF0000"), interpolate ="linear")(9)
names(cols_nodes)=levels(nodes$cut)
nodes$color=cols_nodes[nodes$cut]
###设置边的特征
#edges$width=round(abs(edges$value)*10,0)
edges$width=c(2,2.5,8)[cut(edges$value,3)]
edges$color='#9C9C9C'
###############开始画图
net <- graph_from_data_frame(d=edges, vertices=nodes, directed=F)
V(net)$color <-nodes$color
V(net)$frame.color <- '#4F4F4F'#nodes$frame_color
V(net)$shape='circle'
V(net)$size <- nodes$size*3
V(net)$label.color='black'
V(net)$label.size=nodes$size
E(net)$arrow.mode <- 0 #0:不需要箭头,1:back; 2:forth; 3:both
#E(net)$arrow.size=0.1
#E(net)$arrow.width=0.1
E(net)$edge.color <- edges$color
E(net)$width <- edges$width
plot(net,#vertex.label="",
rescale = T,
seed=1,
layout=layout_nicely,
edge.curved=0
)
##############保存结果，以便用于Cytoscape画图
write.csv(nodes,'nodes.csv')
write.csv(edges,'edges.csv')
从结果来看，这些基因至少可以分为2大簇:免疫反应相关基因(HLA-B, HLA-E, IL7R, IGHG3, IGKC等)基本都聚在了一起，血管生成相关基因(HPGD, SPARC, COL4A1, COL4A2, ANGPT2)也都基本聚在一起。
后续可以通过AI修饰来突出通路。效果如下：
因为这些基因是肺腺癌组织来源的内皮细胞相对于正常肺组织来源内皮细胞的显著差异表达基因，基因表达网络分析结果也高度吻合了研究预期：肿瘤中血管内皮细胞通过上调血管生成通路同时抑制细胞凋亡、免疫反应、脂肪酸等通路来促进肿瘤的发展。
[1]. Single-cell RNA sequencing demonstrates the molecular and cellular reprogramming of metastatic lung adenocarcinoma: Fig3d. Functional association networks between signature genes specific to tumor ECs
[2]. Deciphering cell lineage specification of human lung adenocarcinoma with single-cell RNAs equencing: Fig3d. Gene-gene interaction networks between marker genes in AT2 and AT2-like cell clusters.
[3].Single-cell RNA sequencing reveals distinct tumor microenvironmental patterns in lung adenocarcinoma: Fig5C. Correlation of the proportion of stromal and immune cell clusters, most connected section of correlation network plot shown; Spearman’s correlation statistics, only correlations with rho > 0.7 and p < 0.05 shown.
