这里是佳奥！我们继续学习最重要的三个R包。

1 scater Package(calculateQCMetrics已停用，用addPerCellQC替代)

加载R包

rm(list = ls()) # clear the environment
#load all the necessary libraries
options(warn=-1) # turn off warning message globally
suppressMessages(library(scater))
suppressMessages(library(scRNAseq))

创建测试数据集

library(scRNAseq)
## ----- Load Example Data -----
fluidigm<-ReprocessedFluidigmData()
data(fluidigm)

# Set assay to RSEM estimated counts
assay(fluidigm) <- assays(fluidigm)$rsem_counts
ct <- floor(assays(fluidigm)$rsem_counts)
ct[1:4,1:4]
table(rowSums(ct)==0)
# 这里使用原始表达矩阵，所以有很多基因在所有细胞均无表达量，即表现为没有被检测到，这样的基因是需要过滤掉的。

pheno_data <- as.data.frame(colData(fluidigm))
## 这里需要把Pollen的表达矩阵做成我们的 scater 要求的对象
#data("sc_example_counts")
#data("sc_example_cell_info") 
# 也可以尝试该R包自带的数据集。
# https://bioconductor.org/packages/release/bioc/vignettes/scater/inst/doc/vignette-intro.R
sce <- SingleCellExperiment(
    assays = list(counts = ct), 
    colData = pheno_data
    )
sce
# 后面所有的分析都是基于 sce 这个变量
# 是一个 SingleCellExperiment 对象，被很多单细胞R包采用。
> sce
class: SingleCellExperiment 
dim: 26255 130 
metadata(0):
assays(1): counts
rownames(26255): A1BG A1BG-AS1 ... ZZEF1 ZZZ3
rowData names(0):
colnames(130): SRR1275356 SRR1274090 ... SRR1275366 SRR1275261
colData names(28): NREADS NALIGNED ... Cluster1 Cluster2
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):

探究SingleCellExperiment对象——sce

##sce是S4对象，与getGEO返回的对象相似

> str(sce,max.level=2)

##罗列了所有属性以及下一级菜单的属性
> str(sce)
Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 9 slots
  ..@ int_elementMetadata:Formal class 'DFrame' [package "S4Vectors"] with 6 slots
  .. .. ..@ rownames       : NULL
  .. .. ..@ nrows          : int 26255
  .. .. ..@ elementType    : chr "ANY"
  .. .. ..@ elementMetadata: NULL
  .. .. ..@ metadata       : list()
  .. .. ..@ listData       :List of 1
  .. .. .. ..$ rowPairs:Formal class 'DFrame' [package "S4Vectors"] with 6 slots
  .. .. .. .. .. ..@ rownames       : NULL
  .. .. .. .. .. ..@ nrows          : int 26255
  .. .. .. .. .. ..@ elementType    : chr "ANY"
  .. .. .. .. .. ..@ elementMetadata: NULL
  .. .. .. .. .. ..@ metadata       : list()
  .. .. .. .. .. ..@ listData       : Named list()
  ..@ int_colData        :Formal class 'DFrame' [package "S4Vectors"] with 6 slots
 略

一些质量控制

这里稍作讲解，包括基因层面及细胞层面，但是并不设计实验本身引入的干扰因素。

首先看基因层面的过滤

使用 calculateQCMetrics函数作用于 sce 那个单细胞数据对象后，就可以 用 rowData(object) 可以查看各个基因统计指标：

• mean_counts: 平均表达量
• log10_mean_counts: 归一化 log10-scale.
• pct_dropout_by_counts: 该基因丢失率。
• n_cells_by_counts: 多少个细胞表达了该基因

上面那些指标可以用来过滤，也可以自己重新再次计算一下那些统计学指标。

主要是过滤低表达量基因，还有 线粒体基因 和 ERCC spike-ins 的控制。

exprs(sce) <- log2( calculateCPM( sce ) + 1)
## 只有运行了下面的函数后才有各式各样的过滤指标
genes <- rownames(rowData(sce))
genes[grepl('^MT-',genes)]
genes[grepl('^ERCC-',genes)]
# 比较不幸的是，这个测试数据里面既没有线粒体基因，有没有ERCC序列。
# > genes[grepl('^MT-',genes)]
# character(0)
# > genes[grepl('^ERCC-',genes)]
# character(0)

sce <- calculateQCMetrics(sce, 
        feature_controls = list(ERCC = grep('^ERCC',genes)))

