基于基因集的样品队列分组之PCA

load(file = 'step1-output.Rdata')
cg=c('CD3D','CD3G CD247','IFNG','IL2RG','IRF1','IRF4','LCK','OAS2,STAT1')
cg
cg=cg[cg %in% rownames(dat)]
library("FactoMineR") #画主成分分析图需要加载这两个包
library("factoextra")  
dat.pca <- PCA(as.data.frame(t(dat[cg,]))  )
fviz_pca_ind(dat.pca,
             geom.ind = "point", # show points only (nbut not "text")
             # col.ind =  group_list, # color by groups 
             addEllipses = T,
             legend.title = "Groups"
)
dat.pca$ind$coord

> head(dat.pca$ind$coord[,1:2])
                Dim.1      Dim.2
GSM1419942 -1.5307665  0.1820390
GSM1419943  0.2224343  0.9908672
GSM1419944 -2.1294537  0.3075360
GSM1419945  0.0745490  0.2147092
GSM1419946 -4.8646160 -0.2016098
GSM1419947  0.9149908  0.2421981

group_list=ifelse(dat.pca$ind$coord[,1] < 0 ,'neg','pos')
table(group_list) 
pca_gl = group_list
# 其中 hclust_gl 来自于前面的教程哦
table(pca_gl,hclust_gl)

> table(pca_gl,hclust_gl)
      hclust_gl
pca_gl high low
   neg    1  59
   pos   37  10

library(patchwork)
fviz_pca_ind(dat.pca,
             geom.ind = "point", # show points only (nbut not "text")
             col.ind =  pca_gl, # color by groups 
             addEllipses = T,
             legend.title = "Groups"
)+ 
  fviz_pca_ind(dat.pca,
               geom.ind = "point", # show points only (nbut not "text")
               col.ind =  hclust_gl, # color by groups 
               addEllipses = T,
               legend.title = "Groups"
  )

rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)

library(AnnoProbe)
library(GEOquery) 
gset <- geoChina("GSE58812")
gset
gset[[1]]
a=gset[[1]] #
dat=exprs(a) #a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
dim(dat)#看一下dat这个矩阵的维度
#   GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array
dat[1:4,1:4] #查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
boxplot(dat[,1:4],las=2)  
dat=log2(dat)
boxplot(dat[,1:4],las=2)  
library(limma)
dat=normalizeBetweenArrays(dat)
boxplot(dat[,1:4],las=2)  
pd=pData(a) 
#通过查看说明书知道取对象a里的临床信息用pData
## 挑选一些感兴趣的临床表型。
library(stringr) 
table(pd$`meta:ch1`)
group_list= ifelse(pd$`meta:ch1` == '0' ,'stable','meta')
table(group_list)
 
dat[1:4,1:4]  
dim(dat)
a  
ids=idmap('GPL570','soft')
head(ids)
ids=ids[ids$symbol != '',]
dat=dat[rownames(dat) %in% ids$ID,]
ids=ids[match(rownames(dat),ids$ID),]
head(ids) 
colnames(ids)=c('probe_id','symbol')  
ids$probe_id=as.character(ids$probe_id)
rownames(dat)=ids$probe_id
dat[1:4,1:4] 

ids=ids[ids$probe_id %in%  rownames(dat),]
dat[1:4,1:4]   
dat=dat[ids$probe_id,] 

ids$median=apply(dat,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids
ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s
dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
dat[1:4,1:4]  #保留每个基因ID第一次出现的信息

dat['ACTB',]
dat['GAPDH',]
save(dat,group_list,
     file = 'step1-output.Rdata')

We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.