GEO数据挖掘3

rm(list = ls())
options(stringsAsFactors = F)
load(file = 'step1-output.Rdata')
table(group_list)

## group_list
##     Control Vemurafenib
##           3           3

# 查看数据
dat[1:4,1:4]

##        GSM1052615 GSM1052616 GSM1052617 GSM1052618
## ZZZ3     11.26970   11.12560   11.23260   11.47360
## ZZEF1     9.63412    9.44327    9.67075    9.96091
## ZYX      10.89980   10.93190   10.91850   10.71250
## ZYG11B   10.48080   10.32370   10.51680   10.74500

## 绘制pca
#画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换
dat=t(dat)
#将matrix转换为data.frame
dat=as.data.frame(dat)
#cbind添加分组信息
dat=cbind(dat,group_list)
# 安装需要的包
#设置镜像
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor")

# install.packages(c("FactoMineR", "factoextra"))
library("FactoMineR")
library("factoextra")

# 数据处理，去掉分组信息
dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)
# 绘制pca图
fviz_pca_ind(dat.pca,
             geom.ind = "point", #
             col.ind = dat$group_list, # 根据颜色分组
             addEllipses = TRUE, # 椭圆
             legend.title = "Groups"
)

## 热图绘制
rm(list = ls())
load(file = 'step1-output.Rdata')
dat[1:4,1:4]

##        GSM1052615 GSM1052616 GSM1052617 GSM1052618
## ZZZ3     11.26970   11.12560   11.23260   11.47360
## ZZEF1     9.63412    9.44327    9.67075    9.96091
## ZYX      10.89980   10.93190   10.91850   10.71250
## ZYG11B   10.48080   10.32370   10.51680   10.74500

#apply按行取每一行的方差，从小到大排序,
#也就是取样本间差异较大的基因，取最大的1000个
cg=names(tail(sort(apply(dat,1,sd)),1000))
# install.packages("pheatmap")
library(pheatmap)
# 绘制1000个位点的热图
pheatmap(dat[cg,],# 取差异明显的1000个基因
         show_colnames =F,show_rownames = F)

# 对数据进行归一化
# 因为是按照基因归一化，所以先进行转置，然后再转置回去
n=t(scale(t(dat[cg,])))
# 对绝对值大于2的数取绝对值2
# 使得最后的数据范围控制在2以内
n[n>2]=2 
n[n< -2]= -2
n[1:4,1:4]

##           GSM1052615 GSM1052616 GSM1052617 GSM1052618
## SERPINI1  -0.5864330 -1.3172509 -0.6733724  1.2702722
## SLC44A5   -0.7866286 -0.8619907 -1.0508698  1.1604191
## TIMM8A     0.9119733  0.6284977  1.1376726 -0.7681696
## LINC01530 -0.9223885 -0.9223885 -0.8554885  0.9206052

# 最终的结果显示表达量归一化为2以内
pheatmap(n,show_colnames =F,show_rownames = F)

# 现在的图只是热图，但是没有分组信息
# 添加分组信息
ac=data.frame(g=group_list)
# 添加样本名为行名
rownames(ac)=colnames(n)
# 再次绘制热图
pheatmap(n,show_colnames =T,show_rownames = F,
         annotation_col=ac)

# 可以看出两个分组之间存在不少的差异表达