sce <- addPerCellQC(sce, 
        subsets = list(ERCC = grep('^ERCC',genes)))

# 也没有定义啥细胞是属于control组，所以这里并不需要完全follow教程
# example_sce <- calculateQCMetrics(example_sce, 
# feature_controls = list(ERCC = 1:20, mito = 500:1000),
# cell_controls = list(empty = 1:5, damaged = 31:40))

查看信息

tmp <- as.data.frame(rowData(sce))
colnames(tmp)
head(tmp)

目前只过滤那些在所有细胞都没有表达的基因, 但是这个过滤条件可以自行调整。可以看到基因数量大幅度减少。

keep_feature <- rowSums(exprs(sce) > 0) > 0
table(keep_feature)
sce <- sce[keep_feature,]
sce

然后看细胞层面的过滤

用 colData(object) 可以查看各个样本统计情况

• total_counts: total number of counts for the cell (aka ‘library size’)

• log10_total_counts: total_counts on the log10-scale

• total_features: the number of features for the cell that have expression above the detection limit (default detection limit is zero)

每个数据集都有适合自己的过滤阈值，不要局限于教程

tmp <- as.data.frame(colData(sce))
colnames(tmp) 
# 可以发现细胞质控属性非常多，可以说是实用至极。
## 比如看每个样本测到的基因数量
tf <- sce$total_features_by_counts 
boxplot(tf)
fivenum(tf)

还有很多其它指标，比如哪些细胞的ERCC含量过高，或者某些其高表达量基因占比太高，等等。

一些可视化

可视化函数非常多，所以这个R包才会被引用那么广泛，里面包括如下：

首先是一些细胞距离情况

• plotPCA: produce a principal components plot for the cells.
• plotTSNE: produce a t-distributed stochastic neighbour embedding (reduced dimension) plot for the cells.
• plotDiffusionMap: produce a diffusion map (reduced dimension) plot for the cells.
• plotMDS: produce a multi-dimensional scaling plot for the cells.
• plotReducedDim: plot a reduced-dimension representation of the cells.

然后是一些表达量相关的:

• plotExpression: plot expression levels for a defined set of features.
• plotPlatePosition: plot cells in their position on a plate, coloured by cell metadata and QC metrics or feature expression level.
• plotColData: plot cell metadata and QC metrics.
• plotRowData: plot feature metadata and QC metrics.

首先可视化表达量

可以添加一些细胞属性，展示如下：

# 挑选一些细胞属性（临床信息） 来进行可视化展示。
colnames(tmp)[25:28]
# 主要展现下面这些基因
rownames(sce)[1:6]

展示一些基因在不同细胞分类的表达区别

plotExpression(sce, rownames(sce)[1:6],
               x = "Biological_Condition", 
               exprs_values = "logcounts") 

散点图，展示两个基因的相关性，这里批量展示6对相关性

plotExpression(sce, rownames(sce)[1:6],
               x = rownames(sce)[1])

还可以更复杂

plotExpression(sce, rownames(sce)[1:6],
               colour_by = "Biological_Condition", 
               shape_by = "Cluster1", 
               size_by = rownames(sce)[1])

还可以可视化细胞距离分布

这里可以学习PCA或者tSNE的使用

sce <- runPCA(sce)
# 这里并没有进行任何基因的挑选，就直接进行了PCA，与 Seurat包不一样。
reducedDimNames(sce)

PCA分布图上面添加临床信息，同样的可以看到GW16混在了GW21里面

plotReducedDim(sce, use_dimred = "PCA", 
               colour_by = "Biological_Condition" )

PCA分布图上面添加表达量信息, 但这个基本上不怎么用

plotReducedDim(sce, use_dimred = "PCA", 
               colour_by = rownames(sce)[1], 
               size_by = rownames(sce)[2])

最原始的

plotPCA(sce)

其它降维算法结果的可视化

## ----plot-pca-feature-controls-----
# 这里在真实场景中应用比较多。
sce2 <- runPCA(sce, 
               feature_set = rowData(sce)$is_feature_control)
plotPCA(sce2)

仅仅是选取前20个PC, 可以看到绘图时候细胞分布并没有太大区别

## ----plot-pca-4comp-colby-shapeby-----
sce <- runPCA(sce, ncomponents=20)
plotPCA(sce,  colour_by = "Biological_Condition" )

但是可以挑选指定的PC来可视化, 这里选择的是第4个PC

