Monocle2包学习笔记（一）：Getting Started with Monocle - 图1

Monocle包是一个强大的单细胞转录组数据分析工具，它引入了一个对单细胞进行排序的新方法（拟时序分析），通过利用不同细胞转录调控的异步过程，将它们沿着与生物过程（如细胞分化等）相对应的发育轨迹进行排序。Monocle通过使用先进的机器学习算法（Reversed Graph Embedding反向图嵌入）从单细胞基因组数据中学习显式主图来对不同的单细胞进行排序，该方法对于解决那些复杂的生物学过程具有很高的鲁棒性和准确性。

Monocle2主要执行以下三种类型的分析：

  • 对细胞进行聚类，分类和计数。
  • 构建单细胞发育分化轨迹。
  • 基因差异表达分析。

安装monocle2及其依赖包

  1. source("http://bioconductor.org/biocLite.R")
  2. biocLite("monocle")
  3. biocLite(c("DDRTree", "pheatmap"))
  5. library(monocle)

monocle2基本分析流程

构建CellDataSet对象存储表达数据

  1. # 加载样本的metadata注释信息
  2. pd <- new("AnnotatedDataFrame", data = sample_sheet)
  3. # 加载基因的metadata注释信息
  4. fd <- new("AnnotatedDataFrame", data = gene_annotation)
  5. # 使用newCellDataSet函数构建CellDataSet对象
  6. cds <- newCellDataSet(expr_matrix, phenoData = pd, featureData = fd)

通过marker基因对细胞进行分类

  1. cth <- newCellTypeHierarchy()
  2. # 选择marker基因
  3. MYF5_id <- row.names(subset(fData(cds), gene_short_name == "MYF5"))
  4. ANPEP_id <- row.names(subset(fData(cds), gene_short_name == "ANPEP"))
  6. # 根据marker基因添加细胞类型注释
  7. cth <- addCellType(cth, "Myoblast", classify_func =
  8.     function(x) { x[MYF5_id,] >= 1 })
  9. cth <- addCellType(cth, "Fibroblast", classify_func =
  10.     function(x) { x[MYF5_id,] < 1 & x[ANPEP_id,] > 1 } )
  12. # 使用classifyCells函数根据marker基因对细胞进行分类
  13. cds <- classifyCells(cds, cth, 0.1)

对细胞进行聚类分群

  1. # 使用clusterCells函数对细胞进行聚类分群
  2. cds <- clusterCells(cds)

拟时序分析对细胞进行排序

  1. disp_table <- dispersionTable(cds)
  2. # 根据基因的表达选择排序的基因
  3. ordering_genes <- subset(disp_table, mean_expression >= 0.1)
  4. cds <- setOrderingFilter(cds, ordering_genes)
  6. # 使用reduceDimension函数进行数据降维
  7. cds <- reduceDimension(cds)
  9. # 使用orderCells函数对细胞进行排序
  10. cds <- orderCells(cds)

执行基因表达差异分析

  1. # 使用differentialGeneTest函数进行差异分析
  2. diff_test_res <- differentialGeneTest(cds, fullModelFormulaStr = "~Media")
  3. sig_genes <- subset(diff_test_res, qval < 0.1)

数据导入与CellDataSet对象构建

Monocle2使用CellDataSet对象存储单细胞RNA-seq的基因表达矩阵及其相应的metadata注释信息，该对象衍生自Bioconductor中的ExpressionSet类。

CellDataSet类主要需要以下三个输入文件：

  • exprs，表示数值的基因表达矩阵，其中行是基因，列是细胞。
  • phenoData，一个AnnotatedDataFrame对象，其中行是细胞，列是细胞的属性（例如细胞的类型，培养条件，捕获的天数等）。
  • featureData，一个AnnotatedDataFrame对象，其中行是特征（例如基因），列是基因的属性，例如生物型，gc含量等。

