昨天我在單細(xì)胞天地講解了使用monocle2進(jìn)行擬時(shí)序分析的方法收津，基本上跟著我的代碼走一波就可以學(xué)會(huì)了寂拆，當(dāng)然具體參數(shù)理解需要自行發(fā)力哦誓琼，見(jiàn)：使用monocle做擬時(shí)序分析（單細(xì)胞譜系發(fā)育） 用法只是最基礎(chǔ)的知識(shí)而已，更多的時(shí)候巍扛，我們需要活學(xué)活用领跛，比如課程學(xué)員提到的問(wèn)題，就是因?yàn)樽霾坏交顚W(xué)活用撤奸，他想知道下面的擬時(shí)序分析的熱圖提取基因吠昭，學(xué)員把基因按照發(fā)育順序繪制了熱圖，而這些基因被他分成了3組胧瓜，想拿基因去做GO/KEGG等數(shù)據(jù)庫(kù)進(jìn)行功能注釋矢棚，不知道如何獲取基因名字。
我這里不能拿學(xué)員真實(shí)項(xiàng)目數(shù)據(jù)來(lái)演示府喳，所以還是用我們的老朋友蒲肋，拿scRNAseq包的表達(dá)矩陣測(cè)試，見(jiàn)：使用monocle做擬時(shí)序分析（單細(xì)胞譜系發(fā)育）
首先根據(jù)細(xì)胞發(fā)育譜系來(lái)繪制熱圖
因?yàn)榍懊娴慕坛?使用monocle做擬時(shí)序分析（單細(xì)胞譜系發(fā)育） 我們已經(jīng)把細(xì)胞發(fā)育情況做出來(lái)了钝满，就是NPC細(xì)胞跟另外3種細(xì)胞從生理上就不一樣兜粘，所以是單獨(dú)的發(fā)育軌跡，而 “GW16” and “GW21” 舱沧，“GW21+3” 這種孕期細(xì)胞妹沙，就可以很清晰的看到時(shí)間被反映在我們的擬時(shí)序分析結(jié)果了。
大家可以重溫教程 使用monocle做擬時(shí)序分析（單細(xì)胞譜系發(fā)育） 里面的4個(gè)繪圖代碼：
plot_cell_trajectory(cds, color_by = "Biological_Condition")
# 可以很明顯看到細(xì)胞的發(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)建細(xì)胞發(fā)育譜系熟吏，使用的是不同Biological_Condition的細(xì)胞類型之間使用monocle找到的兩千多個(gè)基因。
簡(jiǎn)單的一個(gè)函數(shù)就可以繪制熱圖：
plot_pseudotime_heatmap(cds[ordering_genes,],
num_clusters = 3,
cores = 1,
show_rownames = T)
可以看到玄窝， 這個(gè)圖就和學(xué)員求助的圖一模一樣啦牵寺，因?yàn)榛驍?shù)量?jī)汕Ф鄠€(gè)，所以畫(huà)出來(lái)肯定是看不清晰的啦恩脂。
既然基因被分成了3組帽氓，想拿基因去做GO/KEGG等數(shù)據(jù)庫(kù)進(jìn)行功能注釋，就需要獲取基因名字俩块。
這個(gè)做不出來(lái)黎休，不怪學(xué)員，因?yàn)檎Ｈ撕茈y想到玉凯，這個(gè)繪圖函數(shù)其實(shí)是可以調(diào)整返回?cái)?shù)據(jù)對(duì)象的势腮，而且官網(wǎng)例子也沒(méi)有提到。
需要看函數(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.
很明顯捎拯，這個(gè)函數(shù)其實(shí)就是pheatmap的一個(gè)包裝罷了，本質(zhì)上也是調(diào)用 hclust 而已盲厌，使用的是ward.D2距離署照。
然后解析熱圖函數(shù)返回對(duì)象
根據(jù)幫助文檔祸泪，我們修改參數(shù)，這樣monocle的plot_pseudotime_heatmap函數(shù)就有返回值了建芙，是一個(gè)對(duì)象没隘。
p=plot_pseudotime_heatmap(cds[ordering_genes,],
num_clusters = 3,
cores = 1,return_heatmap=T,
show_rownames = T)
從pheatmap的對(duì)象里面提取基因名字就很簡(jiǎn)單了，就是在p$tree_row里面
> p$tree_row
Call:
hclust(d = d, method = method)
Cluster method : ward.D2
Number of objects: 2200
就可以拿到基因名對(duì)應(yīng)的cluster啦禁荸，代碼如下：
clusters <- cutree(p$tree_row, k = 3)
clustering <- data.frame(clusters)
clustering[,1] <- as.character(clustering[,1])
colnames(clustering) <- "Gene_Clusters"
table(clustering)
