本文以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/
博主你好，我在GEO上下载了数据，就是cel格式的，GSE27342，但是不知道用R怎么能打开，在网上找了一些教程都是不行，求博主帮忙看一下你那能不能用R打开GSE27342的cel文件，能否发一个命令给我，真的非常感谢!
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
非常感谢博主的热心回复，但是我在调用“GEOquery”包是提示我没有安装，我安装是用以下命令”source(“http://www.bioconductor.org/biocLite.R”)biocLite（”GEOquery”)“，请问这个对吗？为什么我的没有安上”GEOquery“包啊。
请问你的错误提示是什么？同时把你的sessionInfo()提供一下。
博主，我想问一下，我在GEO上下载了cel数据，用R做分析时需要知道affymetrix的平台号，比如HG-U133A，但是我想知道怎么能查出我下载的数据平台号是多少，比如GSE43494（GPL14613和GPL15314），非常感谢博主。
我一般是手动搜一下。如果使用代码的话，通常是下载CEL文件，然后通过annPkgName(data@annotation,type=”db”)来取得。至于如何直接从GPL至平台号，我很少用，代码可能是：
gp14613 <- getGEO('GPL14613') Meta(gp14613)$title 您可以确证一下。
博主，非常感谢您，您帮我看看我这两个的平台是什么，“[miRNA-2_0] Affymetrix Multispecies miRNA-2_0 Array”和“Arraystar Human LncRNA microarray V2.0 (Agilent_033010 Probe Name version)”，我要是从bioconductor上下载注释包平台号应该是什么，麻烦博主帮帮忙。
你中奖了。这两个芯片在bioconductor上都还没有相应的注释文件。所以你需要自己生成注释文件。请参见：http://www.bioconductor.org/packages/release/bioc/html/AnnotationForge.html来自己生成一个吧。注意，第二个不是Affymetrix公司的芯片，是Agilent公司的。
博主，如果一个数据中有两个平台而且还不是一个公司的分析时需要怎么分析啊？比如我这个既有Affymetrix又有Agilent的？
通常是不能这样比较的。但是在一些meta-analysis中，可以先各自分析，在得到normalized结果后可以相互比较。
博主你好，我用getGEOSuppFiles下载的RAW.tar解压居然不是cel文件，而是txt和tabular.txt请问这些文件要怎么用或者怎样生成cel文件啊
那只能说明原本就没有上传cel文件，所以你也无从得到cel文件。你可以直接使用tabular.txt里的值进行分析。
那想请问bioconductor怎么处理GEO上的soft格式文件呢？
这个我倒没有分析过。文本文件，应该不难分析吧。
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", :
cannot open URL 'ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE35nnn/GSE35982/matrix/GSE35982-GPL14767_series_matrix.txt.gz'
In addition: Warning message:
In download.file(sprintf("ftp://ftp.ncbi.nlm.nih.gov/geo/series/%s/%s/matrix/%s", :
InternetOpenUrl failed: 'The FTP session was terminated".
Please help to solve it. Thanks !
我运行的完全没有错误。不能重现你的错误。所以请重新整理你的代码，并且确认你的是最新版本。
我运行的是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]]
非常感谢您的回复。
我也是用GEO2R分析后复制了它的code。但是就是用GEOquery下载不了。。。我把你的命令直接复制了，也运行不了。
把你的sessionInfo信息提供给我看看。
Hi,GEOquery这个问题已经解决了，在自己的电脑上可以下载。可能是国外的电脑firewall有限制。多谢您。
您好，使用 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’
去网页下载pd.hugene.1.0.st.v1至本地安装。有可能是你的网络连接出了问题，或者说bioconductor被墙了。
gset <- gset[[1]]，请问博主这条命令是什么意思，我没看明白。我通过??gset，呈现出很多解释，失去方向？
这个问题很简单，就是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就可以了。
欧博士，在《R语言与Bioconductor生物信息学应用》第五章，分析基因芯片数据，举例采用的是CLL包中的数据集，在5.3.2质量控制中，直接用image(CLLbatch[，1］命令直接作出了灰度图。而我采用的是从GEO下载数据，看到这里我有些衔接不上。Q：（1）您谈及的从GEO下载有下载表达数据及cel文件的两种方法，哪种更适合后续的分析；（2）batch？？在质量控制中及后续分析中都直接采用CLLbatch，那从GEO下载的数据如何操作？
你需要注意的是，CLLbatch是ExpressionSet类的一个实例，而getGEO返回的结果是GDS, GPL, GSM, GSE类，本例中第一次返回就是GSE类，同时返回了GSEMatrix。所以直接出灰度图并不适用。如果想出灰度图，可以从原始文件开始。关于问题1，需要分情况。我一般都从cel文件开始，这样可以避免一些小的错误，比如你需要比对不同来源的样品（meta-analysis），还有你只对其中一小部分感兴趣。很多数据因为都使用rma来分析，而rma每次分析，都会综合考虑所有的样品，所以每次分析都需要重新normalization。所以直接下载他人的GSEMatrix是有风险的。但是，对于一般的数据挖掘，只是希望从他人的已知结果中拿到自己感兴趣的部分，有些时候也是可以直接从matrix下手的。对于初学者，我的经验是从cel走一遍，这样即可以和原作者的分析对比较，从而提高，也可以积累经验。
关于问题2，你可以按照我博客中的教程做。其实都是大同小异。有不明白的再问就好了。
博主，问一个不太相关的问题，illumina下面的nextbio网站做的meta analysis 在R中有什么样的包可以代替吗？nextbio不知道怎样上传自己的实验数据，所以没法和别人的研究进行对比
R中有很多做meta-analysis的包。你可以从网上学习。我对nextbio并不了解，所以帮不到你了。
“The data were normalised using quantile normalisation with Bioconductor (v.2.4) Lumi R-package”，处理后的最终的value代表什么？两个标本比较这个value的时候，怎么算倍数关系？是直接相除，还是相减后再算2的次方?
没有上下文，不好判断。但是通常你看到的值会是一个log2转换后的值，你需要检查一下值的范围，是不是主体落在20以内。
博主你好 我是R的初学者 也是刚读研究生 应老师的课题要求 在GEO上下载了数据GSE2508 但不是.cel格式的raw.data 询问之后下载了series.matrix.txt的文件 但之前接触的是.cel文件 这个txt文件也不知道怎么读入和处理以及通过这个数据得到差异表达基因 能否给我发一个具体详细点的脚本讲解一下整个过程 是20个样本的那个 GPL8300
我是初学者 很多不太了解 望不吝赐教
GSE2508本身没有原始数据，所以你下载后的会是分析过的数据，你得到的就是表达值了，直接使用就好了。
知道了 谢谢帮助
老师您好，我有两个问题想咨询一下~
1）在使用getGEO()之后，下面一般会显示file stored at,…这个是一个临时文件，一般会在C盘里面，如果我下载的多了之后电脑会很卡，但是清除之后，再次使用会再次下载；我想请问的是如果让它store到我想保存的文件夹里面：比如d、e盘等；我知道有一个可以先下载下来再getGEO的选项，但是当我需要的很多时手动下载就会很麻烦~求教
2）第二个问题是，挡我get GSE文件的时候，它一般会自动下载GPL来注释，但是有些GPL会很大很大，所以我想在getGEO的时候使用我已经在电脑中存在的GPL进行注释而不是再次下载GPL，请问这个怎么处理？
谢谢您在百忙之中的回复~
问题1）你可通过调整getGEO中的destdir参数来设置目的目录。
问题2）你可将getGPL设置成FALSE，然后自己对数据进行注释。
博主，你好，我想问一下要是数据来自一个公司不同的平台，要怎么分析啊？
分开来依次分析。最后再使用meta-analysis手段综合。
问题2）你可将getGPL设置成FALSE，然后自己对数据进行注释。
对的老师，我就是想知道怎么自己进行注释呢？如果说我自己下了本地下载了一个GSE，也下载了一个GPL，那么用什么语句进行注释呢？
使用getGEO(filename=”你下载的GPL文件”, AnnotGPL=TRUE)
博主，我最近在研究怎么下载GEO数据，但是R 软件是最新版本，总是提示要用3.0.0以前的版本，不知道这样的问题怎么解决
这种问题比较复杂。如果想得到答案，请把你的错误信息和sessionInfo提供出来。
博主 您好，有一个问题想请教一下，从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)
直接进行分析可以吗？还是需要进行其他的处理？
非常希望得到您的回复。
问题一，不是。你是从CEL文件计算来的eset，并非GEO网站上的表达文本文件。
问题二，你这是做heatmap，而不是二维聚类，二维聚类建议你再读一下其它的博文。做heatmap也不是把所有的值都拿来做图，那样没有太大的意义。
博主 你好 我现在正在学习bioconductor ，我想问一下，从GEO里面下的基因表达谱能整理成数据么，我查了一下能整理的就是cel文件没有找到将表达谱整理成数据的方法的人
不太明白您的意思。
博主 你好 GEO2R中
if (length(gset) > 1) idx <- grep("GPL2507", attr(gset, "names")) else idx <- 1
这里产生idx的目的是？谢谢！
你拿到的gset是一个list，它有可能有多组数据。这里就是猜测应该使用哪一组数据。
你拿到的gset是一个list，它有可能有多组数据。这里就是猜测应该使用哪一组数据。
能否举个多组数据的例子？谢谢！
# 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的目的是什么？谢谢！是否对表达值进行了对数变换，谢谢！
这里就是要猜一下，表达值是否进行了对数变换。如果没有，那就变换一下。这一步不一定是正确的，只是对于大多数来说，应该可以猜个差不多。
博主 您好，您的视频中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”
请您帮忙看一下是什么问题，打扰了！
你从本机上直接访问一下http://www.ncbi.nlm.nih.gov/projects/geo/query/acc.cgi?targ=self&form=xml&view=full&acc=GSE44408看看会得到什么。理论上应该可以访问到一个xml文件。
欧博，你好，看了你的视频觉得比较容易理解。想问一下，你还有哪些地方上传了教学视频吗？
有心无力。实在是事太多，无法集中精力来做视频。
老师，您好！最近用R代码下载GSE数据
getGEOSuppFiles(“GSE3325”, makeDirectory = TRUE, baseDir = getwd())
运行时，下载了一小小部分就自动停止了，请问是网速的问题还是编码或设置哪里出问题了呢？
请检查一下网络，硬盘空间，防火墙，以及长城防火墙。
楼主，怎么实现利用R或其它软件实现GEO数据的下载与处理
（1）实现GEO数据的raw.tar下载
（2）提取样本和表达谱信息
（3）将探针ID对到GeneID
（4）缺失值处理
我想如果你仔细阅读我的博文的话，这些都已经有答案了。
老师，向请问一下，我下载的RAW数据是txt格式的，就是这种：GSM1276464_0007970039_raw.TXT.gz，在跑代码的时候
setwd(‘GSE52802_RAW’)
rawdata <- ReadAffy()
无法读取，这个是因为这个文件是txt格式的吗？应该怎么处理这种文件呢？
你应该先了解清楚你下载的芯片是不是Affy的，如果不是，当然不能使用针对Affy芯片的处理流程了。
谢谢老师，应该是芯片的平台不同，要用illuminate的软件包来处理。
> gset names(gset)
[1] “GSE85341_series_matrix.txt.gz”
> gset_one exprSet dim(exprSet)
[1] 96 3
> par(cex = 0.7)
> n.sample=ncol(gset_one)
> if(n.sample>40) par(cex = 0.5)
> cols boxplot(gset_one, col = cols,main=”expression value”,las=2)
Error in sort.int(x, na.last = na.last, decreasing = decreasing, …) :
‘x’ must be atomic
> data(gset_one)
Warning message:
In data(gset_one) : data set ‘gset_one’ not found
老师，还要麻烦您看一下这个问题，为什么我导入的矩阵文件不能被data()识别呢？而且也不能做箱图呢？我查过help（data），是可以识别txt文件的。而且，我看过其他的做箱图的代码，感觉也并没有对矩阵有特别的处理，很奇怪
首先，你贴上来的代码有些混乱，可能是因为这个博客平台的原因，希望可以重新编辑一下。
然后就是你的问题了。data这个函数是用来调取搜索路径中的数据。这些数据一般是软件包中已经包装好的数据。从你的命令行来看，似乎是你自己生成了一个数据叫get_one，这种情况下是不能使用data的。否则就是理解错了你的代码。
关于boxplot，的确是可以针对矩阵画图的。这里的错误出现在sort.int上。你试着把你的矩阵使用as.matrix转换成matrix，然后看一下这个matrix是什么模式的，是不是数字型。然后再plot.
嗯嗯。谢谢老师
gset <- getGEO("GSE46106", GSEMatrix =TRUE)
https://ftp.ncbi.nlm.nih.gov/geo/series/GSE46nnn/GSE46106/matrix/
Error in function (type, msg, asError = TRUE) :
error:1407742E:SSL routines:SSL23_GET_SERVER_HELLO:tlsv1 alert protocol version
升级你的PHP至7.0以上版本，然后再重装你的package.
老师，您好，我的php是7.2.2版本已经高于7.0，我尝试着重装过包，并没有什么作用还是报一样的错