1 使用簡單的 prcomp 函數(shù)來了解 PCA分析

ng=500 
    nc=20
a1=rnorm(ng*nc);dim(a1)=c(ng,nc) #創(chuàng)建正態(tài)分布隨機矩陣500行，20列
 #dim()檢索或設(shè)置對象的維度
a2=rnorm(ng*nc);dim(a2)=c(ng,nc) #因為是隨機創(chuàng)建，這兩個矩陣不一樣
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)
    dist(a3) #dist()函數(shù)按行分析，按列的話需要轉(zhuǎn)置
    a3=t(a3);dim(a3) ## PCA分析 
pca_dat <- prcomp(a3, scale. = TRUE) #prcomp()主成分分析
p=autoplot(pca_dat) + theme_classic() + ggtitle('PCA plot')
print(p)

屏幕快照 2020-07-05 下午7.29.44.png

可以看到細胞無法被區(qū)分開來。

2 t-SNE分析

set.seed(42)
tsne_out <- Rtsne(a3,pca=FALSE,perplexity=10,theta=0.0) # Run TSNE
    tsnes=tsne_out$Y
    colnames(tsnes) <- c("tSNE1", "tSNE2") #添加列名
ggplot(tsnes, aes(x = tSNE1, y = tSNE2))+ geom_point()

屏幕快照 2020-07-05 下午7.33.23.png

還是無法區(qū)分

3 同樣的正態(tài)分布隨機表達矩陣，但是其中部分細胞+3，PCA可以區(qū)分開來。

 ng=500
    nc=20
    a1=rnorm(ng*nc);dim(a1)=c(ng,nc)
    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)
    a3=t(a3);dim(a3) ## PCA分析 
    
    pca_dat <- prcomp(a3, scale. = TRUE)
    p=autoplot(pca_dat) + theme_classic() + ggtitle('PCA plot')
    print(p)

屏幕快照 2020-07-05 下午7.37.40.png

4 同樣的正態(tài)分布隨機表達矩陣，但是其中部分細胞+3，tSNE可以區(qū)分開來。

set.seed(42)
    tsne_out <- Rtsne(a3,pca=FALSE,perplexity=10,theta=0.0) # Run TSNE
    tsnes=tsne_out$Y
    colnames(tsnes) <- c("tSNE1", "tSNE2")
    ggplot(tsnes, aes(x = tSNE1, y = tSNE2))+ geom_point()

屏幕快照 2020-07-05 下午7.38.46.png

5 真實數(shù)據(jù)的PCA操作

rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
load(file = '../input_rpkm.Rdata')
dat[1:4,1:4]
head(metadata) 
plate=metadata$plate
## 下面是畫PCA的必須操作，需要看不同做PCA的包的說明書。
dat_back=dat

dat=dat_back
dat=t(dat)
dat=as.data.frame(dat)
dat=cbind(dat,plate ) #cbind根據(jù)列進行合并，即疊加所有列 #矩陣添加批次信息
dat[1:4,1:4]
table(dat$plate)
library("FactoMineR")
library("factoextra") 
# The variable plate (index = ) is removed
# before PCA analysis
dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)
fviz_pca_ind(dat.pca,#repel =T,
             geom.ind = "point", # show points only (nbut not "text")
             col.ind = dat$plate, # color by groups
             #palette = c("#00AFBB", "#E7B800"),
             addEllipses = TRUE, # Concentration ellipses
             legend.title = "Groups") 
ggsave('all_cells_PCA_by_plate.png') 

屏幕快照 2020-07-05 下午7.55.52.png