拓端tecdat|R語言繪制圈圖、環形熱圖可視化基因組實戰：展示基因數據比較

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

原文出處：拓端數據部落公眾號

可以使用環狀圖形展示基因數據比較。可以添加多種圖展信息，如熱圖、散點圖等。

本文目標:

• 可視化基因組數據

制作環形熱圖

環形熱圖很漂亮。可以通過R來實現環形熱圖。

首先，讓我們生成一個隨機矩陣，並將其隨機分成五組。

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)

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

輸入數據

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

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

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

heatmap(mat1col_fun1)

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

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

環形布局

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

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

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

1.  
   heatmap(
2.  
   bg.border )

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

1.  
   # 注意，因為在前一個圖中調用了 circos.clear() 。
2.  
   # 現在布局從theta = 0開始（第一個扇區是'e'）。
3.  
   heatmap( levels = c("e", "d", "c", "b", "a))

樹狀圖和行名

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

heatmap(dend.side = "inside")

 

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

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

heatmap(rownames.side = "inside")

 

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

1.  
   dend.side = "inside",
2.  
   rownames.side = "outside"

 

行名的圖形參數可以設置為標量或向量，長度與矩陣中的行數相同。

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

樹狀圖的圖形參數可以通過回調函數直接渲染樹狀圖來設置，這一點將在后面演示。

聚類

默認情況下，數字矩陣是按行聚類的。 cluster參數可以設置為FALSE來關閉聚類。

當然，當cluster被設置為FALSE時，即使dend.side被設置，也不會繪制樹狀圖。

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

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

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

對樹狀圖的回調

聚類產生樹狀圖。回調函數可以在每個樹狀圖生成后應用於相應的類。回調函數可以編輯樹狀圖，例如：1.重新排列樹狀圖，或者2.給樹狀圖着色。

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

• dend: 當前扇區的樹狀圖。
• m: 與當前扇區相對應的子矩陣。
• si: 當前扇區的扇區索引（或扇區名稱）。

默認的回調函數定義如下，它通過對矩陣行均值加權來重新排列樹狀圖。

reorder(dend, rowMeans(m))

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

1.  
   heatmap( col = col_fun1, dend.side = "inside",
2.  
   dendsort(dend)
3.  
   }

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

1.  
    
2.  
   den = function(dend, m, si) {
3.  
   # 當k = 1時，它為整個樹狀圖渲染一種相同的顏色
4.  
   color_branches(dend, k = 1, col = dend_col[si])
5.  
    

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

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

多個熱圖軌跡

如果你制作的環狀圖只包含一個熱圖軌跡，使用heatmap()是非常簡單的。如果你制作一個包含多個軌道的更復雜的環狀圖，你應該了解關於heatmap()的更多細節。

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

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

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

1.  
    
2.  
   mat1[sample(100, 100), ] # 按行隨機排列 mat1
3.  
   heatmap(mat1, split, col_fun1, dend.side = "outside")
4.  
    

 

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

heatmap(mat1, col = col_fun2)

 

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

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

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

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

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

initialize(mat1, split = split)

與其他軌道整合

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

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

• CELL_META$row_dend或簡稱CELL_META$dend：當前扇區的樹狀圖。如果沒有進行聚類，則該值為NULL。
• CELL_META$row_order或簡稱CELL_META$order：聚類后當前扇區中子矩陣的行排序。如果沒有進行聚類，其值為c(1, 2, ..., )。
• CELL_META$subset。原始完整矩陣中指數的子集。這些值的排序是遞增的。

以下是CELL_META$row_dend、CELL_META$row_order和CELL_META$subset在循環熱圖例子中的第一個扇區的輸出。

1.  
   row_order
2.  
    
3.  
   subset

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

1.  
   track(ylim = range(row_mean{
2.  
   lines(CELL_META$cell.xlim, c(0, 0),
3.  
   points(seq_along(y) - 0.5, y, )

同樣地，如果把點數軌道作為第一條軌道，則應事先對布局進行初始化。

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被用來對應矩陣的行。

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))

添加注釋

可以通過設置abels = TRUE來添加區塊的標簽，然而，這並不提供對標簽的任何定制。用戶可以通過定義fun函數來定制自己的標簽，演示如下。

1.  
   heatmap(col = col_fun1)
2.  
   circos.track
3.  
   adj = c(0.5, 0), niceFacing = TRUE)

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

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)

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

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()函數手動生成並添加到圓形圖中。下面是一個添加圖例的簡單例子。在下一節中，你可以找到一個添加許多圖例的更復雜的例子。

1.  
   heatmap(mat1, split = split)
2.  
   clear()
3.  
   grid.draw(lgd)

一個復雜的圓形熱圖的例子

在本節中，我將演示如何制作復雜的圓形熱圖。下圖是正常布局的熱圖，現在我將用圓形布局改變它們。

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

原始熱圖是用隨機數據集生成的。

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

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

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

• mat：一個矩陣，其中各行對應不同的甲基化區域（DMRs）。矩陣中的值是每個樣本中DMR的平均甲基化水平。
• expr：一個矩陣，其中的行對應於與DMR相關的基因（即與DMR最近的基因）。矩陣中的值是每個樣本中每個基因的表達水平。每個基因在不同樣本中的表達量是有比例的。
• direction：甲基化變化的方向（hyper表示腫瘤樣本中甲基化程度較高，hypo表示腫瘤樣本中甲基化程度較低）。
• pvalue：甲基化和相關基因表達之間的相關測試的P值。數值是-log10轉換的。
• type：基因的類型（如蛋白質編碼基因或林肯RNAs）。
• gene：對基因模型的注釋（即基因間、基因內或轉錄起始位置（TSS））。
• dist：從DMRs到被鑒定基因的TSS的距離。
• enhancer：每個DMR中與增強器重疊的部分。

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

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

1.  
   col_meth = colorRamp2(c(0, 0.5, 1), c("blue", "white", "red"))
2.  
   heatmap(mat, split col_meth, track.height = 0.12)

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

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

1.  
   data.frame(
2.  
   from_index
3.  
   to_index

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

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所在的扇區。
7.  
   group1 = km[from_index[i] ]。
8.  
   # DMR2所在的扇區。
9.  
   group2 = km[to_index[i] ]。
10.  
     
11.  
    # 區塊`group1`中的DMRs子集（來自mat_meth的行指數）。
12.  
    data("subset", sector.index = group1)
13.  
    # 扇區`group1`中的行排序。
14.  
    data("row_order", sector.index = group1)
15.  
    # 這是DMR1在`group1`熱圖中的位置。
16.  
    which(subset1[row_order1] == from_index[i])
17.  
     
18.  
    # 區塊`group2`中的DMRs子集（來自mat_meth的行指數）。
19.  
    data("subset", sector.index = group2)
20.  
    # 扇區`group2`中的行排序。
21.  
    ret.data("r sector.indexoup2)
22.  
    # 這是DMR2在`group2`熱圖中的位置。index[i])
23.  
     
24.  
    # 我們取中間的點，在D和DMR2之間畫一個鏈接
25.  
    link(group1, x1 - 0.5, grup2,col = rcoor(1))
26.  
     
27.  
     

我實現了一個函數。現在繪制矩陣行之間的鏈接更簡單。

1.  
   for(i in seq_len(nrow(df_link))) {
2.  
   heatmap.link(from[i],
3.  
   to_index[i],
4.  
   color(1))
5.  
   }

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

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

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

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

現在我們使用gridBase來結合基礎圖形和網格圖形。

1.  
    
2.  
    
3.  
   h = dev.size()[2]
4.  
    
5.  
   draw(lgd_list, x = circle, just = "left")

---

最受歡迎的見解

1.R語言動態圖可視化：如何、創建具有精美動畫的圖

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

3.Python數據可視化-seaborn Iris鳶尾花數據

4.r語言對布豐投針（蒲豐投針）實驗進行模擬和動態

5.R語言生存分析數據分析可視化案例

6.r語言數據可視化分析案例：探索brfss數據數據分析

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

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

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