使用monocle2分析文章数据 / 开普饭

课程笔记

粉丝：有单细胞线上课程吗？

小编：什么？我们的单细胞转录组分析线上课程已经上线好久了，你们竟然都不知道吗，每篇推文后面的课程推荐没人看的吗，小编已哭晕在厕所

好了，戏演完了，下面郑重介绍下我们的单细胞线上课程：（详情戳下方链接）

全网第二个单细胞视频课程预售

这个课程笔记栏目记录了学员们学习单细胞转录组课程的学习笔记

希望大家能有所收获！

前言

第三单元第十一讲：使用monocle2分析文章数据

课程链接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53

载入数据，创建对象

1rm(list = ls()) 2Sys.setenv(R_MAX_NUM_DLLS=999) 3## 首先载入文章的数据 4load(file='../input.Rdata') 5# 原始表达矩阵 6counts=a 7counts[1:4,1:4];dim(counts) # 2万多个基因，768个细胞(需要下一步进行过滤) 8library(stringr) 9# 样本信息 10meta=df 11> head(meta) 12 g plate n_g all 13SS2_15_0048_A3 1 0048 2624 all 14SS2_15_0048_A6 1 0048 2664 all 15SS2_15_0048_A5 2 0048 3319 all 16SS2_15_0048_A4 3 0048 4447 all 17SS2_15_0048_A1 2 0048 4725 all 18SS2_15_0048_A2 3 0048 5263 all 19# 按行/基因检查：每个基因在多少细胞中有表达量 20fivenum(apply(counts,1,function(x) sum(x>0) )) 21boxplot(apply(counts,1,function(x) sum(x>0) )) 22# 按列/样本检查：每个细胞中存在多少表达的基因 23fivenum(apply(counts,2,function(x) sum(x>0) )) 24hist(apply(counts,2,function(x) sum(x>0) ))

按行+按列检查结果

左图显示了：有50%的基因只在低于20个细胞中有表达量，还有许多没有表达量的：

1> table(apply(counts,1,function(x) sum(x>0) )==0) 2 3FALSE TRUE 417429 7153 5# 存在7000多个基因在任何一个细胞中都没表达

右图显示了：大部分细胞都包含2000个以上的有表达的基因

开始使用newCellDataSet()构建monocle对象：

1# ---------首先构建基因的注释信息(feature_data)----------- 2gene_ann <- data.frame( 3 gene_short_name = row.names(counts), 4 row.names = row.names(counts) 5) 6fd <- new("AnnotatedDataFrame", 7 data=gene_ann) 8# ---------然后构建样本的注释信息(sample_data)--------- 9sample_ann=meta 10pd <- new("AnnotatedDataFrame", 11 data=sample_ann) 12 13# ---------开始构建对象--------- 14sc_cds <- newCellDataSet( 15 as.matrix(counts), 16 phenoData = pd, 17 featureData =fd, 18 expressionFamily = negbinomial.size(), 19 lowerDetectionLimit=1) 20sc_cds 21# CellDataSet (storageMode: environment) 22# assayData: 24582 features, 768 samples 23# element names: exprs 24# protocolData: none 25# phenoData 26# sampleNames: SS2_15_0048_A3 SS2_15_0048_A6 ... SS2_15_0049_P24 (768 total) 27# varLabels: g plate ... Size_Factor (5 total) 28# varMetadata: labelDescription 29# featureData 30# featureNames: 0610005C13Rik 0610007P14Rik ... ERCC-00171 (24582 total) 31# fvarLabels: gene_short_name 32# fvarMetadata: labelDescription 33# experimentData: use 'experimentData(object)' 34# Annotation:

接下来质控过滤

=> detectGenes()

1cds=sc_cds 2cds 3## 起初是 24582 features, 768 samples 4#---------首先是对基因的过滤------------- 5cds <- detectGenes(cds, min_expr = 0.1) 6print(head(cds@featureData@data)) 7expressed_genes <- row.names(subset(cds@featureData@data, 8 num_cells_expressed >= 5)) 9length(expressed_genes) 10# 14442 11# 这里需要去掉ERCC基因 12# 去掉ERCC基因 13is.ercc <- grepl("ERCC",expressed_genes) 14length(expressed_genes[!is.ercc]) 15# 14362(看到去掉了80个ERCC) 16cds <- cds[expressed_genes[!is.ercc],] 17cds 18# 过滤基因后是 14362 features, 768 samples 19#---------然后是对细胞的过滤------------- 20# 如果不支持使用pData()函数，可以使用cds@phenoData@data来获得各种细胞注释信息 21cell_anno <- cds@phenoData@data 22> head(cell_anno) 23 g plate n_g all Size_Factor num_genes_expressed 24SS2_15_0048_A3 1 0048 2624 all 0.6693919 3069 25SS2_15_0048_A6 1 0048 2664 all 0.6820532 3040 26SS2_15_0048_A5 2 0048 3319 all 1.0759418 3743 27SS2_15_0048_A4 3 0048 4447 all 0.7294812 5014 28SS2_15_0048_A1 2 0048 4725 all 1.5658507 5128 29SS2_15_0048_A2 3 0048 5263 all 1.6187569 5693 30# 这里简单过滤细胞 31valid_cells <- row.names(cell_anno[cell_anno$num_genes_expressed>2000,] ) 32cds <- cds[,valid_cells] 33cds 34# 最后剩下：14362 features, 693 samples

