对下载好的基因初步分析,进行PCA分析和热图绘制
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)
# 可以看出两个分组之间存在不少的差异表达
此部分的分析较为基础,为GEO分析的初步探索
love&peace