前一阵子练习了几个单细胞测序分析，然而对于所用的R包并没有深入的学习。作为强迫症晚期患者，想知其然，还想知其所以然，所以打算仔细的学习这些R包的使用。这篇文章将深入的学习Monocle2这个包，官方网站http://cole-trapnell-lab.github.io/monocle-release/docs/#introduction，我主要是翻译一下这个包的使用说明和注意事项，一般来说如果你想深入了解一个包，看官网是最靠谱的方法了。里面一些“废话”我就不翻译了，只翻译必要的。以便以后自己分析数据时使用，同时也分享给懒得看英文的童鞋们～
Abstract
单细胞基因表达研究使得我们可以研究在复杂的生物过程和高度异质性细胞群里的转录调节，这些研究有助于基因的发现，鉴定某些细胞的类型，或者标记生物过程的中间态，和两种细胞“命运”的分支。在许多单细胞研究中，细胞个体的基因调节的执行过程并不同步。事实上，每个细胞都是转录过程中的一个“瞬间”。Monocle包提供了一个工具用来分析单细胞测序结果。Monocle引入了一个“伪时间”排序单个细胞的策略，利用细胞个体的异步性把细胞按照一个轨道排列，比如说细胞分化。Monocle利用机器学习的技术(Reversed Graph Embedding)来对细胞进行排序，这个方法能够可靠而准确地解决复杂的生物学过程。Monocle也可以对细胞进行聚类分析。Monocle还可以在聚类后进行差异分析，可以鉴定在不同阶段和不同细胞类型的差异基因。Monocle是专门为单细胞测序分析设计的包，但是也可以用来做其他的分析。Monocle2可以分析非常大的单细胞测序数据（数万个细胞甚至更多）。
Introduction
这篇文章主要提供一个基于Monocle的单细胞测序分析的流程。
Monocle可以帮助你进行三个主要类型的分析：
- 聚类，分类，细胞计数。Monocle可以帮助你鉴定新的细胞类型。
- 构建单细胞轨迹。在发育、疾病以及整个生命过程里，细胞从一个状态转成另一个状态。Monocle可以帮助你发现这些转变。
- 差异表达分析。
在使用Monocle前，先来看一下如何安装Monocle。
安装 Monocle
Monocle运行在R语言环境下。你需要R3.4版本或者更高的版本、Bioconductor V3.5版本， 在这个网址安装Bioconductor:
> source("http://bioconductor.org/biocLite.R")
> biocLite()
安装了Bioconductor，就可以准备安装Monocle和所有需要的依赖包了。
> biocLite("monocle")
检查你的Monocle是否安装正确了，打开你的R，并且打如下代码：
> library(monocle) #实际上就是调用一下
安装最近的版本：
我们推荐你使用Bioconductor安装我们最近的版本。你也可以通过GitHub安装，在R里输入如下代码：
> install.packages("devtools")
> devtools::install_github("cole-trapnell-lab/monocle-release@develop")
有时我们添加了一些features需要你安装额外的包。如果你用上述代码时有报错，你可以安装报错中提示你缺失的包来解决问题。比如说：
> biocLite(c("DDRTree", "pheatmap"))
分析protocol
你不需要知道Monocle包所有的分析功能，有些步骤你可以有很多种方法去做。流程分为多个步骤，我们会标记出来哪些步骤你可以用其他的方法来做。像这样的标记：
流程步骤
Workflow一览（这里先简要的介绍有哪些步骤，下面会具体的把每一步的都详细的讲解）
（一）储存数据到一个CellDataSet对象
第一步就是load你的数据到Monocle的 CellDataSet里:
> pd <- new("AnnotatedDataFrame", data = sample_sheet)
> fd <- new("AnnotatedDataFrame", data = gene_annotation)
> cds <- newCellDataSet(expr_matrix, phenoData = pd, featureData = fd)
（二）用已知marker将细胞分类
利用你对关键marker基因的知识将细胞分类（NOTE:当然有时候你的marker不是很特异或者很难分的很清晰，后面会详细讲解分4种情况来进行细胞分类）：
> cth <- newCellTypeHierarchy()
> MYF5_id <- row.names(subset(fData(cds), gene_short_name == "MYF5"))
> ANPEP_id <- row.names(subset(fData(cds), gene_short_name == "ANPEP"))
> cth <- addCellType(cth, "Myoblast", classify_func =
function(x) { x[MYF5_id,] >= 1 })
> cth <- addCellType(cth, "Fibroblast", classify_func =
function(x) { x[MYF5_id,] < 1 & x[ANPEP_id,] > 1 } )
> cds <- classifyCells(cds, cth, 0.1)
（三）细胞聚类
> cds <- clusterCells(cds)
（四）沿轨迹在伪时间内对细胞进行排序
现在，把你的细胞按照你所研究的生物过程来排序，例如分化，重编程，或者棉衣应答。
> disp_table <- dispersionTable(cds)
> ordering_genes <- subset(disp_table, mean_expression >= 0.1)
> cds <- setOrderingFilter(cds, ordering_genes)
> cds <- reduceDimension(cds)
> cds <- orderCells(cds)
（五）差异分析
这一步就有很多种方法寻找差异基因和处理批次之间的效应。比如说：
> diff_test_res <- differentialGeneTest(cds,
fullModelFormulaStr = "~Media")
> sig_genes <- subset(diff_test_res, qval < 0.1)
开始Monocle（详细讲解）：
Monocle采用Cufflinks或者包来计算基因表达矩阵。Monocle可以处理相对表达值（比如FPKM或者TPM），也可以处理绝对转录数（UMI实验）。Monocle还可以直接使用CellRanger（专门用于10×Genomics单细胞测序）处理得到的转录矩阵。Monocle还可以很好地处理来自其他RNA-Seq工作流程(如sci-RNA-Seq)和Biorad ddSEQ等仪器的数据。尽管Monocle可以用原始counts分析，但除非你按照基因长度对counts进行标准化，否则它们不会与表达值成正比，因此有些Monocle函数可能会产生无意义的结果。如果你没有UMI counts，我们建议你加载FPKM或TPM值进行后续分析，而不是原始counts数。
（一）CellDataSet
Monocle把单细胞表达数据存放在CellDataSet对象里。这个对象源自Bioconductor的ExpressionSet对象，该对象要求三个输入文件：
-
exprs：是一个数值型的表达矩阵，行是基因，列是细胞。 -
phenoData：是一个AnnotatedDataFrame对象,行是细胞，列是细胞属性（比如细胞类型，培养条件，培养天数，等等） 。 -
featureData：是一个AnnotatedDataFrame对象, 行是features（基因），列是基因属性（比如biotype,GC含量，等等）。
表达矩阵必须满足如下要求：
（1）列数与phenoData行数一致。
（2）行数与featureData的data frame行数一致。
另外：
（3）phenoData对象的行名应该与表达矩阵的列名一致。
（4）featureData对象的行名应该与表达矩阵的行名一致。
（5）featureData的其中一列应该命名为"gene_short_name"
你可以按照下面的步骤创建一个CellDataSet对象：
> HSMM_expr_matrix <- read.table("fpkm_matrix.txt")
> HSMM_sample_sheet <- read.delim("cell_sample_sheet.txt")
> HSMM_gene_annotation <- read.delim("gene_annotations.txt")
一旦上面这些table导入之后，你可以创建cellDataSet对象：
> pd <- new("AnnotatedDataFrame", data = HSMM_sample_sheet)
> fd <- new("AnnotatedDataFrame", data = HSMM_gene_annotation)
> HSMM <- newCellDataSet(as.matrix(HSMM_expr_matrix),
phenoData = pd, featureData = fd)
这时将创建一个CellDataSet对象，其中表达值以FPKM表示（Cufflinks方法计算）。默认情况下，Monocle假定你的表达数据以转录本count数为单位，并使用负二项式模型进行下游的差异表达分析步骤。然而，如果你使用的是相对表达值例如TPM或者FPKM数据，请参照下面的内容，如何“告诉”Monocle怎样进行下游的步骤。
如果你的数据里有UMI数据，在创建cellDataSet对象时不应该标准化你的数据。你也不应该用UMI的count去转换相对丰度（比如FPKM/TPM），你也不应该用下面会提到的relative2abs()。monocle内部会做所有需要的标准化步骤。如果你自己标准化可能会对monocle某些关键步骤中断。
（二）从其他包输出文件导入/导出数据
Monocle可以将"Seurat"包的Seurat对象，和"Scater"包的SCESets转换成CellDataSet对象，以便于Monocle的使用。也可以处理Scran的SCESets。要做转换，用importCDS()功能：
# 'data_to_be_imported' 可以是一个Seurat 对象，也可以是一个SCESet
> importCDS(data_to_be_imported)
# 如果我们设置参数'import_all' 为 TRUE，我们就可以把Seurat对象或者SCESet所有的slots导入进来，默认值是FALSE，只导入最小的dataset。
> importCDS(data_to_be_imported, import_all = TRUE)
Monocle还可以导出CellDataSets到Seurat和scater包，用exportCDS()功能：
> lung <- load_lung()
# To convert to Seurat object
> lung_seurat <- exportCDS(lung, 'Seurat')
# To convert to SCESet
> lung_SCESet <- exportCDS(lung, 'Scater')
（三）给你的数据选择一个合适的distribution（分布）（Required！！！）
Monocle可以处理相对表达数据和基于counts的测定(如UMIs)。一般来说，它最适合转录本count数据，特别是UMI数据。无论你的数据类型是什么，为其指定适当的distribution是非常重要的。FPKM/TPM值通常是对数正态分布，而UMIs或counts计数最好使用负二项式来建模。处理counts数据，将负二项分布指定为newCellDataSet的expressionFamily参数:
> HSMM <- newCellDataSet(count_matrix,
phenoData = pd,
featureData = fd,
expressionFamily=negbinomial.size())
expressionFamily参数有几个选项可以选择：
官网说明：请根据你的情况选择正确的expressionFamily参数选项，如果选择错误，可能会导致错误的结果。如果你先用了relative2abs()处理了你的原始数据为FPKM/TPM，你仍然可以使用negative binomial这个选项（但首先你用relative2abs()将相对表达量转换为转录本counts）。会给你更准确的结果。
（四）处理大数据（Recommended）
一些单细胞RNA-Seq实验来自数万个或更多细胞的测量结果。随着仪器的改进和成本的降低，实验将变得越来越庞大，也越来越复杂，包含有许多条件、对照和重复。一个包含50,000个细胞的表达矩阵，并且每个细胞包含有对人类基因组中25,000多个基因的测量就会占用大量的内存。然而，由于目前的protocol通常不会捕获每个细胞中的全部甚至大部分mRNA分子，因此表达矩阵的许多基因为零。使用稀疏矩阵（sparse matrices）可以帮助你在普通的计算机上处理大型数据。我们通常建议大多数用户使用sparseMatrices，因为它可以加速许多计算。
在sparse格式下处理你的数据，可以给Monocle提供一个sparse matrix:
> HSMM <- newCellDataSet(as(umi_matrix, "sparseMatrix"),
phenoData = pd,
featureData = fd,
lowerDetectionLimit = 0.5,
expressionFamily = negbinomial.size())
NOTE：CellRanger输出的一系列RNA-Seq pipelines，已经是sparseMatrix的格式，所以你在做这一步的时候就不应该再用as.matrix()来转换格式了，否则会占用你大量的内存。
如果你的是10X Genomics单细胞测序数据（利用cellrangerRkit）处理的，你可以用下面的代码导入你的数据到Monocle：
> cellranger_pipestance_path <- "/path/to/your/pipeline/output/directory"
> gbm <- load_cellranger_matrix(cellranger_pipestance_path)
> fd <- fData(gbm)
# 这里的2是随便写的，你需要根据你的数据来改相应的数字
# 这个“2”你应该写列数，这个列数与你的featureData的gene short names相一致。
> colnames(fd)[2] <- "gene_short_name"
> gbm_cds <- newCellDataSet(exprs(gbm),
phenoData = new("AnnotatedDataFrame", data = pData(gbm)),
featureData = new("AnnotatedDataFrame", data = fd),
lowerDetectionLimit = 0.5,
expressionFamily = negbinomial.size())
（五）将TPM/FPKM值转换为mRNA counts（Alternative）
如果你的单细胞测序实验里有spike-in control，你可以将这些测量值转换为每个细胞的mRNA的值 (RPC)。RPC值通常比FPKM或TPM值分析起来更容易。事实上，即使没有在实验中加入spike-in control，也有方法可以将FPKM或TPM值转换为RPC值。Monocle 2包含一个名为Census的算法可以执行这种转换。您可以在创建CellDataSet对象之前使用relative2abs()函数将其转换为RPC值，如下所示:
> pd <- new("AnnotatedDataFrame", data = HSMM_sample_sheet)
> fd <- new("AnnotatedDataFrame", data = HSMM_gene_annotation)
# 首先用相对表达水平创建一个CellDataSet
> HSMM <- newCellDataSet(as.matrix(HSMM_expr_matrix),
phenoData = pd,
featureData = fd,
lowerDetectionLimit = 0.1,
expressionFamily = tobit(Lower = 0.1))
# 接下来，转换成RNA的counts
> rpc_matrix <- relative2abs(HSMM, method = "num_genes")
# 现在，可以用RNA counts创建一个新的CellDataSet
> HSMM <- newCellDataSet(as(as.matrix(rpc_matrix), "sparseMatrix"),
phenoData = pd,
featureData = fd,
lowerDetectionLimit = 0.5,
expressionFamily = negbinomial.size())
TPM/FPKM数据和mRNA count数据使用不同的newCellDataSet参数。由于我们使用了Census的mRNA count值，我们已经修改了lowerDetectionLimit这个参数（从0.1改成0.5）来反应新度量的表达量。重要的一点是，我们还修改了expressionFamily参数的选项为negbinomial()，告诉Monocle用负二项分布来进行下游统计分析。如果你不修改这两个参数值，后续可能会产生问题，所以如果你使用Census counts，不要忘记修改这两个参数！！！
（六）估计size factor和分散度（Required！！！）
我们还要用两个函数来pre-calculate一些信息。size facotr帮助我们标准化细胞之间的mRNA的差异。分散度值可以帮助我们进行后续的差异分析。
如果你在上面的步骤里，CellDataSet用的是negbinomial() or negbinomial.size()，那么你只需使用estimateSizeFactors() and estimateDispersions()。
> HSMM <- estimateSizeFactors(HSMM)
> HSMM <- estimateDispersions(HSMM)
（七）过滤低质量的细胞（Recommended）
任何单细胞测序分析的第一步都是要鉴定低质量文库。大部分单细胞处理过程中都会包含一些死细胞，或者空孔。另外去掉doublets也是很重要的：就是有两个甚至更多的细胞的文库。这些细胞会破坏下游步骤，例如伪时间排序和聚类。这部分举例一个比较经典的质控步骤。要看多少细胞里表达某个基因，或者看一个细胞里表达多少基因是非常容易的。Monocle提供了一个非常简单的函数来计算：
> HSMM <- detectGenes(HSMM, min_expr = 0.1)
> print(head(fData(HSMM)))
> HSMM <- detectGenes(HSMM, min_expr = 0.1)
> print(head(fData(HSMM)))
> expressed_genes <- row.names(subset(fData(HSMM),
num_cells_expressed >= 10))
现在这个expressed_genes 里的基因都是至少在10个细胞里有表达的（这里官网写错啦，他写成50了）。我们将在后续生物过程分析里用这个表。你也可以根据自己的需要进行筛选。
你的单细胞测序实验过程中也可以用成像来进行质控（捕获细胞之后，裂解细胞之前）。成像法可以评估你的细胞，确定你的文库里不会包含空孔，或者细胞碎片。一些仪器在捕获细胞时可能会使得你的一个孔里有多于一个的细胞。你应该把这些文库去除掉。特别是空孔和细胞碎片的孔。检查每个细胞的RNA-seq文库的sequence是否达到可接受的程度也是一个过滤的方法。
然而，对于测序的“深度”并没有一个最低标准。但是运用你的判断：如果一个细胞测序只有几千个reads的话，是不太可能产生有意义的结果的。
在CellDataSet对象里，存放每个细胞评分的data在phenoData里。评分属性作为一列。 你可以过滤那些不符合你要求的细胞。你也可以用FastQC来过滤细胞。这类软件可以鉴别RNA-seq文库里严重降解的RNA，也可以鉴别包含有核糖体，线粒体或者其他类型RNA污染的文库。
在我们的例子中，HSMM数据库包含的评分列：
> print(head(pData(HSMM)))
如果你用的是RPC值来测定表达量，也可以根据样本表达量进行过滤：
> pData(HSMM)$Total_mRNAs <- Matrix::colSums(exprs(HSMM))
> HSMM <- HSMM[,pData(HSMM)$Total_mRNAs < 1e6]
> upper_bound <- 10^(mean(log10(pData(HSMM)$Total_mRNAs)) +
2*sd(log10(pData(HSMM)$Total_mRNAs)))
> lower_bound <- 10^(mean(log10(pData(HSMM)$Total_mRNAs)) -
2*sd(log10(pData(HSMM)$Total_mRNAs)))
> qplot(Total_mRNAs, data = pData(HSMM), color = Hours, geom =
"density") +
> geom_vline(xintercept = lower_bound) +
> geom_vline(xintercept = upper_bound)
> HSMM <- HSMM[,pData(HSMM)$Total_mRNAs > lower_bound &
pData(HSMM)$Total_mRNAs < upper_bound]
> HSMM <- detectGenes(HSMM, min_expr = 0.1)
我们还进行了一步过滤，去掉那些mRNA水平非常低或者非常高的细胞。通常来讲，doublets或者triplets是正常细胞mRNA的两倍，所以要把这些细胞也要排除掉。如果你没有条件通过成像法检查你的细胞，那么你可以在这一步来过滤。一般设置的“门槛”是10000和40000 mRNA。你可以根据你的实验需要调整这个门槛。一旦你排除掉不符合你要求的细胞，你应该验证一下储存在CelDataSet里的表达量大致是一个lognormal的分布：
# 先把表达矩阵用log转换一下
> L <- log(exprs(HSMM[expressed_genes,]))
# 标准化每个基因，然后合并你的data
> melted_dens_df <- melt(Matrix::t(scale(Matrix::t(L))))
# 画标准化后的基因表达值的分布
> qplot(value, geom = "density", data = melted_dens_df) +
stat_function(fun = dnorm, size = 0.5, color = 'red') +
xlab("Standardized log(FPKM)") +
ylab("Density")
（二）细胞分类和细胞计数
单细胞实验经常用复杂混合的多细胞类型的样品。从组织里分离出来的细胞往往含有两个、三个、甚至更多不同类型的细胞。在这种情况下，最好是根据细胞类型的特异性marker将细胞分类。在我们的这个例子里，培养的细胞是从肌肉样品来源的（包含成纤维细胞），用于原代细胞培养。成肌细胞表达一些成纤维细胞没有的基因。选择特异性表达的基因，比如说，利用高水平的MYF5可以排除成纤维细胞。另外，成纤维细胞表达高水平的ANPEP(CD13)，而成肌细胞不表达。
有几种细胞分类的方法：
（1）按照细胞类型将细胞分类 （Recommended）
Monocle提供了一个简单的系统来按照你选择的基因来标记细胞。比如说：你可以提供每一个类型细胞的function。这些functions可以是每一个细胞的表达data的input，返回TRUE来告诉Monocle这个细胞满足你定义的这个function。你可以用一个function来返回成肌细胞特异表达基因的细胞TRUR，另一个function为成纤维细胞特异基因，等等。下面是个例子：
# 根据基因名字找到其在表达矩阵的ID
> MYF5_id <- row.names(subset(fData(HSMM), gene_short_name == "MYF5"))
> ANPEP_id <- row.names(subset(fData(HSMM),
gene_short_name == "ANPEP"))
# 这里选取的基因取决于自己的单细胞实验设计
> cth <- newCellTypeHierarchy()
> cth <- addCellType(cth, "Myoblast", classify_func =
function(x) { x[MYF5_id,] >= 1 })
> cth <- addCellType(cth, "Fibroblast", classify_func = function(x)
{ x[MYF5_id,] < 1 & x[ANPEP_id,] > 1 })
> HSMM <- classifyCells(HSMM, cth, 0.1)
classifyCells 功能将所有的gating function去判断每一个细胞是否复合标准，从而进行分类，然后会返回 CellDataSet对象里新的一列,称为 CellType。我们现在可以计数每一个细胞类型里都有多少个细胞。
> table(pData(HSMM)$CellType)
| Fibroblast | Myoblast | Unknown |
| 56 | 85 | 121 |
还可以画图来展示细胞数：
> pie <- ggplot(pData(HSMM),
> aes(x = factor(1), fill = factor(CellType))) + geom_bar(width = 1)
> pie + coord_polar(theta = "y") +
> theme(axis.title.x = element_blank(), axis.title.y = element_blank())
你会发现很多细胞被列为"Unknown"。主要是因为在大部分的单细胞测序实验中，mRNA捕获效率低，比如一个细胞表达很低水平的MYF5，但是我们很不幸没有捕获到它。 当一个细胞不满足你定下的任何一个function标准，就会标记为 "Unknown"。如果它满足多个标准，会被标记为"Ambiguous"。你可以排除这些细胞，但是你会丢掉你的很多data，在这个例子里，你会损失大概一半的细胞。
（2）不根据marker基因将细胞分类 （Alternative）
Monocle提供一种算法，你可以用来划分你的"Unknown"的细胞。这种算法可以通过clusterCells来实现, 根据整体的表达量来把细胞分组。比方说，如果你的细胞表达大量的成肌细胞基因，但是缺少MYF5，我们仍然可以把它识别为成肌细胞。 clustercell可以以无监督的方式使用，也可以以半监督的方式使用。我们先来看看无监督模式：
第一步是确定哪些基因被用来划分细胞群。我们可以用所有的基因，但是很多基因没有足够高的表达量。加上这些基因只会增加背景噪音。我们可以根据平均表达水平过滤基因，我们还可以另外选取那些在细胞之间不经常变动的基因。这些基因能够对细胞状态提供很有效的信息。
> disp_table <- dispersionTable(HSMM)
> unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)
> HSMM <- setOrderingFilter(HSMM, unsup_clustering_genes$gene_id)
> plot_ordering_genes(HSMM)
setOrderingFilter函数会标记基因（在后续clusterCells细胞聚类所用的基因）。 尽管我们也可以提供一个我们自己写的基因列表。plot_ordering_genes功能根据细胞之间的平局表达量来显示一个基因的变异性 。红线显示的是Monocle基于这种关系的预期值。我们把用于细胞分群的基因标记为黑点，其他基因标记为灰点。
下面我们尝试给细胞聚类：
# HSMM@auxClusteringData[["tSNE"]]$variance_explained <- NULL
# 这里看看基因的表达量和基因的变异度之间的关系
> plot_pc_variance_explained(HSMM, return_all = F) # norm_method='log'
> HSMM <- reduceDimension(HSMM, max_components = 2, num_dim = 6,
reduction_method = 'tSNE', verbose = T)
> HSMM <- clusterCells(HSMM, num_clusters = 2)
# 这里先用tSNE的聚类方法处理HSMM数据集，并可视化展示
> plot_cell_clusters(HSMM, 1, 2, color = "CellType",
markers = c("MYF5", "ANPEP"))
Monocle使用t-SNE来给细胞聚类，使用的方法与Seurat包非常相似。
上图里的根据我们提供的gating functions被标记为成肌细胞的用绿色表示，成纤维细胞标记为红色。不表达任何marker的细胞标记为蓝色。在许多实验中，细胞可以被很清楚的划分，但是在这个实验里，细胞聚类并不容易。绿色和红色的细胞并没有清晰的分界线，是由于成肌细胞和混合进来的间质的成纤维细胞表达许多同样的基因（在这个培养条件下）。还有其他的一些因素可以影响分群。实验中的一个影响因素源于实验设计。为了让成肌细胞分化，我们将培养基从高丝裂原生长培养基(GM)转换为低丝裂原分化培养基(DM)。也许这些细胞是基于它们所培养的培养基而聚集在一起的?可以试一下用培养基来分群：
> plot_cell_clusters(HSMM, 1, 2, color = "Media")
Monocle允许我们去掉“不感兴趣”来源变异的影响，从而降低这些因素对聚类的影响。你可以在clusterCells里添加一个参数residualModelFormulaStr，这个参数接受R字符串模型，指定要在聚类之前减去的效果。
# 我们假设只有2种细胞类型，所以在做聚类的时候可以把这个参数添加进去，这样可以去除无关变量的干扰。
> HSMM <- reduceDimension(HSMM, max_components = 2, num_dim = 2,
reduction_method = 'tSNE',
residualModelFormulaStr = "~Media + num_genes_expressed",
verbose = T)
> HSMM <- clusterCells(HSMM, num_clusters = 2)
> plot_cell_clusters(HSMM, 1, 2, color = "CellType")
既然我们已经解释了一些不想要的变异来源，我们准备再尝试一次用无监督聚类来分类细胞：
>HSMM <- clusterCells(HSMM, num_clusters = 2)
>plot_cell_clusters(HSMM, 1, 2, color = "Cluster") +
facet_wrap(~CellType)
现在，大部分的成肌细胞在一个群里，大部分成纤维细胞在另一个群里，然而我们仍然看到有些细胞在两个群里都有，这可能是因为我们的marker基因和CellTypeHierarchy的functions缺乏特异性，但也可能是因为聚类不够理想。为了排除后者，让我们尝试以其半监督模式运行clusterCells。
（3）根据marker基因聚类细胞（Recommended）
首先，我们选择一个不同的基因列表来对细胞聚类。之前我们是选取了那些高表达并且高变异性的基因。现在我们选取那些与我们的marker基因共变的基因。在某种意义上，我们建立一个大的基因列表作为marker，以便于即使一个细胞不表达MYF5，但是它可能表达其他的基因，而被识别为成肌细胞。这种聚类也叫“半监督聚类”。
> marker_diff <- markerDiffTable(HSMM[expressed_genes,],
cth,
residualModelFormulaStr = "~Media + num_genes_expressed",
cores = 1)
markerDiffTable功能可以把所有细胞按照你提供的functions进行分类。 然后它会在识别出不同类型间差异表达的基因之前去掉所有的"Unknown" 和"Ambiguous" functions。该函数返回一个data frame，你可以用这个去挑选用来聚类的基因。通常情况下，最好选择每种细胞类型最具特异性的前10或20个基因。这确保了聚类基因不会被某一种细胞类型的marker所主导。Monocle提供了一个方便的功能，可以根据每种基因的表达受限制程度对基因进行排序。
# 对每个基因增加了pval和qval两列信息，挑选出那些在不同media培养条件下显著差异表达的基因
> candidate_clustering_genes <-
row.names(subset(marker_diff, qval < 0.01))
# 计算这些基因在不同的celltype的specificity值
> marker_spec <-
calculateMarkerSpecificity(HSMM[candidate_clustering_genes,], cth)
> head(selectTopMarkers(marker_spec, 3))
上面最后一行显示了成肌细胞和成纤维细胞的top3的marker gene。 "specificity"评分范围可以从0到1。它越接近1，对细胞类型的限制就越多。你可以用它为已知的细胞类型定义新的marker，也可以挑选出那些你用来纯化新发现的细胞类型的基因。这对于后续的实验非常有价值。
要将这些细胞聚在一起，我们将为每种细胞类型选择前500个标记:
> semisup_clustering_genes <- unique(selectTopMarkers(marker_spec, 500)$gene_id)
> HSMM <- setOrderingFilter(HSMM, semisup_clustering_genes)
> plot_ordering_genes(HSMM)
需要注意的是，我们有一小组基因，它们其中一些在实验中并没有特别高的表达。然而，它们对于区分表达MYF5和表达ANPEP的细胞非常有用。我们已经将它们在聚类中使用，但是即使我们不这样做，我们仍然可以将它们直接提供给clusterCells来使用它们。
# 重新挑选基因，用黑色点基因来进行聚类。
> plot_pc_variance_explained(HSMM, return_all = F)
> HSMM <- reduceDimension(HSMM, max_components = 2, num_dim = 3,
norm_method = 'log',
reduction_method = 'tSNE',
residualModelFormulaStr = "~Media + num_genes_expressed",
verbose = T)
> HSMM <- clusterCells(HSMM, num_clusters = 2)
> plot_cell_clusters(HSMM, 1, 2, color = "CellType")
（4）Imputing cell type （Alternative）（这个title不知道怎么翻译才到位）
需要注意，我们已经减少了成肌细胞群中“污染”的成纤维细胞的数量。但是那些“unknown”的细胞怎么处理?如果你用CellTypeHierarcy来提供clusterCells，Monocle将使用它作为一个整体进行分类，而不仅仅是对单个细胞处理。clustercell计算每个细胞类型在每个聚类里的频率。当一个聚类由某个细胞类型的10%(在本例中为10%)以上组成时，这个聚类中的所有细胞都设置为该类型。如果一个聚类里由多个细胞类型组成，那么整个集群将被标记为“Ambiguous”。如果没有细胞类型符合上述的threshold，则这个聚类被标记为“Unknown”。因此，Monocle可以帮助你归类每个细胞的类型，即使没有现存的marker。
> HSMM <- clusterCells(HSMM,
num_clusters = 2,
frequency_thresh = 0.1,
cell_type_hierarchy = cth)
> plot_cell_clusters(HSMM, 1, 2, color = "CellType",
markers = c("MYF5", "ANPEP"))
可以看出，根据MYF5的表达可以清晰的把细胞聚类。在两个聚类里都有一些细胞表达ANPEP，但是在成肌细胞群里还表达MYF5。这并不令人惊讶，因为ANPEP并不是成纤维细胞的一个特别特异性的marker基因。总之，我们可以成功的把这些细胞分群：
> pie <- ggplot(pData(HSMM), aes(x = factor(1), fill =
factor(CellType))) +
geom_bar(width = 1)
> pie + coord_polar(theta = "y") +
theme(axis.title.x = element_blank(), axis.title.y = element_blank())
构建单细胞轨迹
在发育过程中，细胞应对外界的刺激，从一种“状态”转变为另一种“状态”。细胞在不同的“状态”下表达着不同的基因，产生不同的蛋白行使各自的功能。随着细胞在不同“状态”间转换，其自身也经历着转录水平的重构，一些基因沉默，而一些基因被激活。这些瞬时的“状态”很难被鉴定，因为在更稳定的两个端点状态之间纯化细胞几乎是不可能的。单细胞RNA-Seq可以让你在不需要纯化细胞的情况下看到这些状态。然而，为了做到这一点，我们必须确定每个细胞可能的“状态”范围。
Monocle不是通过实验将细胞纯化成离散的状态，而是使用一种算法来“学习”每个细胞必须经历的基因表达变化的顺序，这是动态生物学过程的一部分。一旦它了解了基因表达变化的整体“轨迹”，Monocle就可以将每个细胞置于其轨迹的适当位置。你就可以使用Monocle的differential analysis toolkit来查找在轨迹过程中基因的调节。
如果该分析过程有多个结果，Monocle将重建一个“分支”轨迹。这些分支对应着细胞的“命运决定”，Monocle提供了强大的工具来识别受影响的基因。后面我们会讲到如何分析分支。Monocle依靠一种称为“反向图嵌入”的机器学习技术来构建单细胞轨迹。
什么是伪时间？
伪时间是衡量一个细胞在生物过程中（比如分化）进程的一个方法。在许多生物学过程中，细胞并不是完全同步进行的。在细胞分化等过程的单细胞表达研究中，捕获的细胞分布在不同的过程里。也就是说，即使在同一时间捕获的细胞群中，有些细胞可能已经进入某个生物过程很久了，而有些细胞甚至还没有开始这个过程。当您想要了解细胞从一种状态过渡到另一种状态时所发生的调节变化的顺序时，这种异步性会产生很大的问题。跟踪同时捕获的细胞的基因表达会产生一个非常压缩的基因动力学，该基因表达的变异性将非常高。Monocle根据每个细胞在轨迹上的进展对其进行排序，从而降低了异步带来的问题。Monocle不是随时间来追踪变化，而是沿着轨迹来追踪的，我们称之为伪时间。伪时间是一个抽象的进程单位：它只是一个细胞到轨迹起点的距离，沿着最短路径测量。轨迹的总长度是根据一个细胞从起始状态移动到结束状态所经历的总转录水平的变化量来定义的。
分析流程：
在我们深入了解沿轨迹排列细胞的细节之前，先来了解一下Monocle在做什么。有三个主要步骤，每个步骤都涉及一个重要的机器学习。
Step 1: 选择定义过程的基因
推断一个单细胞轨迹是一个机器学习的过程。第一步就是选择基因，Monocle将把这些基因作为其机器学习方法的input。这个过程叫做“feature selection”，它对轨迹的形状有很大的影响。在单细胞RNA-Seq中，低水平表达的基因通常非常嘈杂，但某些基因包含着有关细胞状态的重要信息。Monocle通过检查这些基因在细胞群之间的表达模式来对细胞进行排序。Monocle寻找那些以“有趣的”(即不只是嘈杂)方式变化的基因，并利用它们来构建数据。Monocle提供了多种工具来选择基因，并产生一个准确、具有生物学意义的轨迹。你可以用这些工具进行完全的“无监督”分析，或者你也可以用你的专业知识选择基因，从而形成Monocle的轨迹。我们称为这种方式为“半监督”。
Step 2: 数据降维
一旦我们选择了用来给细胞排序的基因，Monocle就会对数据进行降维处理。Monocle使用最近开发的一种称为反向图嵌入的算法来降低数据的维数。
Step 3:按伪时间排序细胞
将表达数据投射到低维空间后，Monocle就准备学习描述细胞如何从一种状态过渡到另一种状态的轨迹。Monocle假设轨迹是“树状”结构的，一端是“根”，另一端是“叶”。Monocle的工作就是尽可能地把最佳的“树状结构”与数据相匹配。这项任务被称为“manifold learning”。一个细胞从生物过程的开始阶段（从根部开始），沿着“主干”前进，直到它到达第一个分支(如果有的话)。然后，这个细胞必须选择一条分支，并且沿着它“走”得越来越远，直到到达一片叶子。一个细胞的伪时间值是它从这片“叶子”返回到“根”的距离。
（下面是对以上三个步骤的详细介绍）
Trajectory step 1: 选择定义细胞过程的基因
比如：我们必须首先要决定用哪些基因来定义细胞的肌生成过程（这里官网开始用他的例子来具体的说明）。我们想要的是一组基因，它们的表达量随着我们研究的生物过程的进展而增加(或减少)。
理想情况下，我们想尽可能少的使用系统生物学的先验知识。我们还想从数据中发现一些重要的排序基因，而不是仅仅依赖于文献和教科书，因为那样的话可能会在排序中引入偏差。我们通常推荐一种更成熟的方法，称为“dpFeature”。
分离排序基因的一种有效方法是比较生物过程开始时和结束时收集的细胞，并寻找差异表达基因。下面的代码就是为了找到从生长培养基到分化培养基的转换中产生差异表达的基因：
> diff_test_res <- differentialGeneTest(HSMM_myo[expressed_genes,],
fullModelFormulaStr = "~Media")
> ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
根据时间点的差异分析来选择基因通常是非常有效的，但是如果我们没有时间排列的数据该怎么办呢?如果这些细胞在我们所研究的生物过程中是不同步的，Monocle通常可以从同时捕获的单个细胞群中重建它们的轨迹。)
一旦我们有了用于排序的基因id列表，我们需要在HSMM对象中设置它们，因为接下来的几个函数都将用到它们：
> HSMM_myo <- setOrderingFilter(HSMM_myo, ordering_genes)
> plot_ordering_genes(HSMM_myo)
Trajectory step 2: 降维
接下来我们把空间降到2维，便于可视化：
> HSMM_myo <- reduceDimension(HSMM_myo, max_components = 2,
method = 'DDRTree')
Trajectory step 3:按照轨迹排序细胞
> HSMM_myo <- orderCells(HSMM_myo)
可视化：
> plot_cell_trajectory(HSMM_myo, color_by = "Hours")
轨迹是树状结构的。我们可以看到，在0时间点收集的细胞位于树的一个尖端附近，而其他细胞分布在两个“分支”上。Monocle并不知道哪条分支作为“开始”，所以我们经常需要使用root_state参数调用ordercell来指定“起点”。首先，我们绘制轨迹，这次按“状态”给细胞上色:
> plot_cell_trajectory(HSMM_myo, color_by = "State")
“State”是Monocle对象的一部分。下面的函数“手动”的定义细胞的“State”，然后将其传递给ordercell：
#实际上就是自己写个函数
> GM_state <- function(cds){
if (length(unique(pData(cds)$State)) > 1){
T0_counts <- table(pData(cds)$State, pData(cds)$Hours)[,"0"]
return(as.numeric(names(T0_counts)[which
(T0_counts == max(T0_counts))]))
} else {
return (1)
}
}
> HSMM_myo <- orderCells(HSMM_myo, root_state = GM_state(HSMM_myo))
> plot_cell_trajectory(HSMM_myo, color_by = "Pseudotime")
如果你画出来的“树”里有大量的"state“，那么很难确定每个”state“落在”树“的位置。你可以将每一个"state"单独画出来，比如说：
> plot_cell_trajectory(HSMM_myo, color_by = "State") +
facet_wrap(~State, nrow = 1)
如果你没有一系列的时间点，你可以根据marker基因的表达来设置”树根“的位置，例如，在这个实验中，一个高度增殖的祖细胞群产生两种有丝分裂后的细胞。所以“树根”所在位置的细胞应该有高水平的增殖marker表达。我们可以使用jitter plot来找出快速增殖所对应的"state":
> blast_genes <- row.names(subset(fData(HSMM_myo),
#比如这里作者就选了3个marker基因
> gene_short_name %in% c("CCNB2", "MYOD1", "MYOG")))
> plot_genes_jitter(HSMM_myo[blast_genes,],
grouping = "State",
min_expr = 0.1)
为了确认排序是正确的，我们可以选择几个肌源性进展的marker来画一下图：
> HSMM_expressed_genes <- row.names(subset(fData(HSMM_myo),
num_cells_expressed >= 10))
> HSMM_filtered <- HSMM_myo[HSMM_expressed_genes,]
> my_genes <- row.names(subset(fData(HSMM_filtered),
gene_short_name %in% c("CDK1", "MEF2C", "MYH3")))
> cds_subset <- HSMM_filtered[my_genes,]
> plot_genes_in_pseudotime(cds_subset, color_by = "Hours")
其他的方法选择排序基因
（1）根据聚类的差异基因排序（Recommended）
我们建议使用“dpFeature”的非监督过程选择基因。
为了使用dpFeature，我们首先要挑选出至少在5%的细胞里都表达的基因：
> HSMM_myo <- detectGenes(HSMM_myo, min_expr = 0.1)
> fData(HSMM_myo)$use_for_ordering <-
fData(HSMM_myo)$num_cells_expressed > 0.05 * ncol(HSMM_myo)
然后，我们将执行PCA分析，以确定每个主成分的variance。我们可以根据散点图来确定需要多少个pca维度。
> plot_pc_variance_explained(HSMM_myo, return_all = F)
然后，我们将使用t-SNE降维（运行上图位于顶部的variance），并将它们进一步投射到二维上。
> HSMM_myo <- reduceDimension(HSMM_myo,
max_components = 2,
norm_method = 'log',
num_dim = 3,
reduction_method = 'tSNE',
verbose = T)
然后在二维t-SNE空间中进行密度peak聚类来鉴定聚类。densityPeak算法根据每个细胞的局部密度(Ρ)和两个细胞之间最近的距离(Δ)进行细胞聚类。默认情况下,clusterCells选择Ρ和Δ定义阈值的95%。我们还可以设置一些想要进行聚类的数量(n)。默认设置通常可以比较好的聚类。
> HSMM_myo <- clusterCells(HSMM_myo, verbose = F)
> plot_cell_clusters(HSMM_myo, color_by = 'as.factor(Cluster)')
> plot_cell_clusters(HSMM_myo, color_by = 'as.factor(Hours)')
我们也可以根据自定义的阈值重新运行聚类：
> HSMM_myo <- clusterCells(HSMM_myo,
rho_threshold = 2,
delta_threshold = 4,
skip_rho_sigma = T,
verbose = F)
> plot_cell_clusters(HSMM_myo, color_by = 'as.factor(Cluster)')
> plot_cell_clusters(HSMM_myo, color_by = 'as.factor(Hours)')
在确定聚类是make sense的情况下，就可以提取差异基因了：
> clustering_DEG_genes <-
differentialGeneTest(HSMM_myo[HSMM_expressed_genes,],
fullModelFormulaStr = '~Cluster',
cores = 1)
我们选取top1000作为排序基因：
> HSMM_ordering_genes <-
row.names(clustering_DEG_genes)[order(clustering_DEG_genes$qval)][1:1000]
> HSMM_myo <-
setOrderingFilter(HSMM_myo,
ordering_genes = HSMM_ordering_genes)
> HSMM_myo <-
reduceDimension(HSMM_myo, method = 'DDRTree')
> HSMM_myo <-
orderCells(HSMM_myo)
> HSMM_myo <-
orderCells(HSMM_myo, root_state = GM_state(HSMM_myo))
> plot_cell_trajectory(HSMM_myo, color_by = "Hours")
（2）选择在细胞间高度变化的基因用来排序（Alternative）
变化很大的基因通常为鉴定细胞亚群或沿轨迹对细胞进行排序提供了大量信息。在RNA-Seq中，一个基因的variance通常取决于它的均值，所以要小心地根据它们的variance选择基因。
> disp_table <- dispersionTable(HSMM_myo)
> ordering_genes <- subset(disp_table,
mean_expression >= 0.5 &
dispersion_empirical >= 1 * dispersion_fit)$gene_id
（3）用已知marker基因进行排序（Alternative）
无监督排序是为了避免在分析中引入偏差。然而，非监督的机器学习有时会focus在一个feature上，而这并不是你实验的重点。例如，当你使用无监督学习时，每个细胞在细胞周期中的位置对轨迹的形状有很大的影响。但是，如果你想研究的生物过程与周期无关该怎么办呢?这时就需要用Monocle的“半监督”排序模式。
以半监督的方式对细胞进行排序非常简单。首先使用CellTypeHierchy定义marker基因，这与我们之前使用它进行细胞类型分类的方式非常相似。然后，用它来选择与这些marker共变的排序基因。最后，根据这些基因对细胞进行排序。因此，无监督和半监督排序之间的唯一区别就是我们用来排序的基因。
正如我们之前看到的，成肌细胞退出细胞周期开始分化，然后通过一系列调控事件使得肌肉收缩所需的一些关键肌肉特异性蛋白的表达。我们可以用cyclin B2 (CCNB2)标记在细胞周期内的细胞，并识别肌管，因为这些细胞表达了高水平的肌球蛋白重链3 (MYH3)。
> CCNB2_id <-
row.names(subset(fData(HSMM_myo), gene_short_name == "CCNB2"))
> MYH3_id <-
row.names(subset(fData(HSMM_myo), gene_short_name == "MYH3"))
> cth <- newCellTypeHierarchy()
> cth <- addCellType(cth,
"Cycling myoblast",
classify_func = function(x) { x[CCNB2_id,] >= 1 })
> cth <- addCellType(cth,
"Myotube",
classify_func = function(x) { x[MYH3_id,] >= 1 })
> cth <- addCellType(cth,
"Reserve cell",
classify_func =
function(x) { x[MYH3_id,] == 0 & x[CCNB2_id,] == 0 })
> HSMM_myo <- classifyCells(HSMM_myo, cth)
现在我们选择与上面的两个基因共变的一些基因：
> marker_diff <- markerDiffTable(HSMM_myo[HSMM_expressed_genes,],
cth,
cores = 1)
#下面这行代码也可以写成这样：semisup_clustering_genes <- row.names(subset(marker_diff, qval < 0.05))
> semisup_clustering_genes <-
row.names(marker_diff)[order(marker_diff$qval)][1:1000]
使用top1000个基因排序产生的轨迹与我们用非监督方法得到的轨迹非常相似，但它更“干净”一些:
> HSMM_myo <- setOrderingFilter(HSMM_myo, semisup_clustering_genes)
> HSMM_myo <- reduceDimension(HSMM_myo, max_components = 2,
method = 'DDRTree', norm_method = 'log')
> HSMM_myo <- orderCells(HSMM_myo)
> HSMM_myo <- orderCells(HSMM_myo, root_state = GM_state(HSMM_myo))
> plot_cell_trajectory(HSMM_myo, color_by = "CellType") +
theme(legend.position = "right")
为了确认排序是正确的，我们可以选择几个肌生成过程的marker。在这个实验中，一个分支对应于成功融合形成肌管的细胞，另一个分支对应于未能完全分化的细胞。现在我们排除后者，只把前者画出来：
> HSMM_filtered <- HSMM_myo[HSMM_expressed_genes,]
> my_genes <- row.names(subset(fData(HSMM_filtered),
gene_short_name %in% c("CDK1", "MEF2C", "MYH3")))
> cds_subset <- HSMM_filtered[my_genes,]
> plot_genes_branched_pseudotime(cds_subset,
branch_point = 1,
color_by = "Hours",
ncol = 1)
差异分析
差异分析是RNA-Seq分析中常见的任务。Monocle可以帮助你找到不同细胞组间差异表达的基因，并评估这些变化的显著性。首先需要把你的细胞分到到两个或更多的组里。这些组由每个CellDataSet的phenoData的列定义。
（一）最基本的差异分析
对人类基因组中的所有基因进行差异表达分析需要相当长的时间。对于我们举的例子这样大的数据集(其中包含数百个细胞)，在单个CPU上进行分析可能需要几个小时。所以先选择一小组我们已知在肌生成中很重要的基因来试一下:
> marker_genes <- row.names(subset(fData(HSMM_myo),
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")))
在成肌细胞数据中，实验开始时收集的细胞在“生长培养基”(GM)中培养，以防止它们分化。在收集细胞后，其余细胞转入“分化培养基”(DM)以促进分化。首先来用Monocle寻找一下哪些基因可能受上述培养基转换的影响：
> diff_test_res <- differentialGeneTest(HSMM_myo[marker_genes,],
fullModelFormulaStr = "~Media")
# 根据 FDR < 10%的基因
> sig_genes <- subset(diff_test_res, qval < 0.1)
> sig_genes[,c("gene_short_name", "pval", "qval")]
Monocle还提供了一些简单的方法来绘制一小组基因的表达，这些基因是根据你在差异分析中使用的因素分组的。这有助于你可视化上述的差异。
> MYOG_ID1 <- HSMM_myo[row.names(subset(fData(HSMM_myo),
gene_short_name %in% c("MYOG", "CCNB2"))),]
> plot_genes_jitter(MYOG_ID1, grouping = "Media", ncol= 2)
这里要注意的是，我们可以控制在上图里怎么样排布你的基因，包括行数和列数。请参阅
plot_genes_jitter手册。
（二）根据细胞类型和状态寻找差异基因
在动态的生物过程中，例如分化，细胞可能呈现出不同的中间态或最终态。我们之前根据几个关键的marker来区分成肌细胞和污染的成纤维细胞。现在我们尝试其他几个也可能区分成纤维细胞和成肌细胞的基因：
> to_be_tested <- row.names(subset(fData(HSMM),
gene_short_name %in% c("UBC", "NCAM1", "ANPEP")))
> cds_subset <- HSMM[to_be_tested,]
#下面这个differentialGeneTest函数是为了验证我们根据细胞类型所寻找的差异基因是不是具有显著差异性
> diff_test_res <- differentialGeneTest(cds_subset,
fullModelFormulaStr = "~CellType")
> diff_test_res[,c("gene_short_name", "pval", "qval")]
> plot_genes_jitter(cds_subset,
grouping = "CellType",
color_by = "CellType",
nrow= 1,
ncol = NULL,
plot_trend = TRUE)
实际上Monocle的差异分析步骤非常的灵活多变，你可以用pData表里任何存在的列来进行分析。比如说，你之前用的是clusterCells来分类你的细胞，那么在差异分析时，你就可以用pData里Cluster这列来寻找差异基因。
（三）根据伪时间功能寻找差异基因
Monocle的主要工作是将细胞按照生物过程(如细胞分化)的顺序进行排列。这样你就可以分析细胞里找到随着生物过程进展而改变的基因。比如说，你可以发现当细胞“成熟”时有哪些基因显著上调。我们看一下一组对肌生成很重要的基因：
> to_be_tested <- row.names(subset(fData(HSMM),
> gene_short_name %in% c("MYH3", "MEF2C", "CCNB2", "TNNT1")))
> cds_subset <- HSMM_myo[to_be_tested,]
> diff_test_res <- differentialGeneTest(cds_subset,
fullModelFormulaStr = "~sm.ns(Pseudotime)")
sm.ns函数是通过表达值来拟合一个曲线，以帮助它将基因表达中的变化描述为一个随生物过程变化的函数。
> diff_test_res[,c("gene_short_name", "pval", "qval")]
> plot_genes_in_pseudotime(cds_subset, color_by = "Hours")
（四）根据伪时间表达pattern聚类基因（热图）
Monocle提供了一个简单的方法可视化所有伪时间依赖的基因。
plot_pseudotime_heatmap接受 CellDataSet 对象（通常包含一小组差异基因）。然后把这些基因聚类，并用pheatmap包画热图。
> diff_test_res <- differentialGeneTest(HSMM_myo[marker_genes,],
fullModelFormulaStr = "~sm.ns(Pseudotime)")
> sig_gene_names <- row.names(subset(diff_test_res, qval < 0.1))
> plot_pseudotime_heatmap(HSMM_myo[sig_gene_names,],
num_clusters = 3,
cores = 1,
show_rownames = T)
多因素差异分析
Monocle可以在多个因素存在的情况下进行差异分析，可以帮助你减去一些不必要的影响因素。在下面的简单例子中，Monocle测试在成肌细胞和成纤维细胞之间的差异表达的3个基因，同时减去Hours的影响（Hours指每个细胞收集的日期）。
> to_be_tested <-
row.names(subset(fData(HSMM),
gene_short_name %in% c("TPM1", "MYH3", "CCNB2", "GAPDH")))
> cds_subset <- HSMM[to_be_tested,]
> diff_test_res <- differentialGeneTest(cds_subset,
fullModelFormulaStr = "~CellType + Hours",
reducedModelFormulaStr = "~Hours")
> diff_test_res[,c("gene_short_name", "pval", "qval")]
> plot_genes_jitter(cds_subset,
grouping = "Hours", color_by = "CellType", plot_trend = TRUE) +
facet_wrap( ~ feature_label, scales= "free_y")
单细胞轨迹的“分支”分析
通常，单细胞的轨迹里有分支。分支的产生是因为细胞有可选择的基因表达模式。比如在发育过程中，当细胞做出“命运”的选择时，分支就会出现在轨迹中:一个发育谱系沿着一条路径前进，而另一个谱系产生第二条路径。Monocle可以分析这些分支事件。用BEAM函数。
> lung <- load_lung()
> plot_cell_trajectory(lung, color_by = "Time")
BEAM的输入对象是排序好的CellDataSet对象。你需要用orderCells 和轨迹分支点的名字来排序你的CellDataSet对象。它会返回一个每个基因的显著值的table。
> BEAM_res <- BEAM(lung, branch_point = 1, cores = 1)
> BEAM_res <- BEAM_res[order(BEAM_res$qval),]
> BEAM_res <- BEAM_res[,c("gene_short_name", "pval", "qval")]
你可以用一种特殊的热图，将所有“分支依赖”的基因可视化。该热图显示的是同一时间点两个谱系的变化。热图的列是伪时间的点，行是基因。从热图中间往右读，是伪时间的一个谱系；往左是另一个谱系。基因是被按照等级聚类的，所以你看到的基因表达模式和谱系的表达模式是非常相似的。
> plot_genes_branched_heatmap(lung[row.names(subset(BEAM_res,
qval < 1e-4)),],
branch_point = 1,
num_clusters = 4,
cores = 1,
use_gene_short_name = T,
show_rownames = T)
我们还可以用plot_genes_branched_pseudotime函数单独画几个基因（细胞命运决定的marker）：
> lung_genes <- row.names(subset(fData(lung),
gene_short_name %in% c("Ccnd2", "Sftpb", "Pdpn")))
> plot_genes_branched_pseudotime(lung[lung_genes,],
branch_point = 1,
color_by = "Time",
ncol = 1)
附录
计算单细胞表达值
用Monocle之前，你必须计算每一个基因在每个细胞里的表达。有很多种方法，我们推荐Cufflinks，你也可以用RSEM, eXpress, Sailfish等其他的包来计算。这里我们展示一个简单的流程利用TopHat和Cufflinks组合来计算
首先你必须有reads矩阵，如果你是双端测序，每个细胞应该是两个文件，像这样：CELL_TXX_YYY.RZ.fastq.gz（XX是你收集细胞的时间点，YY是你在准备文库的时候96孔板，Z一般是1或者2，表示双端测序时left mate或者right mate。所以例如CELL_T24_A01.R1.fastq.gz 意思是细胞在24小时收集的在A01孔里的文库left mate文件）
用TopHat比对reads（这个有点过时了）
tophat -o CELL_T24_A01_thout -G GENCODE.gtf bowtie-hg19-idx CELL_T24_
A01.R1.fastq.gz CELL_T24_A01.R2.fastq.gz
tophat -o CELL_T24_A02_thout -G GENCODE.gtf bowtie-hg19-idx CELL_T24_
A02.R1.fastq.gz CELL_T24_A02.R2.fastq.gz
tophat -o CELL_T24_A03_thout -G GENCODE.gtf bowtie-hg19-idx CELL_T24_
A03.R1.fastq.gz CELL_T24_A03.R2.fastq.gz
用Cufflinks计算基因表达
cuffquant -o CELL_T24_A01_cuffquant_out GENCODE.gtf
CELL_T24_A01_thout/accepted_hits.bam
cuffquant -o CELL_T24_A02_cuffquant_out GENCODE.gtf
CELL_T24_A02_thout/accepted_hits.bam
cuffquant -o CELL_T24_A03_cuffquant_out GENCODE.gtf
CELL_T24_A03_thout/accepted_hits.bam
接下来把所有的表达矩阵merge到一个文件里：
> cuffnorm --use-sample-sheet -o sc_expr_out GENCODE.gtf
sample_sheet.txt
这里的--use-sample-sheet前提是你要有一个表，里面是所有表达矩阵文件的名字，比如这样：
sample_name group
CELL_T24_A01_cuffquant_out/abundances.cxb T24_A01
CELL_T24_A02_cuffquant_out/abundances.cxb T24_A02
CELL_T24_A03_cuffquant_out/abundances.cxb T24_A03
然后就可以快乐的用Monocle分析了～
总结：
后面还有一些数学公式解释Monocle里面一些函数的原理，我就没看了，数学也不好，就不误导别人了。总体来说，过了一遍原始官方说明书，感觉就是觉得“非常细”，我觉得如果有时间的同学在用到每一个非常重要的包，最好还是阅读一下官方说明书。像我这样成天白天做实验，只有晚上才有时间学习生信的人来说，有时为了偷懒，就直接copy别人的代码了。但是还是想多了解一下里面的用法，以及“这一步到底是在干嘛？”总之这次还是有些收获的，之后还会深入学习其他包（有时间的情况下）。生信的东西就是积少成多，我这个小白就感触颇深，今年夏天的时候我连“R包”是什么都不知道，但现在至少知道些东西了～（虽然还是非常的菜。。。）好好学习，天天向上～共勉之！
