在上一篇简书中,我总结了如何下载数据链接如下:[公共数据库挖掘之利用R语言下载并整理GEO数据] https://www.jianshu.com/p/fd5c8389fb53。
数据下载好之后,如何对数据进行处理使其能够用R进行作图,也是关键的一步。
之前我们已经获取了GSE数据,代码如下:
#下载数据,载GSE数据,获得表达量矩阵
options( 'download.file.method.GEOquery' = 'libcurl' )
gset <- getGEO('GSE42872',destdir = ".",
AnnotGPL = F,
getGPL = F)
#保存数据
save(gset,file = 'GSE42872.gset.Rdata')
在代码中,我们没有获得GPL数据,因此,要在下载GPL数据。获得GPL有两种方法,第一种方法就是在上述代码中修改一处,getGPL=T。
ob=gset[[1]]
mat=exprs(ob)
library(GEOquery)
gpl <- getGEO('GPL6244',destdir = ".")
colnames(Table(gpl))
第二种方法,我们要先对该GSE42872的所处的GPL进行查询,找到其GPL值(参考生信菜鸟团,jimmy老师分享的帖子)。
通过该值对应的R包,我们可以下载GPL文件。方法如下:
source("https://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
biocLite("hugene10sttranscriptcluster.db")
suppressPackageStartupMessages(library(hugene10sttranscriptcluster.db))
ls("package:hugene10sttranscriptcluster.db")
#获取对应关系
ids=toTable(hugene10sttranscriptclusterSYMBOL)
length(unique(ids$symbol))
tail(sort(table(ids$symbol)))
table(sort(table(ids$symbol)))
plot(table(sort(table(ids$symbol))))
获得所需要的gpl之后,我们需要对获得的表达矩阵进行过滤。
exprSet=mat
table(rownames(exprSet) %in% ids$probe_id)
dim(exprSet)
exprSet=exprSet[rownames(exprSet) %in% ids$probe_id,]
dim(exprSet)
将表达矩阵和的基因名称与表达量进行整合
new_exprSet = by(exprSet,ids$symbol,function(x) rownames(x)[which.max(rowMeans(x))])
uniprobes = as.character(new_exprSet)
exprSet=exprSet[rownames(exprSet)%in%uniprobes,]
rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]
head(exprSet)
至此,数据整理结束,可以挑选自己感兴趣的基因进行作图了。
例如:
a<- head(exprSet)
a
可以得到下边的数据
image.png
网友评论