  前 · 言 第二單元第六講：聚類算法之PCA與tSNE 還是之前文章附件的圖片，其中b圖是選取兩個(gè)主成分做的PCA圖，c圖是tSNE圖：
關(guān)于PCA的學(xué)習(xí)，之前寫(xiě)過(guò)：
先構(gòu)建一個(gè)非常隨機(jī)的測(cè)試數(shù)據(jù) set.seed(123456789) library(pheatmap) library(Rtsne) library(ggfortify) library(mvtnorm) # 設(shè)置兩個(gè)正態(tài)分布的隨機(jī)矩陣(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 = '') # 先做個(gè)熱圖 pheatmap(a3) 沒(méi)有體現(xiàn)任何的基因差異或者樣本聚類（熱圖中的聚類是自然層次聚類），可以看到樣本名都是無(wú)規(guī)律的交叉顯示 如果做PCA呢？# 先轉(zhuǎn)置一下，讓行為樣本> 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) # 先在原來(lái)數(shù)據(jù)的基礎(chǔ)上添加樣本分組信息(別忘了a3是一個(gè)矩陣，先轉(zhuǎn)換成數(shù)據(jù)框) 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利用了一個(gè)核心函數(shù) tsne_out <- Rtsne(a3,pca=FALSE,perplexity=10,theta=0.0) # 結(jié)果得到一個(gè)列表 > 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中存儲(chǔ)了畫(huà)圖坐標(biāo) > 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") #為坐標(biāo)添加列名 # 基礎(chǔ)作圖代碼 ggplot(tsnes, aes(x = tSNE1, y = tSNE2))+ geom_point() # 在此基礎(chǔ)上添加顏色分組信息，首先還是將tsnes這個(gè)矩陣變成數(shù)據(jù)框，然后增加一列g(shù)roup信息，最后映射在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)) 構(gòu)建一個(gè)有規(guī)律的測(cè)試數(shù)據(jù) nc=20 a1=rnorm(ng*nc);dim(a1)=c(ng,nc) # 和之前的區(qū)別就在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) 熱圖已經(jīng)能看出來(lái)差異了，再看看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)) 真實(shí)數(shù)據(jù)演練 載入RPKM數(shù)據(jù)rm(list = ls())options(stringsAsFactors = F) load(file = '../input_rpkm.Rdata') # 表達(dá)量信息 > 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 #所有數(shù)據(jù)的聚類分組信息 group_list=metadata$g #批次信息 plate=metadata$plate > table(plate) plate 0048 0049 384 384 對(duì)數(shù)據(jù)進(jìn)行PCA# 操作前先備份dat_back=dat # 先對(duì)表達(dá)矩陣進(jìn)行轉(zhuǎn)置，然后轉(zhuǎn)換成數(shù)據(jù)框，就可以添加批次信息了 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", # 只顯示點(diǎn)，不顯示文字 col.ind = dat$plate, # 按分組上色 #palette = c("#00AFBB", "#E7B800"), addEllipses = TRUE, # 添加暈環(huán) legend.title = "Groups"
可以看到兩個(gè)批次之間分不開(kāi)，說(shuō)明沒(méi)有批次效應(yīng) |
|
|
