scRNA包学习Monocle2
前言
第三单元第七讲：使用scRNA包学习Monocle2
课程链接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53
关于monocle2
目前monocle2可以直接利用Bioconductor安装：https://www.bioconductor.org/packages/release/bioc/html/monocle.html
1if (!requireNamespace("BiocManager", quietly = TRUE))
2 install.packages("BiocManager")
3
4BiocManager::install("monocle")
5# 安装的版本是2.12.0
关于这个scRNA包
内容在：https://github.com/jmzeng1314/scRNA_smart_seq2/blob/master/scRNA/study_scRNAseq.html
要使用scRNAseq这个R包，首先要对它进行了解，包中内置了Pollen et al. 2014 的数据集（https://www.nature.com/articles/nbt.2967），到19年8月为止，已经有446引用量了。只不过原文完整的数据是 23730 个基因， 301 个样本【这里只有130个样本文库(高覆盖度、低覆盖度各65个，并且测序深度不同】，这个包中只选取了4种细胞类型：pluripotent stem cells 分化而成的 neural progenitor cells (NPC，神经前体细胞) ，还有 GW16（radial glia，放射状胶质细胞） 、GW21（newborn neuron，新生儿神经元） 、GW21+3(maturing neuron，成熟神经元) ，它们的关系如下图（NPC和其他三类存在较大差别）：
数据大小是50.6 MB，要想知道数据怎么处理的，可以看：https://hemberg-lab.github.io/scRNA.seq.datasets/human/tissues/
加载scRNA包中的数据
1library(scRNAseq)
2data(fluidigm)
3
4> fluidigm
5class: SummarizedExperiment
6dim: 26255 130
7metadata(3): sample_info clusters which_qc
8assays(4): tophat_counts cufflinks_fpkm rsem_counts rsem_tpm
9rownames(26255): A1BG A1BG-AS1 ... ZZEF1 ZZZ3
10rowData names(0):
11colnames(130): SRR1275356 SRR1274090 ... SRR1275366 SRR1275261
12colData names(28): NREADS NALIGNED ... Cluster1 Cluster2
创建CellDataSet对象
=> newCellDataSet()
1# 这个对象需要三个要素
2# 第一个：RSEM表达矩阵（ct = count）
3assay(fluidigm) <- assays(fluidigm)$rsem_counts
4ct <- floor(assays(fluidigm)$rsem_counts)
5ct[1:4,1:4]
6# 第二个：临床信息
7sample_ann <- as.data.frame(colData(fluidigm))
8# 第三个：基因注释信息（必须包含一列是gene_short_name）
9gene_ann <- data.frame(
10 gene_short_name = row.names(ct),
11 row.names = row.names(ct)
12)
13# 然后转换为AnnotatedDataFrame对象
14pd <- new("AnnotatedDataFrame",
15 data=sample_ann)
16fd <- new("AnnotatedDataFrame",
17 data=gene_ann)
18
19# 最后构建CDS对象
20sc_cds <- newCellDataSet(
21 ct,
22 phenoData = pd,
23 featureData =fd,
24 expressionFamily = negbinomial.size(),
25 lowerDetectionLimit=1)
26sc_cds
27# CellDataSet (storageMode: environment)
28# assayData: 26255 features, 130 samples
29# element names: exprs
30# protocolData: none
31# phenoData
32# sampleNames: SRR1275356 SRR1274090 ... SRR1275261 (130 total)
33# varLabels: NREADS NALIGNED ... Size_Factor (29 total)
34# varMetadata: labelDescription
35# featureData
36# featureNames: A1BG A1BG-AS1 ... ZZZ3 (26255 total)
37# fvarLabels: gene_short_name
38# fvarMetadata: labelDescription
39# experimentData: use 'experimentData(object)'
40# Annotation:
注意到构建CDS对象过程中有一个参数是：expressionFamily ，它是选择了一个数据分布，例如FPKM/TPM 值是log-正态分布的；UMIs和原始count值用负二项分布模拟的效果更好。负二项分布有两种方法，这里选用了negbinomial.size，另外一种negbinomial稍微更准确一点，但速度大打折扣，它主要针对非常小的数据集
质控过滤
=> detectGenes()
1cds=sc_cds
2cds ## 原始数据有：26255 features, 130 samples
3# 设置一个基因表达量的过滤阈值，结果会在cds@featureData@data中新增一列num_cells_expressed，记录这个基因在多少细胞中有表达
4cds <- detectGenes(cds, min_expr = 0.1)
5# 结果保存在cds@featureData@data
6print(head(cds@featureData@data))
7# gene_short_name num_cells_expressed
8# A1BG A1BG 10
9# A1BG-AS1 A1BG-AS1 2
10# A1CF A1CF 1
11# A2M A2M 21
12# A2M-AS1 A2M-AS1 3
13# A2ML1 A2ML1 9
在monocle版本2.12.0中，取消了fData函数（此前在2.10版本中还存在），不过在monocle3中又加了回来
如果遇到不能使用fData的情况，就可以采用备选方案：cds@featureData@data
然后可以进行基因过滤 =>subset()
1expressed_genes <- row.names(subset(cds@featureData@data,
2 num_cells_expressed >= 5))
3length(expressed_genes)
4## [1] 13385
5cds <- cds[expressed_genes,]
还可以进行细胞层面的过滤（可选）
1# 依然是：如果不支持使用pData()函数，可以使用cds@phenoData@data来获得各种细胞注释信息
2print(head(cds@phenoData@data))
3# 比如我们看一下细胞注释的第一个NREADS信息
4tmp=pData(cds)
5fivenum(tmp[,1])
6## [1] 91616 232899 892209 8130850 14477100
7
8# 如果要过滤细胞，其实也是利用subset函数，不过这里不会对细胞过滤
9valid_cells <- row.names(cds@phenoData@data)
10cds <- cds[,valid_cells]
11cds
12## CellDataSet (storageMode: environment)
13## assayData: 13385 features, 130 samples
14## element names: exprs
15## protocolData: none
16## phenoData
17## sampleNames: SRR1275356 SRR1274090 ... SRR1275261 (130 total)
18## varLabels: NREADS NALIGNED ... num_genes_expressed (30 total)
19## varMetadata: labelDescription
20## featureData
21## featureNames: A1BG A2M ... ZZZ3 (13385 total)
22## fvarLabels: gene_short_name num_cells_expressed
23## fvarMetadata: labelDescription
24## experimentData: use 'experimentData(object)'
25## Annotation:
聚类
在monocle3中，聚类使用的是
cluster_cells()，利用Louvain community detection的非监督聚类方法，结果保存在cds@clusters$UMAP$clusters
不使用marker基因聚类
使用函数clusterCells()，根据整体的表达量对细胞进行分组。例如，细胞表达了大量的与成肌细胞相关的基因，但就是没有成肌细胞的marker--MYF5 ，我们依然可以判断这个细胞属于成肌细胞。
step1：dispersionTable()
首先就是判断使用哪些基因进行细胞分群。当然，可以使用全部基因，但这会掺杂很多表达量不高而检测不出来的基因，反而会增加噪音。挑有差异的，挑表达量不太低的
1cds <- estimateSizeFactors(cds)
2cds <- estimateDispersions(cds)
3disp_table <- dispersionTable(cds) # 挑有差异的
4unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1) # 挑表达量不太低的
5cds <- setOrderingFilter(cds, unsup_clustering_genes$gene_id) # 准备聚类基因名单
6plot_ordering_genes(cds)
7# 图中黑色的点就是被标记出来一会要进行聚类的基因
step2：plot_pc_variance_explained()
然后选一下主成分
1plot_pc_variance_explained(cds, return_all = F) # norm_method='log'
step3:
根据上面👆的图，选择合适的主成分数量（这个很主观，可以多试几次），这里选前6个成分（大概在第一个拐点处）
1# 进行降维
2cds <- reduceDimension(cds, max_components = 2, num_dim = 6,
3 reduction_method = 'tSNE', verbose = T)
4# 进行聚类
5cds <- clusterCells(cds, num_clusters = 4)
6# Distance cutoff calculated to 0.5225779
7plot_cell_clusters(cds, 1, 2, color = "Biological_Condition")
1> table(cds@phenoData@data$Biological_Condition)
2
3 GW16 GW21 GW21+3 NPC
4 52 16 32 30
需要注意的是，使用的主成分数量会影响结果
前面使用了6个主成分，分的还不错。现在假设使用前16个主成分：
1cds <- reduceDimension(cds, max_components = 2, num_dim = 16,
2 reduction_method = 'tSNE', verbose = T)
3cds <- clusterCells(cds, num_clusters = 4)
4plot_cell_clusters(cds, 1, 2, color = "Biological_Condition")
以下是测试代码！
除此以外，如果有批次效应等干扰因素，也可以在降维（reduceDimension()）的过程中进行排除：
1if(F){
2 cds <- reduceDimension(cds, max_components = 2, num_dim = 6,
3 reduction_method = 'tSNE',
4 residualModelFormulaStr = "~Biological_Condition + num_genes_expressed",
5 verbose = T)
6 cds <- clusterCells(cds, num_clusters = 4)
7 plot_cell_clusters(cds, 1, 2, color = "Biological_Condition")
8}
9# 可以看到，去掉本来的生物学意义后，最后细胞是会被打散的。所以residualModelFormulaStr这个东西的目的就是磨平它参数包含的差异
但是，如果是去除其他的效应：
1# 如果去除生物意义以外的效应
2 cds <- reduceDimension(cds, max_components = 2, num_dim = 6,
3 reduction_method = 'tSNE',
4 residualModelFormulaStr = "~NREADS + num_genes_expressed",
5 verbose = T)
6 cds <- clusterCells(cds, num_clusters = 4)
7 plot_cell_clusters(cds, 1, 2, color = "Biological_Condition")
关于处理批次效应：例如在芯片数据中经常会利用SVA的combat函数。
磨平批次效应实际上就是去掉各个组的前几个主成分
差异分析
=> differentialGeneTest()
1start=Sys.time()
2diff_test_res <- differentialGeneTest(cds,
3 fullModelFormulaStr = "~Biological_Condition")
4end=Sys.time()
5end-start
6# 运行时间在几分钟至十几分钟不等
然后得到差异基因
1sig_genes <- subset(diff_test_res, qval < 0.1)
2> head(sig_genes[,c("gene_short_name", "pval", "qval")] )
3 gene_short_name pval qval
4A1BG A1BG 4.112065e-04 1.460722e-03
5A2M A2M 4.251744e-08 4.266086e-07
6AACS AACS 2.881832e-03 8.275761e-03
7AADAT AADAT 1.069794e-02 2.621123e-02
8AAGAB AAGAB 1.156771e-07 1.021331e-06
9AAMP AAMP 7.626789e-05 3.243869e-04
作图（注意要将基因名变成character）
1cg=as.character(head(sig_genes$gene_short_name))
2# 普通图
3plot_genes_jitter(cds[cg,],
4 grouping = "Biological_Condition", ncol= 2)
5# 还能上色
6plot_genes_jitter(cds[cg,],
7 grouping = "Biological_Condition",
8 color_by = "Biological_Condition",
9 nrow= 3,
10 ncol = NULL )
我们自己也可以根据某个基因的表达量差异和分组信息进行作图（就以A1BG为例）：
1# 以A1BG为例
2boxplot(log10(cds@assayData$exprs["A1BG",]+1) ~ cds@phenoData@data$Biological_Condition)
推断发育轨迹
三步走：从差异分析结果选合适基因=》降维=》细胞排序
step1: 选合适基因
1ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
2cds <- setOrderingFilter(cds, ordering_genes)
3plot_ordering_genes(cds)
step2: 降维
1# 默认使用DDRTree的方法
2cds <- reduceDimension(cds, max_components = 2,
3 method = 'DDRTree')
step3: 细胞排序
1cds <- orderCells(cds)
最后可视化
1plot_cell_trajectory(cds, color_by = "Biological_Condition")
这个图就可以看到细胞的发展过程
另外，plot_genes_in_pseudotime 可以对基因在不同细胞中的表达量变化进行绘图
1plot_genes_in_pseudotime(cds[cg,],
2 color_by = "Biological_Condition")