step2-check.R

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，得到的图差别不大

scale()

1、数据的中心化

所谓数据的中心化是指数据集中的各项数据减去数据集的均值。
例如有数据集1, 2, 3, 6, 3，其均值为3，那么中心化之后的数据集为1-3,2-3,3-3,6-3,3-3，即：-2,-1,0,3,0

2、数据的标准化

所谓数据的标准化是指中心化之后的数据在除以数据集的标准差，即数据集中的各项数据减去数据集的均值再除以数据集的标准差。
例如有数据集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

数据中心化和标准化的意义是一样的，为了消除量纲对数据结构的影响。
在R语言中可以使用scale方法来对数据进行中心化和标准化。

 

参数center表示求数据集各项与均值之差，scale表示数据集各项与均值求差然后再除以数据集的标准差。二者默认为TRUE。

> 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

 

> data
[1] 1 2 3 6 3
> which(data ==6) ##查找，返回data中值为6的下标
[1] 4

 

作者：谢俊飞
链接：https://www.jianshu.com/p/fc82ae05feb9
来源：简书
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。