plotPCA(sce, ncomponents = 4,  
        colour_by = "Biological_Condition" )

意义不大的展示

## ----plot-pca-4comp-colby-sizeby-exprs-----
plotPCA(sce,
        colour_by = rownames(sce)[1], 
        size_by = rownames(sce)[2])

tSNE可视化

## ----plot-tsne-1comp-colby-sizeby-exprs------
# Perplexity of 10 just chosen here arbitrarily. 
set.seed(1000)
# 这里的这个 perplexity 参数很重要
sce <- runTSNE(sce, perplexity=30)
plotTSNE(sce, 
         colour_by = "Biological_Condition" )

## ----plot-tsne-from-pca------
set.seed(1000)
sce <- runTSNE(sce, perplexity=10, use_dimred="PCA", n_dimred=10)
plotTSNE(sce,  
         colour_by = "Biological_Condition" )

这一步会调用destiny

## ----plot-difmap-1comp-colby-sizeby-exprs-------
sce <- runDiffusionMap(sce)
plotDiffusionMap(sce,  
                 colour_by = "Biological_Condition" )

关于SC3包

虽然 SC3包 跟 scater 联系很紧密，但毕竟不属于其知识点，这里就简单运行即可。

预估亚群, 这里显示有24个, 比较多

library(SC3) # BiocManager::install('SC3')
sce <- sc3_estimate_k(sce)
metadata(sce)$sc3$k_estimation
rowData(sce)$feature_symbol <- rownames(rowData(sce))

一步运行sc3的所有分析, 相当耗费时间

这里kn表示的预估聚类数, 考虑到数据集是已知的，我们强行设置为4组, 具体数据要具体考虑。

# 耗费时间
kn <- 4 ## 这里可以选择 3:5 看多种分类结果。
sc3_cluster <- "sc3_4_clusters"
Sys.time()
sce <- sc3(sce, ks = kn, biology = TRUE)
Sys.time()

可视化展示部分, kn就是聚类数

热图: 比较先验分类和SC3的聚类的一致性

sc3_plot_consensus(sce, k = kn, show_pdata = c("Biological_Condition",sc3_cluster))

展示表达量信息

sc3_plot_expression(sce, k = kn, show_pdata =  c("Biological_Condition",sc3_cluster))

展示可能的标记基因

sc3_plot_markers(sce, k = kn, show_pdata =  c("Biological_Condition",sc3_cluster))

在PCA上展示SC3的聚类结果

plotPCA(sce, colour_by =  sc3_cluster )
# sc3_interactive(sce)

显示运行环境

sessionInfo()

2 Seurat Package（最新4.X版本改了很多函数，需要重新查看参考文档）

参考文档

https://satijalab.org/seurat/articles/pbmc3k_tutorial.html

加载R包

rm(list = ls()) # clear the environment
#load all the necessary libraries
options(warn=-1) # turn off warning message globally
suppressMessages(library(Seurat))
# 加载R包，请注意R包版本，可能会有莫名其妙的版本错误
# 单细胞转录组领域发展太快，不同版本的 同一个R包差异很大。

单细胞转录组领域发展太快，不同版本的 同一个R包差异很大，这个Seurat包也不例外，大部分的函数都改了，还专门有一个对照表格供大家学习：

 https://satijalab.org/seurat/essential_commands.html 

也有新的基于3X的教程：（但是现在的版本已经是4.1.1了......）

https://satijalab.org/seurat/pbmc3k_tutorial.html 

创建测试数据集

library(scRNAseq)
## ----- Load Example Data -----
data(fluidigm)
fluidigm<-ReprocessedFluidigmData()
# Set assay to RSEM estimated counts
assay(fluidigm) <- assays(fluidigm)$rsem_counts
ct <- floor(assays(fluidigm)$rsem_counts)
ct[1:4,1:4]

简单看看表达矩阵的性质，主要是基因数量，细胞数量;以及每个细胞表达基因的数量，和每个基因在多少个细胞里面表达。

fivenum(apply(ct,1,function(x) sum(x>0) ))
boxplot(apply(ct,1,function(x) sum(x>0) ))
fivenum(apply(ct,2,function(x) sum(x>0) ))
hist(apply(ct,2,function(x) sum(x>0) ))

names(metadata(fluidigm))
meta <- as.data.frame(colData(fluidigm))
counts <- ct

