非常详细的GEO数据下载和清理，你要看吗-GSE33532
2019-11-21 本文已影响0人
小梦游仙境
来自数据集GSE33532,对于第一步获得表达矩阵，做了还蛮详尽的注解过程，希望能对其他人有帮助。虽然看上去全是黑色代码快，但注解理解示例全在这里面了哦！
rm(list = ls()) ## 魔幻操作，一键清空~
options()$repos
options()$BioC_mirror
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options(stringsAsFactors = F)#在调用as.data.frame的时，将stringsAsFactors设置为FALSE可以避免character类型自动转化为factor类型
################## 1. 下载表达矩阵和临床信息 ###################
################## 得到表达矩阵dat和分组信息group_lis ##################
f='GSE33532_eSet.Rdata'
# https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE33532
library(GEOquery)
# 这个包需要注意两个配置，一般来说自动化的配置是足够的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
gset <- getGEO('GSE33532', destdir=".",
AnnotGPL = F, ## 注释文件
getGPL = F) ## 平台文件
save(gset,file=f) ## 保存到本地
}
load('GSE33532_eSet.Rdata') ## 载入数据
class(gset) #查看数据类型
length(gset) #
class(gset[[1]])
gset
# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的list
a=gset[[1]] #
dat=exprs(a) #a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
dim(dat)
dat[1:4,1:4] #查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
pd=pData(a)
tumor=rownames(pd[grepl('tumor',as.character(pd$title)),])
normal=rownames(pd[grepl('normal',as.character(pd$title)),])
dat=dat[,c(tumor,normal)]
group_list=c(rep('tumor',length(tumor)),
rep('normal',length(normal)))
table(group_list)
dim(dat)
dat[1:4,1:4]
dat <- as.data.frame(dat)
save(dat,group_list,file = 'step1.1_dat_group.Rdata')
################## 2.获得注释信息 ##################
################## 三种方法 ##################
#1.用getGEO下载平台GPL的注释信息
#2.在GEO网站上下载好平台注释信息后倒入R studio
#3.用bioconductor上的包
#方法一 用getGEO下载平台GPL的注释信息
if(F){
library(GEOquery)
gpl<- getGEO('GPL570', destdir=".")
dim(gpl)
colnames(Table(gpl)) #查一下列明
head(Table(gpl)[,c(1,11)])
ids=Table(gpl)[,c(1,11)]
##踩坑放松1:会不会有NA值，检查一下,没问题就不管那些个警告了
tmpids <- na.omit(ids)
dim(tmpids)
head(ids)
save(ids,file='step1.2_ids.Rdata') #三种方法均获得的ids均设置的是'step1.2_ids.Rdata'
}
rm(tmpids)
#方法二 在GEO网站上下载好平台注释信息后导入Rstudio
if(F){
gpl=read.table('GPL570-55999.txt', #提前下载好了平台GPL570的注释信息
sep = '\t',quote = '',
fill=T,header = T)
dim(gpl)
colnames(gpl)
ids=gpl[,c(1,11)]
ids=ids[ids$Gene.Symbol != '',] #去除空的
dim(ids)
head(ids)
save(ids,file='step1.2_ids.Rdata') #三种方法均获得的ids均设置的是'step1.2_ids.Rdata'
}
#方法三 用bioconductor网站上对应平台的基因注释包
#https://mp.weixin.qq.com/s/5_ceRNaRFP3h08Pjd015Wg
if(T){
library(hgu133plus2.db)
ls("package:hgu133plus2.db")
ids<- toTable(hgu133plus2SYMBOL)
dim(ids)
head(ids)
save(ids,file='step1.2_ids.Rdata') #三种方法均获得的ids均设置的是'step1.2_ids.Rdata'
}
################三种获得注释信息的方法，所得到的对应的行数如下：
#> dim(gpl)
#[1] 54675 16
#> dim(ids)
#[1] 45782 2
#> dim(ids)
#[1] 41932 2
################## 3. ID转换 ##################
################## intersect/merge/%in% ##################
################## match函数 ##################
################## （1）用intersect函数 ##################
rm(list = ls())
load('step1.1_dat_group.Rdata')
load('step1.2_ids.Rdata') #此处按照第三种从bioconductor获得的注释包来进行id转换
dim(dat)
dim(ids)
table(rownames(dat) %in% ids$probe_id)
intersect(rownames(dat),ids$probe_id)
length(intersect(rownames(dat),ids$probe_id))
#xt意味着xt(相同)的
xt <- intersect(rownames(dat),ids$probe_id)
dim(xt)
##踩坑放松1:为什么用dim是得到NULL?
length(xt)######xt现在是一个character，所以用length来查看向量的长度
head(xt)
ids_tmp <- ids[ids$probe_id==xt,]#注：想将ids按照取交集后的向量（xt）从ids中取交集
##踩坑放松2:报错-长的对象长度不是短的对象长度的整倍数
length(ids$probe_id)
length(xt)
###tips1：上面会提示报错，因为当向量用'=='进行判断的时候，就是将两个向量中的元素进行一对一地判断是否相同，之所以报‘长的对象长度不是短的对象长度的整倍数’的报错，是因为ids$probe_id有ids的行数那么多（41932个），而xt则是取交集后的23738个。
###tips2：接下来用match函数取子集，就是按照xt的顺序从大子集ids中取，那么取出来的探针的顺序肯定是和xt的顺序是一样的，然后再从dat中同样按照xt的顺序来取子集，因为ids和dat都是按照xt中的探针的顺序取子集，那么二者的dat中的探针名肯定和ids中的探针名是一一对应的了。
#################### 重磅嘉宾-match函数第一次出场 ####################
ids <- ids[match(xt,ids$probe_id),]##ids同样按照probe_id这一列提取子集，是按照坐标位置提取
#用identical函数和setequal函数验证一下(这两个函数都可以)，是不是完全一样，也就是说表达矩阵dat的探针名和交集（xt）是否是一一对应的
identical(ids$probe_id,xt)
setequal(ids$probe_id,xt)
save(ids,file='step1.3_ids.Rdata')
#################### （2）用merge函数 ##################
rm(list = ls())
load('step1.1_dat_group.Rdata')
load('step1.2_ids.Rdata')
ids<- toTable(hgu133plus2SYMBOL)
dat_tmp <- dat #重新赋值一个dat_tmp
ids_tmp <- ids #重新赋值一个ids_tmp
dim(dat_tmp)
dim(ids_tmp)
#用新赋值后的做探索
dat_tmp <- dat_tmp[ids_tmp$probe_id,] #注：想从dat_tmp表达矩阵中按照ids_tmp$probe_id取出的探针名提取子集，dat_tmp的rownames就是探针名，表面看是没问题呀
##踩坑放松3: dat_tmp <- dat_tmp[ids_tmp$probe_id,]无报错，表面看似风平浪静，背后其实翻江倒海，用dim查看大小
dim(dat_tmp) #注：刚才的dat_tmp只有25906个探针，现在竟然和ids_tmp的探针名一样多的个数！表达矩阵竟然变大了！
is.na(dat_tmp) #注：我就想看个行名是不是NA，结果在表达矩阵中返回了TRUE和FALSE，所以又学一了is.na函数
table(is.na(rownames(dat_tmp))) #注：还是不对，改成行名也不对。
dat_tmp[1:4,1:4]
is.na(dat_tmp)#注：观察在dat_tmp中为TRUE也就是为NA值的，行名也都是NA,因此用grepl函数抓取下
table(grepl('NA',rownames(dat_tmp))) #注意区分grep和grepl函数的差别。注意此时的结果
#FALSE TRUE
#23738 18194
#注意观察：上面的FALSE正好为23738个，同前面用intersect函数获得的ids是不是行数相同？
dat_tmp <- na.omit(dat_tmp) #注：用na.omit去掉NA值
dim(dat_tmp)
##为用merge函数作准备-dat新建一列，列名是probe
dat$probe_id <- rownames(dat) #注：dat新建一列，列名是probe，此步是为了后面用merge函数
dim(dat)
####小提问：如何知道dat中第101列的列名
colnames(dat[,101]) #注：dat新建probe这一列后，dat现在有101列，我想看看这第101例的列名
##踩坑放松4:colnames(dat[,101])为什么返回NULL,如何获得第101列的列名？
#分解代码：先看dat[,101]
dat[,101]
##报错原因：把上面的代码分解，dat[,101]，根据“行在左，列在右”，此时101在“，”右边，因此其实取出的是第101这列的所有的探针名，是一个向量，一个长度同行数相等的向量。向量取子集怎么能用colnanmes来提取呢？
#正确操纵如下：
colnames(dat)
colnames(dat)[101]#注，从向量中取子集，可以按照逻辑值和坐标位置，现在我们想取坐标位置为第101位的内容
#关于merge函数，依然用前两步得到的dat和ids
ids_dat_merge <- merge(ids,dat,by.x ='probe_id',by.y='probe_id')
##踩坑放松5:
#(1)merge函数：dat新建的一列，如果列名是probe，也就是说列名不一样，可不可以？
#(2)merge函数：by.x=probe,by.y=probe_id,怎么改？出来的矩阵是什么样的？
#(3)merge函数：如果想要merge的那一列的向量顺序不一样，merge后的矩阵是怎样的,有什么关系？
#下面是merge函数小理解，新建模拟ids_tmp和dat_tmp
if(F){
#第1、2个merge问题
ids_tmp <- data.frame(probe_id = c('1053_at','117_at','121_at','1255_g_at'),
symbols = c('RFC2','HSPA6','PAX8','GUCA1A'))
dat_tmp <- data.frame(row.names = c('1053_at','117_at','121_at','1316_at','1320_at'),
GSM835268=c(1,1,1,1,1),
GSM835269=c(2,2,2,2,2),
probe = c('1053_at','117_at','121_at','1316_at','1320_at'))
ids_tmp
dat_tmp
merge1 <- merge(ids_tmp,dat_tmp,by.x='probe_id',by.y='probe')
merge1
merge2 <- merge(dat_tmp,ids_tmp,by.x='probe',by.y='probe_id')
merge2
#第3个merge问题，如果要merge的探针名顺序不一样
ids_tmp <- data.frame(probe_id = c('1053_at','117_at','121_at','1255_g_at'),
symbols = c('RFC2','HSPA6','PAX8','GUCA1A'))
dat_tmp <- data.frame(row.names = c('121_at','117_at','1053_at','1316_at','1320_at'),
GSM835268=c(1,1,1,1,1),
GSM835269=c(2,2,2,2,2),
probe = c('121_at','117_at','1053_at','1316_at','1320_at'))
ids_tmp
dat_tmp
merge3 <- merge(ids_tmp,dat_tmp,by.x='probe_id',by.y='probe')
merge3
merge4 <- merge(dat_tmp,ids_tmp,by.x='probe',by.y='probe_id')
merge4
#从merge3、merge4相同的那列顺序没有改变，这是最后的矩阵中dat和ids顺序有差别
}
dim(ids_dat_merge) #注：检查数据，dim看大小
ids_dat_merge[1:4,1:4] #注：取子集看数据情况
ids_dat_merge <- na.omit(ids_dat_merge) #注：检查下是否有NA值，注意理解merge函数，通过merge函数获得ids_dat_merge是ids和dat中*都都都*存在的元素。思考，如果想得到merge(x,y)后，在x中有不存在与y中的元素也出现在merge后的矩阵中，用什么参数？
dim(ids_dat_merge)
ids_dat_merge####用$符号看一下，现在的表达矩阵有哪些列（列名）按完$健，按tab健
##接下来进行从ids中按照merge后的ids_dat_merge表达矩阵中按照探针名字，从ids中取子集
ids_tmp <- ids[ids_dat_merge$probe_id,] #注：我已经得到了dat和ids共有列探针（probe）merge后的表达矩阵了，想按照merge后的ids_dat_merge中的probe_id这一列提取探针名字，然后将ids按照这个探针名字提取子集从而想生成一个ids_tmp，从而获得和表达表达矩阵的探针名相同的ids注释！
dim(ids_tmp)
ids_tmp[1:4,1:2]#注：一看数据具体情况，竟然都是NA!
#踩坑放松6：从数据框中提取子集，可以按rownames也就是行名，也可以按照坐标，可是为什么都是NA呢？
dim(ids)
ids[1:4,1:2]
#答：ids是行名为1、2、3、4...，并不是探针名！从ids中取子集时，要么按照探针名(行名)，要么按照坐标位置,要么按照逻辑值，就这三种，三种，三种！重要的事情说三遍哦！而ids_dat_merge$probe_id后得到的是探针名probe_id，而ids的行名又不是探针名probe_id，而是1、2、3、4...,所以此处取出的都是NA!
#既然不是，那我把ids的行名变成探针名probe_id不就好了嘛！
rownames(ids_tmp) <- ids_dat_merge$probe_id
##看下大小
dim(ids_tmp)
ids_tmp[1:4,1:2] #注：虽然重新复制了探针名，但是需要注意此处的ids_tmp由于上一步取子集的错误操作，现在都是NA,所以需要重新获得一个ids_tmp.
#################### 重磅嘉宾-match函数第二次出场 ####################
ids_tmp <- ids[match(ids_dat_merge$probe_id,ids$probe_id),]
#注：上面的match函数，从ids$probe_id的探针名中按照ids_dat_merge$probe_id在ids$probe_id中的排名位置，也就是坐标位置，从ids中按照这个坐标位置取子集，得到ids_tmp
######下面是match函数小理解
if(T){
use_id <- c('A','B','C','D')
use_id
m<-data.frame(u=c('B','C','D','A'))
match(use_id,m$u)
a<-m[match(use_id,m$u),]
a
}
dim(ids_tmp)
ids_tmp[1:4,1:2]
ids_tmp <- na.omit(ids_tmp)#注：被NA弄怕了，再用na.omit检查下，match函数用熟练或用对了，就可省略
ids <- ids_tmp
save(ids,file='step1.3_ids.Rdata')
###################### 此处小贴士 ###################
#其实前面的到的ids_dat_merge矩阵中已经含有基因symbol这一列了，是不是可以就把矩阵中的symbol变成行名就可以了呢？不就已经完成id转换了嘛？那试一下
ids_dat_merge[1:4,1:4]
rownames(ids_dat_merge) <- ids_dat_merge$symbol
##踩坑放松8:rownames(ids_dat_merge) <- ids_dat_merge$symbol 报错，不允许有重复的‘row.names’,因此基因名symbol需要去除重复！
#################### （3）用%in%函数 ##################
rm(list = ls())
load('step1.1_dat_group.Rdata')
load('step1.2_ids.Rdata')
dim(dat)
dim(ids)
ids <- ids[ids$probe_id %in% rownames(dat),]#注：想从dat的rownames中提取出ids$symbol后存在于rownames(dat)中的元素，%in%是进行判断，逻辑值为TRUE的，从ids中提取出来。
dim(ids)
ids[1:4,1:2]#注：此时为什么是ids[1:4,1:2]，如果用ids[1:4,1:4]会报错undefined columns selected，是因为ids本身就只有两列。
##踩坑放松9:此处是温馨小提示：同前面踩坑放松6是类似的道理，虽然现在ids的行名是1、2、3、4，但是%in%是进行判断，得到的是逻辑值，所以不用管行名是什么，我们就是按照TRUE的提取！三种提取子集的方法：名字，坐标，逻辑值！
#下面是%in%函数小理解,新建模拟ids_tmp和dat_tmp.
if(F){
ids_tmp <- data.frame(probe_id = c('1053_at','117_at','121_at','1255_g_at'),
symbols = c('RFC2','HSPA6','PAX8','GUCA1A'))
dat_tmp <- data.frame(row.names = c('1053_at','117_at','121_at','1316_at','1320_at'),
GSM835268=c(1,1,1,1,1),
GSM835269=c(2,2,2,2,2))
ids_tmp
dat_tmp
m <- ids_tmp[ids_tmp$probe_id %in% rownames(dat_tmp),]
m
n <- dat_tmp[rownames(dat_tmp) %in% ids_tmp$probe_id,]
n
#如果探针名顺序不一样
ids_tmp <- data.frame(probe_id = c('1053_at','117_at','121_at','1255_g_at'),
symbols = c('RFC2','HSPA6','PAX8','GUCA1A'))
dat_tmp <- data.frame(row.names = c('121_at','1316_at','1053_at','117_at','1320_at'),
GSM835268=c(1,1,1,1,1),
GSM835269=c(2,2,2,2,2))
#注意上面dat_tmp中的row.names的顺序。
ids_tmp
dat_tmp
m <- ids_tmp[ids_tmp$probe_id %in% rownames(dat_tmp),]
m
n <- dat_tmp[rownames(dat_tmp) %in% ids_tmp$probe_id,]
n
#上面探针名顺序不一样，可以帮助理解取子集，从哪个矩阵取，就按哪个矩阵中的逻辑值取
}
identical(ids$probe_id,rownames(dat))
ids[1:4,1:2]
dat[1:4,1:4]
#################### 重磅嘉宾-match函数第三次出场 ####################
dat <- dat[match(ids$probe_id,rownames(dat)),]#注：ids的probe_id是全部存在于dat的rownames里的,因此用match函数从dat中按照ids的顺序提取子集，一定是能得到探针完全对应的。此处就是按照坐标在提取子集，这个坐标指的就是a在b中的顺序，注意，如果数据框的行名是1、2、3、4，与这个1、2、3、4无关，仅仅就是a中的元素一个一个试验在b中排第几位，比如a中的元素依次在b中的元素进行一一对应，排第1位记个1，排第5位记个5，排第10位记个10，然后从b中按照这个1、5、10提取b中的元素。所以从b中提取的元素就是和a完全相同的。
save(ids,file = 'step1.3_ids.Rdata')#注，上面的dat取子集是按照坐标位置提取的，此处并没有保存dat，但是要知道此出dat用match函数提取子集的方法。
####################### 4.将探针名转换为基因名 ################
###################### 基因名symbol去重复 ######################
rm(list = ls())
load('step1.1_dat_group.Rdata')
load(file='step1.3_ids.Rdata')
dim(dat)
dat[1:4,1:4]
dim(ids)
ids[1:4,1:2]
dat <- dat[ids$probe_id,]#注：从dat表达矩阵（其实是数据框，就是这么叫一下表达矩阵）中按照ids$probe_id也就是探针名字，提取子集。
######下面是数据框提取子集的小理解，按照rownames的名字或者坐标位置均可以提取子集，示例如下：
if(F){
x <- data.frame(row.names = c('a','b','c'),xm = c('xh','xl','xh'),nl=c(15,18,21))
x
x['a',]
x[1,]
identical(x['a',],x[1,])
}
dim(dat)
dim(ids)
ids[1:4,1:2]
#同理，看看dat和ids能不能一一对应
identical(rownames(dat),ids$probe_id)
setequal(rownames(dat),ids$probe_id)
rownames(dat) <- ids$symbol
#再踩坑一次！这里rownames(dat) <- ids$symbol #这个时候由于rownames的symbol有重复的，会报错，所以先需要进行去重复,还是和上面相同的报错，再一次提醒我们：不允许有重复的'row.names'，因此需要基因名去重复
##################基因名symbol去重复如下：
if(T){
ids$median=apply(dat,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids
ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s
dim(ids)
ids[1:3,1:3]
dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中probe_id所在的每一行组成一个新的dat，这时获得的dat的探针名就是和ids的probe_id这一列是完全对应的
rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat，作为dat的rowname
dim(dat)
dat[1:4,1:4]
#identical(rownames(dat),ids$symbol) #可验证，也可不验证，只要能确保 dat=dat[ids$probe_id,]这步是正确的
boxplot(dat)
dat=log(dat+1)
save(dat,ids,group_list,file = 'step1.4_dat_ids_group_list.Rdata')
}
