根据表达矩阵进行分群-2
上面我们使用了RPKM矩阵，下面的Seurat将会使用原始表达矩阵。当然也是推荐使用原始矩阵进行分析的
3.1 下载原始表达矩阵
链接在：https://raw.githubusercontent.com/IStevant/XX-XY-mouse-gonad-scRNA-seq/master/data/female_count.Robj
load(file="../female_count.Robj")
load('../female_rpkm.Rdata')
# 直接对细胞和基因过滤
female_count <- female_count[rownames(female_count) %in% rownames(females),!colnames(female_count) %in% grep("rep",colnames(female_count), value=TRUE)]
> female_count[1:3,1:3]
E10.5_XX_20140505_C01_150331_1 E10.5_XX_20140505_C02_150331_1
eGFP 19582 526
Gnai3 2218 122
Pbsn 0 0
E10.5_XX_20140505_C03_150331_1
eGFP 4786
Gnai3 4
Pbsn 0
save(female_count,file = '../female_count.Rdata')
3.2 对细胞操作=》细胞发育时期的获取
load('../female_count.Rdata')
female_stages <- sapply(strsplit(colnames(female_count), "_"), `[`, 1)
names(female_stages) <- colnames(female_count)
> table(female_stages)
female_stages
E10.5 E11.5 E12.5 E13.5 E16.5 P6
68 100 103 99 85 108
3.3 使用Seurat V3
构建对象
sce_female <- CreateSeuratObject(counts = female_count,
project = "sce_female",
min.cells = 1, min.features = 0)
> sce_female
An object of class Seurat
16765 features across 563 samples within 1 assay
Active assay: RNA (16765 features)
添加样本注释信息
sce_female <- AddMetaData(object = sce_female,
metadata = apply(female_count, 2, sum),
col.name = 'nUMI_raw')
sce_female <- AddMetaData(object = sce_female,
metadata = female_stages,
col.name = 'female_stages')
数据归一化
sce_female <- NormalizeData(sce_female)
sce_female[["RNA"]]@data[1:3,1:3]
找差异基因HVGs
sce_female <- FindVariableFeatures(sce_female,
selection.method = "vst",
nfeatures = 2000)
# HVGs可视化
VariableFeaturePlot(sce_female)
seurat3_HVGs <- VariableFeatures(sce_female)
# 检查与之前得到的HVGs重合度
load('females_hvg_matrix.Rdata')
load('seurat3_HVGs.Rdata')
length(intersect(rownames(females_data),seurat3_HVGs))
# 结果和之前822个HVGs有434个重合
数据标准化
# 默认只对FindVariableFeatures得到的HVGs进行操作
sce_female <- ScaleData(object = sce_female,
vars.to.regress = c('nUMI_raw'),
model.use = 'linear',
use.umi = FALSE)
PCA降维
sce_female <- RunPCA(sce_female,
features = VariableFeatures(object = sce_female))
降维后聚类
# 这里可以多选一些PCs
sce_female <- FindNeighbors(sce_female, dims = 1:20)
sce_female <- FindClusters(sce_female, resolution = 0.3)
进行tSNE
ElbowPlot(sce_female)
sce_female_tsne <- RunTSNE(sce_female, dims = 1:9)
tSNE结果可视化
# 6个发育时间
DimPlot(object = sce_female_tsne, reduction = "tsne",
group.by = 'female_stages')
# 4个cluster
DimPlot(sce_female_tsne, reduction = "tsne")
比较两次的聚类结果
cluster1 <- read.csv('female_clustering.csv')
cluster2 <- as.data.frame(Idents(sce_female_tsne))
# 把它们放在一起比较，前提条件是它们的行名相同
> identical(cluster1[,1],rownames(cluster2))
[1] TRUE
> table(cluster1[,2],cluster2[,1])
0 1 2 3
C1 224 3 13 0
C2 6 0 84 0
C3 12 177 0 1
C4 0 0 0 43
这也说明了，不同方法虽然选择的HVGs数量不同，也不完全一样，聚类的参数也不同，但最后真正的生物学意义是不会去掉的。只能说，最后选多少群是根据分析的人根据自己的理解去解释，只要参数变化，就会有各种不同的结果
4 使用更简单的函数去分群
rm(list = ls())
options(warn=-1)
options(stringsAsFactors = F)
load('../female_rpkm.Rdata')
# 根据分群获得颜色
cluster <- read.csv('female_clustering.csv')
color <- rainbow(4)[as.factor(cluster[,2])]
> table(color)
color
#00FFFFFF #8000FFFF #80FF00FF #FF0000FF
190 43 90 240
# 取前1000个sd最大的基因作为HVGs
choosed_count <- females
# 表达矩阵过滤
choosed_count <- choosed_count[apply(choosed_count, 1, sd)>0,]
choosed_count <- choosed_count[names(head(sort(apply(choosed_count, 1, sd),decreasing = T),1000)),]
进行PCA分析
pca_out <- prcomp(t(choosed_count),scale. = T)
> pca_out$x[1:3,1:3]
PC1 PC2 PC3
E10.5_XX_20140505_C01_150331_1 13.21660 -4.1600782 1.5287334
E10.5_XX_20140505_C02_150331_1 13.73109 -0.2848806 -0.8443587
E10.5_XX_20140505_C03_150331_1 10.89558 -0.2720221 -3.3839651
library(ggfortify)
autoplot(pca_out, col=color) +theme_classic()+ggtitle('PCA plot')
进行tSNE
library(Rtsne)
# 依旧选前9个
tsne_out <- Rtsne(pca_out$x[,1:9], perplexity = 10,
pca = F, max_iter = 2000,
verbose = T)
tsnes_cord <- tsne_out$Y
colnames(tsnes_cord) <- c('tSNE1','tSNE2')
ggplot(tsnes_cord, aes(x=tSNE1, y = tSNE2)) + geom_point(col=color) + theme_classic()+ggtitle('tSNE plot')
除了之前的HCPC和seurat分群，还可以利用DBSCAN、kmeans分群
# 这个运行会非常慢！
if(T){
library(Rtsne)
N_tsne <- 50
tsne_out <- list(length = N_tsne)
KL <- vector(length = N_tsne)
set.seed(1234)
for(k in 1:N_tsne)
{
tsne_out[[k]]<-Rtsne(t(log2(females+1)),initial_dims=30,verbose=FALSE,check_duplicates=FALSE,
perplexity=27, dims=2,max_iter=5000)
KL[k]<-tail(tsne_out[[k]]$itercosts,1)
print(paste0("FINISHED ",k," TSNE ITERATION"))
}
names(KL) <- c(1:N_tsne)
opt_tsne <- tsne_out[[as.numeric(names(KL)[KL==min(KL)])]]$Y
}
# DBSCAN结果
library(dbscan)
plot(opt_tsne, col=dbscan(opt_tsne,eps=3.1)$cluster,
pch=19, xlab="tSNE dim 1", ylab="tSNE dim 2")
# kmeans结果
plot(opt_tsne, col=kmeans(opt_tsne,centers = 4)$clust,
pch=19, xlab="tSNE dim 1", ylab="tSNE dim 2")
比较它们的差异
# 其中kmeans是4群
> table(kmeans(opt_tsne,centers = 4)$clust,dbscan(opt_tsne,eps=3.5)$cluster)
0 1 2 3 4
1 2 0 0 206 0
2 1 106 0 0 0
3 0 93 10 0 0
4 1 138 0 1 5本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。
原始发表：2019-11-15，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
目录