实验记录9：Monocle包的使用方法
梗概
参照官网教程：http://cole-trapnell-lab.github.io/monocle-release/docs/#getting-started-with-monocle
或在R中运行以下命令：
browseVignettes("monocle")
关于Monocle
Monocle introduced the strategy of ordering single cells in pseudotime,
placing them along a trajectory corresponding to a biological process such as cell differentiation. Monocle learns this trajectory directly from the data, in either a fully unsupervised or a semi-supervised manner. It also performs differential gene expression and clustering to identify important genes and cell states.
Monocle是用于分析single-cell RNA-seq的R包，能够将细胞按照模拟的时间顺序进行排列，显示它们的发展轨迹如细胞分化等生物学过程。Monocle从数据中用无监督或半监督（？）的方式直接获得这个轨迹。除此之外，它还能够做差异表达分析和聚类来揭示重要的基因和细胞。
在发育过程中，为了响应刺激和生命，细胞从一种功能性“状态”转变为另一种功能性“状态”。处于不同状态的细胞表达不同的基因组，产生蛋白质和代谢物的动态重复序列，从而完成它们的工作。当细胞在状态之间移动时，经历转录重组的过程，一些基因被沉默，另一些被新激活。这些瞬态通常难以表征，因为在更稳定的端点状态之间纯化细胞可能是困难的或不可能的。单细胞RNA-Seq可以让您在不需要纯化的情况下查看这些状态。但是，为此，我们必须确定每个单元格的可能状态范围。
Monocle介绍了使用RNA-Seq进行单细胞轨迹分析的策略。Monocle不是通过实验将细胞纯化为离散状态，而是使用算法来学习每个细胞必须经历的基因表达变化的序列，作为动态生物过程的一部分。一旦它了解了基因表达变化的整体“轨迹”，Monocle就可以将每个细胞放置在轨迹中的适当位置。然后，您可以使用Monocle的差异分析工具包来查找在轨迹过程中受到调控的基因，如“ 查找作为假性函数的变化的基因 ”一节中所述。如果该过程有多个结果，Monocle将重建“分支”轨迹。这些分支对应于细胞“决策”，Monocle提供了强大的工具来识别受其影响的基因并参与制作它们。您可以在“分析单细胞轨迹中的分支”一节中查看如何分析分支。Monocle依赖于称为反向图嵌入的机器学习技术来构建单细胞轨迹。您可以 在Monocle理论背后的部分中阅读更多关于Monocle方法的理论基础，或参考小插图末尾显示的参考文献
什么是伪时间？
Pseudotime是衡量单个细胞通过细胞分化等过程取得多少进展的量度。在许多生物过程中，细胞不能完美同步。在诸如细胞分化的过程的单细胞表达研究中，捕获的细胞可以在进展方面广泛分布。也就是说，在完全同时捕获的细胞群中，一些细胞可能很远，而其他细胞甚至可能尚未开始该过程。当您想要了解细胞从一个状态转换到另一个状态时发生的调节变化的顺序时，这种不同步会产生重大问题。跟踪同时捕获的细胞的表达产生非常压缩的基因动力学感觉，以及该基因的明显变异性' 的表达会非常高。通过沿着学习轨迹根据其进度对每个单元进行排序，Monocle减轻了由于不同步引起的问题。Monocle不是根据时间跟踪表达式的变化，而是根据轨迹上的进度跟踪变化，我们称之为pseudotime'。Pseudotime是一个抽象的进展单元：它只是一个单元格与轨迹起点之间的距离，沿着最短路径测量。轨迹的总长度是根据细胞从起始状态移动到最终状态时经历的转录变化的总量来定义的。
实验记录
伪时间分析实验技术流程主要包括三个步骤，而每一步都是一个机器学习任务：
step1：选择输入基因用于机器学习
这个过程称为feature selection（特征选择），这些基因对轨迹的形状有着最重要的影响。我们应该要选择的是最能反映细胞状态的基因。
如果直接通过Seurat输出的一些重要的基因（比如每个cluster中的高FC值基因）作为输入对象的话就能够实现一个“无监督”分析。或者也可以利用生物学知识手动选择一些重要的基因进行“半监督”分析。
step2：数据降维
利用Reversed graph embedding算法将数据降维。没有太懂。
相对于PCA来说这个算法更能够反应数据的内部结构（据说是这样）
step3：将细胞按照伪时间排序
这个过程称为manifold learning（流形学习）。Monocle利用轨迹来描述细胞如何从一个状态转换到另一个状态。得到的是一个树状图，树的根部或树干表示的是细胞最初的状态（如果有的话），随着细胞的变化就沿着分叉树枝发展。一个细胞的伪时间值（pseudotime value）为它的位置沿着树枝到根部的距离。
以下是正式实验记录：
step0：数据输入
下载与安装
实验要求——
R version 3.4 or higher, Bioconductor version 3.5, and monocle 2.4.0 or higher to have access to the latest features.
查看R版本：
sessionInfo()
R version 3.4.4 (2018-03-15)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: OS X El Capitan 10.11.6
安装bioconductor以及monocle
source("http://bioconductor.org/biocLite.R")
biocLite("monocle")
加载程辑包：
library(monocle)
将Seurat Object导入
Monocle包可以利用importCDS()命令将Seurat包中的数据导入，同时也能用exportCDS()将数据导出到其他包中。
由于之前的Seurat变量已经被覆盖，故重新创建一个。
spleen_monocle <- CreateSeuratObject(raw.data = spleen.data, min.cells = 3, min.genes = 200, project = "10X_spleen")
spleen_monocle <- importCDS(spleen_monocle, import_all = TRUE)
查看数据：
spleen_monocle
CellDataSet (storageMode: environment)
assayData: 15655 features, 1959 samples
element names: exprs
protocolData: none
phenoData
sampleNames: AAACCTGAGGTGATAT AAACCTGAGTCATCCA ... TTTGTCAGTCGTCTTC (1959 total)
varLabels: nGene nUMI orig.ident Size_Factor
varMetadata: labelDescription
featureData
featureNames: RP11-34P13.7 FO538757.2 ... AC240274.1 (15655 total)
fvarLabels: gene_short_name
fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
Annotation:
head(fData(spleen_monocle))
gene_short_name
RP11-34P13.7 RP11-34P13.7
FO538757.2 FO538757.2
AP006222.2 AP006222.2
RP4-669L17.10 RP4-669L17.10
RP11-206L10.9 RP11-206L10.9
LINC00115 LINC00115
head(pData(spleen_monocle))
nGene nUMI orig.ident Size_Factor
AAACCTGAGGTGATAT 1072 3999 10X_spleen NA
AAACCTGAGTCATCCA 1562 4640 10X_spleen NA
AAACCTGCAGTGAGTG 995 3489 10X_spleen NA
AAACGGGAGACTGGGT 1704 7397 10X_spleen NA
AAACGGGAGACTTGAA 976 3102 10X_spleen NA
AAACGGGAGATAGGAG 1236 4548 10X_spleen NA
15655个基因，1959个细胞，与之前创建的的Seurat对象相符。
（细胞名是一些barcode序列）
step1：选择输入基因用于机器学习
①脾脏数据subset的获取
因为脾脏是在4个时间点经过缺血处理的，
spleen_monocle <- detectGenes(spleen_monocle)
查看数据，添加了一列num_cells_expressed
head(fData(spleen_monocle))
gene_short_name num_cells_expressed
RP11-34P13.7 RP11-34P13.7 3
FO538757.2 FO538757.2 214
AP006222.2 AP006222.2 157
RP4-669L17.10 RP4-669L17.10 10
RP11-206L10.9 RP11-206L10.9 25
LINC00115 LINC00115 35
print(pData(spleen_monocle))
nGene nUMI orig.ident percent.mito res.0.6 Size_Factor num_genes_expressed
AAACCTGAGGTGATAT 1072 3999 10X_spleen 0.0115086315 3 NA 1070
AAACCTGAGTCATCCA 1562 4640 10X_spleen 0.0370929480 5 NA 1560
AAACCTGCAGTGAGTG 995 3489 10X_spleen 0.0395528805 3 NA 995
AAACGGGAGACTGGGT 1704 7397 10X_spleen 0.0285250777 0 NA 1704
AAACGGGAGACTTGAA 976 3102 10X_spleen 0.0493230174 2 NA 976
AAACGGGAGATAGGAG 1236 4548 10X_spleen 0.0268249780 1 NA 1236
AAACGGGGTCCGAACC 2430 6863 10X_spleen 0.0379119277 6 NA 2425
AAACGGGTCAGGTAAA 1256 4750 10X_spleen 0.0381052632 1 NA 1256
AAACGGGTCGTCACGG 1887 8095 10X_spleen 0.0302730755 5 NA 1885
AAAGATGAGCGATGAC 1473 5729 10X_spleen 0.0185088179 4 NA 1471
…………
step2：数据降维
step3：将细胞按照伪时间排序
(此段内容舍弃)将Seurat Object导入
Monocle包可以利用importCDS()命令将Seurat包中的数据导入，同时也能用exportCDS()将数据导出到其他包中。
spleen_monocle <- importCDS(spleen, import_all = TRUE)
查看数据：
spleen_monocle
CellDataSet (storageMode: environment)
assayData: 15655 features, 1940 samples
element names: exprs
protocolData: none
phenoData
sampleNames: AAACCTGAGGTGATAT AAACCTGAGTCATCCA ... TTTGTCAGTCGTCTTC (1940 total)
varLabels: nGene nUMI ... Size_Factor (6 total)
varMetadata: labelDescription
featureData
featureNames: RP11-34P13.7 FO538757.2 ... AC240274.1 (15655 total)
fvarLabels: gene_short_name
fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
Annotation:
15655个基因，1940个细胞，与之前的Seurat对象相符。
选择数据分布类型
需要根据实验的数据类型选择数据的分布，比如是相对表达量的数据/基于计数的数据（如UMI），官方认为：“FPKM/TPM values are generally log-normally distributed, while UMIs or read counts are better modeled with the negative binomial.”即利用负二项式的模型会更加适合UMI计数的实验数据。
参数描述.png
spleen_M <- newCellDataSet(spleen_monocle, expressionFamily=negbinomial.size())
报错：spleen_monocle不是矩阵.
Error in newCellDataSet(spleen_monocle, expressionFamily = negbinomial.size()) :
Error: argument cellData must be a matrix (either sparse from the Matrix package or dense)
In addition: Warning message:
In newCellDataSet(spleen_monocle, expressionFamily = negbinomial.size()) :
Warning: featureData must contain a column verbatim named 'gene_short_name' for certain functions
预估size factors和分散
spleen_monocle <- estimateSizeFactors(spleen_monocle)
spleen_monocle <- estimateDispersions(spleen_monocle)
removing 9 outliers.
比较奇怪，在输入这两个命令之后显示去除了9个异常值，但是查看spleen_monocle变量时返回：
CellDataSet (storageMode: environment)
assayData: 15655 features, 1940 samples
element names: exprs
protocolData: none
phenoData
sampleNames: AAACCTGAGGTGATAT AAACCTGAGTCATCCA ... TTTGTCAGTCGTCTTC (1940 total)
varLabels: nGene nUMI ... Size_Factor (6 total)
varMetadata: labelDescription
featureData
featureNames: RP11-34P13.7 FO538757.2 ... AC240274.1 (15655 total)
fvarLabels: gene_short_name
fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
Annotation:
> dim(spleen_monocle)
Features Samples
15655 1940
仍然是15655个基因，1940个细胞，与之前无异。因此采用了另一种输入方式。
remove(spleen_monocle)