检测了 counts 和 meta 两个变量，后面需要使用

identical(rownames(meta),colnames(counts))

这里需要把Pollen的表达矩阵做成我们的Seurat要求的对象

Pollen <- CreateSeuratObject(counts = counts, 
                             meta.data =meta,
                             min.cells = 3, 
                             min.genes = 200, 
                             project = "Pollen")
Pollen
## 后续所有的分析都基于这个 Pollen 变量，是一个对象
# An object of class Seurat in project Pollen 

> Pollen
An object of class Seurat 
14656 features across 130 samples within 1 assay 
Active assay: RNA (14656 features, 0 variable features)

检查表达矩阵

将Pollen赋值给sce，目的是代码复用

sce <- Pollen

为原信息增加线粒体基因的比例，如果线粒体基因所占比例过高，意味着这可能是死细胞

由于包的更新，这里运行会提示sce没有名为raw.data的对象，需要进行修改。

mito.genes <- grep(pattern = "^MT-", x = rownames(x = sce@meta.data), value = TRUE)
# 恰好这个例子的表达矩阵里面没有线粒体基因
percent.mito <- Matrix::colSums(sce@raw.data[mito.genes, ]) / Matrix::colSums(sce@raw.data)
## 也可以加入很多其它属性，比如 ERCC 等。

# AddMetaData adds columns to object@meta.data, and is a great place to stash QC stats
sce <- AddMetaData(object = sce, metadata = percent.mito,
                   col.name = "percent.mito")

这里绘图，可以指定分组，前提是这个分组变量存在于meta信息里面

这里的例子是：'Biological_Condition'

VlnPlot(object = sce, features.plot = c("nGene", "nUMI", "percent.mito"), group.by = 'Biological_Condition', nCol = 3)
GenePlot(object = sce, gene1 = "nUMI", gene2 = "nGene")

可以看看高表达量基因是哪些

tail(sort(Matrix::rowSums(sce@raw.data)))
## 散点图可视化任意两个基因的一些属性（通常是细胞的度量）
GenePlot(object = sce, gene1 = "SOX11", gene2 = "EEF1A1")
# 散点图可视化任意两个细胞的一些属性（通常是基因的度量）

CellPlot(sce,sce@cell.names[3],sce@cell.names[4],do.ident = FALSE)

表达矩阵的归一化

起初sce对象里面的data就是原始表达矩阵

identical(sce@raw.data,sce@data)
sce <- NormalizeData(object = sce, 
                     normalization.method = "LogNormalize", 
                     scale.factor = 10000,
                     display.progress = F)

经过了归一化，sce对象里面的data被改变。

identical(sce@raw.data,sce@data)

寻找波动比较明显的基因，后续用这些基因而非全部基因进行分析，主要为了降低计算量。

FindVariableGenes函数当前版本无法调用

sce <- FindVariableGenes(object = sce, mean.function = ExpMean, dispersion.function = LogVMR, 
                         x.low.cutoff = 0.0125, 
                         x.high.cutoff = 3, 
                         y.cutoff = 0.5)
# 通过调整参数可以得到不同数量的 var.genes
length(sce@var.genes)

对归一化后的矩阵进行去除混杂因素和降维

对矩阵进行回归建模，以及scale,主要是为了去除一些文库大小，线粒体基因含量，ERCC含量等因素。

sce <- ScaleData(object = sce, 
                 vars.to.regress = c("nUMI"),##去除文库大小影响
                 display.progress = F)

现在sce对象的 sce@scale.data 也有了数值

运行PCA进行线性降维，这里仅仅是挑选高变化的基因组成的表达矩阵进行PCA分析。

sce <- RunPCA(object = sce, 
              pc.genes = sce@var.genes, 
              do.print = TRUE, 
              pcs.print = 1:5, 
              genes.print = 5)
sce@dr

这样就能拿到PC的基因的重要性占比情况。

tmp <- sce@dr$pca@gene.loadings
VizPCA( sce, pcs.use = 1:2)

PCAPlot(sce, dim.1 = 1, dim.2 = 2,
        group.by = 'Biological_Condition')
sce <- ProjectPCA(sce, do.print = FALSE)

因为细胞数量不多，所以可以全部画出来

PCHeatmap(object = sce, 
          pc.use = 1, 
          cells.use = ncol(sce@data), 
          do.balanced = TRUE, 
          label.columns = FALSE)

