拟时间分析R包Monocle2使用说明

在单细胞RNA测序项目中，细胞谱系分析是常规的分析内容之一。所谓细胞谱系分析，主要是通过时间轴的细胞轨迹构建、基因表达变化情况来揭示细胞从一种类型转变为另一种类型、或者一种状态转变为另一种状态。其中最常用的R包为Monocle2。

一、分析对象

分析对象主要是基于CellDataSet的Class，此Class需要三个输入文件构成，分别为:

1. exprs, a numeric matrix of expression values, where rows are genes, and columns are cells;即表达矩阵。对于表达矩阵，如果数据量较大，建议采用sparseMatrices 将之转化稀疏矩阵。但是，其实上如CellRanger等大部分的上游工具输出的已经是稀疏矩阵之格式，则可以直接输入。

2. phenoData, an AnnotatedDataFrame object, where rows are cells, and columns are cell attributes (such as cell type, culture condition, day captured, etc.)  即表型矩阵，每一行为一个细胞，每一列为细胞属性信息，如细胞类型、培养条件和捕获天数等。

3. featureData, an AnnotatedDataFrame object, where rows are features (e.g. genes), and columns are gene attributes, such as biotype, gc content, etc. 。特征矩阵，主要是基因的一些情况，如生物类型，gc含量等。

我们首先分别读入这三个文件，然后采用newCellDataSet来创造一个CellDataSet对象。命令如下：

默认情况下，Monocle2认为输入的为表达量的绝对值，符合负二项分布，采用的family function = negbinomial.size()；如果输入的表达量为TPM或者FPKM，则可以采用函数relative2abs将相对表达量转化为绝对表达量，即mRNAs per cell (RPC)，这可能比直接使用tobit()分布函数效果更好。不同类型的输入变量有不同的分布函数，需要指定，具体如下：

二、轨迹构建

在发育过程中，应对外界刺激和在生命周期中，细胞从一种功能状态转化为另一种。在这种转化中，不同状态的细胞表达不同的基因集，这种表达的动态变化是受到转录调控的，在这个过程中，一些基因表达沉默，同时其他基因的表达重新激活。Monocle2采用机器学习算法reversed graph embedding来学习基因表达变化的顺序，从而构建轨迹。在构建轨迹之前，我们需要先理解何为“Pseudotime”：其实在许多生物过程中，细胞间的进程都不是完美同步的，有的细胞可能在某个过程的前端，而有的细胞还没开始，此称之为“异步”；这种异步会让细胞从一个状态转变为另一个状态的调控变化的理解变得困难。通过将细胞根据这个过程进行排序而形成一个轨迹，从而追踪伴随轨迹的过程变化函数，称之为“pseudotime”。

轨迹构建主要有三个步骤：

1）Choosing genes for ordering：选择基因，即是机器学习中选择特征的过程。选择的基因是影响轨迹形状构建的主要因素。单细胞RNA测序中，虽然低表达的基因往往为噪音，但是也可能包含重要的信息。基因的选择有半监督和非监督两种方式。

2) Reducing the dimensionality of the data：数据降维。使用的函数为reduceDimension。可以选择对函数reduceDimension选择两种算法：

（1）Independent Component Analysis，为Monocle最初使用的线性降维算法；

（2）DDRTree，为非线性降维技术，是目前Monocle2默认采用的降维算法，建议采用。

3) Ordering the cells in pseudotime

轨迹构建。Monocle假设轨迹是一个树状结构，一端为“根（root）”,另一端为“枝叶（leaaves）”Monocle 采用命令orderCells函数来对细胞进行排序，该函数有个重要的参数可以用来设置root of the tree，作为轨迹的起点。一旦轨迹构建完成，则采用函数plot_cell_trajectory来可视化轨迹图。

Monocle2构建的轨迹图

在轨迹构建中，最关键的是特征基因的选择，不同的选择将会带来不同轨迹结果。其中，挑选基因有几种方法：

1）无监督。即主要是依靠数据本身而不是先验知识来确定特征基因。具体而言：（1）如果有时间节点信息，可以基于时间的差异分析；（2）Selecting genes with high dispersion across cells：选择在细胞间基因表达差异大的基因（3）Selecting genes based on PCA loading：通过irlba R包进行PCA分析，选择a good set of ordering genes.（4）Unsupervised feature selection based on density peak clustering：基本思想是先对细胞进行聚类，然后基于不同的group进行差异分析，发现的差异表达基因作为marker基因进行轨迹的构建。

2）半监督：You first define genes that mark progress using the CellTypeHierchy system, very similar to how we used it for cell type classification. Then, you use it to select ordering genes that co-vary with these markers. Finally, you order the cell based on these genes just as we do in unsupervised ordering. So the only difference between unsupervised and semi-supervised ordering is in which genes we use for ordering. ---先确定已知的marker基因（可能是基于先验的生物知识得来），再寻找与这些marker基因共表达变化的基因，最后基于这些基因进行细胞的轨迹构建。

三、拟时间函数的基因表达变化表示

通常需要展现一些有意义的基因的表达情况，过程如下：

1）首先确定要展现的基因集：

to_be_tested <- row.names(subset(fData,

gene_short_name %in% c("MYH3", "MEF2C", "CCNB2", "TNNT1")))

2）然后从之前构建的class 中选择这些基因的信息

cds_subset <- HSMM_myo[to_be_tested,]

3）最后使用命令plot_genes_in_pseudotime(cds_subset, color_by="Hours")来展现基因的表达情况：

参考文献：Monocle: Cell counting, differential expression, and trajectory analysis for single-cell RNA-Seq experiments

多生信分析需求，请加微信13621202201