Gene Expression Omnibus database (GEO)是由NCBI负责维护的一个数据库,设计初衷是为了收集整理各种表达芯片数据,但是后来也加入了甲基化芯片,甚至高通量测序数据!
1、GEO Platform (GPL) 芯片平台
2、GEO Sample (GSM) 样本ID号
3、GEO Series (GSE) study的ID号
4、GEO Dataset (GDS) 数据集的ID号
1、用GEOquery从GEO数据库下载数据
#这应该是所有启动R写代码时的操作
> rm(list = ls())#清空所有变量
> options(stringsAsFactors = F)#字符串不设置成因子
正式代码
> library(GEOquery)#加载包
#插播:[[1]]列表取子集,最里面1只能是位置
#用GEOquery包里面的函数getGEO,并用GSE42872该芯片
> eSet <- getGEO("GSE42872",
destdir = '.', #目标目录,下载并读取在当前路径;默认为工作目录
getGPL = F)#不同时下载临床信息#此时的临床信息大;
#getGEO(GEO = NULL, filename = NULL, destdir = tempdir(), SElimits=NULL,GSEMatrix=TRUE,AnnotGPL=FALSE)
-
View(eSet)
-
class(eSet)#显示eSet是一个列表
> class(eSet)
[1] "list"
从GEO下载GSE的数据,读取,赋值给了eSet,当本地目录下有该数据集时,不会下载了;
当信息下载不完整,重新下载或自行网页下载并放在工作目录下。
此时的eSet是一个大的列表,可以用View,一个是实验数据
既然是列表,那么从列表中提取信息就是 eSet[[1]]
(1)提取表达矩阵exp
> exp <- exprs(eSet[[1]])#eSet的第一项就是实验数据,提取表达矩阵。
> exp[1:4,1:4]
> GSM1052615 GSM1052616 GSM1052617 GSM1052618
7892501 7.24559 6.80686 7.73301 6.18961
7892502 6.82711 6.70157 7.02471 6.20493
7892503 4.39977 4.50781 4.88250 4.36295
7892504 9.48025 9.67952 9.63074 9.69200
#View(exp)样本名(竖行名)芯片ID(横行)
#这个时候需要肉眼检查下数据的情况,看是否取对数了
#exp = log2(exp+1)#有表达量exp为0的数据没有意义,所以才加0
exp.png
(2)提取临床信息——将样本表型信息从数据框中提取出
##将样本表型信息从数据框中提取出来【取出来的是表型、样本的数据框】
> pd <- pData(eSet[[1]])
> View(pd)#有汉字乱码
> colnames(pd)
[1] "title" "geo_accession"
[3] "status" "submission_date"
[5] "last_update_date" "type"
[7] "channel_count" "source_name_ch1"
[9] "organism_ch1" "characteristics_ch1"
[11] "characteristics_ch1.1" "treatment_protocol_ch1"
[13] "growth_protocol_ch1" "molecule_ch1"
[15] "extract_protocol_ch1" "label_ch1"
[17] "label_protocol_ch1" "taxid_ch1"
[19] "hyb_protocol" "scan_protocol"
[21] "description" "data_processing"
[23] "platform_id" "contact_name"
[25] "contact_email" "contact_institute"
[27] "contact_address" "contact_city"
[29] "contact_zip/postal_code" "contact_country"
[31] "supplementary_file" "data_row_count"
[33] "cell line:ch1" "cell type:ch1"
pd.png
(3)提取芯片平台编号
> gpl <- eSet[[1]]@annotation#对应平台是GPL6244;
> gpl
[1] "GPL6244"
#View(gpl)#接下来就需要知道GPL6244对应的注释包
> save(pd,exp,gpl,file = "step1output.Rdata")#保存了数据
#总结,下载表达数据,临床信息
网友评论