rm(list = ls()) options(stringsAsFactors = F)
# 注意查看下载文件的大小，检查数据
f='GSE17215_eSet.Rdata'
# 使用GEOquery包下载GEO数据（需提前下载安装此包）
library(GEOquery)
# 下载可能会出错，可以试试更改镜像、翻墙，重启，rm（list=ls（），据说早上的网好容易下载...看运气了，下2为更改镜像options(repos<- c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/") )
options("BioC_mirror"<- "https://mirrors.ustc.edu.cn/bioc/")
if(!file.exists(f)){  
gset <- getGEO('GSE17215', destdir=".",                 
AnnotGPL = F,     ## 注释文件                 
getGPL = F)       ## 平台文件
save(gset,file=f)   ## 保存到本地
}
# 先保存后载入，养成好习惯
load('GSE17215_eSet.Rdata')
# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的
listclass(gset)
length(gset)
class(gset[[1]])
# not list易于操作
a=gset[[1]]
# 获取表达矩阵
dat=exprs(a)dim(dat)
# 加载注释包
library(hgu133a.db)
# 看一下包里有什么
ls("package:hgu133a.db")
# 获取探针与基因的对应关系
ids=toTable(hgu133aSYMBOL)
head(ids)


#这块代码报错
#Error in dat[ids$probe_id, ] : subscript out of bounds
# 筛选ids中包含的probe_id，将两个表行名设为同序
dat=dat[ids$probe_id,]




dat[1:4,1:4]
# 取dat表达量的中位数，生成ids的median列
ids$median=apply(dat,1,median)
# 将ids按照symbol和median排序（降序），可以看一下前后对比关系进行理解ids=ids[order(ids$symbol,ids$median,decreasing = T),]
# 去除重复的symbol值，保留第一位，即表达量最大的一位（dim一下看看前后对比）
ids=ids[!duplicated(ids$symbol),]
# 同上
dat=dat[ids$probe_id,]
# 先同序，再直接改为symbol（dim一下瞧瞧两者行数）
rownames(dat)=ids$symboldat[1:4,1:4]  
dim(dat)
# 好，进入正题
ng='ACTR3B ANLN BAG1 BCL2 BIRC5 BLVRA CCNB1 CCNE1 CDC20 CDC6 CDCA1 CDH3 CENPF CEP55 CXXC5 EGFR ERBB2 ESR1 EXO1 FGFR4 FOXA1 FOXC1 GPR160 GRB7 KIF2C KNTC2 KRT14 KRT17 KRT5 MAPT MDM2 MELK MIA MKI67 MLPH MMP11 MYBL2 MYC NAT1 ORC6L PGR PHGDH PTTG1 RRM2 SFRP1 SLC39A6 TMEM45B TYMS UBE2C UBE2T'
# 把基因加上“”，使之与rownames(dat)格式一致
ng=strsplit(ng,' ')[[1]]
# 看一下，多少基因存在于dat中，可得41个TRUE&9个FALSE
table(ng %in%  rownames(dat))
# 清洗掉不存在的ng，注意这一步存在排序（连同下一步理解）
ng=ng[ng %in%  rownames(dat)]
# 再改行名
dat=dat[ng,]
# 画图
dat=log2(dat)pheatmap::pheatmap(dat,scale = 'row')




红色标注的地方求解为啥会报错，怎么解决，求高手指导！