原文鏈接：http://tecdat.cn/?p=23891

原文出處：拓端數(shù)據(jù)部落公眾號

可以使用環(huán)狀圖形展示基因數(shù)據(jù)比較?？梢蕴砑佣喾N圖展信息，如熱圖、散點(diǎn)圖等。

本文目標(biāo):

• 可視化基因組數(shù)據(jù)

制作環(huán)形熱圖

環(huán)形熱圖很漂亮。可以通過R來實現(xiàn)環(huán)形熱圖。

首先，讓我們生成一個隨機(jī)矩陣，并將其隨機(jī)分成五組。

1. 


2. mat1 = rbind(cbind(matrix(rnorm(50*5, mean = 1), nr = 50),

3. matrix(rnorm(50*5, mean = -1), nr = 50)),

4. cbind(matrix(rnorm(50*5, mean = -1), nr = 50),

5. matrix(rnorm(50*5, mean = 1), nr = 50))

6. )

下面的圖是熱圖的正常布局。
?

Heatmap(mat1, row_split = split)

在接下來的章節(jié)中，我將演示如何將其可視化。

輸入數(shù)據(jù)

heatmap()的輸入應(yīng)該是一個矩陣（或者一個將被轉(zhuǎn)換為單列矩陣的向量）。如果矩陣被分割成組，必須用split參數(shù)指定一個分類變量。注意spilt的值應(yīng)該是一個字符向量或一個因子。如果它是一個數(shù)字向量，它將被轉(zhuǎn)換為字符。

顏色是矩陣中數(shù)值的重要美學(xué)映射。用戶必須用用戶定義的顏色模式指定col參數(shù)。如果矩陣是連續(xù)數(shù)字，如果矩陣是字符，col的值應(yīng)該是一個命名的顏色向量。

下面的圖是之前熱圖的圓形版本。請注意，矩陣的行沿圓形方向分布，矩陣的列沿徑向方向分布。在下面的圖中，圓形被分成五個部分，每個部分對應(yīng)一個行組。

heatmap(mat1col_fun1)

有一件事非常重要，那就是在創(chuàng)建圓形熱圖之后，你必須完全刪除布局。

如果沒有指定split，就只有一個大的扇區(qū)包含完整的熱圖。

環(huán)形布局

與生成的其他圓形圖類似，環(huán)形布局可以在制作圖之前由par()控制。

熱圖軌道的參數(shù)可以在circos()函數(shù)中控制，如track.height（軌道的高度）和bg.border（軌道的邊界）。

在下面的例子中，通過設(shè)置show.sector.labels參數(shù)，增加了扇區(qū)的標(biāo)簽。扇區(qū)的順序是c("a", "b", "c", "d", "e")，按時鐘方向排列。你可以在下面的圖中看到，a扇區(qū)從?θ=90°θ=90°開始。

1. heatmap(

2. bg.border )

如果split參數(shù)的值是一個因子，那么因子水平的順序控制熱圖的順序。如果split是一個簡單的向量，熱圖的順序是unique(split)。?

1. # 注意，因為在前一個圖中調(diào)用了 circos.clear() 。

2. # 現(xiàn)在布局從theta = 0開始（第一個扇區(qū)是'e'）。

