我們在數(shù)據(jù)分析過程中击碗，有了一群重點關(guān)注的基因（比如差異表達基因鸡岗、預(yù)后相關(guān)基因鸳劳、特定通路或signature成員基因）婚瓜，接下來就可以分析他們之間的關(guān)聯(lián)宝鼓。gene-gene相互作用網(wǎng)絡(luò)除了可以用常見的string工具畫PPI網(wǎng)絡(luò)圖，還可以納入更多的信息闰渔，從而找到重點基因或突出差異席函。
比如下面這兩幅圖：
這種圖有其他差異分析結(jié)果展示圖铐望，比如熱圖冈涧、火山圖、箱線圖等無法帶來的好處：我們不僅可以看出哪些基因高表達或低表達了以及每個基因差異表達水平的高低正蛙，還以知道差異表達基因之間的關(guān)聯(lián)--諸如他們是否處在同一通路上督弓、表達模式上的相關(guān)性及其程度，甚至還能幫助我們找到其中的關(guān)鍵基因乒验。那么我們就可以針對這些關(guān)鍵基因或通路進行后續(xù)分析愚隧。
類似的還可以畫下面這種圖：
上圖表示的是細胞類型之間的關(guān)聯(lián)，可以一眼看出在腫瘤中最強烈相關(guān)的幾種細胞類型锻全。
這種網(wǎng)絡(luò)圖的重點不外乎vertices(點)的大小和顏色狂塘，edges(邊)的粗細和顏色录煤。抓住這點，就可以很容易的根據(jù)自己的數(shù)據(jù)復(fù)現(xiàn)出來類似的圖荞胡。關(guān)鍵是要明白自己的數(shù)據(jù)里應(yīng)該應(yīng)用什么值來進行對應(yīng)的設(shè)置妈踊。上面圖1和圖2兩個圖中都是用的差異表達倍數(shù)（log FoldChange）來設(shè)置每個基因節(jié)點的顏色（和大小）泪漂。那edges(邊)是什么呢廊营？可以用基因功能相似性系數(shù)也可以用表達值相關(guān)性系數(shù)。
畫網(wǎng)絡(luò)圖的工具除了大名鼎鼎的Cytoscape外萝勤，還有很多其他工具露筒，其中就包括R的igraph包。接下來我們就分別用這兩個來復(fù)現(xiàn)這個圖敌卓。
###安裝包
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)
###############################數(shù)據(jù)準備
###參考文獻中用的GO數(shù)據(jù)庫中BP分類下的通路慎式，當然也可使使用其他通路
hsGO <- msigdbr(species = "Homo sapiens", category = "C5",subcategory='BP')# %>% dplyr::select(gs_exact_source,gs_name, gene_symbol)
###差異分析結(jié)果
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中很多通路是高度相似的，直接用會導(dǎo)致很多基因功能相似性過高假哎，導(dǎo)致網(wǎng)絡(luò)圖中邊的數(shù)量過多瞬捕，所以先做富集得到重點通路
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系數(shù)計算基因之間的功能相似性
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,]
############### 設(shè)置網(wǎng)絡(luò)圖的節(jié)點和邊
nodes=tumor_ec_markers0
edges=jaccard_index_df0
###設(shè)置節(jié)點的特征
#用pvalue控制節(jié)點的大小
nodes$size <- c(5,6,8)[cut(abs(log10(nodes$Bonferroni)),3)]
#用FoldChange控制節(jié)點的顏色和深淺
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]
###設(shè)置邊的特征
#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
)
##############保存結(jié)果，以便用于Cytoscape畫圖
write.csv(nodes,'nodes.csv')
write.csv(edges,'edges.csv')
從結(jié)果來看舵抹，這些基因至少可以分為2大簇:免疫反應(yīng)相關(guān)基因(HLA-B, HLA-E, IL7R, IGHG3, IGKC等)基本都聚在了一起肪虎，血管生成相關(guān)基因(HPGD, SPARC, COL4A1, COL4A2, ANGPT2)也都基本聚在一起。
后續(xù)可以通過AI修飾來突出通路惧蛹。效果如下：
因為這些基因是肺腺癌組織來源的內(nèi)皮細胞相對于正常肺組織來源內(nèi)皮細胞的顯著差異表達基因扇救，基因表達網(wǎng)絡(luò)分析結(jié)果也高度吻合了研究預(yù)期：腫瘤中血管內(nèi)皮細胞通過上調(diào)血管生成通路同時抑制細胞凋亡、免疫反應(yīng)香嗓、脂肪酸等通路來促進腫瘤的發(fā)展迅腔。
[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.
