寫在前面
我8月26號踏著歡快的腳步回家度假，原計劃9月5號回珠海，結(jié)果返程飛機被連環(huán)取消了三次，我滯留在家5天了，已經(jīng)搭好戲臺子在家講課，打算等沒課的周天再走。。。

1.從一個表達矩陣開始

相關性網(wǎng)絡圖一般都是經(jīng)過若干步驟選出來幾個感興趣的基因展示一下，下面的這個例子用我的小R包tinyarray獲取表達矩陣和做探針注釋、差異分析了，相當于對芯片常規(guī)分析的一個簡化操作。

#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就是經(jīng)過整理的差異分析結(jié)果表格。

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

把表達矩陣轉(zhuǎn)換一下，并選了top10差異基因用于做后續(xù)的網(wǎng)絡圖

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.計算相關性

根據(jù)相關系數(shù)和p值做了個分組，用于后面的配色了，我用的是0.3，這個閾值可以調(diào)整。

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就是網(wǎng)絡圖的輸入數(shù)據(jù)了

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.專用的相關性網(wǎng)絡圖R包

發(fā)現(xiàn)一個哭笑不得的現(xiàn)象，如果這些基因相關性太強了，這個專用的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)

奇妙的體驗。