關(guān)于什么是“擬時(shí)序分析”这吻，可以參考本期推送的另一篇推文吊档。這一篇直接演示代碼
monocle2這個(gè)軟件用得太多了，很多文章都是monocle2的圖橘原。因?yàn)橹皇褂帽磉_(dá)矩陣作為輸入籍铁，相比于其他軟件涡上，已經(jīng)很方便了趾断。不過話說回來，我總感覺這種軌跡推斷像玄學(xué)吩愧，大半的結(jié)果是調(diào)整出來的/事先已知的芋酌，比如擬時(shí)序里面的起始狀態(tài)經(jīng)常需要自己指定。
在做擬時(shí)序之前雁佳，最好先跑完seurat標(biāo)準(zhǔn)流程脐帝，并注釋好細(xì)胞類型同云。我這里使用的數(shù)據(jù)都已經(jīng)注釋好細(xì)胞類型，并且事先已經(jīng)知道哪一個(gè)細(xì)胞類型可能是初始狀態(tài)堵腹。
1. 導(dǎo)入數(shù)據(jù)炸站，創(chuàng)建對(duì)象，預(yù)處理
library(monocle)
library(tidyverse)
expr_matrix=read.table("count_mat.txt",header = T,row.names = 1,sep = "\t",stringsAsFactors = F)
#10X的數(shù)據(jù)使用UMI count矩陣
cell_anno=read.table("cell_anno.txt",header = T,row.names = 1,sep = "\t",stringsAsFactors = F)
gene_anno=read.table("gene_anno.txt",header = T,row.names = 1,sep = "\t",stringsAsFactors = F)
#“gene_short_name”為列名
pd=new("AnnotatedDataFrame",data = cell_anno)
fd=new("AnnotatedDataFrame",data = gene_anno)
test=newCellDataSet(as(as.matrix(expr_matrix),'sparseMatrix'),phenoData = pd,featureData = fd)
#大數(shù)據(jù)集使用稀疏矩陣疚顷，節(jié)省內(nèi)存旱易，加快運(yùn)算
test <- estimateSizeFactors(test)
test <- estimateDispersions(test)
test=detectGenes(test,min_expr = 0.1) #計(jì)算每個(gè)基因在多少細(xì)胞中表達(dá)
2. 選擇基因
選擇研究的生物學(xué)過程涉及到的基因集，這一步對(duì)于軌跡形狀的影響很大腿堤。
可以選擇數(shù)據(jù)集中的高變基因阀坏，或者是在seurat中分析得到的marker基因列表。如果是時(shí)間序列數(shù)據(jù)笆檀，可以直接比較時(shí)間點(diǎn)選差異基因忌堂。總之選擇的基因要含有盡可能多的信息酗洒。
我這里直接用的各種亞群差異基因的集合
marker_gene=read.table("seurat_marker_gene.txt",header=T,sep="\t",stringsAsFactors=F)
test_ordering_genes=unique(marker_gene$gene)
test=setOrderingFilter(test,ordering_genes = test_ordering_genes)
#指明哪些基因用于后續(xù)的聚類/排序
3. 推斷軌跡士修，并按照擬時(shí)序給細(xì)胞排序
test=reduceDimension(test,reduction_method = "DDRTree",max_components = 2, norm_method = 'log',residualModelFormulaStr = "~num_genes_expressed")
#residualModelFormulaStr減少其他因素的影響，比如不同樣本樱衷、不同批次
test=orderCells(test)
4. 繪圖
plot_cell_trajectory(test,color_by = "celltype")
ggsave("celltype.pdf",device = "pdf",width = 8,height = 9,units = c("cm"))
plot_cell_trajectory(test,color_by = "State")
ggsave("State.pdf",device = "pdf",width = 8,height = 9,units = c("cm"))
plot_cell_trajectory(test,color_by = "Pseudotime")
ggsave("Pseudotime.pdf",device = "pdf",width = 8,height = 9,units = c("cm"))
plot_cell_trajectory(test,color_by = "celltype")+facet_wrap(~celltype,nrow=1)
ggsave("celltypeb.pdf",device = "pdf",width = 21,height = 9,units = c("cm"))
有時(shí)候（大多數(shù)時(shí)候）李命，擬時(shí)序的方向或是根節(jié)點(diǎn)弄錯(cuò)了，還需要手動(dòng)更改
test=orderCells(test,root_state = 3)
plot_cell_trajectory(test,color_by = "Pseudotime")
ggsave("Pseudotime.pdf",device = "pdf",width = 8,height = 9,units = c("cm"))
5. 找差異基因
這里是指找隨擬時(shí)序變化的差異基因箫老，以及不同state之間的差異基因楚午。這兩個(gè)都是monocle里面的概念尿背，與seurat里面找的差異基因不同。
如果在做monocle2的時(shí)候，想展示這種差異基因赡麦，就需要做這一步。
expressed_genes=row.names(subset(fData(test),num_cells_expressed>=10)) #在部分基因里面找
pseudotime_de <- differentialGeneTest(test[expressed_genes,],
fullModelFormulaStr = "~sm.ns(Pseudotime)")
pseudotime_de <- pseudotime_de[order(pseudotime_de$qval), ]
states_de <- differentialGeneTest(test[expressed_genes,],
fullModelFormulaStr = "~State")
states_de <- states_de[order(states_de$qval), ]
saveRDS(test, file = "test_monocle.rds")
write.table(pseudotime_de, file = "pseudotime_de.rds", quote = FALSE, sep = '\t', row.names = FALSE, col.names = TRUE)
write.table(states_de, file = "states_de.rds", quote = FALSE, sep = '\t', row.names = FALSE, col.names = TRUE)
6. 分支點(diǎn)的分析
解決的問題：沿著擬時(shí)序的方向仰禽，經(jīng)過不同的branch除秀，發(fā)生了哪些基因的變化？
經(jīng)常在文獻(xiàn)里面看到的monocle2熱圖就是這種分析涣达。
BEAM_res=BEAM(test,branch_point = 1,cores = 1)
#會(huì)返回每個(gè)基因的顯著性在辆，顯著的基因就是那些隨不同branch變化的基因
#這一步很慢
BEAM_res=BEAM_res[,c("gene_short_name","pval","qval")]
saveRDS(BEAM_res, file = "BEAM_res.rds")
畫圖
tmp1=plot_genes_branched_heatmap(test[row.names(subset(BEAM_res,qval<1e-4)),],
branch_point = 1,
num_clusters = 4, #這些基因被分成幾個(gè)group
cores = 1,
branch_labels = c("Cell fate 1", "Cell fate 2"),
#hmcols = NULL, #默認(rèn)值
hmcols = colorRampPalette(rev(brewer.pal(9, "PRGn")))(62),
branch_colors = c("#979797", "#F05662", "#7990C8"), #pre-branch, Cell fate 1, Cell fate 2分別用什么顏色
use_gene_short_name = T,
show_rownames = F,
return_heatmap = T #是否返回一些重要信息
)
pdf("branched_heatmap.pdf",width = 5,height = 6)
tmp1$ph_res
dev.off()
monocle2的熱圖怎么看
- 從你想研究的節(jié)點(diǎn)（第4步圖中的1節(jié)點(diǎn)）到根（一般是人為指定，圖中是Ucell對(duì)應(yīng)的狀態(tài)）都是pre-branch度苔。
- 然后沿著擬時(shí)序的方向匆篓，State 1對(duì)應(yīng)的枝是fate 1（對(duì)應(yīng)圖中的Gcell），State 2對(duì)應(yīng)的枝是fate 2（對(duì)應(yīng)Ccell）寇窑。軟件作者是這么說的：Cell fate 1 corresponds to the state with small id while cell fate 2 corresponds to state with bigger id
- 熱圖中的基因是BEAM函數(shù)找到的branch特異的基因鸦概。從熱圖中間向兩邊看，能看到不同的fate/branch基因的動(dòng)態(tài)變化甩骏。熱圖中列是擬時(shí)序窗市，行是基因先慷，熱圖中間是擬時(shí)序的開始。
返回的列表tmp1包含熱圖對(duì)象咨察，熱圖的數(shù)值矩陣论熙，熱圖主題顏色，每一行的注釋（基因?qū)儆谀囊粋€(gè)group）等信息摄狱。
根據(jù)行注釋提取出基因group赴肚，可以去做富集分析，后期加到熱圖的旁邊二蓝。像這樣：
gene_group=tmp1$annotation_row
gene_group$gene=rownames(gene_group)
library(clusterProfiler)
library(org.Hs.eg.db)
allcluster_go=data.frame()
for (i in unique(gene_group$Cluster)) {
small_gene_group=filter(gene_group,gene_group$Cluster==i)
df_name=bitr(small_gene_group$gene, fromType="SYMBOL", toType=c("ENTREZID"), OrgDb="org.Hs.eg.db")
go <- enrichGO(gene = unique(df_name$ENTREZID),
OrgDb = org.Hs.eg.db,
keyType = 'ENTREZID',
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.2,
readable = TRUE)
go_res=go@result
if (dim(go_res)[1] != 0) {
go_res$cluster=i
allcluster_go=rbind(allcluster_go,go_res)
}
}
head(allcluster_go[,c("ID","Description","qvalue","cluster")])
也可以換一種方式來展示具體的基因誉券，這些基因可以是：
- 隨擬時(shí)序、state而改變的基因（第5步）
- 細(xì)胞亞群的marker基因（seurat得到的差異基因）
- 分支點(diǎn)分析得到的分支特異的基因（第6步BEAM函數(shù)得到的基因）
test_genes=c("TFF3","GUCA2B")
pdf("genes_branched_pseudotime.pdf",width = 9,height = 4)
plot_genes_branched_pseudotime(test[test_genes,],
branch_point = 1,
color_by = "celltype",
cell_size=2,
ncol = 2)
dev.off()
因水平有限刊愚，有錯(cuò)誤的地方踊跟，歡迎批評(píng)指正！
在公眾號(hào)可以獲取本文的測試數(shù)據(jù)和全部代碼鸥诽。
