根据表达矩阵进行分群-1

Original 刘小泽单细胞天地 2022-06-07

收录于合集 #第二期单细胞视频笔记汇总 16个

课程笔记

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

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

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

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

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

希望大家能有所收获！

作者 | 单细胞天地小编刘小泽

课程链接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=55
这次会介绍如何针对表达矩阵进行分群，不一定需要包装好的函数。对应视频第三单元5-6讲

前言

目的是得到这个图

将会用到来自作者的包装好的analysis_functions.R代码：

https://github.com/IStevant/XX-XY-mouse-gonad-scRNA-seq/blob/master/scripts/analysis_functions.R

这个代码有1800多行，将会贯穿整个分析，正是这些DIY的代码，才让文章的图显得与众不同

1 首先创造表达矩阵

先下载作者上游定量处理好的数据：female_rpkm.Robj https://github.com/IStevant/XX-XY-mouse-gonad-scRNA-seq/raw/master/data/female_rpkm.Robj
一会要用的基因列表：https://github.com/IStevant/XX-XY-mouse-gonad-scRNA-seq/blob/master/data/prot_coding.csv

 1load(file="female_rpkm.Robj")
 2
 3## 去掉重复细胞
 4#（例如：同一个细胞建库两次，这里作者用“rep”进行了标记） 
 5grep("rep",colnames(female_rpkm))
 6colnames(female_rpkm)[256:257]
 7female_rpkm <- female_rpkm[,!colnames(female_rpkm) %in% grep("rep",colnames(female_rpkm), value=TRUE)]
 8
 9## 只保留编码基因(去掉类似：X5430419D17Rik、BC003331等)
10prot_coding_genes <- read.csv(file="prot_coding.csv", row.names=1)
11females <- female_rpkm[rownames(female_rpkm) %in% as.vector(prot_coding_genes$x),]
12save(females,file = 'female_rpkm.Rdata')

2 然后使用包装好的代码进行tSNE

2.1 对细胞操作=》细胞发育时期的获取

细胞是从6个时间点取出的，于是先找到这6个时间点

 1load('../female_rpkm.Rdata')
 2> dim(females)
 3[1] 21083   563
 4
 5> head(colnames(females))
 6[1] "E10.5_XX_20140505_C01_150331_1" "E10.5_XX_20140505_C02_150331_1"
 7[3] "E10.5_XX_20140505_C03_150331_1" "E10.5_XX_20140505_C04_150331_2"
 8[5] "E10.5_XX_20140505_C06_150331_2" "E10.5_XX_20140505_C07_150331_3"
 9
10## 取下划线分隔的第一部分
11female_stages <- sapply(strsplit(colnames(females), "_"), `[`, 1)
12# 或者
13female_stages <- sapply(strsplit(colnames(females), "_"), 
14                        function(x)x[1])
15# 再或者
16female_stages <- stringr::str_split(colnames(females),'_', simplify = T)[,1]
17
18names(female_stages) <- colnames(females)
19> table(female_stages)
20female_stages
21E10.5 E11.5 E12.5 E13.5 E16.5    P6 
22   68   100   103    99    85   108

2.2 对基因操作=》基因过滤与统计

去掉在所有细胞都不表达的基因

1> (dim(females))
2[1] 21083   563
3> females <- females[rowSums(females)>0,]
4> (dim(females))
5[1] 16765   563

可以看到去掉了4000多个

计算各种统计指标

 1# 利用apply函数对每行（每个基因）进行统计
 2mean_per_gene <- apply(females, 1, mean, na.rm = TRUE) 
 3sd_per_gene <- apply(females, 1, sd, na.rm = TRUE) 
 4mad_per_gene <- apply(females, 1, mad, na.rm = TRUE) 
 5cv = sd_per_gene/mean_per_gene
 6library(matrixStats)
 7var_per_gene <- rowVars(as.matrix(females))
 8cv2=var_per_gene/mean_per_gene^2
 9# 存储统计结果
10cv_per_gene <- data.frame(mean = mean_per_gene,
11                          sd = sd_per_gene,
12                          mad=mad_per_gene,
13                          var=var_per_gene,
14                          cv=cv,
15                          cv2=cv2)
16rownames(cv_per_gene) <- rownames(females)
17head(cv_per_gene)
18# 根据表达量过滤统计结果
19cv_per_gene=cv_per_gene[cv_per_gene$mean>1,]
20# 简易的可视化
21with(cv_per_gene,plot(log10(mean),log10(cv2)))

CV值，它表示变异系数(coefficient of variation)。变异系数又称离散系数或相对偏差 ，我们肯定都知道标准偏差，也就是sd值，sd描述了数据值偏离算术平均值的程度。这个相对偏差CV描述的是标准偏差与平均值之比。

sd值，它和均值mean、方差var一样，都是对一维数据进行的分析，如果出现两组数据测量尺度差别太大或数据量纲存在差异的话，直接用标准差就不合适了
CV变异系数就可以解决这个问题，它利用原始数据标准差和原始数据平均值的比值来各自消除尺度与量纲的差异。

复杂一点的统计可视化：

