結合MATLAB、Python、R語言，在求得顯著差異的邊（節點對）之後，怎麼畫circle圖
阿新 • • 發佈：2021-01-12
先來看看成果圖：
OK，開始畫圖：
實驗背景宣告：在腦影像分析中，我們首先構建腦網路，然後使用雙樣本t對比兩組人的連線差異，然後使用以上的圖進行視覺化，一般紅色連線代表顯著升高，綠色代表顯著下降。（非必須，根據實際需求設計，如上圖中紅色代表相應的連線差異與HAMD抑鬱量表評分顯著相關，綠色表示不相關）。這裡呢，我們研究了一組病人以及年齡性別匹配的健康被試的fMRI的資料，首先進行fMRI預處理，然後基於中科院246模板提取ROISignals，然後構建腦網路（細節不做介紹）。最後進行連線差異對比。然後得到病人組相對於正常被試組有哪些節點之間的連線權重出現了顯著降低或顯著升高。
Step1: R語言的畫圖程式碼MyCircle.R，這裡我們先給出程式碼，根據程式碼的檔案需求，準備材料。在以下程式碼中，開頭即指明瞭需要的材料，有
LobeID.txt；
GyrusID.txt；
ROIID.txt；
Links.txt
LobeID <- read.table("./Circle edges/LobeID.txt", header = TRUE, sep = "\t")
head(LobeID)
GeneID <- read.table("./Circle edges/GyrusID.txt", header = TRUE, sep = "\t")
head(GeneID)
GeneFe <- read.table("./Circle edges/ROIID.txt", header = TRUE, sep = "\t")
head(GeneFe)
Links <- read.table("./Circle edges/Links.txt", header = TRUE, sep = "\t")
head(Links)
LobeID$LobeID <- factor(LobeID$LobeID, levels = LobeID$LobeID)
library(stringi)
library(circlize)
pdf("circlize1.pdf",width = 17.5,height = 10)
par(mar=rep(0,4))
circos.clear()
circos.par(start.degree = 90, #從哪裡開始畫，沿著逆時針順序
gap.degree = 2, #每兩個sectors間的間隔大小，此處是每兩個組間的間隔
track.margin = c(0,0.08), #gap.degree規定了左右，track.margin規定上下
cell.padding = c(0,0,0,0) #？
)
#### 根據"LobeID.txt"畫組
#圖是多個組圍成的一個圈，其中每個組及其亞組都畫在一個單獨的座標系裡
circos.initialize(factors = LobeID$LobeID,
xlim = cbind(LobeID$LobeStart, LobeID$LobeEnd))
circos.trackPlotRegion(ylim = c(0, 1), factors = LobeID$LobeID, track.height=0.1,
#panel.fun for each sector
panel.fun = function(x, y) {
#select details of current sector
name <- get.cell.meta.data("sector.index")
i <- get.cell.meta.data("sector.numeric.index")
xlim <- get.cell.meta.data("xlim")
ylim <- get.cell.meta.data("ylim")
#組ID文字的設定
circos.text(x = mean(xlim),
y = 1,
labels = name,
cex = 1.7, #組ID文字大小
#col = "darkgrey", #組ID的顏色
adj = c(0.5,-0.8)
#如果想讓組ID呈放射狀排列，就用下面三行代替上面那行
#niceFacing = TRUE, #組ID頭朝上
#facing = "reverse.clockwise",
#adj = c(1.5, -0.8) #組ID的位置
)
#畫組
circos.rect(xleft = xlim[1],
ybottom = ylim[1],
xright = xlim[2],
ytop = ylim[2],
col = "blue",
border = "black")
})
#### 根據"GyrusID.txt"畫亞組
library(graphics)
#bgcol <- rainbow(sum(LobeID$LobeEnd), s = 1, v = 1, start = 0, end = max(1, sum(LobeID$LobeEnd) - 1)/sum(LobeID$LobeEnd), alpha = 1)
for (i in seq_len(nrow(GeneID))){
ylim <- c(0, 1)
circos.rect(sector.index = GeneID$LobeID[i],
track.index = 1,
xleft = GeneID$featureStart[i], # + 0.01, #讓feature之間流出空隙
ybottom = ylim[1],
xright = GeneID$featureEnd[i], # - 0.01, #讓feature之間流出空隙
ytop = ylim[2],
#如果想讓每個方塊顏色不同，用彩虹色畫每一個Gene feature就執行下面這行
#col = bgcol[i],
#如果想讓每個組用同一個顏色，呼叫輸入檔案中寫好的顏色，就執行下面這行
col = "white",
border = paste("#", GeneID$barColor[i], sep = ""), lwd = 8)
}
# 亞組的文字
for (i in seq_len(nrow(GeneID))){
ylim <- c(0, 1)
circos.text((GeneID$featureStart[i] + GeneID$featureEnd[i])/2,
ylim[1] + 0.5,
GeneID$GyrusID[i],
sector.index = GeneID$LobeID[i],
facing = "inside",
cex = 0.9)
}
circos.trackPlotRegion(ylim = c(0, 1), factors = LobeID$LobeID, track.height=0.1,
#panel.fun for each sector
panel.fun = function(x, y) {
#select details of current sector
name <- get.cell.meta.data("sector.index")
i <- get.cell.meta.data("sector.numeric.index")
xlim <- get.cell.meta.data("xlim")
ylim <- get.cell.meta.data("ylim")
#畫組
circos.rect(xleft = xlim[1],
ybottom = ylim[1],
xright = xlim[2],
ytop = ylim[2],
col = "blue",
border = "white",lwd = 6)
})
#提需求的小夥伴還想讓圓圈圖中每個方塊顏色不同，整體顯示彩虹色，就像circlize_rainbow.pdf那樣的效果，就把下面這個for迴圈裡的“col = ...”行換成“col = bgcol[i], ”
#執行下面兩行，設定彩虹色
library(graphics)
bgcol <- rainbow(sum(LobeID$LobeEnd), s = 1, v = 1, start = 0, end = max(1, sum(LobeID$LobeEnd) - 1)/sum(LobeID$LobeEnd), alpha = 1)
for (i in seq_len(nrow(GeneFe))){
ylim <- c(0, 1)
circos.rect(sector.index = GeneFe$LobeID[i],
track.index = 2,
xleft = GeneFe$featureStart[i], # + 0.01, #讓feature之間流出空隙
ybottom = ylim[1],
xright = GeneFe$featureEnd[i], # - 0.01, #讓feature之間流出空隙
ytop = ylim[2],
#如果想讓每個方塊顏色不同，用彩虹色畫每一個Gene feature就執行下面這行
#col = bgcol[i],
#如果想讓每個組用同一個顏色，呼叫輸入檔案中寫好的顏色，就執行下面這行
col = paste("#", GeneFe$barColor[i], sep = ""),
border = "gray")
}
# 亞組的文字
for (i in seq_len(nrow(GeneFe))){
ylim <- c(0, 1)
circos.text(GeneFe$featureStart[i] + 0.5,
ylim[1] + 0.5,
GeneFe$ROIID[i],
sector.index = GeneFe$LobeID[i],
facing = "inside",
cex = 0.8)
}
#### 根據"Links.txt"畫連線
for(i in seq_len(nrow(Links))){
circos.link(sector.index1 = Links$ROI_1[i],
point1 = Links[i, 2],
sector.index2 = Links$ROI_2[i],
point2 = Links[i, 4],
#在輸入檔案中已經寫好了連線的顏色
col = paste("#", Links$link_color[i], sep = ""),
#如果不想在輸入檔案中寫顏色，還可以執行下面這行，在這裡指定顏色
#col = "red",
border = FALSE,
lwd = 2, #連線的粗細
rou = 0.64, #連線起始和終止點的y值（圓半徑的百分比）
h.ratio = 0.7 #連線的拱形高度，數值越大，拱形越低
#h.ratio = Links$link_height[i] #在輸入檔案中指定拱形高度
)
}
#Legend of brain regions
GeneCol<-data.frame(legendID=GeneID$GyrusID,legendCol=GeneID$barColor)
GeneCol_uniq<-base::unique(GeneCol)
legend(1,0.2,
bty = "n",#不要邊框
pch = 16,#畫成正方形，圓點是16
cex = 1.2,
legend = c("SFG:Superior_Frontal_Gyrus","PrG:Precentral_Gyrus","PCL:Paracentral_Lobule",
"STG:Superior_Temporal_Gyrus","MTG:Middle_Temporal_Gyrus","FuG:Fusiform_Gyrus",
"SPL:Superior_Parietal_Lobule","IPL:Inferior_Parietal_Lobuleb","PoG:Postcentral_Gyrus",
"MVOcC:Medio_Ventral_Occipital_Cortex","LOcC:lateral_Occipital_Cortex"),
col = paste0("#",GeneCol_uniq$legendCol),
horiz = FALSE)
#Legend of links
LinksCol<-data.frame(Links$link_color)
LinksCol_uniq<-base::unique(LinksCol)
legend(1,0.4,
bty = "n",
lty = 1,#畫線
lwd = 5,#線的粗細
legend = c("Not Significantly correlated with HAMD","Significantly correlated with HAMD"),
col = paste0("#",LinksCol_uniq$Links.link_col),
horiz = FALSE)
dev.off()
Step2: 準備材料
以上提到的幾個檔案，內容是這個樣子的：
LobeID.txt；
GyrusID.txt；
ROIID.txt；
Links.txt
根據以上的內容就可以正常使用畫圖了，如果你不會寫程式碼，就自己一點一點的敲文字就好了。
如果熟悉基本的文字處理，接下來用程式碼輔助製作材料就可以啦，這裡我使用MATLAB+python輔助實現，都是很簡單的指令碼，
可以提高準確性，提高效率，給大家參考。一下程式碼不是最好的，跟隨著作者的意識流生成的，可優化空間很大，但是一定能達到目標：
Base材料：
edge_for_circles: 代表節點對，及具有顯著差異的邊
name_246：代表中科院246模板的子節點名字（label）---->http://atlas.brainnetome.org/
衍生材料：
根據以上矩陣，使用以下matlab程式碼製作衍生材料，然後將edge_all_s儲存為mat檔案，命名為"edge_num.mat"，如下圖：
edge_all = unique(edge_for_cirles); %得到獨一無二的腦區代號 edge_all_s = sort(edge_all); %對以上腦區編號排個序 edge_s = sortrows(edge_for_cirles); %對節點對進行擴充套件排序，按照起點的排序進行排序 name_all = name_246(edge_all_s,2); %得到節點名的mat矩陣
（1）LobeID.txt製作過程：
我們需要知道涉及的Lobe有哪些，而上面的程式碼可以知道所有的子區，根據子區的number，定位Lobe，這裡使用了一個python的指令碼把相應的名字讀出來了，程式碼中出現的
full_nameOfROIs.csv檔案，以及所有的程式碼檔案皆可以在本人github網站下載：https://github.com/chaochaobar/circles_for_brain
import scipy.io as sio
import pandas as pd
datafile = './edge_num.mat'
mat = sio.loadmat(datafile)
full_name = pd.read_csv("./full_nameOfROIs.csv")
index_global = mat['edge_all_s']
save_file_name = "for_circle.node"
i = 0
for line in open('./node246.node'):
i = i + 1
a = i // 2
b = i % 2
if b == 0:
app = full_name["ROI_R"][a - 1]
elif b == 1:
app = full_name["ROI_L"][a]
if i in index_global:
if i in range(0, 69):
Lobe = "Frontal_Lobe"
if i in range(0, 15):
Gyrus = "Superior_Frontal_Gyrus"
elif i in range(15, 29):
Gyrus = "Middle_Frontal_Gyrus"
elif i in range(29, 41):
Gyrus = "Inferier_Frontal_Gyrus"
elif i in range(41, 53):
Gyrus = "Orbital_Gyrus"
elif i in range(53, 65):
Gyrus = "Precentral_Gyrus"
elif i in range(65, 69):
Gyrus = "Paracentral_Lobule"
elif i in range(69, 125):
Lobe = "Temporal_Lobe"
if i in range(69, 81):
Gyrus = "Superior_Temporal_Gyrus"
elif i in range(81, 89):
Gyrus = "Middle_Temporal_Gyrus"
elif i in range(89, 103):
Gyrus = "Inferior_Temporal_Gyrus"
elif i in range(103, 109):
Gyrus = "Fusiform_Gyrus"
elif i in range(109, 121):
Gyrus = "Parahippocampal_Gyrus"
elif i in range(121, 125):
Gyrus = "posterior_Superior_Temporal_Sulcus"
elif i in range(125, 163):
Lobe = "Parietal_Lobe"
if i in range(125, 135):
Gyrus = "Superior_Parietal_Lobule"
elif i in range(135, 147):
Gyrus = "Inferior_Parietal_Lobule"
elif i in range(147, 155):
Gyrus = "Precuneus"
elif i in range(155, 163):
Gyrus = "Postcentral_Gyrus"
elif i in range(163, 175):
Lobe = "Insular_Lobe"
Gyrus = "Insular_Gyrus"
elif i in range(175, 189):
Lobe = "Limbic_Lobe"
Gyrus = "Cingulate_Gyrus"
elif i in range(189, 211):
Lobe = "Occipital_Lobe"
if i in range(189, 199):
Gyrus = "Medio_Ventral_Occipital_Cortex"
elif i in range(199, 211):
Gyrus = "lateral_Occipital_Cortex"
elif i in range(211, 247):
Lobe = "Subcortical_Nuclei"
if i in range(211, 215):
Gyrus = "Amygdala"
elif i in range(215, 219):
Gyrus = "Hippocampus"
elif i in range(219, 231):
Gyrus = "Basal_Ganglia"
elif i in range(213, 247):
Gyrus = "Thalamus"
print(Lobe, ",", Gyrus, ",", line)
with open(save_file_name, "a+") as f:
f.write("{Lobe} {Gyrus} {line}".format(Lobe = Lobe, Gyrus = Gyrus, line = line))
f.close()
生成了“for_circle.node”檔案之後，使用matlab讀取，在製作資料材料ROIID.txtyiji Links.txt檔案的時候，先用以下matlab程式碼將每一列的內容生成出來，然後貼上到excel裡面，最後貼上到txt檔案裡面，有點繞是吧，這就是意識流，用最快的方法，即便有點蠻力，但是省去了過多的思考，半自動化：
%製作ROIID:先使用python讀lobe和gyrus的名字，然後執行下面的語句
[lobe,gyrus,x,y,z,t,name_sub] = textread('for_circle.node','%s %s %d %d %d %f %s');
%%% 製作LobeID, GyrusID, ROIID, LinkID,使用excel
for i = 1:len_sig
index_start = find(edge_all==edge_s(i, 1));
index_end = find(edge_all==edge_s(i, 2));
edges_all_name(i,1:2) = [name_all(index_start), name_all(index_end)];
end
%%% start lobe , end lobe
for i = 1:len_sig
cur_start = edges_all_name{i, 1};
index_start = find(strcmp(C,cur_start));
cur_end = edges_all_name{i,2};
index_end = find(strcmp(C, cur_end));
edges_all_name_2(i, 1) = A(index_start, 1);
edges_all_name_2(i, 2) = A(index_end, 1);
end
%%% start lobe , end lobe
for i = 1:len_sig
cur_start = edges_all_name{i, 1};
index_start = find(strcmp(C,cur_start));
cur_end = edges_all_name{i,2};
index_end = find(strcmp(C, cur_end));
edges_all_name_2(i, 1) = B(index_start, 1);
edges_all_name_2(i, 2) = B(index_end, 1);
end
%%% 製作LinkID的程式碼
point1 = zeros(len_sig, 1);
point2 = zeros(len_sig, 1);
for i = 1:len_sig
cur_start = edges_all_name{i, 1};
index_start = find(strcmp(C,cur_start));
cur_end = edges_all_name{i,2};
index_end = find(strcmp(C, cur_end));
point1(i, 1) = D(index_start, 1) + 0.5;
point2(i,1) = D(index_end, 1) + 0.5;
end
程式碼以及資料檔案製作過程都有了，動手實踐起來吧小夥伴們，有什麼問題，歡迎聯絡：[email protected]
&n