(四) 03_PCA降维和heatmap热图

> rm(list = ls())  #清除所有变量
> load(file = "step2output.Rdata")#加载之前的数据
#输入数据：exp表达数据和group_list因子化分组
#Principal Component Analysis
#http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials

PCA降维——第一步，将exp数据转置，降维为甚需转置呢？

> exp[1:4,1:4]#探针id和样本名
        GSM1052615 GSM1052616 GSM1052617 GSM1052618
7892501    7.24559    6.80686    7.73301    6.18961
7892502    6.82711    6.70157    7.02471    6.20493
7892503    4.39977    4.50781    4.88250    4.36295
7892504    9.48025    9.67952    9.63074    9.69200
> t(exp[1:4,1:4])#转置
           7892501 7892502 7892503 7892504
GSM1052615 7.24559 6.82711 4.39977 9.48025
GSM1052616 6.80686 6.70157 4.50781 9.67952
GSM1052617 7.73301 7.02471 4.88250 9.63074
GSM1052618 6.18961 6.20493 4.36295 9.69200
#View(t(exp))

第二步，PCA降维

> {
  dat=as.data.frame(t(exp))#exp表达矩阵数据框
  library(FactoMineR)#FactoMineR（负责分析标准化数据）
  library(factoextra)#factoextra（负责可视化)
  dat.pca <- PCA(dat, graph = FALSE) 
  #dat从exp得来;graph逻辑值，TRUE的话自动画图
 #FactoMineR中的PCA()自动帮助dat进行标准化
  #下面用factoextra包里面的函数进行可视化操作
  fviz_pca_ind(dat.pca,#fviz_pca_ind对单个变量进行画图
               geom.ind = "point", # 只显示点
               col.ind = group_list, # color by groups用颜色指示值
               addEllipses = TRUE, # Concentration ellipses
               legend.title = "Groups"
  )
  ggsave('all_samples_PCA.png')#保存在工作目录下
> }

all_samples_PCA.png

dim1和dim2分别代表什么？

• 热图——第一步，数据准备

数据exp.png

> cg=names(tail(sort(apply(exp,1,sd)),1000))#sd是标准差
#apply(X, MARGIN, FUN, ...)#X 阵列
#MARGIN 1表示矩阵行，2表示矩阵列，也可以是c(1,2)
#对数组或者矩阵的一个维度使用函数生成值得列表或者数组、向量。
#sort（）默认升序
#view(cg)
#head(cg)
#head(exp)
> n=exp[cg,]#筛选行,exp中取1000子集
> head(n)
        GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619 GSM1052620
7893963    5.79288    5.27317    6.08883    4.52216    5.10400    4.48187
7895624    7.71782    8.01318    8.40404    7.21848    8.01405    6.59685
8109830   11.47330   11.30540   11.42210   10.37890   10.24600   10.05490
7953844   10.10830    9.97290   10.09840    8.98978    8.82521    8.81154
8156126    9.18541    8.87702    9.23965    7.97783    7.97500    7.84363
8152119    7.09276    6.98717    7.13838    8.26911    8.25215    8.26818

第二步，将分组数据转换成数据框，如下图所示

> annotation_col=data.frame(group=group_list)#转成数据框方便增加,列名
            group
GSM1052615 control
GSM1052616 control
GSM1052617 control
GSM1052618   treat
GSM1052619   treat
GSM1052620   treat

分组数据框.png

第三步，画热图

> rownames(annotation_col)=colnames(n) #行名
[1] "GSM1052615" "GSM1052616" "GSM1052617" "GSM1052618" "GSM1052619" "GSM1052620"
#View(annotation_col)
> library(pheatmap)
> pheatmap(n,
         show_colnames =F,#是否显示列名
         show_rownames = F,
         annotation_col=annotation_col,#增加分组的栏
         scale = "row")#横的比较
> dev.off()

Rplot.png