其中，这三个输入文件必须满足以下要求：

  • 基因表达矩阵的列数与phenoData的行数相同。
  • 基因表达矩阵的行数与featureData的行数相同。
  • phenoData对象的行名称应与表达矩阵的列名称相匹配。
  • featureData对象的行名应与表达矩阵的行名相匹配。
  • featureData中必须有一个名为“gene_short_name”的列。
  1. library(monocle)
  2. # 读取所需的三个输入文件
  3. # 基因表达矩阵
  4. HSMM_expr_matrix <- read.table("fpkm_matrix.txt")
  5. # 细胞注释信息
  6. HSMM_sample_sheet <- read.delim("cell_sample_sheet.txt")
  7. # 基因注释信息
  8. HSMM_gene_annotation <- read.delim("gene_annotations.txt")
  10. # 构建CellDataSet对象
  11. pd <- new("AnnotatedDataFrame", data = HSMM_sample_sheet)
  12. fd <- new("AnnotatedDataFrame", data = HSMM_gene_annotation)
  14. HSMM <- newCellDataSet(as.matrix(HSMM_expr_matrix),
  15.                        phenoData = pd, featureData = fd)

默认情况下，Monocle2假定我们读取的表达数据是以转录计数的原始count值，并使用负二项分布模型来进行下游的基因表达差异分析。

注意：如果我们有原始的UMI counts计数矩阵，则在创建CellDataSet对象之前不应该对其进行规归一化处理，也不应尝试将UMI计数转换为相对丰度（如将其转换为FPKM / TPM数据）。Monocle2将在内部执行所有必需的标准化步骤，自己提前进行归一化处理可能会破坏Monocle的一些关键步骤。

Monocle2也可以通过importCDS和exportCDS函数将其他格式的数据导入为CellDataSet对象，或将CellDataSet对象导出为其他格式。

  1. # Where 'data_to_be_imported' can either be a Seurat object or an SCESet.
  2. # 将其他格式的数据导入为CellDataSet对象
  3. importCDS(data_to_be_imported)
  5. # We can set the parameter 'import_all' to TRUE if we'd like to import all the slots from our Seurat object or SCESet. (Default is FALSE or only keep minimal dataset)
  6. importCDS(data_to_be_imported, import_all = TRUE)
  8. lung <- load_lung()
  9. # To convert to Seurat object
  10. # 将CellDataSet对象导出为Seurat对象
  11. lung_seurat <- exportCDS(lung, 'Seurat')
  13. # To convert to SCESet
  14. # 将CellDataSet对象导出为SCESet对象
  15. lung_SCESet <- exportCDS(lung, 'Scater')

Choosing a distribution for your data

Monocle2既可以使用基因的相对表达丰度，也可以使用基于计数的原始count值作为输入的表达矩阵。通常，它最适合使用转录本的原始count值，尤其是UMI数据。无论我们使用哪种类型的数据，都必须为其指定适当的分布。FPKM / TPM值通常呈对数正态分布，而UMI或read counts可更好地使用负二项分布建模。对于原始的count计数数据，请将负二项式分布指定为newCellDataSet的expressionFamily参数。

  1. #Do not run
  2. HSMM <- newCellDataSet(count_matrix,
  3.                 phenoData = pd,
  4.                 featureData = fd,
  5.                 expressionFamily=negbinomial.size())

以下是一些常用的expressionFamily值及其适用条件：
Monocle2包学习笔记（一）：Getting Started with Monocle - 图2

Notes：Using the wrong expressionFamily for your data will lead to bad results, errors from Monocle, or both. However, if you have FPKM/TPM data, you can still use negative binomial if you first convert your relative expression values to transcript counts using relative2abs(). This often leads to much more accurate results than using tobit(). See the section on Converting TPM to mRNA Counts for details.

Working with large data sets

