1. 文献怎么用 monocle2 ?
(1) ex1
01.png 02.png作者处理2种细胞使用的策略还不一样!
- 中性粒细胞 neutrophil 是先subset出来,去除污染细胞,然后Seurat聚类,找类之间 top20 DEG 做 cell order(就是monocle2的排序);
根据前体细胞的marker基因的表达,选择那个state为起点。
- 中性粒细胞 neutrophil 是先subset出来,去除污染细胞,然后Seurat聚类,找类之间 top20 DEG 做 cell order(就是monocle2的排序);
- 而 monocyte 是subset出cluster1,去除污染细胞,然后像上文一样cluster,然后DEG是 Monocle 的 differentialGeneTest() 找的,过滤条件为 qval<1e-10 and cell>20 这样的基因 order cell;
起点是包含最多WT库的分支。
- 而 monocyte 是subset出cluster1,去除污染细胞,然后像上文一样cluster,然后DEG是 Monocle 的 differentialGeneTest() 找的,过滤条件为 qval<1e-10 and cell>20 这样的基因 order cell;
(2) ex2
https://www.sciencedirect.com/science/article/pii/S0092867417305962
p1.pngFigure 6. The CD8+ and CD4+ T Cell State Transition Analysis Based on Integrated Expression and TCR Clonality
(A) The ordering of CD8+ T cells along pseudotime in a two-dimensional state-space defined by Monocle2. Cell orders are inferred from the expression of most dispersed genes across CD8+ T cell populations sans MAIT. Each point corresponds to a single cell, and each color represents a T cell cluster.
// 使用方差最大的基因 dispersed genes,没说多少个。
(B) The same pseudotime plot as in (A) for four clusters of CD4+ T helper cells.
(3) ex3
https://www.nature.com/articles/s41586-018-0694-x
To characterize the potential process of T cell functional changes and determine the potential lineage differentiation between diverse T cells, we applied the Monocle (version 2) algorithm50 with the top 700 signature genes of CD8+ T cells excluding MAIT cells, based on the rank of F statistic generated by ANOVA (Supplementary Table 5).
// 使用ANOVA的F统计量,top 700 个基因。
(4) ex4
https://www.nature.com/articles/s41467-020-20019-0
First, the top 1000 HVGs were selected using the function ‘FindVariableFeatures’ in Seurat v3.
(10) tips
- 建议 过滤、分群、差异基因都在 Seurat 中完成,然后只用 monocle2 做 拟时序分析。
- 不同分组间的细胞尽量不要放在一起做轨迹分析,同一组的生物学重复可以一起分析。明显没有生物学相关性的细胞也不要放在一起做轨迹分析。
2. 从 Seurat 对象新建 Monocle 对象
#' get Monocle2 obj from Seurat obj
#'
#' version 0.2
#'
#' @param sce Seurat obj
#'
#' @return Monocle2 obj
#' @export
#'
#' @examples
#' cds = Seurat2Monocle2(sce)
Seurat2Monocle2=function(sce, do.init=T){
# get data
#expr_matrix=as.matrix(sce@assays$RNA@counts)
expr_matrix=GetAssayData(sce, assay = 'RNA', slot = 'counts')
cell_data=sce@meta.data
gene_data=data.frame(
gene_short_name=row.names(sce), #must have this column
row.names = row.names(sce)
)
# new obj
pd <- new("AnnotatedDataFrame", data = cell_data)
fd <- new("AnnotatedDataFrame", data = gene_data)
cds <- newCellDataSet(expr_matrix,
phenoData = pd,
featureData = fd,
lowerDetectionLimit = 0.5,
expressionFamily=negbinomial.size())
# init processing
if(do.init){
cds <- estimateSizeFactors(cds) #add column: Size_Factor
cds <- estimateDispersions(cds) #Removing 323 outliers
cds=detectGenes(cds, min_expr = 0.1) #add column: num_genes_expressed
}
print( head(pData(cds)) )
#
return(cds)
}
示例:
# subset of Seurat obj
sce=subset(scObj_nue, idents = c(9,10), invert=T,
downsample=50)
sce=subset(sce, subset= origin=="APC")
sce #298
DimPlot(sce, label=T, split.by = "origin")
# 1.load data to monocle2
cds = Seurat2Monocle2(sce)
class(cds)
3. 拟时序分析
The ordering workflow
Step 1: choosing genes that define progress
Step 2: reducing the dimensionality of the data
Step 3: ordering the cells in pseudotime
除了第一步比较灵活外,剩下2步是固定的。最后是各种可视化。
(1) step1 获取差异基因列表 DEG list
这个gene list 对拟时序结构具有决定作用。如果结果不理想,这里是优化的重点。
## Trajectory Step 1: choosing genes that define progress----
#第1步) 获取 输入用于排序的基因 symbol list
Monocle官网教程提供了4个选择方法:
* 选择发育差异表达基因
* 选择clusters差异表达基因
* 选择离散程度高的基因
* 自定义发育marker基因
前三种都是无监督分析方法,细胞发育轨迹生成完全不受人工干预;
最后一种是半监督分析方法,可以使用先验知识辅助分析。
##1)使用clusters差异表达基因
deg.cluster <- FindAllMarkers(scRNA.Osteoclastic)
diff.genes <- subset(deg.cluster,p_val_adj<0.05)$gene
HSMM <- setOrderingFilter(HSMM, diff.genes)
plot_ordering_genes(HSMM)
##2)使用seurat选择的高变基因
var.genes <- VariableFeatures(scRNA.Osteoclastic)
HSMM <- setOrderingFilter(HSMM, var.genes)
plot_ordering_genes(HSMM)
##3)使用monocle选择的高变基因
disp_table <- dispersionTable(HSMM)
disp.genes <- subset(disp_table, mean_expression >= 0.1 & dispersion_empirical >= 1 * dispersion_fit)$gene_id
HSMM <- setOrderingFilter(HSMM, disp.genes)
plot_ordering_genes(HSMM)
##4) Monocle2 的 differentialGeneTest() 计算的差异基因,比如保留 qval<1e-10 且 cell number >10 的基因
diff_test_res <- differentialGeneTest(cds[expressed_genes,],
fullModelFormulaStr = "~percent.mt")
# 20:35 --> 20:37 特别慢,仅仅1k个细胞就需要10min。
ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
cds <- setOrderingFilter(cds, ordering_genes)
plot_ordering_genes(cds)
(2). step2 降维
## Trajectory step 2: reduce data dimensionality----
# cds <- setOrderingFilter(cds, ordering_genes)
cds <- reduceDimension(cds, max_components = 2, method = 'DDRTree')
(3). 按照轨迹排序细胞
## Trajectory step 3: order cells along the trajectory ----
cds <- orderCells(cds)
(4). 可视化
plot_cell_trajectory(cds, color_by = "Pseudotime")
plot_cell_trajectory(cds, color_by = "seurat_clusters")
plot_cell_trajectory(cds, color_by = "State")
# "State" is just Monocle's term for the segment of the tree.
## 如何重叠太厉害,还可以分面
plot_cell_trajectory(HSMM, color_by="XXX") +
facet_wrap(~XXX, nrow=1)
如果自动指定的不对,还可以手动指定哪个 state 是起点。
如果没有时间序列,就需要根据生物学知识找到root。
比如: 高度增殖的祖细胞群产生两种有丝分裂后的细胞。所以根应该有表达高水平增殖标记的细胞。
cds<- orderCells(cds, root_state = 1)
plot_cell_trajectory(cds, color_by = "Pseudotime")
更多可视化函数:
plot_cell_trajectory(HSMM_myo, color_by = "State",cell_size = 0.75) ####以state进行着色
plot_cell_trajectory(HSMM_myo, color_by = "State",cell_size = 0.75)+facet_wrap(~State, nrow = Cluster_num) ##绘制state的分面图
plot_cell_trajectory(HSMM_myo, color_by = "Pseudotime",cell_size = 0.75) ##根据拟时间值着色
plot_cell_trajectory(HSMM_myo, color_by = "Cluster",cell_size = 0.75) ##根据细胞聚类的进行着色
plot_cell_trajectory(HSMM_myo, color_by = "Cluster",cell_size = 0.75)+facet_wrap(~Cluster, nrow =Cluster_num) ###绘制clusster的分面图
plot_cell_trajectory(HSMM_myo, color_by = "seurat_clusters",cell_size = 0.75) +scale_color_manual(values=col) ##如果有Seurat生的rds文件的话,按照seurat中分的群进行着色,如果不想用ggplot的默认色,可以提供颜色列表col list。
plot_cell_trajectory(HSMM_myo, color_by = "seurat_clusters", cell_size = 0.75) + facet_wrap (~seurat_clusters, nrow =Cluster_num) ###按照seurat中分的群绘制分面图。
plot_cell_trajectory(HSMM_myo, color_by = "stim",cell_size = 0.75) ###按照样本进行着色
plot_cell_trajectory(HSMM_myo, color_by = "stim",cell_size = 0.75)+facet_wrap(~stim, ncol = 2)
##绘制样本着色的分面图。
(5) 查看核心基因轨迹
# 根据 cds子集对象 画伪时间上的基因表达
cds_subset <- HSMM_filtered[my_genes,]
plot_genes_in_pseudotime(cds_subset, color_by = "seurat_clusters")
plot_genes_in_pseudotime(cds_subset, color_by = "State")
(6) 查看核心基因表达热图
diff_test_res = differentialGeneTest(cds[marker_genes,], fullModelFormulaStr="~sm.ns(Pseudotime)")
sig_gene_names = row.names(subset(diff_test_rtes, qval<0.01)) #原文例子是 0.1,是不是太大了?
# 主要是这个基因列表
plot_pseudotime_heatmap(cds[sig_gene_names,], num_clusters=6, cores=1, show_rownames=T)
热图可以采用top50的基因进行展示,当然,也可以自己筛选。
一行为每个基因对应的表达量。一列为一个拟时间点,一个拟时间值可能代表多个细胞。
4. BEAM 分析: 单细胞轨迹的“分支”分析
# 查看分支点:就是黑色圆圈中的白色数字
plot_cell_trajectory(cds, color_by = "State")
# BEAM进行统计分析:我们分析branch_point = 1这个分支处的细胞命名分叉是如何进行的。
BEAM_res = BEAM(cds[expressed_genes, ], branch_point=1, cores=4)
# 对影响分析的基因根据qvalue进行排序
BEAM_res=BEAM_res[order(BEAM_res$qval),]
BEAM_res=BEAM_res[, c("gene_short_name", "pval", "qval")]
# 绘制分支热图
my_branched_heatmap = plot_genes_branched_heatmap(
cds[row.names(subset(BEAM_res, qval<1e-4)),],
branch_point=1,
num_clusters=4, cores=2,
use_gene_short_name=T,
show_rownames=T)
软件作者是这么说的:Cell fate 1 corresponds to the state with small id while cell fate 2 corresponds to state with bigger id
能自定义颜色吗?如果基因过多,就别不显示了,看不清。
tmp1=plot_genes_branched_heatmap(test[row.names(subset(BEAM_res,qval<1e-4)),],
branch_point = 1,
num_clusters = 4, #这些基因被分成几个group
cores = 1,
branch_labels = c("Cell fate 1", "Cell fate 2"), #规定标签
#hmcols = NULL, #默认值
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()
返回的列表tmp1包含热图对象,热图的数值矩阵,热图主题颜色,每一行的注释(基因属于哪一个group)等信息。
根据行注释提取出基因group,可以去做富集分析,后期加到热图的旁边。像这样:
https://www.cnblogs.com/TOP-Bio/p/14856979.html
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")])
(2) 单个基因的伪时间曲线,类似 plot_genes_in_pseudotime(),不同的是该函数显示2个发育路径。
展示具体的基因,这些基因可以是:
- 随拟时序、state而改变的基因(第5步)
- 细胞亚群的marker基因(seurat得到的差异基因)
- 分支点分析得到的分支特异的基因(第6步BEAM函数得到的基因)
We can plot a couple of these genes, such as Pdpn and Sftpb (both known markers of fate in this system), using the plot_genes_branched_pseudotime function, which works a lot like the plot_genes_in_pseudotime function, except it shows two kinetic trends, one for each lineage, instead of one.
另外,可以绘制基因在所有细胞中的表达量在不同分支上的表达量情况。我们可以根据分叉的state位置来判定这个基因对分支的影响。在拟时间分支树构建的过程中,其实是有很多细小的分支,最终被拟合成一条直线,便于可视化的友好性。下图中的Branch表示挑选了两个不同的未被拟合的实际分支。
plot_genes_branched_pseudotime(
HSMM[rownames(subset(BEAM_res, qval<1e-8)),],
branch_point=1,
color_by="orig.ident",
ncol=3)
> str(HSMM)
# 示例2
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()
网友评论