3. heatmap( levels = c("e", "d", "c", "b", "a))

樹狀圖和行名

默認(rèn)情況下，數(shù)字矩陣是按行聚類的，因此，有聚類產(chǎn)生的樹狀圖。side參數(shù)控制樹狀圖相對于熱圖軌道的位置。注意，樹枝圖是在一個分離的軌道上。

heatmap(dend.side = "inside")

?

樹狀圖的高度是由dend.track.height參數(shù)控制的。

矩陣的行名可以通過設(shè)置rownames.side參數(shù)來繪制。行名也會被繪制在一個分離的軌道中。

heatmap(rownames.side = "inside")

?

矩陣的行名和樹狀圖可以同時繪制。當(dāng)然，它們不能在熱圖軌道的同一側(cè)。

1. dend.side = "inside",

2. rownames.side = "outside"

?

行名的圖形參數(shù)可以設(shè)置為標(biāo)量或向量，長度與矩陣中的行數(shù)相同。

heatmap(col = col_fun1, rownames.side = "outside")

樹狀圖的圖形參數(shù)可以通過回調(diào)函數(shù)直接渲染樹狀圖來設(shè)置，這一點(diǎn)將在后面演示。

聚類

默認(rèn)情況下，數(shù)字矩陣是按行聚類的。 cluster參數(shù)可以設(shè)置為FALSE來關(guān)閉聚類。

當(dāng)然，當(dāng)cluster被設(shè)置為FALSE時，即使dend.side被設(shè)置，也不會繪制樹狀圖。

聚類方法和距離方法由clustering.method和distance.method參數(shù)控制。

請注意heatmap()不直接支持對矩陣列的聚類。 你應(yīng)該在使用heatmap()之前應(yīng)用列的重新排序，例如。

hclust(dist(t(mat1)))$order

對樹狀圖的回調(diào)

聚類產(chǎn)生樹狀圖?；卣{(diào)函數(shù)可以在每個樹狀圖生成后應(yīng)用于相應(yīng)的類。回調(diào)函數(shù)可以編輯樹狀圖，例如：1.重新排列樹狀圖，或者2.給樹狀圖著色。

在circos.heatmap()中，一個用戶定義的函數(shù)應(yīng)該被設(shè)置為callback參數(shù)。該用戶定義的函數(shù)應(yīng)該有三個參數(shù)。

• dend: 當(dāng)前扇區(qū)的樹狀圖。

• m: 與當(dāng)前扇區(qū)相對應(yīng)的子矩陣。

• si: 當(dāng)前扇區(qū)的扇區(qū)索引（或扇區(qū)名稱）。

默認(rèn)的回調(diào)函數(shù)定義如下，它通過對矩陣行均值加權(quán)來重新排列樹狀圖。

reorder(dend, rowMeans(m))

下面的例子通過dendsort()對每個扇區(qū)的樹狀圖重新排序。?

1. heatmap( col = col_fun1, dend.side = "inside",

2. dendsort(dend)

3. }

我們可以使用color()來渲染樹狀圖的邊緣。 例如，為五個區(qū)的樹枝圖分配不同的顏色。這里，樹枝圖軌道的高度由height參數(shù)增加。

1. 


2. den = function(dend, m, si) {

3. # 當(dāng)k = 1時，它為整個樹狀圖渲染一種相同的顏色

4. color_branches(dend, k = 1, col = dend_col[si])

5. 


或者如果矩陣沒有被分割，我們可以給子樹狀圖分配不同的顏色。?

1. 


2. color_branches(dend, k = 4, col = 2:5)

多個熱圖軌跡

如果你制作的環(huán)狀圖只包含一個熱圖軌跡，使用heatmap()是非常簡單的。如果你制作一個包含多個軌道的更復(fù)雜的環(huán)狀圖，你應(yīng)該了解關(guān)于heatmap()的更多細(xì)節(jié)。

heatmap()的第一次調(diào)用實際上是初始化布局，即應(yīng)用聚類和拆分矩陣。樹狀圖和分割變量是內(nèi)部存儲的。這就是為什么你應(yīng)該明確地調(diào)用clear()來刪除所有的內(nèi)部變量，這樣可以確保當(dāng)你制作一個新的圓形熱圖時，heatmap()的第一次調(diào)用是在一個新的環(huán)境中。

heatmap()的第一次調(diào)用決定了所有軌道的行順序（循環(huán)方向的順序），因此，接下來的軌道中的矩陣共享與第一個軌道中相同的行順序。另外，后面軌道中的矩陣也會根據(jù)第一個heatmap軌道中的分割情況進(jìn)行分割。

如果在第一個熱圖軌道中沒有應(yīng)用聚類，則使用行的自然排序（即c（1，2，...，n））。

1. 


2. mat1[sample(100, 100), ] # 按行隨機(jī)排列 mat1

3. heatmap(mat1, ?split, col_fun1, dend.side = "outside")

4. 


?

如果我切換兩個軌道，你可以看到現(xiàn)在的聚類是由第一個熱圖軌道控制的，也就是藍(lán)-紅熱圖軌道。

heatmap(mat1, col = col_fun2)

?

你可能想問，如果我不希望聚類是由第一個軌道決定的，而第二個或第三個軌道呢？解決辦法很簡單。實際上，初始化可以通過明確調(diào)用initialize()函數(shù)來手動完成。

在initialize()中，你指定你想應(yīng)用聚類的任何矩陣以及分割變量，然后，下面的heatmap()調(diào)用都共享這個布局。

在下面的例子中，全局布局是由mat1決定的，它在第二個軌道中被可視化。我在第一個軌道中設(shè)置了side = "outside"，實際上你可以發(fā)現(xiàn)樹狀圖實際上是根據(jù)第二個軌道中的矩陣生成的。

circos.heatmap.initialize(mat1, split = split)

在下一個例子中，熱圖布局是由mat1生成的，而兩個熱圖軌道分別只包含五列。?

initialize(mat1, split = split)

與其他軌道整合

其他非熱圖軌道整合。在環(huán)形布局中，x軸和y軸上的值只是數(shù)字索引。假設(shè)在一個扇形區(qū)域內(nèi)有nr行和nc列的熱圖，熱圖行的繪制間隔為（0，1），c（1，2），...，c（nr-1，nr），熱圖列也類似。同時，原始矩陣也被重新排序。如果增加更多的軌道，需要考慮所有這些影響，以確保與熱圖軌道有正確的對應(yīng)關(guān)系。

熱圖布局完成后，軌道/扇區(qū)/單元的額外信息可以通過特殊變量CELL_META來檢索。單元/扇區(qū)的附加元數(shù)據(jù)列舉如下，它們對于正確對應(yīng)熱圖軌道非常重要。

• CELL_META$row_dend或簡稱CELL_META$dend：當(dāng)前扇區(qū)的樹狀圖。如果沒有進(jìn)行聚類，則該值為NULL。

• CELL_META$row_order或簡稱CELL_META$order：聚類后當(dāng)前扇區(qū)中子矩陣的行排序。如果沒有進(jìn)行聚類，其值為c(1, 2, ..., )。

• CELL_META$subset。原始完整矩陣中指數(shù)的子集。這些值的排序是遞增的。

以下是CELL_META$row_dend、CELL_META$row_order和CELL_META$subset在循環(huán)熱圖例子中的第一個扇區(qū)的輸出。

1. row_order

2. 


3. subset

在下面的例子中，我添加了一個軌跡，它將mat1中前五列的行平均值可視化。我添加了cell.padding = c(0.02, 0, 0.02, 0)，這樣最大和最小的點(diǎn)就不會與單元格的上下邊界重疊了。

1. track(ylim = range(row_mean{

2. lines(CELL_META$cell.xlim, c(0, 0),

3. points(seq_along(y) - 0.5, y, )

同樣地，如果把點(diǎn)數(shù)軌道作為第一條軌道，則應(yīng)事先對布局進(jìn)行初始化。

1. initialize(mat1, split = split)

2. # 這與前面的例子相同

3. track(ylim = range(row_mean), panel.fun = function(x, y) {

4. circose(y > 0, "red", "blue"))

5. }, cell.padding = c(0.02, 0, 0.02, 0))

6. # 不需要在這里指定 "分割"。

7. 


Boxplots被用來對應(yīng)矩陣的行。

1. circos.heatmap(mat1, split = split, col = col_fun1)

2. circos.track(ylim = range(mat1), panel.fun = function(x, y) {

3. m = mat1[CELL_META$subset, 1:5, drop = FALSE]

4. m = m[CELL_META$row_order, , drop = FALSE]

5. n = nrow(m)

6. # circos.boxplot is applied on matrix columns, so here we transpose it.

7. circos.boxplot(t(m), pos = 1:n - 0.5, pch = 16, cex = 0.3)

8. circos.lines(CELL_META$cell.xlim, c(0, 0), lty = 2, col = "grey")

9. }, cell.padding = c(0.02, 0, 0.02, 0))

添加注釋

可以通過設(shè)置abels = TRUE來添加區(qū)塊的標(biāo)簽，然而，這并不提供對標(biāo)簽的任何定制。用戶可以通過定義fun函數(shù)來定制自己的標(biāo)簽，演示如下。

1. heatmap(col = col_fun1)

2. circos.track

3. adj = c(0.5, 0), niceFacing = TRUE)

heatmap()不直接支持矩陣的列名，但也可以通過自行定義fun函數(shù)輕松添加。在下面的例子中，我通過par()中的after參數(shù)在最后一個扇區(qū)（第五扇區(qū)）后設(shè)置了較大的空間（10度，用戶通常需要嘗試幾個值來獲得最佳空間），之后我在fun中繪制了最后一個扇區(qū)中的列名。?

1. par(gap.after = c(2, 2, 2, 2, 10))

2. heatmap(mat1, split , col = col_fun1)

3. track(track.index = 1

4. 


5. }, bg.border = NA)

下一個例子添加了矩形和標(biāo)簽來顯示矩陣中的兩組列。fun里面的代碼很簡單。convert_x()將x方向上的單位轉(zhuǎn)換為環(huán)形坐標(biāo)系中測量的適當(dāng)數(shù)值。?

1. par(gap.after = c(2, 2, 2, 2, 10))

2. heatmap(mat1, split , col = col_fun1)

3. track(track.index = 1

4. circos.text(cell.xlim[2] + convert_x(3, "mm"), 7.5,

5. }, bg.border = NA)

circlize不生成圖例，但是圖例可以由Legend()函數(shù)手動生成并添加到圓形圖中。下面是一個添加圖例的簡單例子。在下一節(jié)中，你可以找到一個添加許多圖例的更復(fù)雜的例子。

1. heatmap(mat1, split = split)

2. clear()

3. grid.draw(lgd)

一個復(fù)雜的圓形熱圖的例子

在本節(jié)中，我將演示如何制作復(fù)雜的圓形熱圖。下圖是正常布局的熱圖，現(xiàn)在我將用圓形布局改變它們。

熱圖直觀地顯示了DNA甲基化、基因表達(dá)和其他基因組水平信息之間的相關(guān)性。

原始熱圖是用隨機(jī)數(shù)據(jù)集生成的。

與原始熱圖類似，通過對甲基化矩陣（mat_meth）的行進(jìn)行k-means聚類，將所有熱圖的行分成5組。

km = kmeans(mat_meth, centers = 5)$cluster

現(xiàn)在有以下矩陣/向量需要被可視化為熱圖。

• mat：一個矩陣，其中各行對應(yīng)不同的甲基化區(qū)域（DMRs）。矩陣中的值是每個樣本中DMR的平均甲基化水平。

• expr：一個矩陣，其中的行對應(yīng)于與DMR相關(guān)的基因（即與DMR最近的基因）。矩陣中的值是每個樣本中每個基因的表達(dá)水平。每個基因在不同樣本中的表達(dá)量是有比例的。

• direction：甲基化變化的方向（hyper表示腫瘤樣本中甲基化程度較高，hypo表示腫瘤樣本中甲基化程度較低）。

• pvalue：甲基化和相關(guān)基因表達(dá)之間的相關(guān)測試的P值。數(shù)值是-log10轉(zhuǎn)換的。

• type：基因的類型（如蛋白質(zhì)編碼基因或林肯RNAs）。

• gene：對基因模型的注釋（即基因間、基因內(nèi)或轉(zhuǎn)錄起始位置（TSS））。

• dist：從DMRs到被鑒定基因的TSS的距離。

• enhancer：每個DMR中與增強(qiáng)器重疊的部分。

在這些變量中，mat_meth、mat_expr、cor_pvalue、dist和anno_enhancer是數(shù)字變量，我為它們設(shè)置了顏色映射函數(shù)。對于其他變量，我設(shè)置了命名的顏色向量。

在下面的代碼中，我在heatmap()的第一次調(diào)用中指定了分裂，這是甲基化熱圖。軌道的高度是手動調(diào)整的。

1. col_meth = colorRamp2(c(0, 0.5, 1), c("blue", "white", "red"))

2. heatmap(mat, split col_meth, track.height = 0.12)

循環(huán)熱圖看起來很美! 由于矩陣中的行是基因組區(qū)域（差異甲基化區(qū)域），如果我們能在一些區(qū)域之間建立聯(lián)系，例如三維染色體結(jié)構(gòu)中的物理相互作用，那么這個圖就會更漂亮、更有用。

在下面的代碼中，我在DMRs之間生成一些隨機(jī)的相互作用。df_link中的每一行意味著有一個從第i個DMR到第j個DMR的互動。

1. data.frame(

2. from_index

3. to_index

在圓形熱圖上找到這些DMR的位置是很棘手的。請看下面代碼中的注釋。注意這里的子集和行序元數(shù)據(jù)是通過get.data()函數(shù)明確指定扇形索引來檢索的。

1. 


2. for(i in seq_len(nrow)) {

3. # 讓我們把索引為df_link$from_index[i]的DMR稱為DMR1。

4. # 另一個索引為df_link$to_index[i]的為DMR2。

5. 


6. # DMR1所在的扇區(qū)。

7. group1 = km[from_index[i] ]。

8. # DMR2所在的扇區(qū)。

9. group2 = km[to_index[i] ]。

10. 


11. # 區(qū)塊`group1`中的DMRs子集（來自mat_meth的行指數(shù)）。

12. data("subset", sector.index = group1)

13. # 扇區(qū)`group1`中的行排序。

14. data("row_order", sector.index = group1)

15. # 這是DMR1在`group1`熱圖中的位置。

16. which(subset1[row_order1] == from_index[i])

17. 


18. # 區(qū)塊`group2`中的DMRs子集（來自mat_meth的行指數(shù)）。

19. data("subset", sector.index = group2)

20. # 扇區(qū)`group2`中的行排序。

21. ret.data("r sector.indexoup2)

22. # 這是DMR2在`group2`熱圖中的位置。index[i])

24. # 我們?nèi)≈虚g的點(diǎn)，在D和DMR2之間畫一個鏈接

25. link(group1, x1 - 0.5, grup2,col = rcoor(1))

我實現(xiàn)了一個函數(shù)?，F(xiàn)在繪制矩陣行之間的鏈接更簡單。

1. for(i in seq_len(nrow(df_link))) {

2. heatmap.link(from[i],

3. to_index[i],

4. color(1))

5. }

添加鏈接后，繪圖看起來更漂亮了!

圖例對于理解熱圖非常重要。

繪制圓形圖的函數(shù)只是前面代碼的一個封裝，沒有任何修改。

圖例對于理解熱圖非常重要。按照該鏈接的說明，我們需要一個繪制圓形圖的函數(shù)和一個Legends對象。

現(xiàn)在我們使用gridBase來結(jié)合基礎(chǔ)圖形和網(wǎng)格圖形。

1. 


2. 


3. h = dev.size()[2]

4. 


5. draw(lgd_list, x = circle, just = "left")

最受歡迎的見解

1.R語言動態(tài)圖可視化：如何、創(chuàng)建具有精美動畫的圖

2.R語言生存分析可視化分析

3.Python數(shù)據(jù)可視化-seaborn Iris鳶尾花數(shù)據(jù)

4.r語言對布豐投針（蒲豐投針）實驗進(jìn)行模擬和動態(tài)

5.R語言生存分析數(shù)據(jù)分析可視化案例

6.r語言數(shù)據(jù)可視化分析案例：探索brfss數(shù)據(jù)數(shù)據(jù)分析

7.R語言動態(tài)可視化：制作歷史全球平均溫度的累積動態(tài)折線圖動畫gif視頻圖

8.R語言高維數(shù)據(jù)的主成分pca、 t-SNE算法降維與可視化分析案例報告

9.python主題LDA建模和t-SNE可視化