前 · 言 第二单元第六讲:聚类算法之PCA与tSNE 还是之前文章附件的图片,其中b图是选取两个主成分做的PCA图,c图是tSNE图:
关于PCA的学习,之前写过:
先构建一个非常随机的测试数据 set.seed(123456789) library(pheatmap) library(Rtsne) library(ggfortify) library(mvtnorm) # 设置两个正态分布的随机矩阵(500*20) ng=500 nc=20 a1=rnorm(ng*nc);dim(a1)=c(ng,nc) a2=rnorm(ng*nc);dim(a2)=c(ng,nc) a3=cbind(a1,a2) > dim(a3) [1] 500 40 # 添加列名 colnames(a3)=c(paste0('cell_01_',1:nc), paste0('cell_02_',1:nc)) # 添加行名 rownames(a3)=paste('gene_',1:ng,sep = '') # 先做个热图 pheatmap(a3) 没有体现任何的基因差异或者样本聚类(热图中的聚类是自然层次聚类),可以看到样本名都是无规律的交叉显示 如果做PCA呢?# 先转置一下,让行为样本> a3=t(a3);dim(a3) [1] 40 500 # prcomp()主成分分析 pca_dat <- prcomp(a3, scale. = TRUE) p=autoplot(pca_dat) + theme_classic() + ggtitle('PCA plot') print(p) # 先在原来数据的基础上添加样本分组信息(别忘了a3是一个矩阵,先转换成数据框) df=cbind(as.data.frame(a3),group=c(rep('b1',20),rep('b2',20))) autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = 'group')+theme_bw() 另外看下tsne利用了一个核心函数 tsne_out <- Rtsne(a3,pca=FALSE,perplexity=10,theta=0.0) # 结果得到一个列表 > str(tsne_out) List of 14 $ N : int 40 $ Y : num [1:40, 1:2] -36.7 -28 -168 -33.4 22.4 ... $ costs : num [1:40] 0.00348 -0.00252 0.01496 0.01646 0.00951 ... # 其中在Y中存储了画图坐标 > head(tsne_out$Y,3) [,1] [,2] [1,] -36.72621 -78.03709 [2,] -28.00151 33.30229 [3,] -167.98560 -80.26850 tsnes=tsne_out$Y colnames(tsnes) <- c("tSNE1", "tSNE2") #为坐标添加列名 # 基础作图代码 ggplot(tsnes, aes(x = tSNE1, y = tSNE2))+ geom_point() # 在此基础上添加颜色分组信息,首先还是将tsnes这个矩阵变成数据框,然后增加一列group信息,最后映射在geom_point中 tsnes=as.data.frame(tsnes) group=c(rep('b1',20),rep('b2',20)) tsnes$group=group ggplot(tsnes, aes(x = tSNE1, y = tSNE2))+ geom_point(aes(col=group)) 构建一个有规律的测试数据 nc=20 a1=rnorm(ng*nc);dim(a1)=c(ng,nc) # 和之前的区别就在a2这里,都加了3 a2=rnorm(ng*nc)+3;dim(a2)=c(ng,nc) a3=cbind(a1,a2) colnames(a3)=c(paste0('cell_01_',1:nc),paste0('cell_02_',1:nc)) rownames(a3)=paste('gene_',1:ng,sep = '') pheatmap(a3) 热图已经能看出来差异了,再看看PCA a3=t(a3);dim(a3)df=cbind(as.data.frame(a3),group=c(rep('b1',20),rep('b2',20))) autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = 'group')+theme_bw() tsne也是如此 set.seed(42)tsne_out <- Rtsne(a3,pca=FALSE,perplexity=10,theta=0.0) tsnes=tsne_out$Y colnames(tsnes) <- c("tSNE1", "tSNE2") tsnes=as.data.frame(tsnes) group=c(rep('b1',20),rep('b2',20)) tsnes$group=group ggplot(tsnes, aes(x = tSNE1, y = tSNE2))+ geom_point(aes(col=group)) 真实数据演练 载入RPKM数据rm(list = ls())options(stringsAsFactors = F) load(file = '../input_rpkm.Rdata') # 表达量信息 > dat[1:2,1:3] SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5 0610007P14Rik 0 0 74.95064 0610009B22Rik 0 0 0.00000 # 样本属性 > head(metadata,3) g plate n_g all SS2_15_0048_A3 1 0048 3065 all SS2_15_0048_A6 2 0048 3036 all SS2_15_0048_A5 1 0048 3742 all #所有数据的聚类分组信息 group_list=metadata$g #批次信息 plate=metadata$plate > table(plate) plate 0048 0049 384 384 对数据进行PCA# 操作前先备份dat_back=dat # 先对表达矩阵进行转置,然后转换成数据框,就可以添加批次信息了 dat=dat_back dat=t(dat) dat=as.data.frame(dat) dat=cbind(dat,plate ) > dim(dat_back) [1] 12689 768 > dim(dat) [1] 768 12690 library("FactoMineR") library("factoextra") dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE) fviz_pca_ind(dat.pca, # repel =T, geom.ind = "point", # 只显示点,不显示文字 col.ind = dat$plate, # 按分组上色 #palette = c("#00AFBB", "#E7B800"), addEllipses = TRUE, # 添加晕环 legend.title = "Groups" 可以看到两个批次之间分不开,说明没有批次效应 |
|