跟着官网学习单细胞Monocle V3包（上）
大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~
就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~
这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！
豆豆写于19.8.21+28
这次花费的时间较长，个人感觉这个包比Seurat、Scater难学，说明文档写的实在是太长了，于是分为两章发送！
前言
依旧主打全人工理解（不是翻译）
对这款软件的了解主要是因为它做的发育轨迹、拟时序分析很漂亮，这也是它的主打
Monocle的官网在：
版本2：https://cole-trapnell-lab.github.io/monocle-release/docs/#installing-monocle
版本3：https://cole-trapnell-lab.github.io/monocle3/monocle3_docs/
首先看新版本的特性和升级之处：
处理的细胞数增加很多（millions of cells）
针对发育生物学领域，做了一些重大改进：
发育轨迹的研究流程优化
支持UMAP推断发育轨迹，无缝衔接到
reduceDimension函数，或者使用独立的函数UMAP支持多个祖源（toots）的发育推断
利用这个算法approximate graph abstraction 理解各条发育轨迹的分离及平行发育的轨迹
3D构建发育轨迹
除了发育轨迹以外，官方还介绍它在亚型发现、差异分析方面表现不错
另外，版本3在持续不断优化，官方称每几个周就会增加一些新功能
其中的算法细节看2019年发表的Cao & Spielmann et al.
它的主要分析流程是：
注意：以下是介绍版本3（版本2的使用和3存在一些不同之处，会在”跟着官网学习“这一部分的最后标注）
安装版本3
要求R版本3.5以上，Bioconductor版本3.5以上
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
# 一些依赖包
bioc_pkgs <- c('BiocGenerics', 'DelayedArray', 'DelayedMatrixStats',
'limma', 'S4Vectors', 'SingleCellExperiment',
'SummarizedExperiment')
for (bpkg in bioc_pkgs){
if (! require(bpkg,character.only=T) ) {
BiocManager::install(bpkg,ask = F,update = F)
require(bpkg,character.only=T)
}
}
# (可选)如果要在细胞聚类时设定分辨率参数(resolution)，就要安装一个python包
install.packages("reticulate")
reticulate::py_install("louvain")
# 现在安装3版本，还是要通过github
devtools::install_github('cole-trapnell-lab/monocle3')
# 重启Rstudio，检测是否成功
library(monocle3)
Monocle示例1-细胞聚类及鉴定亚群
创建对象
整体代码就是这样（先了解概况然后下面有实际练习）：
cds <- new_cell_data_set(expression_matrix,
cell_metadata = cell_metadata,
gene_metadata = gene_annotation)
但是这个对象cell_data_set的设置需要好好理解：
这个对象的灵感来自SingleCellExperiment，因此它的操作也是类似的
expression_matrix是一个行为基因，列为细胞样本的表达矩阵cell_metadata是一个数据框，行为细胞，列是细胞的属性（比如细胞类型、培养环境、培养时间等）gene_metadata是一个数据框，行为feature信息（比如基因），列是基因属性（比如GC含量）【这个很新鲜，因为其他的基于SingleCellExperiment的包如Scater只需要表达矩阵和样本信息即可】。更重要的是，这个数据框中必须有这么一列：`gene_short_name`，其中保存基因名（简称或Symbol标准名均可，用于作图）并且要满足：
expression_matrix的行与gene_metadata的行对应；expression_matrix的列与cell_metadata的行对应
然后测试一下（下载速度会很慢！）
小提示：使用服务器
axel下载，速度会加快（不过这要看服务器的网络配置），如果给力的话，下载均速能到120kb/s，一两分钟就下载好
如果网络真的不给力也没关系，我帮大家下载下来了，在豆豆和花花的公众号”生信星球“后台回复”monocle“就可以得到这三个文件啦：
# 下载数据(如果已经有了，可以跳过)
expression_matrix <- readRDS(url("http://staff.washington.edu/hpliner/data/cao_l2_expression.rds"))
cell_metadata <- readRDS(url("http://staff.washington.edu/hpliner/data/cao_l2_colData.rds"))
gene_annotation <- readRDS(url("http://staff.washington.edu/hpliner/data/cao_l2_rowData.rds"))
# 加载、探索数据
expression_matrix <- readRDS(file = 'cao_l2_expression.rds')
cell_metadata <- readRDS(file = 'cao_l2_colData.rds')
gene_annotation <- readRDS(file = 'cao_l2_rowData.rds')
dim(expression_matrix)
expression_matrix[1:3,1:3]
cell_metadata[1:3,1:3]
head(gene_annotation)
# 创建CDS(Cell Data Set)对象
cds <- new_cell_data_set(expression_matrix,
cell_metadata = cell_metadata,
gene_metadata = gene_annotation)
cds
数据预处理
这一步就是告诉Monocle，你想怎么对数据进行归一化、标准化，初步降维是使用PCA（针对标准的RNA-seq，需要指定计算的主成分数）还是LSA（Latent semantic analysis，针对ATAC-seq）。这些都是为后面的聚类操作（tSNE、UMAP做准备）
# 比较耗时
cds <- preprocess_cds(cds, num_dim = 100) #大约2 mins
# 这里设置这么多主成分，可能是因为细胞数太多(4万多个)，成分太少不足以代表整体
这个preprocess_cds看下帮助文档就能得知：
降维方法
method = c("PCA", "LSI")，默认是PCA默认维度
num_dim是50初步降维前归一化
norm_method默认是log处理设置降维方法是PCA时，自动进行标准化处理
检查一下使用的主成分（PCs）是否能够抓取最主要的基因表达变化信息
plot_pc_variance_explained(cds)
# 很像Seurat的ElbowPlot()
可以看到，超过100个主成分后对整体变化的贡献值就没有太多了，而且每多一个PC就会对下游数据分析造成压力
继续降维及可视化
关于维度有必要好好再了解下：
（详见
?reduce_dimension） 那我就编一个”降维“小故事吧
每个细胞中都有2万多基因，可以被当做高维空间上的一个点，而这里说的每个维度就是就不同基因的表达量，也就是说，我们现在有2万个维度上的4万个点；处理的维度越多，后面的轨迹推断也就是越难
好在基因之间不是毫无瓜葛的，各种基因之间总会有错综复杂的联系，通过排除这种”冥冥中的相关性“就会提取出它们之间最不同的地方，更能说明整体的一个特征。于是线性降维就能将2万个维度降到这里的100个，那么之后再处理的话，可能线性降维也有压力。”专人专事：“~于是找来了更专业的非线性降维工具tSNE、UMAP，它们最终能实现再将这100个成分继续压缩，最终得到我们能理解的二维或者3维空间（那么其中包含的样本特征一定是精华中的精华）
英文的解释就是：Each cell can be viewed as a point in a high-dimensional space, where each dimension describes the expression of a different gene.
这里主要采用非线性降维的方法：一种是目前很流行的tSNE，另一种是2018年发布的基于Python的UMAP（UMAP处理大量细胞的算法更优秀；另外与t-SNE相比，UMAP可以保留更多的细胞亚群连续性）
cds <- reduce_dimension(cds, preprocess_method = "PCA",reduction_method = c("UMAP"))
# 大约运行1.5 min
关于这个函数：
多种继续降维的算法：如
"UMAP", "tSNE", "PCA" and "LSI"，不过函数默认使用UMAP关于运行加速：
基于BiocParallel支持多线程运行，使用cores定义；或者使用umap.fast_sgd=TRUE参数可以加速降维
但是自己加速的话，它会友情提示一下：# Note: reduce_dimension will produce slightly different output each time you run it unless you set 'umap.fast_sgd = FALSE' and 'cores = 1'默认得到两个维度
需要定义上一步
preprocess_method使用的算法，默认是LSI
进行可视化：
plot_cells(cds)
图中的每个点都代表cds对象中的不同细胞，可以看到不同的细胞分成了不同的群，有的群有几千个细胞，有的群只有少数几个，我们这里直接使用测试数据中做好的细胞标记
# 数据中一共提供了30类细胞(细胞类型及数量见：)
table(colData(cds)$"cao_cell_type")
# 根据这个上色
plot_cells(cds, color_cells_by="cao_cell_type")
# (另外除了用cao_cell_type细胞类型，还能用下面的任意一个)
> colnames(colData(cds))
[1] "plate" "cao_cluster" "cao_cell_type" "cao_tissue" "Size_Factor"
可以看到许多细胞类型再UMAP结果中靠的很近
根据基因表达量进行映射：
plot_cells(cds, genes=c("cpna-2", "egl-21", "ram-2", "inos-1"))
如果想使用tSNE方法降维的话：
cds <- reduce_dimension(cds, reduction_method="tSNE")
# 然后对tSNE结果可视化
plot_cells(cds, reduction_method="tSNE", color_cells_by="cao_cell_type")
其实可以看到，tSNE结果的各个细胞亚群内部的关联性不如UMAP（比如这个stem cells）
检查、移除批次效应 => residual_model_formula_str参数
为了更加真实地反映差异基因的来源是生物学因素，而不是技术因素（比如细胞板批次、细胞培养差别、测序差别等），就首先要检查有没有所谓的”技术批次“对数据差异产生影响。
当进行降维分析时，应该经常检查批次效应，最好在colData(cds) 添加一列，其中注明细胞的批次信息（比如来自哪个细胞板），然后就能根据这个对细胞的批次进行可视化，结果最好是各个细胞亚群的不同批次混在一起，这样我们就能认为：降维所采用的”特异性“的特征并不来自批次
plot_cells(cds, color_cells_by="plate", label_cell_groups=FALSE)
这里看到还不错，但是如果真的担心不同批次会产生其他的影响，可以在降维之前对批次进行校正，保证后面的各个细胞群只来自一个批次，那么它们之间的分群差异就更可能是由于生物因素导致
# 假入要对批次进行校正
cds = preprocess_cds(cds, num_dim = 100, residual_model_formula_str = "~ plate")
# 校正后再降维
cds = reduce_dimension(cds)
# 再次检查批次效应
plot_cells(cds, color_cells_by="plate", label_cell_groups=FALSE)
细胞聚类 => cluster_cells()
这是细胞分群过程中非常重要的一步。Monocle利用Louvain community detection 这个非监督聚类方法（它是phenoGraph算法的一部分），它和我们常见的Kmeans、hclust层次聚类等还不太一样，这个方法更倾向于找到一个网络中联系紧密的部分，而经常忽略节点的特性；而我们常见的聚类呢，它一般会忽视整体中各个部分的联系，而通过计算两个节点(目标)之间的距离（如欧式距离、曼哈顿距离、余弦相似度等） 找到相似的特性
这一步依赖一个Python模块louvain
关于Python模块在R中的安装： 可以看我写过的Linux/Mac/Windows的Rstudio安装Python模块总报错，怎么破？
做的话很简单，一个函数，结果保存在了cds@clusters$UMAP$clusters中：
cds = cluster_cells(cds, resolution=c(10^seq(-6,-1)))
plot_cells(cds)
# 这个可视化函数如果不加任何参数，它默认对不同的cluster上色
# 那么有哪些cluster呢？用下面函数查看：
cds@clusters$UMAP$clusters #或者
clusters(cds)
#一共有85个cluster
不过这样细分看的太混杂，于是 <a href="https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1663-x" "="" style="font-size: inherit; line-height: inherit; color: rgb(30, 107, 184); overflow-wrap: break-word;">Alex Wolf et al 利用他们开发的 PAGA 算法将一些小的cluster聚合成一个大的partition，结果也是保存在了cluster_cells的计算结果中cds@clusters$UMAP$partitions
plot_cells(cds, color_cells_by="partition", group_cells_by="partition")
# 检查有多少partition
cds@clusters$UMAP$partitions
# 整合成了36个大的partition
除了默认按照cluster编号进行标记不同的亚群，还可以按照每个cluster对应的细胞类型进行可视化：
plot_cells(cds, color_cells_by="cao_cell_type")
不过这样又很多重复的细胞类型，我们可以对label进行简化：
plot_cells(cds, color_cells_by="cao_cell_type", label_groups_by_cluster=FALSE)
归根结底，都是对细胞的类型、坐标、颜色映射。就是说：
cds@colData$cao_cell_type中存储了每个细胞对应的细胞类型；cds@reducedDims$UMAP中存储了每个细胞UMAP降维后的坐标；cds@clusters$UMAP$partitions或cds@clusters$UMAP$clusters中存储了每个细胞聚类后对应的细胞分群；接下来就是根据坐标对每个细胞进行分群上色、标记类型
找marker基因 => top_markers()
确定了各种亚群以后，我们就想知道什么导致了它们的差异，这就是寻找marker基因的过程，利用top_markers()函数
marker_test_res = top_markers(cds, group_cells_by="partition", reference_cells=1000, cores=8)
# 设置reference_cells是随机挑出来这些数量的细胞作为参照，然后让top_markers和参照集中的基因进行显著性检验；另外reference_cells还可以是来自colnames(cds)的细胞名
marker_test_res[1:4,1:4]
结果得到一个数据框，其中包含了每个partition （因为上面我们设置了group_cells_by）的特异表达基因。
然后我们可以自行过滤：
top_specific_markers = marker_test_res %>%
filter(fraction_expressing >= 0.10) %>%
group_by(cell_group) %>%
top_n(1, pseudo_R2)
# 基因id去重
top_specific_marker_ids = unique(top_specific_markers %>% pull(gene_id))
然后可以对每组的marker基因可视化=> plot_genes_by_group()：
plot_genes_by_group(cds,
top_specific_marker_ids,
group_cells_by="partition",
ordering_type="maximal_on_diag",
max.size=3)
如果要多看几个marker，只需要修改一下top_n
top_specific_markers = marker_test_res %>%
filter(fraction_expressing >= 0.10) %>%
group_by(cell_group) %>%
top_n(3, pseudo_R2)
top_specific_marker_ids = unique(top_specific_markers %>% pull(gene_id))
plot_genes_by_group(cds,
top_specific_marker_ids,
group_cells_by="partition",
ordering_type="cluster_row_col",
max.size=3)
对细胞类型进行注释
因为我们上面的测试数据中给出了细胞类型，所以可以进行聚类后的可视化，但实际上，我们通常只能看到不同的cluster数字和它们的分布（也就是类似于as.character(partitions(cds))的结果），因此利用marker基因去重新定义细胞类型是非常关键的一步。
# 先将partitions的分组由因子型转为字符型
colData(cds)$assigned_cell_type = as.character(partitions(cds))
# 再对字符型重新定义
colData(cds)$assigned_cell_type = dplyr::recode(colData(cds)$assigned_cell_type,
"1"="Body wall muscle",
"2"="Germline",
"3"="Unclassified neurons",
"4"="Seam cells",
"5"="Coelomocytes",
"6"="Pharyngeal epithelia",
"7"="Vulval precursors",
"8"="Non-seam hypodermis",
"9"="Intestinal/rectal muscle",
"10"="Touch receptor neurons",
"11"="Pharyngeal neurons",
"12"="Am/PH sheath cells",
"13"="NA",
"14"="Unclassified neurons",
"15"="flp-1(+) interneurons",
"16"="Canal associated neurons",
"17"="Pharyngeal gland",
"18"="Other interneurons",
"19"="Ciliated sensory neurons",
"20"="Ciliated sensory neurons",
"21"="Ciliated sensory neurons",
"22"="Ciliated sensory neurons",
"23"="Ciliated sensory neurons",
"24"="Ciliated sensory neurons",
"25"="Oxygen sensory neurons",
"26"="Ciliated sensory neurons",
"27"="Unclassified neurons",
"28"="Pharyngeal gland",
"29"="Ciliated sensory neurons",
"30"="Ciliated sensory neurons",
"31"="Ciliated sensory neurons",
"32"="Ciliated sensory neurons",
"33"="Pharyngeal muscle",
"34"="Failed QC")
plot_cells(cds, group_cells_by="partition", color_cells_by="assigned_cell_type")
# Seurat中用RenameIdents 进行细胞类型重定义
# 另外，想从中取子集、过滤的话
cds[,colData(cds)$assigned_cell_type != "Failed QC"]
基于Garnett的自动化注释
下面这部分先了解下方法就行，日后有需要再细学
前面对细胞类型进行手动注释有点麻烦，而且cluster编号一旦更改，又要手动去注释一遍。
利用Monocle团队开发的 Garnett ，它可以根据marker基因对细胞进行分类
预处理
# step-1：首先根据上面得到的细胞类型assigned_cell_type，找top_marker
assigned_type_marker_test_res = top_markers(cds,
group_cells_by="assigned_cell_type",
reference_cells=1000,
cores=8)
# step-2：过滤（阈值自定义）
garnett_markers = assigned_type_marker_test_res %>%
filter(marker_test_q_value < 0.01 & specificity >= 0.5) %>%
group_by(cell_group) %>%
top_n(5, marker_score)
# step-3：去重复
garnett_markers = garnett_markers %>% group_by(gene_short_name) %>%
filter(n() == 1)
# step-4：生成marker文件
generate_garnett_marker_file(garnett_markers, file="./marker_file.txt")
得到这个文件，其实也只是自动化实现的第一步，我们自己可以根据已有的知识在其中增加或删除一些marker信息；另外如果发现多个细胞类型中的marker基因有大量重合，就要考虑这几种细胞类型可能都是某一种的细胞的亚型。还有很多关于marker文件修改的帮助信息：Garnett documentation
方法一：自己使用Garnett
# 安装Garnett
devtools::install_github("cole-trapnell-lab/garnett", ref="monocle3")
library(garnett)
# 安装物种注释库(这里的物种是C. elegans秀丽隐杆线虫)，目的是为了做基因ID转换
BiocManager::install("org.Ce.eg.db")
colData(cds)$garnett_cluster = clusters(cds)
worm_classifier <- train_cell_classifier(cds = cds,
marker_file = "./marker_file.txt",
db=org.Ce.eg.db::org.Ce.eg.db,
cds_gene_id_type = "ENSEMBL",
num_unknown = 50,
marker_file_gene_id_type = "SYMBOL",
cores=8)
# 这样自己就做了一个分组数据集，官方也建议将高质量的训练集汇总到：https://github.com/cole-trapnell-lab/garnett/issues
#
cds = classify_cells(cds, worm_classifier,
db = org.Ce.eg.db::org.Ce.eg.db,
cluster_extend = TRUE,
cds_gene_id_type = "ENSEMBL")
plot_cells(cds,
group_cells_by="partition",
color_cells_by="cluster_ext_type")
方法二：使用别人提供的分组数据集
例如作者提供了已经做好的数据集（https://cole-trapnell-lab.github.io/garnett/classifiers/ceWhole），我们直接下载就能使用：
load(ceWhole)
cds = classify_cells(cds, ceWhole,
db = org.Ce.eg.db,
cluster_extend = TRUE,
cds_gene_id_type = "ENSEMBL")
点击底部的“阅读原文”，获得更好的阅读体验哦😻
初学生信，很荣幸带你迈出第一步。
我们是生信星球，一个不拽术语、通俗易懂的生信知识平台。由于是2018年新号，竟然没有留言功能。需要帮助或提出意见请后台留言、联系微信或发送邮件到jieandze1314@gmail.com，每一条都会看到的哦~