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

生信技能树jimmy

发布于 2020-03-30 14:40:18

1.1K0

文章被收录于专栏：单细胞天地单细胞天地

目的是得到这个图

将会用到来自作者的包装好的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

load(file="female_rpkm.Robj")

## 去掉重复细胞
#（例如：同一个细胞建库两次，这里作者用“rep”进行了标记） 
grep("rep",colnames(female_rpkm))
colnames(female_rpkm)[256:257]
female_rpkm <- female_rpkm[,!colnames(female_rpkm) %in% grep("rep",colnames(female_rpkm), value=TRUE)]

## 只保留编码基因(去掉类似：X5430419D17Rik、BC003331等)
prot_coding_genes <- read.csv(file="prot_coding.csv", row.names=1)
females <- female_rpkm[rownames(female_rpkm) %in% as.vector(prot_coding_genes$x),]
save(females,file = 'female_rpkm.Rdata')

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

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

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

load('../female_rpkm.Rdata')
> dim(females)
[1] 21083   563

> head(colnames(females))
[1] "E10.5_XX_20140505_C01_150331_1" "E10.5_XX_20140505_C02_150331_1"
[3] "E10.5_XX_20140505_C03_150331_1" "E10.5_XX_20140505_C04_150331_2"
[5] "E10.5_XX_20140505_C06_150331_2" "E10.5_XX_20140505_C07_150331_3"

## 取下划线分隔的第一部分
female_stages <- sapply(strsplit(colnames(females), "_"), `[`, 1)
# 或者
female_stages <- sapply(strsplit(colnames(females), "_"), 
                        function(x)x[1])
# 再或者
female_stages <- stringr::str_split(colnames(females),'_', simplify = T)[,1]

names(female_stages) <- colnames(females)
> table(female_stages)
female_stages
E10.5 E11.5 E12.5 E13.5 E16.5    P6 
   68   100   103    99    85   108

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

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

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

可以看到去掉了4000多个

计算各种统计指标

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

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

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

复杂一点的统计可视化：

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

library(psych)
pairs.panels(cv_per_gene, 
             method = "pearson", # correlation method
             hist.col = "#00AFBB",
             density = TRUE,  # show density plots
             ellipses = TRUE # show correlation ellipses
)

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

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

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

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

females_data <- log(females_data+1)
> females_data[1:4,1:4]
         E10.5_XX_20140505_C01_150331_1 E10.5_XX_20140505_C02_150331_1
Ngfr                                  0                              0
Slc22a18                              0                              0
Tspan32                               0                              0
Gmpr                                  0                              0
         E10.5_XX_20140505_C03_150331_1 E10.5_XX_20140505_C04_150331_2
Ngfr                          0.4204863                       3.619946
Slc22a18                      0.0000000                       0.000000
Tspan32                       0.0000000                       0.000000
Gmpr                          0.0000000                       0.000000

save(females_data,file = 'females_hvg_matrix.Rdata')

2.3 6个发育时期RtSNE分析

先是PCA

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

female_sub_pca <- FactoMineR::PCA(
  t(females_data), 
  ncp = ncol(females_data), 
  graph=FALSE
)

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

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

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

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

# 6个时期给定6个颜色
female_stagePalette <-     c(
  "#2754b5", 
  "#8a00b0", 
  "#d20e0f", 
  "#f77f05", 
  "#f9db21",
  "#43f14b"
)
female_t_sne <- run_plot_tSNE(
  pca=female_sub_pca,
  pc=significant_pcs,
  iter=5000,
  conditions=female_stages,
  colours=female_stagePalette
)

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

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

# 使用9个显著主成分重新跑PCA
res.pca <- FactoMineR::PCA(
  t(females_data), 
  ncp = significant_pcs, 
  graph=FALSE
)
# 作者根据经验认为分成4群比较好解释，于是设置4
res.hcpc <- FactoMineR::HCPC(
  res.pca, 
  graph = FALSE,
  min=4
)
# 得到分群结果
female_clustering <- res.hcpc$data.clust$clust
> table(female_clustering)
female_clustering
  1   2   3   4 
 90 240 190  43 
# 重新命名
female_clustering <- paste("C", female_clustering, sep="")
names(female_clustering) <- rownames(res.hcpc$data.clust)
# 将C1和C2调换位置
female_clustering[female_clustering=="C1"] <- "C11"
female_clustering[female_clustering=="C2"] <- "C22"
female_clustering[female_clustering=="C22"] <- "C1"
female_clustering[female_clustering=="C11"] <- "C2"
> table(female_clustering)
female_clustering
 C1  C2  C3  C4 
240  90 190  43 

write.csv(female_clustering, file="female_clustering.csv")

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

# 为4种cluster设置颜色
female_clusterPalette <- c(
  "#560047", 
  "#a53bad", 
  "#eb6bac", 
  "#ffa8a0"
)

> head(female_t_sne)
                                  tSNE_1    tSNE_2  cond
E10.5_XX_20140505_C01_150331_1 -2.714291 -24.47912 E10.5
E10.5_XX_20140505_C02_150331_1 -1.580757 -26.45072 E10.5
E10.5_XX_20140505_C03_150331_1 -1.577123 -25.36753 E10.5
E10.5_XX_20140505_C04_150331_2 -6.677577 -20.00208 E10.5
E10.5_XX_20140505_C06_150331_2  3.442235 -23.32570 E10.5
E10.5_XX_20140505_C07_150331_3  3.793953 -23.33955 E10.5

# 作者包装的函数
female_t_sne_new_clusters <- plot_tSNE(
  tsne=female_t_sne, 
  conditions=female_clustering, 
  colours= female_clusterPalette
)
ggsave('tSNE_cluster.pdf')