对于一些大型的单细胞RNA-seq数据集，如具有50,000个细胞和25,000+个基因组成的基因表达矩阵会占用大量的内存。但是，由于当前的测序技术通常无法捕获每个细胞中的全部或大多数mRNA分子，因此许多基因的表达均会被检测为零。我们通常建议大多数用户使用稀疏矩阵sparseMatrices的表达数据，它会加快许多计算的速度。
我们可以使用Matrix包将基因表达矩阵转换为稀疏矩阵提供给Monocle：

  1. HSMM <- newCellDataSet(as(umi_matrix, "sparseMatrix"),
  2.                 phenoData = pd,
  3.                 featureData = fd,
  4.                 lowerDetectionLimit = 0.5,
  5.                 expressionFamily = negbinomial.size())

如果我们有10X Genomics产生的counts数据，可以使用cellrangerRkit包中的load_cellranger_matrix函数直接加载基因表达矩阵。

  1. cellranger_pipestance_path <- "/path/to/your/pipeline/output/directory"
  2. gbm <- load_cellranger_matrix(cellranger_pipestance_path)
  3. # 获取基因metadata注释信息
  4. fd <- fData(gbm)
  6. # The number 2 is picked arbitrarily in the line below.
  7. # Where "2" is placed you should place the column number that corresponds to your featureData's gene short names.
  8. # 将基因注释信息的第二列命名为gene_short_name
  9. colnames(fd)[2] <- "gene_short_name"
  11. # 构建CellDataSet对象
  12. gbm_cds <- newCellDataSet(exprs(gbm),
  13.                   phenoData = new("AnnotatedDataFrame", data = pData(gbm)),
  14.                   featureData = new("AnnotatedDataFrame", data = fd),
  15.                   lowerDetectionLimit = 0.5,
  16.                   expressionFamily = negbinomial.size())

Converting TPM/FPKM values into mRNA counts

如果我们得到的基因表达矩阵是已经进行了归一化处理的数据（如FPKM或TPM值），我们可以使用relative2abs函数在创建CellDataSet对象之前将其转换为RPC值，如下所示：

  1. pd <- new("AnnotatedDataFrame", data = HSMM_sample_sheet)
  2. fd <- new("AnnotatedDataFrame", data = HSMM_gene_annotation)
  4. # First create a CellDataSet from the relative expression levels
  5. HSMM <- newCellDataSet(as.matrix(HSMM_expr_matrix),
  6.                 phenoData = pd,
  7.                 featureData = fd,
  8.                 lowerDetectionLimit = 0.1,
  9.                 expressionFamily = tobit(Lower = 0.1))
  11. # Next, use it to estimate RNA counts
  12. rpc_matrix <- relative2abs(HSMM, method = "num_genes")
  14. # Now, make a new CellDataSet using the RNA counts
  15. HSMM <- newCellDataSet(as(as.matrix(rpc_matrix), "sparseMatrix"),
  16.                 phenoData = pd,
  17.                 featureData = fd,
  18.                 lowerDetectionLimit = 0.5,
  19.                 expressionFamily = negbinomial.size())

变异系数计算和数据过滤预处理

Estimate size factors and dispersions

Monocle2使用estimateSizeFactors和estimateDispersions两个函数计算Size factors和”dispersion” values，Size factors可以帮助我们对不同细胞中mRNA表达的差异进行归一化处理，”dispersion” values将有助于后续的基因差异表达分析。

  1. HSMM <- estimateSizeFactors(HSMM)
  2. HSMM <- estimateDispersions(HSMM)

Notes: estimateSizeFactors() and estimateDispersions() will only work, and are only needed, if you are working with a CellDataSet with a negbinomial() or negbinomial.size() expression family.

Filtering low-quality cells

