1.利用已有的Seurat对象转换为CellDataSet对象
加载之前seutat标准流程跑完的数据，如果没有可以走一遍流程，教程点我
#加载需要的包
library(Seurat)
# devtools::install_github('satijalab/seurat-data')
library(SeuratData)
library(monocle)
library(ggplot2)
library(patchwork)
library(dplyr)
#加载之前seutat标准流程跑完的数据，如果没有可以走一遍流程
load(file = '../section-01-cluster/basic.sce.pbmc.Rdata')
pbmc
table(Idents(pbmc))
DimPlot(pbmc, reduction = 'umap', label = TRUE, pt.size = 0.5) + NoLegend()
sce=pbmc
table( Idents(sce ))
table(sce@meta.data$seurat_clusters)
table(sce@meta.data$orig.ident)
将之前的数据调出来后，我们看一下monocle包的newCellDataSet函数。
?newCellDataSet
#Creates a new CellDateSet object.
#cellData expression data matrix for an experiment
#phenoData data frame containing attributes of individual cells
#featureData data frame containing attributes of features (e.g. genes)
可以看到这三个数据在sce中都可以调出来
head(sce@meta.data )
sample_ann <- sce@meta.data
sample_ann$celltype=Idents(sce)
head(sample_ann)
# rownames(sample_ann)=sample_ann[,1]
gene_ann <- data.frame(gene_short_name = rownames(sce@assays$RNA) , row.names = rownames(sce@assays$RNA)
)
head(gene_ann)
#newCellDataSet要求featureData和phenoData格式为AnnotatedDataFrame
pd <- new("AnnotatedDataFrame",data=sample_ann)
fd <- new("AnnotatedDataFrame",data=gene_ann)
ct=as.data.frame(sce@assays$RNA@counts)
ct[1:4,1:4]
创建CellDataSet对象
sc_cds <- newCellDataSet(as.matrix(ct), phenoData = pd,featureData =fd,expressionFamily = negbinomial.size(),lowerDetectionLimit=1)
sc_cds
convert the Seurat object to a SingleCellExperiment object
pbmc_sce <- as.SingleCellExperiment(pbmc)
convert singlecellexperment to Seurat
as.Seurat(对象名)
Seurat::as.CellDataSet()函数可以直接将Seurat对象转化为monocle2的对象，进行monocle2的拟时分析
sc_cds<-Seurat::as.CellDataSet(sce)
2.monocle的标准流程
2.1 数据过滤
library(monocle)
sc_cds
colnames(fData(sc_cds))
colnames(pData(sc_cds))
sc_cds <- detectGenes(sc_cds, min_expr = 1)
colnames(fData(sc_cds))
colnames(pData(sc_cds))
#detectGenes设置一个阈值，在在此阈值之上的 CellDataSet对象中，计算在每个细胞中基因表达的数量pdata(cds)$num_genes_expressed,计算每个基因在多少细胞中表达fdata(cds)$num_cells_expressed。
sc_cds <- sc_cds[fDat(sc_cds$num_cells_expressed > 10, ]
# 数值可以自行摸索 featureData返回一个包含变量值和变量元数据信息的对象。 fvarLabels 返回测量变量名称的字符向量。 fData 返回一个数据框，其中基因为行，变量为列。 fvarMetadata 返回一个数据框，其中变量名称作为行，描述标签（例如，测量单位）作为列。
sc_cds
cds <- sc_cds
2.2 标准化和归一化
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
2.3 挑选表达量不太低的基因用于后续分析
# dispersionTable差异基因，并不是所有的基因都有作用，所以先进行挑选，合适的基因用来进行聚类。
disp_table <- dispersionTable(cds)
## 挑表达量不太低的基因，在各细胞中平均表达量大于0.1
unsup_clustering_genes <- subset(disp_table,mean_expression >= 0.1)
unsup_clustering_genes
## 准备聚类基因名单
cds <- setOrderingFilter(cds, unsup_clustering_genes$gene_id)
#按均值与离散度绘制基因，突出显示选择排序的基因，每个灰点都是一个基因。黑点是包含在最后一次调用 setOrderingFilter中的那些。红色曲线显示通过estimateDispersions() 学习的均方差模型。
plot_ordering_genes(cds)
#从归一化后的数据中绘制基于 PCA 的每个主成分的方差百分比
plot_pc_variance_explained(cds,
2.4 降维聚类分群
# 降维，其中 num_dim 参数选择基于上面的PCA图
cds <- reduceDimension(cds, max_components = 2, num_dim = 6,reduction_method = 'tSNE', verbose = T)
#聚类
cds <- clusterCells(cds, num_clusters = 6)
plot_cell_clusters(cds, 1, 2 )
2.5 查看monocle和seurat聚类分群的差别
#phenoData 返回一个包含变量值和变量元数据信息的对象。 varLabels 返回测量变量的字符向量。 pData 返回一个数据框，其中样本为行，变量为列。 varMetadata 返回一个数据框，其中变量名称为行，描述标签（例如，测量单位）为列。
head(rownames(pData(cds)))
colnames(pData(cds))
table(pData(cds)$Cluster)
table(pData(cds)$seurat_clusters)
table(pData(cds)$Cluster,pData(cds)$seurat_clusters)
table(pData(cds)$Cluster,pData(cds)$celltype)
可以看到 monocle 给细胞重新定义了亚群，亚群数量可以自己选择，整体来说，monocle和seurat 各自独立流程定义的亚群的一致性还不错。
