放弃Venn-Upset-花瓣图,拥抱二分网络
生物信息學習的正確姿勢
NGS系列文章包括NGS基礎、在線繪圖、轉錄組分析?(Nature重磅綜述|關于RNA-seq你想知道的全在這)、ChIP-seq分析?(ChIP-seq基本分析流程)、單細胞測序分析?(重磅綜述:三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程)、DNA甲基化分析、重測序分析、GEO數據挖掘(典型醫學設計實驗GEO數據分析 (step-by-step))、批次效應處理等內容。
寫在前面
讓點隨機排布在一個區域,保證點之間不重疊,并且將點的圖層放到最上層,保證節點最清晰,然后邊可以進行透明化,更加突出節點的位置。這里我新構建了布局函數 PolyRdmNotdCirG 來做這個隨機排布。調用的是packcircles包的算法。使用和其他相似函數一樣,這里我們重點介紹一下使用這種算法構造的二分網絡布局。
微生物網絡
ggClusterNet 安裝
ggClusterNet包依賴的R包均在cran或者biocductor中,所以未能成功安裝,需要檢查依賴是否都順利安裝。如果網路問題,無法下載R包,可以在github中手動下載安裝:具體安裝方法參考:玩轉R包
#---ggClusterNet devtools::install_github("taowenmicro/ggClusterNet") #--如果無法安裝請檢查網絡或者換個時間導入R包和輸入文件
#--導入所需R包#------- library(ggplot2) library(ggrepel) library(ggClusterNet) library(phyloseq) library(dplyr)# 數據內置 #-----導入數據#------- data(ps)#--可選 #-----導入數據#------- ps = readRDS("../ori_data/ps_liu.rds")這里我們提取一部分OTU,節省出圖時間。
# ps data(ps)ps_sub = filter_taxa(ps, function(x) sum(x ) > 20 , TRUE) ps_sub = filter_taxa(ps_sub, function(x) sum(x ) < 30 , TRUE) ps_subdiv_network函數 用于計算共有和特有關系
這個函數是之前我寫的專門用于從OTU表格整理成Gephi的輸入文件,所以大家直接用這個函數即可轉到gephi進行操作。這次為了配合二分網絡,我設置了參數flour = TRUE,代表是否僅僅提取共有部分和特有部分。
# ?div_network result = div_network(ps_sub,num = 6)edge = result[[1]] head(edge)# levels(edge$target) # node = result[[2]] # head(node) # # tail(node) data = result[[3]] dim(data)#----計算節點坐標 # flour參數,設置是否僅僅展示共有和特有的二分網絡div_culculate函數 核心算法,用于計算二分網絡的節點和邊的表格
參數解釋:
distance = 1.1:
中心一團點到樣本點距離
distance2 = 1.5:
中心點模塊到獨有OTU點之間距離
distance3 = 1.3:
樣本點和獨有OTU之間的距離
order = FALSE :
節點是否需要隨機擾動效果
對OTU進行注釋,方便添加到圖形上
為了讓節點更加豐富,這里我對節點文件添加了注釋信息。
# table(plotdata$elements) node = plotdata[plotdata$elements == unique(plotdata$elements), ]otu_table = as.data.frame(t(vegan_otu(ps_sub))) tax_table = as.data.frame(vegan_tax(ps_sub)) res = merge(node,tax_table,by = "row.names",all = F) dim(res) head(res) row.names(res) = res$Row.names res$Row.names = NULL plotcord = resxx = data.frame(mean =rowMeans(otu_table)) head(xx) plotcord = merge(plotcord,xx,by = "row.names",all = FALSE) head(plotcord) # plotcord$Phylum row.names(plotcord) = plotcord$Row.names plotcord$Row.names = NULL head(plotcord)p = ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2),data = edge, size = 0.3,color = "yellow") +geom_point(aes(X1, X2,fill = Phylum,size =mean ),pch = 21, data = plotcord) +geom_point(aes(X1, X2),pch = 21, data = groupdata,size = 5,fill = "blue",color = "black") +geom_text_repel(aes(X1, X2,label = elements ), data = groupdata) +theme_void()pggsave("4.png",p,width = 12,height = 8)map = as.data.frame(sample_data(ps_sub))map$Group2 <- rep(c("A1","A2","A3","A4","A5","A6"),3)sample_data(ps_sub) <- map# ?div_network result = div_network(ps_sub,num = 3,group = "Group2",flour = TRUE)edge = result[[1]] head(edge)# levels(edge$target) # node = result[[2]] # head(node) # # tail(node)data = result[[3]] dim(data)#----計算節點坐標 # flour參數,設置是否僅僅展示共有和特有的二分網絡result <- div_culculate(table = result[[3]],distance = 1.1,distance2 = 1.5,distance3 = 1.3,order = FALSE)edge = result[[1]] head(edge)plotdata = result[[2]] head(plotdata)groupdata <- result[[3]]# table(plotdata$elements) node = plotdata[plotdata$elements == unique(plotdata$elements), ]otu_table = as.data.frame(t(vegan_otu(ps_sub))) tax_table = as.data.frame(vegan_tax(ps_sub)) res = merge(node,tax_table,by = "row.names",all = F) dim(res) head(res) row.names(res) = res$Row.names res$Row.names = NULL plotcord = resxx = data.frame(mean =rowMeans(otu_table)) head(xx) plotcord = merge(plotcord,xx,by = "row.names",all = FALSE) head(plotcord) # plotcord$Phylum row.names(plotcord) = plotcord$Row.names plotcord$Row.names = NULL head(plotcord)p = ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2),data = edge, size = 0.3,color = "yellow") +geom_point(aes(X1, X2,fill = Phylum,size =mean ),pch = 21, data = plotcord) +geom_point(aes(X1, X2),pch = 21, data = groupdata,size = 5,fill = "blue",color = "black") +geom_text_repel(aes(X1, X2,label = elements ), data = groupdata) +theme_void() p ggsave("4.png",p,width = 12,height = 8)map = as.data.frame(sample_data(ps_sub))map = map[1:12,]# map$Group2 <- rep(c("A1","A2","A3","A4","A5","A6"),2) sample_data(ps_sub) <- mapresult = div_network(ps_sub,num = 3,group = "Group",flour = TRUE)edge = result[[1]] head(edge)# levels(edge$target) # node = result[[2]] # head(node) # # tail(node)data = result[[3]] dim(data)result <- div_culculate(table = result[[3]],distance = 1.1,distance2 = 1.5,distance3 = 1.3,order = FALSE)edge = result[[1]] head(edge)plotdata = result[[2]] head(plotdata)groupdata <- result[[3]]# table(plotdata$elements) node = plotdata[plotdata$elements == unique(plotdata$elements), ]otu_table = as.data.frame(t(vegan_otu(ps_sub))) tax_table = as.data.frame(vegan_tax(ps_sub)) res = merge(node,tax_table,by = "row.names",all = F) dim(res) head(res) row.names(res) = res$Row.names res$Row.names = NULL plotcord = resxx = data.frame(mean =rowMeans(otu_table)) head(xx) plotcord = merge(plotcord,xx,by = "row.names",all = FALSE) head(plotcord) # plotcord$Phylum row.names(plotcord) = plotcord$Row.names plotcord$Row.names = NULL head(plotcord)p = ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2),data = edge, size = 0.3,color = "yellow") +geom_point(aes(X1, X2,fill = Phylum,size =mean ),pch = 21, data = plotcord) +geom_point(aes(X1, X2),pch = 21, data = groupdata,size = 5,fill = "blue",color = "black") +geom_text_repel(aes(X1, X2,label = elements ), data = groupdata) +theme_void()p# ggsave("4.png",p,width = 12,height = 22)在線繪制 Venn 圖和二分網絡圖,點擊閱讀原文或掃描二維碼訪問吧!
往期精品(點擊圖片直達文字對應教程)
后臺回復“生信寶典福利第一波”或點擊閱讀原文獲取教程合集
總結
以上是生活随笔為你收集整理的放弃Venn-Upset-花瓣图,拥抱二分网络的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 中科院动物所团队揭示了空间种子传播集合网
- 下一篇: 使用Mfuzz包进行基因表达的时间趋势分