使用Bioconductor下载GEO(Gene Expression Omnibus)上的数据 57

本文以GSE46106数据为例,讲述如何从GEO上下载数据。分为两个方面,第一,直接从GEO上下载表达数据,第二,直接从GEO上下载CEL文件。

首先我们来看第一种。

> library(GEOquery)
> gset <- getGEO("GSE46106", GSEMatrix =TRUE)
Found 1 file(s)
GSE46106_series_matrix.txt.gz
trying URL 'ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE46nnn/GSE46106/matrix/GSE46106_series_matrix.txt.gz'
ftp data connection made, file length 4110183 bytes
opened URL
==================================================
downloaded 3.9 Mb
 
File stored at: 
/var/folders/n4/11sc2xz13k56hl85z_h2rgq00000gn/T//RtmpBFQdpL/GPL570.soft
> length(gset)
[1] 1
> gset <- gset[[1]]
> head(pData(gset)[,1:5])
> # load NCBI platform annotation
> gpl <- annotation(gset)
> platf <- getGEO(gpl, AnnotGPL=TRUE)
File stored at: 
/var/folders/n4/11sc2xz13k56hl85z_h2rgq00000gn/T//RtmpBFQdpL/GPL570.annot.gz
There were 30 warnings (use warnings() to see them)
> ncbifd <- data.frame(attr(dataTable(platf), "table"))
> eset <- exprs(gset)
> head(eset[,1:5])
          GSM1123782 GSM1123783 GSM1123784 GSM1123785 GSM1123786
1007_s_at    10.1689    10.5247    10.8179    10.3539    10.4964
1053_at       9.6002     7.9436     9.8653     9.9733    10.1960
117_at        5.6808     5.0301     3.7654     2.7751     2.8179
121_at        4.2268     4.6148     4.6147     4.4977     4.6147
1255_g_at     2.1869     2.1869     2.1869     2.1869     2.1869
1294_at       2.1874     2.1874     2.1874     2.1874     2.1874
> head(ncbifd[,1:5])
         ID                                  Gene.title Gene.symbol Gene.ID UniGene.title
