第三单元第十一讲：使用monocle2分析文章数据
课程链接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53
载入数据，创建对象
代码语言：javascript复制rm(list = ls())
Sys.setenv(R_MAX_NUM_DLLS=999)
## 首先载入文章的数据
load(file='../input.Rdata')
# 原始表达矩阵
counts=a
counts[1:4,1:4];dim(counts) # 2万多个基因，768个细胞(需要下一步进行过滤)
library(stringr)
# 样本信息
meta=df
> head(meta)
g plate n_g all
SS2_15_0048_A3 1 0048 2624 all
SS2_15_0048_A6 1 0048 2664 all
SS2_15_0048_A5 2 0048 3319 all
SS2_15_0048_A4 3 0048 4447 all
SS2_15_0048_A1 2 0048 4725 all
SS2_15_0048_A2 3 0048 5263 all
# 按行/基因检查：每个基因在多少细胞中有表达量
fivenum(apply(counts,1,function(x) sum(x>0) ))
boxplot(apply(counts,1,function(x) sum(x>0) ))
# 按列/样本检查：每个细胞中存在多少表达的基因
fivenum(apply(counts,2,function(x) sum(x>0) ))
hist(apply(counts,2,function(x) sum(x>0) ))
按行 按列检查结果
左图显示了：有50%的基因只在低于20个细胞中有表达量，还有许多没有表达量的：
代码语言：javascript复制> table(apply(counts,1,function(x) sum(x>0) )==0)
FALSE TRUE
17429 7153
# 存在7000多个基因在任何一个细胞中都没表达
右图显示了：大部分细胞都包含2000个以上的有表达的基因
开始使用newCellDataSet()构建monocle对象：
# ---------首先构建基因的注释信息(feature_data)-----------
gene_ann <- data.frame(
gene_short_name = row.names(counts),
row.names = row.names(counts)
)
fd <- new("AnnotatedDataFrame",
data=gene_ann)
# ---------然后构建样本的注释信息(sample_data)---------
sample_ann=meta
pd <- new("AnnotatedDataFrame",
data=sample_ann)
# ---------开始构建对象---------
sc_cds <- newCellDataSet(
as.matrix(counts),
phenoData = pd,
featureData =fd,
expressionFamily = negbinomial.size(),
lowerDetectionLimit=1)
sc_cds
# CellDataSet (storageMode: environment)
# assayData: 24582 features, 768 samples
# element names: exprs
# protocolData: none
# phenoData
# sampleNames: SS2_15_0048_A3 SS2_15_0048_A6 ... SS2_15_0049_P24 (768 total)
# varLabels: g plate ... Size_Factor (5 total)
# varMetadata: labelDescription
# featureData
# featureNames: 0610005C13Rik 0610007P14Rik ... ERCC-00171 (24582 total)
# fvarLabels: gene_short_name
# fvarMetadata: labelDescription
# experimentData: use 'experimentData(object)'
# Annotation:
接下来质控过滤
=> detectGenes()
代码语言：javascript复制cds=sc_cds
cds
## 起初是 24582 features, 768 samples
#---------首先是对基因的过滤-------------
cds <- detectGenes(cds, min_expr = 0.1)
print(head(cds@featureData@data))
expressed_genes <- row.names(subset(cds@featureData@data,
num_cells_expressed >= 5))
length(expressed_genes)
# 14442
# 这里需要去掉ERCC基因
# 去掉ERCC基因
is.ercc <- grepl("ERCC",expressed_genes)
length(expressed_genes[!is.ercc])
# 14362(看到去掉了80个ERCC)
cds <- cds[expressed_genes[!is.ercc],]
cds
# 过滤基因后是 14362 features, 768 samples
#---------然后是对细胞的过滤-------------
# 如果不支持使用pData()函数，可以使用cds@phenoData@data来获得各种细胞注释信息
cell_anno <- cds@phenoData@data
> head(cell_anno)
g plate n_g all Size_Factor num_genes_expressed
SS2_15_0048_A3 1 0048 2624 all 0.6693919 3069
SS2_15_0048_A6 1 0048 2664 all 0.6820532 3040
SS2_15_0048_A5 2 0048 3319 all 1.0759418 3743
SS2_15_0048_A4 3 0048 4447 all 0.7294812 5014
SS2_15_0048_A1 2 0048 4725 all 1.5658507 5128
SS2_15_0048_A2 3 0048 5263 all 1.6187569 5693
# 这里简单过滤细胞
valid_cells <- row.names(cell_anno[cell_anno$num_genes_expressed>2000,] )
cds <- cds[,valid_cells]
cds
# 最后剩下：14362 features, 693 samples
然后进行归一化
=> estimateSizeFactors()
代码语言：javascript复制library(dplyr)
colnames(phenoData(cds)@data)
## 必要的归一化
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
然后降维聚类
step1：dispersionTable()
代码语言：javascript复制disp_table <- dispersionTable(cds)
unsup_clustering_genes <- subset(disp_table,
mean_expression >= 0.1)
cds <- setOrderingFilter(cds, unsup_clustering_genes$gene_id)
cds
plot_ordering_genes(cds)
# 图中黑色的点就是被标记出来一会要进行聚类的基因
step2：plot_pc_variance_explained()
代码语言：javascript复制plot_pc_variance_explained(cds, return_all = F) # norm_method='log'
step3:
关于聚类：monocle2一共有三种方法：densityPeak", "louvain", "DDRTree"
默认使用density peak的方法，但是对于大型数据（例如有5万细胞）推荐louvain方法
# ---------进行降维---------
cds <- reduceDimension(cds, max_components = 2, num_dim = 6,
reduction_method = 'tSNE', verbose = T)
# ---------进行聚类---------
# (这里先设置num_clusters，不一定最后真就分成5群；我们这里设置5，最后只能得到4群；如果设成4，结果就得到3群)
cds <- clusterCells(cds, num_clusters = 5)
# Distance cutoff calculated to 1.818667
# color使用的这些数据就在：cds$Cluster
plot_cell_clusters(cds, 1, 2, color = "Cluster")
因为之前还做过层次聚类，所以还可以对比一下：
代码语言：javascript复制plot_cell_clusters(cds, 1, 2, color = "g")
看到monocle使用其他的聚类算法确实不如使用自己的结果得到的效果好
再次对比不同分群结果的基因数量差异：
代码语言：javascript复制boxplot(cds@phenoData@data$num_genes_expressed~cds@phenoData@data$Cluster)
boxplot(cds@phenoData@data$num_genes_expressed~cds@phenoData@data$g)
去除一些影响因素
因为这几群的细胞中基因表达数量是有差别的，因此我们可以在聚类之前先去掉这部分影响，从而关注它们真正的生物学影响
代码语言：javascript复制## 去除检测到基因数量效应
cds <- reduceDimension(cds, max_components = 2, num_dim = 2,
reduction_method = 'tSNE',
residualModelFormulaStr = "~num_genes_expressed",
verbose = T)
cds <- clusterCells(cds, num_clusters = 5)
plot_cell_clusters(cds, 1, 2, color = "Cluster")
差异分析
=> differentialGeneTest()
代码语言：javascript复制start=Sys.time()
diff_test_res <- differentialGeneTest(cds,
fullModelFormulaStr = "~Cluster")
end=Sys.time()
end-start
# Time difference of 4.724045 mins
挑差异基因
代码语言：javascript复制# 选择FDR-adjusted p-value(也就是q值) < 10%的基因作为差异基因
sig_genes <- subset(diff_test_res, qval < 0.1)
dim(sig_genes)
> head(sig_genes[,c("gene_short_name", "pval", "qval")] )
gene_short_name pval qval
0610007P14Rik 0610007P14Rik 3.377244e-03 1.277429e-02
0610010F05Rik 0610010F05Rik 1.243943e-02 3.924761e-02
0610011F06Rik 0610011F06Rik 2.338530e-03 9.285587e-03
1110004E09Rik 1110004E09Rik 2.487600e-02 7.003903e-02
1110020A21Rik 1110020A21Rik 2.318327e-02 6.618129e-02
1110059E24Rik 1110059E24Rik 5.193533e-09 4.856089e-08
作图
代码语言：javascript复制cg=as.character(head(sig_genes$gene_short_name))
# 还能上色
plot_genes_jitter(cds[cg,],
grouping = "Cluster",
color_by = "Cluster",
nrow= 3,
ncol = NULL )
推断发育轨迹
step1: 选合适基因 => setOrderingFilter()
代码语言：javascript复制ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
cds <- setOrderingFilter(cds, ordering_genes)
plot_ordering_genes(cds)
step2: 降维 => reduceDimension()
代码语言：javascript复制# 默认使用DDRTree的方法
cds <- reduceDimension(cds, max_components = 2,
method = 'DDRTree')
step3: 细胞排序 => orderCells()
代码语言：javascript复制cds <- orderCells(cds)
最后可视化
代码语言：javascript复制plot_cell_trajectory(cds, color_by = "Cluster") 
