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个：：tail之后顺序还是从小到大的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[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的列名，即对每一个探针标记上分组信息
> ac g 1 Control 2 Control 3 Control 4 Vemurafenib 5 Vemurafenib 6 Vemurafenib > rownames(ac) = colnames(n) > ac g GSM1052615 Control GSM1052616 Control GSM1052617 Control GSM1052618 Vemurafenib GSM1052619 Vemurafenib GSM1052620 Vemurafenib
pheatmap(n,show_colnames =F,show_rownames = F, annotation_col=ac,filename = ‘heatmap_top1000_sd.png‘)
> pheatmap(n,show_colnames = F,show_rownames = F, + annotation_col = ac,file = "heatmap_top1000_sd.png") ##所以可以直接在代码中保存文件， > pheatmap(n,show_colnames = F,show_rownames = F, + annotation_col = ac) ##然后我想直接看看这张图长啥样子，输入这行命令之后没啥反应，plot窗口没有图 > pheatmap(n,show_colnames = F,show_rownames = F) > dev.off() ##之前不记得看什么的时候记过，画图没反应就用这个函数dev.off() null device 1 > pheatmap(n,show_colnames = F,show_rownames = F) > pheatmap(n,show_colnames = F,show_rownames = F,annotation = ac) ##不过没太懂annotation这个参数是啥意思
> pheatmap(n,show_colnames = F,show_rownames = F,annotation_col = ac) ##并且感觉这里加不加_col，得到的图差别不大
所谓数据的中心化是指数据集中的各项数据减去数据集的均值。
例如有数据集1, 2, 3, 6, 3，其均值为3，那么中心化之后的数据集为1-3,2-3,3-3,6-3,3-3，即：-2,-1,0,3,0
所谓数据的标准化是指中心化之后的数据在除以数据集的标准差，即数据集中的各项数据减去数据集的均值再除以数据集的标准差。
例如有数据集1, 2, 3, 6, 3，其均值为3，其标准差为1.87，那么标准化之后的数据集为(1-3)/1.87,(2-3)/1.87,(3-3)/1.87,(6-3)/1.87,(3-3)/1.87，即：-1.069,-0.535,0,1.604,0
> options(digits = 3) ##这里表示系统输出小数点后三位 > data <- c(1,2,3,6,3) > data [1] 1 2 3 6 3 > scale(data,center = T,scale = F) ##scale()中心化：各项减均值 [,1] [1,] -2 [2,] -1 [3,] 0 [4,] 3 [5,] 0 attr(,"scaled:center") [1] 3 > scale(data,center = T,scale = T) ##scale()中心化之后标准化：各项减均值之后再除以标准差 [,1] [1,] -1.069 [2,] -0.535 [3,] 0.000 [4,] 1.604 [5,] 0.000 attr(,"scaled:center") [1] 3 attr(,"scaled:scale") [1] 1.87 > scale(data,center = F,scale = F) ##都为F的话就没动静？？ [,1] [1,] 1 [2,] 2 [3,] 3 [4,] 6 [5,] 3 > scale(data,center = F,scale = T) ##这里没有中心化，只有scale为T,作何解？ [,1] [1,] 0.260 [2,] 0.521 [3,] 0.781 [4,] 1.562 [5,] 0.781 attr(,"scaled:scale") [1] 3.84
原文：https://www.cnblogs.com/SWTwanzhu/p/13083999.html