Monocle2使用detectGenes函数查看基因的表达情况，并根据基因的表达预测其是否用于后续的细胞排序分析。

  1. HSMM <- detectGenes(HSMM, min_expr = 0.1)
  2. print(head(fData(HSMM)))
  3.     gene_short_name    biotype    num_cells_expressed    use_for_ordering
  4. ENSG00000000003.10    TSPAN6    protein_coding    184    FALSE
  5. ENSG00000000005.5    TNMD    protein_coding    0    FALSE
  6. ENSG00000000419.8    DPM1    protein_coding    211    FALSE
  7. ENSG00000000457.8    SCYL3    protein_coding    18    FALSE
  8. ENSG00000000460.12    C1orf112    protein_coding    47    TRUE
  9. ENSG00000000938.8    FGR    protein_coding    0    FALSE
  11. print(head(pData(HSMM)))
  12. Library    Well    Hours    Media    Mapped.Fragments    Pseudotime    State    Size_Factor    num_genes_expressed
  13. T0_CT_A01    SCC10013_A01    A01    0    GM    1958074    23.916673    1    1.392811    6850
  14. T0_CT_A03    SCC10013_A03    A03    0    GM    1930722    9.022265    1    1.311607    6947
  15. T0_CT_A05    SCC10013_A05    A05    0    GM    1452623    7.546608    1    1.218922    7019
  16. T0_CT_A06    SCC10013_A06    A06    0    GM    2566325    21.463948    1    1.013981    5560
  17. T0_CT_A07    SCC10013_A07    A07    0    GM    2383438    11.299806    1    1.085580    5998
  18. T0_CT_A08    SCC10013_A08    A08    0    GM    1472238    67.436042    2    1.099878    6055
  20. # 选取特定的基因
  21. expressed_genes <- row.names(subset(fData(HSMM),
  22.     num_cells_expressed >= 50))

现在，expressed_genes向量中存储着至少在50个细胞中表达的基因，我们将用其对后续的细胞按照生物学过程的顺序进行排序的分析。

  1. # 对细胞进行过滤
  2. valid_cells <- row.names(subset(pData(HSMM),
  3.             Cells.in.Well == 1 &
  4.             Control == FALSE &
  5.             Clump == FALSE &
  6.             Debris == FALSE &
  7.             Mapped.Fragments > 1000000))
  8. HSMM <- HSMM[,valid_cells]
  10. # 根据基因的表达对基因进行过滤
  11. pData(HSMM)$Total_mRNAs <- Matrix::colSums(exprs(HSMM))
  12. HSMM <- HSMM[,pData(HSMM)$Total_mRNAs < 1e6]
  14. upper_bound <- 10^(mean(log10(pData(HSMM)$Total_mRNAs)) +
  15.             2*sd(log10(pData(HSMM)$Total_mRNAs)))
  16. lower_bound <- 10^(mean(log10(pData(HSMM)$Total_mRNAs)) -
  17.             2*sd(log10(pData(HSMM)$Total_mRNAs)))
  18. # 查看基因的表达分布
  19. qplot(Total_mRNAs, data = pData(HSMM), color = Hours, geom =
  20. "density") +
  21. geom_vline(xintercept = lower_bound) +
  22. geom_vline(xintercept = upper_bound)

Monocle2包学习笔记（一）：Getting Started with Monocle - 图3

  1. # 过滤掉那些极高表达和极低表达的细胞
  2. HSMM <- HSMM[,pData(HSMM)$Total_mRNAs > lower_bound &
  3.       pData(HSMM)$Total_mRNAs < upper_bound]
  4. HSMM <- detectGenes(HSMM, min_expr = 0.1)
  1. # Log-transform each value in the expression matrix.
  2. L <- log(exprs(HSMM[expressed_genes,]))
  4. # Standardize each gene, so that they are all on the same scale, Then melt the data with plyr so we can plot it easily
  5. melted_dens_df <- melt(Matrix::t(scale(Matrix::t(L))))
  7. # Plot the distribution of the standardized gene expression values.
  8. qplot(value, geom = "density", data = melted_dens_df) +
  9. stat_function(fun = dnorm, size = 0.5, color = 'red') +
  10. xlab("Standardized log(FPKM)") +
  11. ylab("Density")

Monocle2包学习笔记（一）：Getting Started with Monocle - 图4

参考来源：http://cole-trapnell-lab.github.io/monocle-release/docs/

Monocle2包学习笔记（一）：Getting Started with Monocle - 图5
Monocle2包学习笔记（一）：Getting Started with Monocle - 图6