使用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 )

推断发育轨迹

step1: 选合适基因 => setOrderingFilter()

1ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
2cds <- setOrderingFilter(cds, ordering_genes)
3plot_ordering_genes(cds)

step2: 降维 => reduceDimension()

1# 默认使用DDRTree的方法 
2cds <- reduceDimension(cds, max_components = 2,
3                            method = 'DDRTree')

step3: 细胞排序 => orderCells()

1cds <- orderCells(cds)

最后可视化

1plot_cell_trajectory(cds, color_by = "Cluster")  

“

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

”

---