1. 流程
Monocle主要可以进行以下三种分析：
- 细胞的聚类、分类和计数。 Single-cell RNA-Seq experiments allow you to discover new (and possibly rare) subtypes of cells. Monocle helps you identify them.
- 重建单细胞轨迹。 In development, disease, and throughout life, cells transition from one state to another. Monocle helps you discover these transitions.
- 差异表达分析。 Characterizing new cell types and states begins with comparing them to other, better understood cells. Monocle includes a sophisticated but easy to use system for differential expression.
2. 实操
2.1 安装
install.packages("devtools")
devtools::install_github("cole-trapnell-lab/monocle-release@develop")
或者
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("monocle")
2.2 创建CellDataSet--CDS对象
# 从Seurat对象中提取构建CDS对象所需要的3个输入文件：表达矩阵信息、基因信息和表型信息
library(SeuratData) #加载seurat数据集
getOption('timeout')
options(timeout=10000)
#InstallData("pbmc3k")
data("pbmc3k")
sce <- pbmc3k.final
library(Seurat)
table(Idents(sce))
sce$celltype = Idents(sce)#确定细胞类型栏
sce= sce[,sce$celltype %in% c('CD14+ Mono','FCGR3A+ Mono')]#这个表示找sce数据集中这两种单核细胞，并将它们提取出来
table(Idents(sce))
library(monocle)
sample_ann <- sce@meta.data #准备meta.data
head(sample_ann)
#rownames(sample_ann)=sample_ann[,1]
gene_ann <- data.frame(
gene_short_name = rownames(sce@assays$RNA) , ##提取seurat对象中的基因-rownames
row.names = rownames(sce@assays$RNA)
) #准备基因名字
head(gene_ann)
pd <- new("AnnotatedDataFrame", #将准备的meta.data和基因名字转化为AnnotatedDataFrame文件
data=sample_ann)
fd <- new("AnnotatedDataFrame",
data=gene_ann)
ct=as.data.frame(sce@assays$RNA@counts)#准备表达矩阵--这里用的count
ct[1:4,1:4]
##创建CellDataSet对象
sc_cds <- newCellDataSet(
as.matrix(ct), #表达矩阵
phenoData = pd,#metadata，即每个细胞的信息
featureData =fd,#基因名称
expressionFamily = negbinomial.size(),#负二项分布有两种方法，这里选用了negbinomial.size
lowerDetectionLimit=1)
sc_cds
FPKM/TPM值通常是对数正态分布的，而UMIs或读计数使用负二项更好地建模。要处理计数数据，需要将负二项分布指定为newCellDataSet的expressionFamily参数,negbinomial.size()和negbinomial()：输入的表达矩阵为UMI,一般适用于10x的数据；negbinomial()的结果更准确，但是计算更耗时；一般建议采用negbinomial.size()。稀疏矩阵用negbinomial.size() tobit():适用于输入的表达矩阵为FPKM或者TPM, 构建monocle2的class时会自动进行log化计算 gaussianff():输入为log化后的FPKM或者TPM。(目前在单细胞数据中，FPKM已不多用，smart-seq2平台数据一般采用TPM)
可以使用不同的数据来源来创建CDS对象
2.2.1 直接读取表达矩阵来创建
library(data.table)
##读取数据
data <- fread("fpkm.txt",data.table = F,header = T)
pd <- fread("metadata.txt",data.table = F,header = T)
fd <- fread("gene_annotations.txt",data.table = F,header = T)
##创建
pd <- new("AnnotatedDataFrame", data = pd)
fd <- new("AnnotatedDataFrame", data = fd)
HSMM <- newCellDataSet(as.matrix(data),
phenoData = pd, featureData = fd,
expressionFamily = tobit())
###如果数据量大，建议转化为稀疏矩阵
HSMM <- newCellDataSet(as(as.matrix(data), "sparseMatrix"),
phenoData = pd,
featureData = fd,
expressionFamily = tobit())
2.2.2 将Seurat对象直接转化为CellDataSet对象--不同R版本操作有差异
importCDS(pbmc)
如果我们想要从Seurat对象或SCESet中导入所有的插槽，我们可以设置参数'import_all'为TRUE。#(默认为FALSE或只保留最小数据集)
2.3 估计size factor和离散度
size facotr帮助我们标准化细胞之间的mRNA的差异。 离散度值可以帮助我们进行后续的差异分析。 （类似于seurat的数据归一化处理）
cds <- estimateSizeFactors(cds)
#这一操作会在pData(cds)中添加一列SizeFactors
cds <- estimateDispersions(cds)
#获得负二项式分布数据的离散估计
所有中间值（基因方面的离散估计、趋势拟合的拟合离散估计等）都存储在 中
mcols(dds)，有关这些列的信息在mcols(mcols(dds))?没找到
与seurat把标准化后的表达矩阵保存在对象中不同，monocle只保存一些中间结果在对象中，需要用时再用这些中间结果转化。经过上面三个函数的计算，mycds对象中多了SizeFactors、Dipersions、num_cells_expressed和num_genes_expressed等信息。
2.4 质控-- 基于基因的过滤（基于表达量的过滤？）
sc_cds <- detectGenes(sc_cds, min_expr = 1) #对细胞表达的基因数量进行筛选，至少表达一个基因
#这一操作会在pData(cds)中添加一列num_genes_expressed detectGenes: Sets the global expression detection threshold to be used with this CellDataSet. Counts how many cells each feature in a CellDataSet object that are detectably expressed above a minimum threshold. Also counts the number of genes above this threshold are detectable in each cell.
sc_cds <- sc_cds[fData(sc_cds)$num_cells_expressed > 10, ]#对每个基因表达的细胞量进行筛选，至少有10个细胞表达这个基因 Create a functional data object
#这一操作会在fData(cds)中添加一列num_genes_expressed
sc_cds
2.5 细胞分类
Monocle官网教程提供了4个分类方法：
- Classifying cells by type
- Clustering cells without marker genes
- Clustering cells using marker genes
- Imputing cell type
建议先将细胞注释好再进行Monocle分析，不建议使用monocle做细胞分类。
2.6 轨迹定义基因选择及可视化和构建轨迹
The ordering workflow
- Step 1: 选择基因
- Step 2: 基因降维
- Step 3: 拟时序排列
2.6.1 选择基因
Monocle官网教程提供了4个选择方法：
- 选择发育差异表达基因?
- 选择clusters差异表达基因
- 选择离散程度高的基因❤我
- 自定义发育marker基因 前三种都是无监督分析方法，细胞发育轨迹生成完全不受人工干预；最后一种是半监督分析方法，可以使用先验知识辅助分析。
选择基因--基因筛选--进行排序--排序后可视化
#使用seurat选择的高变基因⚠️--
express_genes <- VariableFeatures(pbmc)
cds <- setOrderingFilter(cds, express_genes)
plot_ordering_genes(cds)
#使用clusters差异表达基因--cluster差异表达基因
deg.cluster <- FindAllMarkers(pbmc)
express_genes <- subset(deg.cluster,p_val_adj<0.05)$gene
cds <- setOrderingFilter(cds, express_genes)
plot_ordering_genes(cds)
#使用monocle选择的高变基因⚠️--离散程度高的基因
disp_table <- dispersionTable(cds)
disp.genes <- subset(disp_table, mean_expression >= 0.1 & dispersion_empirical >= 1 * dispersion_fit)$gene_id
cds <- setOrderingFilter(cds, disp.genes)
plot_ordering_genes(cds)
我用的下面这个--无监督聚类 离散程度比较高的
disp_table <- dispersionTable(cds) #计算均值和标准差
unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)#在上面计算的结果中筛选表达均值大于0.1的基因
cds <- setOrderingFilter(cds, unsup_clustering_genes$gene_id)
##这一步是很重要的，在我们得到想要的基因列表后，我们需要使用setOrderingFilter将它嵌入cds对象，后续的一系列操作都要依赖于这个list。
#setOrderingFilter之后，这些基因被储存在cds@featureData@data[["use_for_ordering"]]，可以通过table(cds@featureData@data[["use_for_ordering"]])查看
plot_ordering_genes(cds) #根据表达量排序，X轴是mean_expression,Y轴是dispersion_empirical
第一步是决定用哪些基因来进行聚类（特征选择）。我们可以使用所有的基因，但是我们将包括很多没有表达到足够高水平来提供有意义的信号的基因，它们只会给系统增加噪音。我们可以根据平均表达水平筛选基因，我们还可以选择细胞间异常变异的基因。这些基因往往对细胞状态有很高的信息含量。
师姐用的下面这个--半监督聚类（数据替换掉了，她用的下面的数据）--作为下面带心心部分的补充
pData(cds)$Cluster=pData(cds)$celltype
table(pData(cds)$Cluster)
#Test genes for differential expression
Sys.time()
diff_test_res <- differentialGeneTest(cds,
fullModelFormulaStr = "~Cluster")
Sys.time()
# Select genes that are significant at an FDR < 10%
sig_genes <- subset(diff_test_res, qval < 0.1)
sig_genes=sig_genes[order(sig_genes$pval),]
head(sig_genes[,c("gene_short_name", "pval", "qval")] )
cg=as.character(head(sig_genes$gene_short_name))
# 挑选差异最显著的基因可视化
# Plots expression for one or more genes as a jittered, grouped...
plot_genes_jitter(cds[cg,],
grouping = "Cluster",
color_by = "Cluster",
nrow= 3,
ncol = NULL )
师姐是下面这个---❤❤❤
#也可输入seurat筛选出的高变基因：
# expressed_genes <- VariableFeatures(pbmc)
使用monocle中按照细胞类型寻找差异基因的方式，有点像seurat中根据细胞类型进行FindAllmarkers
diff <-differentialGeneTest(cds[expressed_genes,],fullModelFormulaStr="~cell_type",cores=1)
#~后面是表示对谁做差异分析的变量，理论上可以为p_data的任意列名
head(diff)
##差异表达基因作为轨迹构建的基因,差异基因的选择标准是qval<0.01,decreasing=F表示按数值增加排序
deg <- subset(diff, qval < 0.01) #选出2829个基因
deg <- deg[order(deg$qval,decreasing=F),]
head(deg)
##差异基因的结果文件保存
write.table(deg,file="train.monocle.DEG.xls",col.names=T,row.names=F,sep="\t",quote=F)
#轨迹构建基因可视化
ordergene <- rownames(deg)
cds <- setOrderingFilter(cds, ordergene)
pdf("train.ordergenes.pdf")
plot_ordering_genes(cds)
dev.off()
#图黑色的点表示用来构建轨迹的差异基因，灰色表示背景基因。红色的线是根据基因表达大小和离散度分布的趋势(可以看到，找到的基因属于离散度比较高的基因) ⚠️选择的用于排序的基因数目一般在2000左右比较合适
gene数太多的话也可以选择top基因
ordergene <- row.names(deg)[order(deg$qval)][1:400]
2.6.2 降维
cds <- reduceDimension(cds, max_components = 2,
method = 'DDRTree')
2.6.3 拟时间轴轨迹构建
cds <- orderCells(cds)
#⚠️使用root_state参数可以设置拟时间轴的根，如下面的拟时间着色图中可以看出，左边是根。根据state图可以看出，根是State1，若要想把另一端设为根，可以按如下操作
#cds <- orderCells(cds, root_state = 5) #把State5设成拟时间轴的起始点
2.6.4 可视化
根据cds@phenoData@data中的表型信息(metadata)对细胞上色
#以pseudotime值上色 (Pseudotime是monocle2基于细胞基因表达信息计算的概率，表示时间的先后。)
pdf("train.monocle.pseudotime.pdf",width = 7,height = 7)
plot_cell_trajectory(cds,color_by="Pseudotime", size=1,show_backbone=TRUE)
dev.off()
图中123是分叉点，不代表特别意义，也不代表时间先后。这张图按照sudotime，时间先后是从左往右，左边是起点（root)。起点的设置可以使用orderCells的root_state参数，将右侧设置为起点。
以细胞类型上色
pdf("train.monocle.celltype.pdf",width = 7,height = 7)
plot_cell_trajectory(cds,color_by="cell_type", size=1,show_backbone=TRUE)
dev.off()
以细胞状态上色
pdf("train.monocle.state.pdf",width = 7,height = 7)
plot_cell_trajectory(cds, color_by = "State",size=1,show_backbone=TRUE)
dev.off()
state的多少是monocle算出来的，不能调整，与输入的用于轨迹学习的基因有关。分叉和顶点之间或者顶点和顶点之间为一个state，与发育轨迹时间先后没有关系，与细胞类型也不完全相关。
按照seurat分群排序细胞
pdf("seurat.clusters.pdf",width = 7,height = 7)
plot_cell_trajectory(cds, color_by = "seurat_clusters")
dev.off()
以细胞状态上色（拆分）“分面”轨迹图，以便更容易地查看每个状态的位置。
pdf("train.monocle.state.faceted.pdf",width = 10,height = 7)
plot_cell_trajectory(cds, color_by = "State") + facet_wrap("~State", nrow = 1)
dev.off()
改颜色
library(ggsci)
p1=plot_cell_trajectory(cds, color_by = "cell_type") + scale_color_npg()
p2=plot_cell_trajectory(cds, color_by = "State") + scale_color_nejm()
colour=c("#DC143C","#0000FF","#20B2AA","#FFA500","#9370DB","#98FB98","#F08080")
p3=plot_cell_trajectory(cds, color_by = "State") + scale_color_manual(values = colour)
p1|p2|p3
树形图
p1 <- plot_cell_trajectory(cds, x = 1, y = 2, color_by = "celltype") +
theme(legend.position='none',panel.border = element_blank()) + #去掉第一个的legend
scale_color_manual(values = colour)
p2 <- plot_complex_cell_trajectory(cds, x = 1, y = 2,
color_by = "celltype")+
scale_color_manual(values = colour) +
theme(legend.title = element_blank())
p1|p2
提取感兴趣的细胞（进行后续分析）
#比如对State7的细胞感兴趣
pdata <- Biobase::pData(cds)
s.cells <- subset(pdata, State=="7") %>% rownames()
save(s.cells, file = "Monocle_state7.rda")
保存结果
write.csv(pData(cds), "pseudotime.csv")
save(cds, file = "cds.rda")
2.7. 指定基因的可视化
2.7.1 直接查看一些基因随细胞状态等的表达变化
top基因
##选择前4个top基因并将其对象取出
keygenes <- head(ordergene,4)
cds_subset <- cds[keygenes,]
##可视化：以state/celltype/pseudotime进行
p1 <- plot_genes_in_pseudotime(cds_subset, color_by = "State")
p2 <- plot_genes_in_pseudotime(cds_subset, color_by = "cell_type")
p3 <- plot_genes_in_pseudotime(cds_subset, color_by = "Pseudotime")
plotc <- p1|p2|p3
ggsave("Genes_pseudotimeplot.pdf", plot = plotc, width = 16, height = 8)
指定基因
s.genes <- c("SELL","CCR7","IL7R", "CD84","CCL5","S100A4")
p1 <- plot_genes_jitter(cds[s.genes,], grouping = "State", color_by = "State")
p2 <- plot_genes_violin(cds[s.genes,], grouping = "State", color_by = "State")
p3 <- plot_genes_in_pseudotime(cds[s.genes,], color_by = "State")
plotc <- p1|p2|p3
ggsave("Genes_Jitterplot.pdf", plot = plotc, width = 16, height = 8)
由于state1-7并不代表分化时间先后，因此这张图还要结合其他的图来看
拟时序展示单个基因表达量
colnames(pData(cds))
pData(cds)$CCL5 = log2( exprs(cds)['CCL5',]+1)
p1=plot_cell_trajectory(cds, color_by = "CCL5") + scale_color_gsea()
pData(cds)$S100A4 = log2(exprs(cds)['S100A4',]+1)
p2=plot_cell_trajectory(cds, color_by = "S100A4") + scale_color_gsea()
library(patchwork)
p1+p2
2.8 寻找拟时相关的基因（拟时差异基因）
官方给出的差异分析有三大方法：
- 1、Basic Differential Analysis
- 2、Finding Genes that Distinguish Cell Type or State
- 3、Finding Genes that Change as a Function of Pseudotime
我们重点关注第三个：根据伪时间功能寻找差异基因
- sm.ns函数指出Monocle应该通过表达式值拟合自然样条曲线，以帮助它将表达式的变化描述为进程的函数。
寻找拟时差异基因（qvalue体现基因与拟时的密切程度）绘制热图
#这里是把排序基因（ordergene）提取出来做回归分析，来找它们是否跟拟时间有显著的关系
#如果不设置，就会用所有基因来做它们与拟时间的相关性
Time_diff <- differentialGeneTest(cds[ordergene,], cores = 1,
fullModelFormulaStr = "~sm.ns(Pseudotime)")
Time_diff <- Time_diff[,c(5,2,3,4,1,6,7)] #把gene放前面，也可以不改
write.csv(Time_diff, "Time_diff_all.csv", row.names = F)
Time_genes <- Time_diff %>% pull(gene_short_name) %>% as.character()
p=plot_pseudotime_heatmap(cds[Time_genes,], num_clusters=4, show_rownames=T, return_heatmap=T)
ggsave("Time_heatmapAll.pdf", p, width = 5, height = 10)
把每个cluster的基因单独提出来做分析
# 根据帮助文档，我们修改参数，这样monocle的plot_pseudotime_heatmap函数就有返回值了，是一个对象。
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)
write.csv(clustering, "Time_clustering_all.csv", row.names = F)
#作者：Seurat_Satija
#链接：https://www.jianshu.com/p/972e65671f45
#来源：简书
#著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。
这样就把每个基因属于哪个cluster提取出来了，后续可以做每个cluster的富集分析。 绘出的热图可以让我们观测到假时间依赖性基因中的不同基因模块在不同的时间内共同变化，能比较好的回答时间序列基因表达中“哪些基因遵循相似的动力学趋势”这一常见问题。
显著差异基因按热图结果排序并保存
hp.genes <- p$tree_row$labels[p$tree_row$order]
Time_diff_sig <- Time_diff[hp.genes, c("gene_short_name", "pval", "qval")]
write.csv(Time_diff_sig, "Time_diff_sig.csv", row.names = F)
手动选择基因来绘制热图，查看其表达模式
marker_genes <- row.names(subset(fData(cds),
gene_short_name %in% c("MEF2C", "MEF2D", "MYF5",
"ANPEP", "PDGFRA","MYOG",
"TPM1", "TPM2", "MYH2",
"MYH3", "NCAM1", "TNNT1",
"TNNT2", "TNNC1", "CDK1",
"CDK2", "CCNB1", "CCNB2",
"CCND1", "CCNA1", "ID1")))
diff_test_res <- differentialGeneTest(cds[marker_genes,],
fullModelFormulaStr = "~sm.ns(Pseudotime)")
sig_gene_names <- row.names(subset(diff_test_res, qval < 0.1))
plot_pseudotime_heatmap(cds[sig_gene_names,],
num_clusters = 6,
cores = 1,
show_rownames = T)
2.9 单细胞轨迹的“分支”分析--没搞太懂，后面用到了再看看？
Monocle提供了一个特殊的统计测试：分支表达式分析建模，或BEAM。 BEAM(Branched expression analysis modeling)是一种统计方法，用于寻找以依赖于分支的方式调控的基因。
