text: ## 
  ### ---------------
  ###
  ### Create: Jianming Zeng
  ### Date: 2018-12-20 15:43:52
  ### Email: jmzeng1314@163.com
  ### Blog: http://www.bio-info-trainee.com/
  ### Forum:  http://www.biotrainee.com/thread-1376-1-1.html
  ### CAFS/SUSTC/Eli Lilly/University of Macau
  ### Update Log: 2018-12-20  First version
  ### Update Log: 2019-09-10 基于R version 3.5.1 (2018-07-02)
  ### ---------------
  
  
  rm(list = ls())  ## 魔幻操作，一键清空~
  options(stringsAsFactors = F)
  load(file = 'step1-output.Rdata')
  table(group_list)
  # 每次都要检测数据
  dat[1:4,1:4]
  ## 下面是画PCA的必须操作，需要看说明书。
  dat=t(dat)#画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换
  dat=as.data.frame(dat)#将matrix转换为data.frame
  dat=cbind(dat,group_list) #cbind横向追加，即将分组信息追加到最后一列
  library("FactoMineR")#画主成分分析图需要加载这两个包
  library("factoextra") 
  # The variable group_list (index = 54676) is removed
  # before PCA analysis
  dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)#现在dat最后一列是group_list，需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
  fviz_pca_ind(dat.pca,
               geom.ind = "point", # show points only (nbut not "text")
               col.ind = dat$group_list, # color by groups
               # palette = c("#00AFBB", "#E7B800"),
               addEllipses = TRUE, # Concentration ellipses
               legend.title = "Groups"
  )
  ggsave('all_samples_PCA.png')
  
  
  
  
  rm(list = ls())  ## 魔幻操作，一键清空~
  load(file = 'step1-output.Rdata') #此步为一个小插曲，即计算一下从第一行开是计算每一行的sd值，知道最后一行所需要的时间
  dat[1:4,1:4] 
  
  cg=names(tail(sort(apply(dat,1,sd)),1000))#apply按行（'1'是按行取，'2'是按列取）取每一行的方差，从小到大排序，取最大的1000个
  library(pheatmap)
  pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #对那些提取出来的1000个基因所在的每一行取出，组合起来为一个新的表达矩阵
  n=t(scale(t(dat[cg,]))) # 'scale'可以对log-ratio数值进行归一化
  n[n>2]=2 
  n[n< -2]= -2
  n[1:4,1:4]
  pheatmap(n,show_colnames =F,show_rownames = F)
  ac=data.frame(g=group_list)
  rownames(ac)=colnames(n) #把ac的行名给到n的列名，即对每一个探针标记上分组信息（是'noTNBC'还是'TNBC'）
  ## 可以看到TNBC具有一定的异质性，拿它来区分乳腺癌亚型指导临床治疗还是略显粗糙。
  pheatmap(n,show_colnames =F,show_rownames = F,
           annotation_col=ac,filename = 'heatmap_top1000_sd.png')
  
  
  
  
  