PCHeatmap(object = sce, 
          pc.use = 1:10, 
          cells.use = ncol(sce@data), 
          do.balanced = TRUE, 
          label.columns = FALSE)

基于PCA结果看看细胞如何分群

重点: 需要搞懂这里的 resolution 参数，而且降维算法可以选PCA或者ICA ， 分群算法也可以选择。

sce <- FindClusters(object = sce, 
                    reduction.type = "pca", 
                    dims.use = 1:10, force.recalc = T,
                    resolution = 0.9, print.output = 0,
                    save.SNN = TRUE)
PrintFindClustersParams(sce)
table(sce@meta.data$res.0.9)

看单细胞分群后的tSNE图

跟前面的 RunPCA 函数功能差不多，都是为了降维。这里为了节省计算量，首先使用PCA的线性降维结果，再进行tSNE

sce <- RunTSNE(object = sce, 
               dims.use = 1:10, 
               do.fast = TRUE, 
               perplexity=10)
# note that you can set do.label=T to help label individual clusters
TSNEPlot(object = sce)

可以看到，虽然说有4类细胞，但是 GW16和GW21没有区分开来，需要探索参数。

table(meta$Biological_Condition)
table(meta$Biological_Condition,sce@meta.data$res.0.9)
TSNEPlot(object = sce,group.by = 'Biological_Condition')

对每个分类都寻找其marker基因

# 下面的代码是需要适应性修改，因为不同的数据集分组不一样，本次是3组，所以演示3组后的代码。
markers_df <- FindMarkers(object = sce, ident.1 = 0, min.pct = 0.25)
print(x = head(markers_df))
markers_genes =  rownames(head(x = markers_df, n = 5))
VlnPlot(object = sce, features.plot =markers_genes, use.raw = TRUE, y.log = TRUE)
FeaturePlot(object = sce, 
            features.plot =markers_genes, 
            cols.use = c("grey", "blue"), 
            reduction.use = "tsne")

markers_df <- FindMarkers(object = sce, ident.1 = 1, min.pct = 0.25)
print(x = head(markers_df))
markers_genes =  rownames(head(x = markers_df, n = 5))
VlnPlot(object = sce, features.plot =markers_genes, use.raw = TRUE, y.log = TRUE)
FeaturePlot(object = sce, 
            features.plot =markers_genes, 
            cols.use = c("grey", "blue"), 
            reduction.use = "tsne")

markers_df <- FindMarkers(object = sce, ident.1 = 2, min.pct = 0.25)
print(x = head(markers_df))
markers_genes =  rownames(head(x = markers_df, n = 5))
VlnPlot(object = sce, features.plot =markers_genes, use.raw = TRUE, y.log = TRUE)
FeaturePlot(object = sce, 
            features.plot =markers_genes, 
            cols.use = c("grey", "blue"), 
            reduction.use = "tsne") 

展现各个分类的marker基因的表达情况

sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE, min.pct = 0.25, 
                              thresh.use = 0.25)
DT::datatable(sce.markers)
library(dplyr)
sce.markers %>% group_by(cluster) %>% top_n(2, avg_logFC)
top10 <- sce.markers %>% group_by(cluster) %>% top_n(10, avg_logFC)

# setting slim.col.label to TRUE will print just the cluster IDS instead of# every cell name
DoHeatmap(object = sce, genes.use = top10$gene, slim.col.label = TRUE, remove.key = TRUE)

FeaturePlot(object = sce, 
            features.plot =top10$gene, 
            cols.use = c("grey", "blue"), 
            reduction.use = "tsne")

Seurat总结

counts矩阵进来后被包装为对象，方便操作。

然后一定要经过 NormalizeData 和 ScaleData 的操作

函数 FindVariableGenes 可以挑选适合进行下游分析的基因集。

函数 RunPCA 和 RunTSNE 进行降维

函数 FindClusters 直接就分群了，非常方便
函数 FindAllMarkers 可以对分群后各个亚群找标志基因。

函数 FeaturePlot 可以展示不同基因在所有细胞的表达量
函数 VlnPlot 可以展示不同基因在不同分群的表达量差异情况
函数 DoHeatmap 可以选定基因集后绘制热图

显示运行环境

sessionInfo()

3 monocle Package

引言

教程，当然是以官网为主：

http://cole-trapnell-lab.github.io/monocle-release/docs/
http://cole-trapnell-lab.github.io/monocle-release/tutorials/

