劉小澤寫于19.7.5-第二單元第六講：聚類算法之PCA與tSNE
筆記目的：根據(jù)生信技能樹的單細(xì)胞轉(zhuǎn)錄組課程探索smart-seq2技術(shù)相關(guān)的分析技術(shù)
課程鏈接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53
還是之前文章附件的圖片，其中b圖是選取兩個(gè)主成分做的PCA圖崖疤，c圖是tSNE圖：
幾個(gè)常用函數(shù)的轉(zhuǎn)置t(transpose)秘车，傻傻分不清？： 計(jì)算距離介紹過
dist()函數(shù)戳晌，它是按行為操作對象鲫尊，而聚類是要對樣本聚類，因此要先將我們平時(shí)見到的表達(dá)矩陣(行為基因沦偎，列為樣本)轉(zhuǎn)置疫向；同樣PCA也是對行/樣本進(jìn)行操作，也是需要先轉(zhuǎn)置豪嚎；另外歸一化的scale()函數(shù)雖然是對列進(jìn)行操作搔驼，但它的對象是基因，因此也需要轉(zhuǎn)置
關(guān)于PCA的學(xué)習(xí)侈询，之前寫過：
- StatQuest-PCA學(xué)習(xí)：http://www.reibang.com/p/b83ac8f7f5a7
- StatQuest--在R中拆解PCA：http://www.reibang.com/p/8a74508c3737
先構(gòu)建一個(gè)非常隨機(jī)的測試數(shù)據(jù)
# 設(shè)置隨機(jī)種子舌涨，可以重復(fù)別人使用的隨機(jī)數(shù)
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)
沒有體現(xiàn)任何的基因差異或者樣本聚類（熱圖中的聚類是自然層次聚類），可以看到樣本名都是無規(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)
可以看到每組的20個(gè)細(xì)胞都分不開温技，但每組具體有哪些樣本還是看不出來，因此這里為每組加上顏色來表示
# 先在原來數(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ù)Rtsne()
set.seed(42)
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中存儲了畫圖坐標(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ī)律的測試數(shù)據(jù)
ng=500
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)能看出來差異了，再看看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
對數(shù)據(jù)進(jìn)行PCA
# 操作前先備份
dat_back=dat
# 先對表達(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è)批次之間分不開慕淡，說明沒有批次效應(yīng)
