【原】擬時序分析的熱圖提取基因問題

健明 2021-07-15

展開全文

昨天我在單細胞天地講解了使用monocle2進行擬時序分析的方法，基本上跟著我的代碼走一波就可以學會了,，當然具體參數(shù)理解需要自行發(fā)力哦,，見：使用monocle做擬時序分析（單細胞譜系發(fā)育）用法只是最基礎(chǔ)的知識而已，更多的時候,，我們需要活學活用,，比如課程學員提到的問題，就是因為做不到活學活用,，他想知道下面的擬時序分析的熱圖提取基因,，學員把基因按照發(fā)育順序繪制了熱圖，而這些基因被他分成了3組,，想拿基因去做GO/KEGG等數(shù)據(jù)庫進行功能注釋,，不知道如何獲取基因名字。

我這里不能拿學員真實項目數(shù)據(jù)來演示,，所以還是用我們的老朋友,，拿scRNAseq包的表達矩陣測試，見：使用monocle做擬時序分析（單細胞譜系發(fā)育）

首先根據(jù)細胞發(fā)育譜系來繪制熱圖

因為前面的教程使用monocle做擬時序分析（單細胞譜系發(fā)育）我們已經(jīng)把細胞發(fā)育情況做出來了,，就是NPC細胞跟另外3種細胞從生理上就不一樣,，所以是單獨的發(fā)育軌跡，而 “GW16” and “GW21” ,，“GW21+3” 這種孕期細胞,，就可以很清晰的看到時間被反映在我們的擬時序分析結(jié)果了。

大家可以重溫教程使用monocle做擬時序分析（單細胞譜系發(fā)育）里面的4個繪圖代碼：

plot_cell_trajectory(cds, color_by = "Biological_Condition")  
# 可以很明顯看到細胞的發(fā)育軌跡 
plot_cell_trajectory(cds, color_by = "State")
plot_cell_trajectory(cds, color_by = "Pseudotime")
plot_cell_trajectory(cds, color_by = "State") +
  facet_wrap(~State, nrow = 1)

我們前面構(gòu)建細胞發(fā)育譜系,，使用的是不同Biological_Condition的細胞類型之間使用monocle找到的兩千多個基因,。

簡單的一個函數(shù)就可以繪制熱圖：

plot_pseudotime_heatmap(cds[ordering_genes,],
                        num_clusters = 3,
                        cores = 1,
                        show_rownames = T)

可以看到，這個圖就和學員求助的圖一模一樣啦,，因為基因數(shù)量兩千多個,，所以畫出來肯定是看不清晰的啦。

既然基因被分成了3組,，想拿基因去做GO/KEGG等數(shù)據(jù)庫進行功能注釋,，就需要獲取基因名字。

這個做不出來,，不怪學員,，因為正常人很難想到，這個繪圖函數(shù)其實是可以調(diào)整返回數(shù)據(jù)對象的,，而且官網(wǎng)例子也沒有提到,。

需要看函數(shù)的幫助文檔，如下：

hclust_method    
The method used by pheatmap to perform hirearchical clustering of the rows.
hclust_method = "ward.D2"
return_heatmap    
Whether to return the pheatmap object to the user.

很明顯,，這個函數(shù)其實就是pheatmap的一個包裝罷了,，本質(zhì)上也是調(diào)用 hclust 而已，使用的是ward.D2距離,。

然后解析熱圖函數(shù)返回對象

根據(jù)幫助文檔,，我們修改參數(shù)，這樣monocle的plot_pseudotime_heatmap函數(shù)就有返回值了,，是一個對象,。

p=plot_pseudotime_heatmap(cds[ordering_genes,],
                        num_clusters = 3,
                        cores = 1,return_heatmap=T,
                        show_rownames = T)

從pheatmap的對象里面提取基因名字就很簡單了，就是在p$tree_row里面

> p$tree_row

Call:
hclust(d = d, method = method)

Cluster method   : ward.D2 
Number of objects: 2200

就可以拿到基因名對應的cluster啦,，代碼如下：

clusters <- cutree(p$tree_row, k = 3)
clustering <- data.frame(clusters)
clustering[,1] <- as.character(clustering[,1])
colnames(clustering) <- "Gene_Clusters"
table(clustering)