读取geo表达矩阵

getwd()
setwd("G:/r/r language introduction_2020_12_05/R in action/2020_1205_geo/GEO-master/task7-GSE44076-colon")
getwd()
rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)

library(GEOquery)
# 这个包需要注意两个配置，一般来说自动化的配置是足够的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
  gset <- getGEO('GSE44076', destdir=".",
                 AnnotGPL = F,     ## 注释文件
                 getGPL = F)       ## 平台文件
  save(gset,file=f)   ## 保存到本地
}

getwd()
load('GSE44076_eSet.Rdata')  ## 载入数据
class(gset)  #查看数据类型
length(gset)  #
class(gset[[1]])
gset
exprs(gset[[1]])
gset[[1]]

a=gset[[1]] #
dat=exprs(a) #a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
dim(dat)#看一下dat这个矩阵的维度
# GPL13667
dat[1:4,1:4] #查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
boxplot(dat[,1:10],las=2)

pd=pData(a) #通过查看说明书知道取对象a里的临床信息用pData
head(pd)
head(pd)[1:6,1:3]
dim(pd)
table(pd$title)

library(stringr)
g=str_split(pd$title,' ',simplify = T)[,1]
table(g)
length(g)
dat=dat[,g!='Mucosa']
dim(dat)

group_list=g[g!='Mucosa']
table(group_list)

dat[1:4,1:4] 

############

load(file = 'step1-output.Rdata')
table(group_list)
# 每次都要检测数据
dat[1:4,1:4]
## 下面是画PCA的必须操作，需要看说明书。

## 下面是画PCA的必须操作，需要看说明书。
dat=t(dat)#画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换
dat=as.data.frame(dat)#将matrix转换为data.frame
head(dat)[1:2,1:3]
dat=cbind(dat,group_list) #cbind横向追加，即将分组信息追加到最后一列
head(dat)[1:2,(ncol(dat)-2):ncol(dat)]

1

# 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')