if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")
if (!requireNamespace("monocle"))
    BiocManager::install("monocle")

加载R包

rm(list = ls()) # clear the environment
#load all the necessary libraries
options(warn=-1) # turn off warning message globally
suppressMessages(library(monocle))

创建数据集

后续分析的前提就是将数据构建成monocle需要的对象。

因此这里先介绍一下monocle需要的用来构建 CellDataSet 对象的三个数据集

• 表达量矩阵exprs:数值矩阵 行名是基因, 列名是细胞编号.
• 细胞的表型信息phenoData: 第一列是细胞编号，其他列是细胞的相关信息
• 基因注释featureData: 第一列是基因编号, 其他列是基因对应的信息

并且这三个数据集要满足如下要求:

表达量矩阵必须：

• 保证它的列数等于phenoData的行数
• 保证它的行数等于featureData的行数

而且

• phenoData的行名需要和表达矩阵的列名匹配
• featureData和表达矩阵的行名要匹配
• featureData至少要有一列"gene_short_name", 就是基因的symbol

如下是几个例子

测试数据集

该R包团队给出的例子有点可怕。

需要使用网络公共数据，取决于网速，是浙江大学郭老师的40万小鼠单细胞转录组数据。

第一个rds文件大小接近3G，太难下载了，即使侥幸下载，一般人的电脑也没办法处理它。

下载地址:

http://trapnell-lab.gs.washington.edu/public_share/MCA_merged_mat.rds
http://trapnell-lab.gs.washington.edu/public_share/MCA_All-batch-removed-assignments.csv

第二个csv文件也有 26.5M

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE72857

MCA <- readRDS("./MCA_merged_mat.rds")
cell.meta.data <- read.csv("./MCA_All-batch-removed-assignments.csv", row.names = 1)

39855 x 405191 sparse Matrix of class "dgCMatrix"

class(MCA)
## 这是一个对象，需要仔细理解，这里略

overlapping_cells <- intersect(row.names(cell.meta.data), colnames(MCA))
gene_ann <- data.frame(gene_short_name = row.names(MCA), row.names = row.names(MCA))

pd <- new("AnnotatedDataFrame",
          data=cell.meta.data[overlapping_cells, ])
fd <- new("AnnotatedDataFrame",
          data=gene_ann)

下面构造了一个MCA数据集的子集，基于：overlapping_cells

MCA_cds <- newCellDataSet(
  MCA[, overlapping_cells], 
  phenoData = pd,
  featureData =fd,
  expressionFamily = negbinomial.size(),
  lowerDetectionLimit=1)

save(MCA_cds,file = './MCA_cds_monocle_example.Rdata')   

但是数据集实在是太大，我这里只能载入前面步骤保存的小数据集。

## 这个文件很大，我放在了Rmd路径下面。
load(file = './MCA_cds_monocle_example.Rdata')

MCA_cds
# 是一个 CellDataSet 对象， 有着24万的单细胞，演示软件用法实在是太难

这里需要学习 CellDataSet 对象，就跟之前GEO视频教程教大家的 ExpressionSet 对象类似，可以看到这里的数据集仍然是24万细胞。

scRNAseq R包中的数据集

我们选择 scRNAseq 这个R包。

这个包内置的是 Pollen et al. 2014 数据集，人类单细胞细胞，分成4类，分别是 pluripotent stem cells 分化而成的 neural progenitor cells (“NPC”) ，还有 “GW16” and “GW21” ，“GW21+3” 这种孕期细胞。

library(scRNAseq)
## ----- Load Example Data -----
data(fluidigm)
# Set assay to RSEM estimated counts
assay(fluidigm)  <-  assays(fluidigm)$rsem_counts
ct <- floor(assays(fluidigm)$rsem_counts)
ct[1:4,1:4] 
sample_ann <- as.data.frame(colData(fluidigm))

准备Monocle对象需要的phenotype data 和 feature data 以及表达矩阵，从 scRNAseq 这个R包里面提取这三种数据

gene_ann <- data.frame(
  gene_short_name = row.names(ct), 
  row.names = row.names(ct)
)

pd <- new("AnnotatedDataFrame",
          data=sample_ann)
fd <- new("AnnotatedDataFrame",
          data=gene_ann)

构建Monocle后续分析的所用对象，主要是根据包的说明书，仔细探索其需要的构建对象的必备元素

注意点: 因为表达矩阵是counts值，所以注意 expressionFamily 参数

