根据表达矩阵进行分群-1 / 四六文摘

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

课程链接在：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')

“

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

前言