然后进行归一化

=> estimateSizeFactors()

1library(dplyr) 2colnames(phenoData(cds)@data) 3## 必要的归一化 4cds <- estimateSizeFactors(cds) 5cds <- estimateDispersions(cds)

然后降维聚类

step1：dispersionTable()

1disp_table <- dispersionTable(cds) 2unsup_clustering_genes <- subset(disp_table, 3 mean_expression >= 0.1) 4cds <- setOrderingFilter(cds, unsup_clustering_genes$gene_id) 5cds 6plot_ordering_genes(cds) 7# 图中黑色的点就是被标记出来一会要进行聚类的基因

step2：plot_pc_variance_explained()

1plot_pc_variance_explained(cds, return_all = F) # norm_method='log'

step3:

关于聚类：monocle2一共有三种方法：densityPeak", "louvain", "DDRTree"

默认使用density peak的方法，但是对于大型数据（例如有5万细胞）推荐louvain方法

1# ---------进行降维--------- 2cds <- reduceDimension(cds, max_components = 2, num_dim = 6, 3 reduction_method = 'tSNE', verbose = T) 4# ---------进行聚类--------- 5# (这里先设置num_clusters，不一定最后真就分成5群；我们这里设置5，最后只能得到4群；如果设成4，结果就得到3群) 6cds <- clusterCells(cds, num_clusters = 5) 7# Distance cutoff calculated to 1.818667 8 9# color使用的这些数据就在：cds$Cluster 10plot_cell_clusters(cds, 1, 2, color = "Cluster")

因为之前还做过层次聚类，所以还可以对比一下：

1plot_cell_clusters(cds, 1, 2, color = "g")

看到monocle使用其他的聚类算法确实不如使用自己的结果得到的效果好

再次对比不同分群结果的基因数量差异：

1boxplot(cds@phenoData@data$num_genes_expressed~cds@phenoData@data$Cluster) 2boxplot(cds@phenoData@data$num_genes_expressed~cds@phenoData@data$g)

去除一些影响因素

因为这几群的细胞中基因表达数量是有差别的，因此我们可以在聚类之前先去掉这部分影响，从而关注它们真正的生物学影响

1## 去除检测到基因数量效应 2cds <- reduceDimension(cds, max_components = 2, num_dim = 2, 3 reduction_method = 'tSNE', 4 residualModelFormulaStr = "~num_genes_expressed", 5 verbose = T) 6cds <- clusterCells(cds, num_clusters = 5) 7plot_cell_clusters(cds, 1, 2, color = "Cluster")

差异分析

=> differentialGeneTest()

1start=Sys.time() 2diff_test_res <- differentialGeneTest(cds, 3 fullModelFormulaStr = "~Cluster") 4end=Sys.time() 5end-start 6# Time difference of 4.724045 mins

挑差异基因

1# 选择FDR-adjusted p-value(也就是q值) < 10%的基因作为差异基因 2sig_genes <- subset(diff_test_res, qval < 0.1) 3dim(sig_genes) 4> head(sig_genes[,c("gene_short_name", "pval", "qval")] ) 5 gene_short_name pval qval 60610007P14Rik 0610007P14Rik 3.377244e-03 1.277429e-02 70610010F05Rik 0610010F05Rik 1.243943e-02 3.924761e-02 80610011F06Rik 0610011F06Rik 2.338530e-03 9.285587e-03 91110004E09Rik 1110004E09Rik 2.487600e-02 7.003903e-02 101110020A21Rik 1110020A21Rik 2.318327e-02 6.618129e-02 111110059E24Rik 1110059E24Rik 5.193533e-09 4.856089e-08

作图

1cg=as.character(head(sig_genes$gene_short_name)) 2# 还能上色 3plot_genes_jitter(cds[cg,], 4 grouping = "Cluster", 5 color_by = "Cluster", 6 nrow= 3, 7 ncol = NULL )