1 1007_s_at discoidin domain receptor tyrosine kinase 1        DDR1     780          <NA>
2   1053_at replication factor C (activator 1) 2, 40kDa        RFC2    5982          <NA>
3    117_at        heat shock 70kDa protein 6 (HSP70B')       HSPA6    3310          <NA>
4    121_at                                paired box 8        PAX8    7849          <NA>
5 1255_g_at     guanylate cyclase activator 1A (retina)      GUCA1A    2978          <NA>
6   1294_at ubiquitin-like modifier activating enzyme 7        UBA7    7318          <NA>

接着是第二种。

> library(GEOquery)
> getGEOSuppFiles("GSE46106")
[1] "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE46nnn/GSE46106/suppl/"
trying URL 'ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE46nnn/GSE46106/suppl//GSE46106_RAW.tar'
ftp data connection made, file length 399247360 bytes
opened URL
===============================================
downloaded 380.8 Mb
 
trying URL 'ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE46nnn/GSE46106/suppl//filelist.txt'
ftp data connection made, file length 3308 bytes
opened URL
==================================================
downloaded 3308 bytes
> setwd("GSE46106/")
> dir()
[1] "filelist.txt"     "GSE46106_RAW.tar"
> untar("GSE46106_RAW.tar")
> files <- dir(pattern="gz$")
> sapply(files, gunzip)
GSM1123782_2147TG11_U133plus2.CEL.gz GSM1123783_2147TG15_U133plus2.CEL.gz  GSM1123784_2147TG1_U133plus2.CEL.gz 
                            33380229                             32529207                             33390958 
 GSM1123785_2147TG6_U133plus2.CEL.gz  GSM1123786_2277TG1_U133plus2.CEL.gz  GSM1123787_2277TG9_U133plus2.CEL.gz 
                            33210730                             33082826                             33531908 
GSM1123788_2665TG15_U133plus2.CEL.gz  GSM1123789_2665TG1_U133plus2.CEL.gz  GSM1123790_2665TG6_U133plus2.CEL.gz 
                            32232774                             33380782                             33432878 
 GSM1123791_3104TG1_U133plus2.CEL.gz  GSM1123792_3104TG3_U133plus2.CEL.gz  GSM1123793_3107TG1_U133plus2.CEL.gz 
                            32850644                             32588606                             32584342 
 GSM1123794_3107TG5_U133plus2.CEL.gz  GSM1123795_3143TG1_U133plus2.CEL.gz  GSM1123796_3143TG5_U133plus2.CEL.gz 
                            33358932                             32369391                             32806108 
 GSM1123797_3204TG1_U133plus2.CEL.gz  GSM1123798_3204TG5_U133plus2.CEL.gz  GSM1123799_3561TG1_U133plus2.CEL.gz 
                            32557782                             33016324                             32883029 
 GSM1123800_3561TG5_U133plus2.CEL.gz  GSM1123801_3611TG1_U133plus2.CEL.gz  GSM1123802_3611TG5_U133plus2.CEL.gz 
                            33454274                             32657648                             33126603 
 GSM1123803_3613TG1_U133plus2.CEL.gz  GSM1123804_3613TG5_U133plus2.CEL.gz  GSM1123805_3807TG1_U133plus2.CEL.gz 
                            32581902                             33193826                             32791917 
 GSM1123806_3807TG5_U133plus2.CEL.gz  GSM1123807_3824TG2_U133plus2.CEL.gz  GSM1123808_3887TG1_U133plus2.CEL.gz 
                            33266865                             33078422                             32479022 
 GSM1123809_3887TG5_U133plus2.CEL.gz  GSM1123810_3904TG1_U133plus2.CEL.gz  GSM1123811_3904TG5_U133plus2.CEL.gz 
                            33291287                             32645537                             33425796 
 GSM1123812_3936TG1_U133plus2.CEL.gz  GSM1123813_3936TG5_U133plus2.CEL.gz  GSM1123814_3963TG1_U133plus2.CEL.gz 
                            32939923                             33317196                             32725410 
 GSM1123815_3963TG5_U133plus2.CEL.gz  GSM1123816_4013TG1_U133plus2.CEL.gz  GSM1123817_4169TG2_U133plus2.CEL.gz 
                            33276082                             33151406                             33167461 
 GSM1123818_4175TG1_U133plus2.CEL.gz  GSM1123819_4272TG1_U133plus2.CEL.gz  GSM1123820_4400TG1_U133plus2.CEL.gz 
                            33139829                             33111635                             33206021 
 GSM1123821_4400TG5_U133plus2.CEL.gz  GSM1123822_4664TG1_U133plus2.CEL.gz  GSM1123823_4849TG1_U133plus2.CEL.gz 
                            33026806                             33364390                             33300210 
 GSM1123824_4888TG1_U133plus2.CEL.gz  GSM1123825_4913TG1_U133plus2.CEL.gz 
                            33299396                             33207528 
> filelist <- dir(pattern="CEL$")
> library(affy)
> library(annotate)
> data <- ReadAffy(filenames=filelist)
> affydb<-annPkgName(data@annotation,type="db")
> require(affydb, character.only=TRUE)
> eset<-rma(data,verbose=FALSE)
> eset.e <- exprs(eset)
> library(annaffy)
> symbols<-as.character(aafSymbol(as.character(rownames(eset)),affydb))
> genes<-as.character(aafUniGene(as.character(rownames(eset)),affydb))

视频链接:
youtube

土豆网

更多:http://www2.warwick.ac.uk/fac/sci/moac/people/students/peter_cock/r/geo/

57 thoughts on “使用Bioconductor下载GEO(Gene Expression Omnibus)上的数据

  1. Reply R初学者 12月 27,2013 3:28 上午

    博主你好,我在GEO上下载了数据,就是cel格式的,GSE27342,但是不知道用R怎么能打开,在网上找了一些教程都是不行,求博主帮忙看一下你那能不能用R打开GSE27342的cel文件,能否发一个命令给我,真的非常感谢!

    • Reply admin 12月 27,2013 10:24 上午


      library(GEOquery)
      getGEOSuppFiles("GSE27342")
      setwd("GSE27342/")
      untar("GSE27342_RAW.tar")
      files < - dir(pattern="gz$") filelist <- sapply(files, gunzip) filelist <- dir(pattern="CEL$") head(filelist) ##library(XML) ##url <- "http://www.ncbi.nlm.nih.gov/geo/geo2r/?acc=GSE27342" ##url <- "GEO2R - GEO - NCBI.html" ##tables <- readHTMLTable(url) ##length(tables) ##tables <- tables[[1]] ##head(tables) library(oligo) affyGeneFS <- read.celfiles(filelist)

      更多资料,请参考:使用xps和oligo分析Affymetrix Exon/Gene ST Arrays

  2. Reply R初学者 12月 28,2013 1:38 上午

    非常感谢博主的热心回复,但是我在调用“GEOquery”包是提示我没有安装,我安装是用以下命令”source(“http://www.bioconductor.org/biocLite.R”)biocLite(”GEOquery”)“,请问这个对吗?为什么我的没有安上”GEOquery“包啊。

  3. Reply R初学者 1月 13,2014 8:22 下午

    博主,我想问一下,我在GEO上下载了cel数据,用R做分析时需要知道affymetrix的平台号,比如HG-U133A,但是我想知道怎么能查出我下载的数据平台号是多少,比如GSE43494(GPL14613和GPL15314),非常感谢博主。

    • Reply admin 1月 15,2014 9:46 上午

      我一般是手动搜一下。如果使用代码的话,通常是下载CEL文件,然后通过annPkgName(data@annotation,type=”db”)来取得。至于如何直接从GPL至平台号,我很少用,代码可能是:
      gp14613 <- getGEO('GPL14613') Meta(gp14613)$title 您可以确证一下。

  4. Reply R初学者 1月 15,2014 10:29 上午

    博主,非常感谢您,您帮我看看我这两个的平台是什么,“[miRNA-2_0] Affymetrix Multispecies miRNA-2_0 Array”和“Arraystar Human LncRNA microarray V2.0 (Agilent_033010 Probe Name version)”,我要是从bioconductor上下载注释包平台号应该是什么,麻烦博主帮帮忙。

    • Reply admin 1月 16,2014 4:04 下午

      你中奖了。这两个芯片在bioconductor上都还没有相应的注释文件。所以你需要自己生成注释文件。请参见:http://www.bioconductor.org/packages/release/bioc/html/AnnotationForge.html来自己生成一个吧。注意,第二个不是Affymetrix公司的芯片,是Agilent公司的。

  5. Reply R初学者 1月 17,2014 12:22 上午

    博主,如果一个数据中有两个平台而且还不是一个公司的分析时需要怎么分析啊?比如我这个既有Affymetrix又有Agilent的?

  6. Reply kaji331 5月 7,2014 7:39 上午

    博主你好,我用getGEOSuppFiles下载的RAW.tar解压居然不是cel文件,而是txt和tabular.txt请问这些文件要怎么用或者怎样生成cel文件啊

  7. Reply Jieyun 6月 6,2014 2:59 下午

    Hi, I tried to use GEOquery to download GSE35982_series_matrix.txt.
    Below is my code.
    library(Biobase)
    library(GEOquery)
    library(limma)

    # load series and platform data from GEO

    gset 1) idx <- grep("GPL14767", attr(gset, "names")) else idx <- 1
    gset <- gset[[idx]]

    However, the system showed "Error in download.file(sprintf("ftp://ftp.ncbi.nlm.nih.gov/geo/series/%s/%s/matrix/%s&quot;, :
    cannot open URL 'ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE35nnn/GSE35982/matrix/GSE35982-GPL14767_series_matrix.txt.gz&#039;
    In addition: Warning message:
    In download.file(sprintf("ftp://ftp.ncbi.nlm.nih.gov/geo/series/%s/%s/matrix/%s&quot;, :
    InternetOpenUrl failed: 'The FTP session was terminated".
    Please help to solve it. Thanks !

    • Reply admin 6月 6,2014 8:33 下午

      我运行的完全没有错误。不能重现你的错误。所以请重新整理你的代码,并且确认你的是最新版本。
      我运行的是GEO2R的标准代码。
      # Differential expression analysis with limma
      library(Biobase)
      library(GEOquery)
      library(limma)

      # load series and platform data from GEO

      gset < - getGEO("GSE35982", GSEMatrix =TRUE) if (length(gset) > 1) idx <- grep("GPL14767", attr(gset, "names")) else idx <- 1 gset <- gset[[idx]]

  8. Reply Jieyun 6月 12,2014 9:45 下午

    非常感谢您的回复。
    我也是用GEO2R分析后复制了它的code。但是就是用GEOquery下载不了。。。我把你的命令直接复制了,也运行不了。

  9. Reply R 菜鸟 10月 20,2014 3:52 上午

    您好,使用 affyGeneFS <- read.celfiles(geneCELs)时出现错误,如何解决?

    Loading required package: pd.hugene.1.0.st.v1
    Attempting to obtain 'pd.hugene.1.0.st.v1' from BioConductor website.
    Checking to see if your internet connection works…
    Package 'pd.hugene.1.0.st.v1' was not found in the BioConductor repository.
    The 'pdInfoBuilder' package can often be used in situations like this.
    错误于read.celfiles(geneCELs) :
    The annotation package, pd.hugene.1.0.st.v1, could not be loaded.
    此外: 警告信息:
    In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE, :
    there is no package called ‘pd.hugene.1.0.st.v1’

  10. Reply amiee2406 10月 21,2014 10:09 上午

    gset <- gset[[1]],请问博主这条命令是什么意思,我没看明白。我通过??gset,呈现出很多解释,失去方向?

    • Reply admin 10月 24,2014 8:34 上午

      这个问题很简单,就是gset只是从geo从抓回的数据,它可能是多个数据,所以返回结果保存在了一个list当中。为了以后操作的方便,对于返回长度为1的list,我将其结果从list中抽取出来。就是类似:
      a <- list(a=c(1,2,3)) 以后每次访问c(1,2,3),我都要写成a[[1]]这样,感觉不方便,于是 a <- a[[1]] 这样以后访问c(1,2,3)就只需要写成a就可以了。

  11. Reply amiee2406 10月 21,2014 10:34 上午

    欧博士,在《R语言与Bioconductor生物信息学应用》第五章,分析基因芯片数据,举例采用的是CLL包中的数据集,在5.3.2质量控制中,直接用image(CLLbatch[,1]命令直接作出了灰度图。而我采用的是从GEO下载数据,看到这里我有些衔接不上。Q:(1)您谈及的从GEO下载有下载表达数据及cel文件的两种方法,哪种更适合后续的分析;(2)batch??在质量控制中及后续分析中都直接采用CLLbatch,那从GEO下载的数据如何操作?

    • Reply admin 10月 24,2014 8:44 上午

      你需要注意的是,CLLbatch是ExpressionSet类的一个实例,而getGEO返回的结果是GDS, GPL, GSM, GSE类,本例中第一次返回就是GSE类,同时返回了GSEMatrix。所以直接出灰度图并不适用。如果想出灰度图,可以从原始文件开始。关于问题1,需要分情况。我一般都从cel文件开始,这样可以避免一些小的错误,比如你需要比对不同来源的样品(meta-analysis),还有你只对其中一小部分感兴趣。很多数据因为都使用rma来分析,而rma每次分析,都会综合考虑所有的样品,所以每次分析都需要重新normalization。所以直接下载他人的GSEMatrix是有风险的。但是,对于一般的数据挖掘,只是希望从他人的已知结果中拿到自己感兴趣的部分,有些时候也是可以直接从matrix下手的。对于初学者,我的经验是从cel走一遍,这样即可以和原作者的分析对比较,从而提高,也可以积累经验。
      关于问题2,你可以按照我博客中的教程做。其实都是大同小异。有不明白的再问就好了。

  12. Reply Ying 10月 23,2014 4:07 上午

    博主,问一个不太相关的问题,illumina下面的nextbio网站做的meta analysis 在R中有什么样的包可以代替吗?nextbio不知道怎样上传自己的实验数据,所以没法和别人的研究进行对比

  13. Reply sunjian 10月 28,2014 7:41 下午

    “The data were normalised using quantile normalisation with Bioconductor (v.2.4) Lumi R-package”,处理后的最终的value代表什么?两个标本比较这个value的时候,怎么算倍数关系?是直接相除,还是相减后再算2的次方?

    • Reply admin 10月 30,2014 2:12 下午

      没有上下文,不好判断。但是通常你看到的值会是一个log2转换后的值,你需要检查一下值的范围,是不是主体落在20以内。

  14. Reply JBLXuan 10月 29,2014 12:02 下午

    博主你好 我是R的初学者 也是刚读研究生 应老师的课题要求 在GEO上下载了数据GSE2508 但不是.cel格式的raw.data 询问之后下载了series.matrix.txt的文件 但之前接触的是.cel文件 这个txt文件也不知道怎么读入和处理以及通过这个数据得到差异表达基因 能否给我发一个具体详细点的脚本讲解一下整个过程 是20个样本的那个 GPL8300

    我是初学者 很多不太了解 望不吝赐教

  15. Reply JBLXuan 10月 31,2014 6:49 上午

    知道了 谢谢帮助

  16. Reply Robin Li 11月 7,2014 1:21 上午

    老师您好,我有两个问题想咨询一下~
    1)在使用getGEO()之后,下面一般会显示file stored at,…这个是一个临时文件,一般会在C盘里面,如果我下载的多了之后电脑会很卡,但是清除之后,再次使用会再次下载;我想请问的是如果让它store到我想保存的文件夹里面:比如d、e盘等;我知道有一个可以先下载下来再getGEO的选项,但是当我需要的很多时手动下载就会很麻烦~求教
    2)第二个问题是,挡我get GSE文件的时候,它一般会自动下载GPL来注释,但是有些GPL会很大很大,所以我想在getGEO的时候使用我已经在电脑中存在的GPL进行注释而不是再次下载GPL,请问这个怎么处理?
    谢谢您在百忙之中的回复~

    • Reply admin 11月 12,2014 1:48 下午

      问题1)你可通过调整getGEO中的destdir参数来设置目的目录。
      问题2)你可将getGPL设置成FALSE,然后自己对数据进行注释。

  17. Reply 菜鸟 11月 7,2014 9:40 下午

    博主,你好,我想问一下要是数据来自一个公司不同的平台,要怎么分析啊?

  18. Reply Robin Li 11月 22,2014 6:02 上午

    问题2)你可将getGPL设置成FALSE,然后自己对数据进行注释。

    对的老师,我就是想知道怎么自己进行注释呢?如果说我自己下了本地下载了一个GSE,也下载了一个GPL,那么用什么语句进行注释呢?

  19. Reply 初学者 1月 3,2015 7:25 下午

    博主,我最近在研究怎么下载GEO数据,但是R 软件是最新版本,总是提示要用3.0.0以前的版本,不知道这样的问题怎么解决

  20. Reply washionton 1月 5,2015 2:20 上午

    博主 您好,有一个问题想请教一下,从GEO上下载了*.CEL文件后,用下面的R脚本处理
    mydata <- ReadAffy()
    eset <- rma(mydata)
    write.exprs(eset, file="mydata.txt")
    得到mydata.txt文件。
    问题一:我所得到的mydata.txt文件中数据是GEO网站处理过的吗?还是直接从原始实验中得到的数据,没有经过任何处理?
    问题二:得到mydata.txt文件后,做表达谱芯片的二维聚类分析,
    mydata2 <- exprs(eset)
    mydatascale <- t(scale(t(mydata2)))
    heatmap.2(as.matrix(mydatascale)
    直接进行分析可以吗?还是需要进行其他的处理?
    非常希望得到您的回复。

    • Reply admin 1月 6,2015 3:58 下午

      问题一,不是。你是从CEL文件计算来的eset,并非GEO网站上的表达文本文件。
      问题二,你这是做heatmap,而不是二维聚类,二维聚类建议你再读一下其它的博文。做heatmap也不是把所有的值都拿来做图,那样没有太大的意义。

  21. Reply lingling 1月 7,2015 6:50 上午

    博主 你好 我现在正在学习bioconductor ,我想问一下,从GEO里面下的基因表达谱能整理成数据么,我查了一下能整理的就是cel文件没有找到将表达谱整理成数据的方法的人

  22. Reply biostata 2月 20,2016 4:14 上午

    博主 你好 GEO2R中
    if (length(gset) > 1) idx <- grep("GPL2507", attr(gset, "names")) else idx <- 1
    这里产生idx的目的是?谢谢!

  23. Reply biostata 2月 20,2016 4:16 上午

    # log2 transform
    ex <- exprs(gset)

    qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))

    LogC 100) ||
    (qx[6]-qx[1] > 50 && qx[2] > 0) ||
    (qx[2] > 0 && qx[2] 1 && qx[4] < 2)
    if (LogC) { ex[which(ex <= 0)] <- NaN
    exprs(gset) <- log2(ex) }

    这里qx和logC的目的是什么?谢谢!是否对表达值进行了对数变换,谢谢!

    • Reply admin 2月 22,2016 9:54 上午

      这里就是要猜一下,表达值是否进行了对数变换。如果没有,那就变换一下。这一步不一定是正确的,只是对于大多数来说,应该可以猜个差不多。

  24. Reply keyzhi 4月 5,2016 3:02 上午

    博主 您好,您的视频中pheno=getGSEDataTables(“GSE44408”),这个命令我一运行就出现error : Unknown IO error
    I/O warning : failed to load external entity “http://www.ncbi.nlm.nih.gov/projects/geo/query/acc.cgi?targ=self&form=xml&view=full&acc=GSE44408”
    Error in UseMethod(“xmlRoot”) :
    no applicable method for ‘xmlRoot’ applied to an object of class “NULL”

    请您帮忙看一下是什么问题,打扰了!

    • Reply admin 4月 5,2016 8:29 上午

      你从本机上直接访问一下http://www.ncbi.nlm.nih.gov/projects/geo/query/acc.cgi?targ=self&form=xml&view=full&acc=GSE44408看看会得到什么。理论上应该可以访问到一个xml文件。

  25. Reply Liang L 5月 15,2016 10:00 下午

    欧博,你好,看了你的视频觉得比较容易理解。想问一下,你还有哪些地方上传了教学视频吗?

  26. Reply 杨建课 7月 7,2016 11:02 上午

    老师,您好!最近用R代码下载GSE数据
    getGEOSuppFiles(“GSE3325”, makeDirectory = TRUE, baseDir = getwd())

    运行时,下载了一小小部分就自动停止了,请问是网速的问题还是编码或设置哪里出问题了呢?

  27. Reply 小江 2月 28,2017 1:54 上午

    楼主,怎么实现利用R或其它软件实现GEO数据的下载与处理
    (1)实现GEO数据的raw.tar下载
    (2)提取样本和表达谱信息
    (3)将探针ID对到GeneID
    (4)缺失值处理

Leave a Reply

  

  

  

%d 博主赞过: