写在前面
我8月26号踏着欢快的脚步回家度假，原计划9月5号回珠海，结果返程飞机被连环取消了三次，我滞留在家5天了，已经搭好戏台子在家讲课，打算等没课的周天再走。。。
1.从一个表达矩阵开始
相关性网络图一般都是经过若干步骤选出来几个感兴趣的基因展示一下，下面的这个例子用我的小R包tinyarray获取表达矩阵和做探针注释、差异分析了，相当于对芯片常规分析的一个简化操作。
#devtools::install_github("xjsun1221/tinyarray")
library(tinyarray)
gse = "GSE42872"
geo = geo_download(gse,destdir=tempdir(),by_annopbrobe = FALSE)
Group = rep(c("control","treat"),each = 3)
Group = factor(Group)
find_anno(geo$gpl)
## [1] "`library(hugene10sttranscriptcluster.db);ids <- toTable(hugene10sttranscriptclusterSYMBOL)` and `ids <- AnnoProbe::idmap('GPL6244')` are both avaliable"
ids <- AnnoProbe::idmap(geo$gpl,destdir = tempdir())
deg = get_deg(geo$exp,Group,ids)
这个deg就是经过整理的差异分析结果表格。
head(deg)
## logFC AveExpr t P.Value adj.P.Val B probe_id
## 1 5.780170 7.370282 82.94833 3.495205e-12 1.163798e-07 16.32898 8133876
## 2 -4.212683 9.106625 -68.40113 1.437468e-11 2.393169e-07 15.71739 7965335
## 3 5.633027 8.763220 57.61985 5.053466e-11 4.431880e-07 15.04752 7972259
## 4 -3.801663 9.726468 -57.21112 5.324059e-11 4.431880e-07 15.01709 7972217
## 5 3.263063 10.171635 50.51733 1.324638e-10 8.821294e-07 14.45166 8129573
## 6 -3.843247 9.667077 -45.87910 2.681063e-10 1.487856e-06 13.97123 8015806
## symbol change ENTREZID
## 1 CD36 up 948
## 2 DUSP6 down 1848
## 3 DCT up 1638
## 4 SPRY2 down 10253
## 5 MOXD1 up 26002
## 6 ETV4 down 2118
把表达矩阵转换一下，并选了top10差异基因用于做后续的网络图
exp = trans_array(geo$exp,ids)
cg = deg$symbol[deg$change!="stable"]
set.seed(10086)
exp = exp[head(cg,10),]
exp[1:4,1:4]
## GSM1052615 GSM1052616 GSM1052617 GSM1052618
## CD36 4.54610 4.40210 4.49239 10.25060
## DUSP6 11.25310 11.20760 11.17820 6.98791
## DCT 5.81479 5.91209 6.11324 11.54910
## SPRY2 11.64830 11.63730 11.59630 7.90508
2.计算相关性
根据相关系数和p值做了个分组，用于后面的配色了，我用的是0.3，这个阈值可以调整。
library(corrplot)
M = cor(t(exp))
testRes = cor.mtest(t(exp), conf.level = 0.95)$p
library(tidyverse)
g = pivot_longer(rownames_to_column(as.data.frame(M),var = "from"),
cols = 2:(ncol(M)+1),
names_to = "to",
values_to = "cor")
gp = pivot_longer(rownames_to_column(as.data.frame(testRes)),
cols = 2:(ncol(M)+1),
names_to = "gene",
values_to = "p")
g$p = gp$p
g = g[g$from!=g$to,]
g$group = case_when(g$cor>0.3 & g$p<0.05 ~ "positive",
g$cor< -0.3 & g$p<0.05 ~ "negative",
T~"not" )
head(g)
## # A tibble: 6 x 5
## from to cor p group
## <chr> <chr> <dbl> <dbl> <chr>
## 1 CD36 DUSP6 -1.00 0.0000000933 negative
## 2 CD36 DCT 0.999 0.00000196 positive
## 3 CD36 SPRY2 -1.00 0.000000226 negative
## 4 CD36 MOXD1 1.00 0.0000000840 positive
## 5 CD36 ETV4 -0.999 0.00000208 negative
## 6 CD36 DTL -0.999 0.00000168 negative
上面这个g就是网络图的输入数据了
3.画图
正相关和负相关分别用红色和蓝色表示咯。
library(igraph)
network = graph_from_data_frame(d=g[g$group!="not",c(1,2,3,5)], directed=F)
my_color = c("#2874C5","#f87669")[as.numeric(as.factor(E(network)$group))]
par(bg="white", mar=c(0,0,0,0))
plot(network,
vertex.size=20,
layout=layout.circle,
vertex.label.cex=0.7,
vertex.frame.color="transparent",
edge.width=abs(E(network)$cor)*3,
edge.color=my_color,edge.curved = 0.2)
4.专用的相关性网络图R包
发现一个哭笑不得的现象，如果这些基因相关性太强了，这个专用的R包画图还叠到一起了，不是个例，我换电脑试过了哈哈。
library(tidyverse)
t(exp) %>% corrr::correlate() %>%
corrr::network_plot(colours = c("#2874C5", "white", "#f87669"),repel = F,min_cor = .5)
反而是他们的相关性不特别强时，才比较好看，错落有致了咧。
exp = trans_array(geo$exp,ids)
set.seed(100)
cg = sample(1:nrow(exp),10)
exp = exp[cg,]
t(exp) %>% corrr::correlate() %>%
corrr::network_plot(colours = c("#2874C5", "white", "#f87669"),repel = F,min_cor = .5)
奇妙的体验。