sc_cds <- newCellDataSet(
  ct, 
  phenoData = pd,
  featureData =fd,
  expressionFamily = negbinomial.size(),
  lowerDetectionLimit=1)
sc_cds

下面是monocle对新构建的CellDataSet 对象的标准操作

library(dplyr)
colnames(phenoData(sc_cds)@data)
## 为了电脑的健康，我这里选择小数据集。 
sc_cds <- estimateSizeFactors(sc_cds)
sc_cds <- estimateDispersions(sc_cds)

本地加载RPKM数据

你也可以从本地RPKM值文件读入R语言后构造 CellDataSet 对象，下面是简单的例子：

#do not run
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")
# Once these tables are loaded, you can create the CellDataSet object like this:

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)

值得注意的是，因为monocle和前面我们讲解的scater，还有Seurat，它们基于的对象都不一样，所以monocle干脆提供了转换函数：

# 加入你把上面的 HSMM 赋值给 lung ，然后使用函数进行转换：
lung  <-  HSMM
# To convert to Seurat object
lung_seurat <- exportCDS(lung, 'Seurat')

# To convert to SCESet
lung_SCESet <- exportCDS(lung, 'Scater')

直接读取10X结果

因为10X实在是太流行了，所以monocle跟Seurat一样，也提供了直接读取10X上游分析结果的接口函数（其实是使用另外一个R包），因为本文数据来源于smart-seq2，所以并不演示下面的代码：

cellranger_pipestance_path <- "/path/to/your/pipeline/output/directory"
gbm <- load_cellranger_matrix(cellranger_pipestance_path)

fd <- fData(gbm)

# The number 2 is picked arbitrarily in the line below.
# Where "2" is placed you should place the column number that corresponds to your
# featureData's 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())

假设使用monocle3版本

具体教程以官网为主：

http://cole-trapnell-lab.github.io/monocle-release/monocle3/  

这里仍然是演示monocle2，所以下面代码不要运行，仅仅作为演示而已。

devtools::install_github("cole-trapnell-lab/monocle-release", ref="monocle3_alpha")
library(monocle)
cds <- sc_cds
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
cds <- preprocessCDS(cds, num_dim = 20)
cds <- reduceDimension(cds, reduction_method = 'UMAP')
cds <- partitionCells(cds)
cds <- learnGraph(cds,  RGE_method = 'SimplePPT')
plot_cell_trajectory(cds,
                     color_by = "cell_type2") +
                     scale_color_manual(values = cell_type_color)

首先是质控

这里通常也是对基因 和 细胞进行质控，质控指标需要根据项目来进行具体探索，这里只是演示一下用法。

cds=sc_cds
cds
## 起初是： 26255 features, 130 samples 
cds <- detectGenes(cds, min_expr = 0.1)
print(head(fData(cds)))
expressed_genes <- row.names(subset(fData(cds),
                                    num_cells_expressed >= 5))
length(expressed_genes)
cds <- cds[expressed_genes,]
cds
# 过滤基因后是：assayData: 13385 features, 130 samples 
print(head(pData(cds)))
tmp=pData(cds)
fivenum(tmp[,1])
fivenum(tmp[,30])
# 这里并不需要过滤细胞，如果想过滤，就自己摸索阈值，然后挑选细胞即可。
valid_cells <- row.names(pData(cds) )
cds <- cds[,valid_cells]
cds 

聚类

单细胞转录组最重要的就是把细胞分群啦，这里可供选择的算法非常多，我们首先演示PCA结果。

# 并不是所有的基因都有作用，所以先进行挑选，合适的基因用来进行聚类。
disp_table <- dispersionTable(cds)
unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)
cds <- setOrderingFilter(cds, unsup_clustering_genes$gene_id)
plot_ordering_genes(cds) 
plot_pc_variance_explained(cds, return_all = F) # norm_method='log'

# 其中 num_dim 参数选择基于上面的PCA图
cds <- reduceDimension(cds, max_components = 2, num_dim = 6,
                        reduction_method = 'tSNE', verbose = T)
cds <- clusterCells(cds, num_clusters = 4) 
plot_cell_clusters(cds, 1, 2, color = "Biological_Condition")
table(pData(cds)$Biological_Condition)

可以看到，GW21 也是被打散在其它分组里面。

值得注意的是：这里选择不同的PC个数进行后续分析，是会影响聚类结果的。具体可以参考：

https://davetang.org/muse/2017/10/01/getting-started-monocle/

排除干扰因素后聚类

跟前面的质控步骤一样，所谓的干扰因素，也是看自己对数据集的认识情况来自己摸索的，比如我们这里

tmp=pData(cds)
fivenum(tmp[,1])
fivenum(tmp[,30])
colnames(tmp)
# 放在 residualModelFormulaStr 里面的是需要去除的
cds <- reduceDimension(cds, max_components = 2, num_dim = 2,
                        reduction_method = 'tSNE',
                        residualModelFormulaStr = "~Biological_Condition + num_genes_expressed",
                        verbose = T)
cds <- clusterCells(cds, num_clusters = 2)
plot_cell_clusters(cds, 1, 2, color = "Biological_Condition")

## 上面去除了Biological_Condition，所以接下来聚类它们就被打散了。

cds <- reduceDimension(cds, max_components = 2, num_dim = 2,
                        reduction_method = 'tSNE',
                        residualModelFormulaStr = "~NREADS + num_genes_expressed",
                        verbose = T)
cds <- clusterCells(cds, num_clusters = 2)
plot_cell_clusters(cds, 1, 2, color = "Biological_Condition")

差异分析

这个是转录组数据的常规分析了，在单细胞转录组领域也是如此，monocle这个包提供 differentialGeneTest 函数来做差异分析，作用就是挑选那些在某些类别细胞里面高表达的基因，假设其为那一组细胞的marker基因。

在我们的例子里面可以是已知的细胞分类，或者是自己推断的聚类结果。

这一步时间会比较久，也可以不给定分组信息，使用reduced model来计算各个基因的差异与否，是另外一种算法了。

Sys.time()
diff_test_res <- differentialGeneTest(cds,
                                      fullModelFormulaStr = "~Biological_Condition")
Sys.time()
# 可以看到运行耗时

# Select genes that are significant at an FDR < 10%
sig_genes <- subset(diff_test_res, qval < 0.1)
head(sig_genes[,c("gene_short_name", "pval", "qval")] )

cg=as.character(head(sig_genes$gene_short_name))
plot_genes_jitter(cds[cg,], 
                  grouping = "Biological_Condition", ncol= 2)

下面是不同的可视化参数的效果

plot_genes_jitter(cds[cg,],
                  grouping = "Biological_Condition",
                  color_by = "Biological_Condition",
                  nrow= 3,
                  ncol = NULL )

这个时候可以考虑去选用多个差异分析R包来进行不同的比较。

寻找marker基因

使用起来有点复杂，需要预先给定好分组情况，具体自行看说明书。3个R包其实是共通的。

推断发育轨迹

前面介绍的monocle的功能都只能说是中规中矩，而这个推断发育轨迹才是monocle的拿手好戏，也是它荣升为3大R包的核心竞争力。

第一步: 挑选合适的基因。有多个方法，例如提供已知的基因集，这里选取统计学显著的差异基因列表

ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
cds <- setOrderingFilter(cds, ordering_genes)
plot_ordering_genes(cds)
cds <- reduceDimension(cds, max_components = 2,
                            method = 'DDRTree')

第二步: 降维。降维的目的是为了更好的展示数据。函数里提供了很多种方法, 不同方法的最后展示的图都不太一样, 其中“DDRTree”是Monocle2使用的默认方法

cds <- reduceDimension(cds, max_components = 2,
                            method = 'DDRTree')

第三步: 对细胞进行排序

cds <- orderCells(cds)

最后两个可视化函数，对结果进行可视化

plot_cell_trajectory(cds, color_by = "Biological_Condition")  

可以很明显看到细胞的发育轨迹，正好对应 pluripotent stem cells 分化而成的 neural progenitor cells (“NPC”) ，还有 “GW16” and “GW21” ，“GW21+3” 这种孕期细胞。

plot_genes_in_pseudotime可以展现marker基因，本例子随便选取了6个差异表达基因。

plot_genes_in_pseudotime(cds[cg,],
                         color_by = "Biological_Condition")

最开始挑选合适基因，除了我们演示的找统计学显著的差异表达基因这个方法外，还可以于已知的标记基因，主要是基于生物学背景知识。

如果是已知基因列表，就需要自己读取外包文件，导入R里面来分析。

显示运行环境

sessionInfo()

三个包介绍结束。

下一篇我们回到文章数据，使用上述介绍的R包进行分析。

我们下一篇再见！