其实就是求每列之间的相关性

1library(psych)
2pairs.panels(cv_per_gene, 
3             method = "pearson", # correlation method
4             hist.col = "#00AFBB",
5             density = TRUE,  # show density plots
6             ellipses = TRUE # show correlation ellipses
7)

可以得到不同统计指标的关系

再用作者包装的函数：getMostVarGenes()

1females_data <- getMostVarGenes(females, fitThr=2)
2> dim(females_data)
3[1] 822 563

这个函数也找了822个变化比较大的基因，用于下游分析，这其实也很像Seurat的FindVariableFeatures()做的事情

 1females_data <- log(females_data+1)
 2> females_data[1:4,1:4]
 3         E10.5_XX_20140505_C01_150331_1 E10.5_XX_20140505_C02_150331_1
 4Ngfr                                  0                              0
 5Slc22a18                              0                              0
 6Tspan32                               0                              0
 7Gmpr                                  0                              0
 8         E10.5_XX_20140505_C03_150331_1 E10.5_XX_20140505_C04_150331_2
 9Ngfr                          0.4204863                       3.619946
10Slc22a18                      0.0000000                       0.000000
11Tspan32                       0.0000000                       0.000000
12Gmpr                          0.0000000                       0.000000
13
14save(females_data,file = 'females_hvg_matrix.Rdata')

2.3 6个发育时期RtSNE分析

先是PCA

针对上面的822个HVGs进行操作

1female_sub_pca <- FactoMineR::PCA(
2  t(females_data), 
3  ncp = ncol(females_data), 
4  graph=FALSE
5)

然后挑选最显著的主成分，作为tSNE的输入

记得在Seurat中是使用ElbowPlot() 关注肘部的PC，这里不需要观察，直接返回最优解

1significant_pcs <- jackstraw::permutationPA(
2  female_sub_pca$ind$coord, 
3  B = 100, 
4  threshold = 0.05, 
5  verbose = TRUE, 
6  seed = NULL
7)$r
8> significant_pcs
9[1] 9

然后使用上面`jackstraw`挑出的显著主成分进行tSNE

 1# 6个时期给定6个颜色
 2female_stagePalette <-     c(
 3  "#2754b5", 
 4  "#8a00b0", 
 5  "#d20e0f", 
 6  "#f77f05", 
 7  "#f9db21",
 8  "#43f14b"
 9)
10female_t_sne <- run_plot_tSNE(
11  pca=female_sub_pca,
12  pc=significant_pcs,
13  iter=5000,
14  conditions=female_stages,
15  colours=female_stagePalette
16)

2.4 根据PCA结果进行层次聚类

采用的方法是：Hierarchical Clustering On Principle Components (HCPC)

 1# 使用9个显著主成分重新跑PCA
 2res.pca <- FactoMineR::PCA(
 3  t(females_data), 
 4  ncp = significant_pcs, 
 5  graph=FALSE
 6)
 7# 作者根据经验认为分成4群比较好解释，于是设置4
 8res.hcpc <- FactoMineR::HCPC(
 9  res.pca, 
10  graph = FALSE,
11  min=4
12)
13# 得到分群结果
14female_clustering <- res.hcpc$data.clust$clust
15> table(female_clustering)
16female_clustering
17  1   2   3   4 
18 90 240 190  43 
19# 重新命名
20female_clustering <- paste("C", female_clustering, sep="")
21names(female_clustering) <- rownames(res.hcpc$data.clust)
22# 将C1和C2调换位置
23female_clustering[female_clustering=="C1"] <- "C11"
24female_clustering[female_clustering=="C2"] <- "C22"
25female_clustering[female_clustering=="C22"] <- "C1"
26female_clustering[female_clustering=="C11"] <- "C2"
27> table(female_clustering)
28female_clustering
29 C1  C2  C3  C4 
30240  90 190  43 
31
32write.csv(female_clustering, file="female_clustering.csv")

还是基于之前tSNE坐标，对聚类得到的4个cluster可视化

 1# 为4种cluster设置颜色
 2female_clusterPalette <- c(
 3  "#560047", 
 4  "#a53bad", 
 5  "#eb6bac", 
 6  "#ffa8a0"
 7)
 8
 9> head(female_t_sne)
10                                  tSNE_1    tSNE_2  cond
11E10.5_XX_20140505_C01_150331_1 -2.714291 -24.47912 E10.5
12E10.5_XX_20140505_C02_150331_1 -1.580757 -26.45072 E10.5
13E10.5_XX_20140505_C03_150331_1 -1.577123 -25.36753 E10.5
14E10.5_XX_20140505_C04_150331_2 -6.677577 -20.00208 E10.5
15E10.5_XX_20140505_C06_150331_2  3.442235 -23.32570 E10.5
16E10.5_XX_20140505_C07_150331_3  3.793953 -23.33955 E10.5
17
18# 作者包装的函数
19female_t_sne_new_clusters <- plot_tSNE(
20  tsne=female_t_sne, 
21  conditions=female_clustering, 
22  colours= female_clusterPalette
23)
24ggsave('tSNE_cluster.pdf')