GEO数据库是NCBI开发的基因表达数据库,主要接收通过高通量测序、基因芯片等方法获得的基因表达数据——这就方便大家利用他人数据发文章了。
而进行GEO数据挖掘的第一步就是进行数据下载,但是进入网站点点来下载,查询搜索工作就不少,下载下来的数据还不一定能看懂。有没有什么方法可以解决这个问题呢?当然有——R包GEOquery!下面就针对芯片数据,教大家用GEOquery包完成下载工作。
GEO数据
在下载之前要先了解GEO数据库具体存放的四类数据:GSE、GDS、GSM、和GPL。
一个GSE Accession对应的是整个研究项目的系列的数据,可能涉及不同平台;
一个GDS Accession对应的一个同一平台的数据集;
一个GSM Accession对应单一样品的数据信息,它只能是单一平台的数据,往往,GSE 和GDS中会包含多个GSM数据;
一个GPL Accession,则对应一个platform信息。
R包安装与加载
GEOquery
## try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
biocLite("GEOquery")
Biobase
## try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
biocLite("Biobase")
正确加载
library('Biobase')
library('GEOquery')
setwd("F:/GEO") ############有需要可以设置路径
利用GSE Accession
通过阅读文献查找感兴趣的GSE Accession,下载对应的表达数据和平台信息等,可以利用GEOquery中的getGEO()函数下载series_matrix.txt。例如GSE70213:
> gse = getGEO("GSE70213", GSEMatrix =TRUE, destdir = ".", getGPL = T, AnnotGPL = T)
###destdir设置当前目录,getGPL 和AnnotGPL都设置TRUE,可以下载和获得平台的注释文件
gse为列表数据,对应的GSM是单平台,则length为1,之后分别利用Biobase包中的exprs()、pData()和fData()获得表达数据、样品处理分组等信息、芯片平台的设计注释信息,还可以利用annotation()函数了解对应GPL Accession,譬如exprs()函数:
> exprSet=exprs(gse[[1]])
> head(exprSet,2)
GSM1720833 GSM1720834 GSM1720835 GSM1720836 GSM1720837 GSM1720838 GSM1720839 GSM1720840 GSM1720841 GSM1720842
10338001 2041.40800 2200.86100 2323.7600 3216.26300 2362.77500 2195.31800 2013.35900 2146.25800 1785.9460 2067.04100
10338002 63.78059 65.08438 58.3082 75.86145 66.95605 43.81526 49.11361 51.29279 48.9604 42.14286
GSM1720843 GSM1720844 GSM1720845 GSM1720846 GSM1720847 GSM1720848 GSM1720849 GSM1720850 GSM1720851 GSM1720852
10338001 1769.1150 1720.77400 1847.42900 2214.69800 2279.51500 2530.45600 2303.26400 2358.83400 1701.40000 1970.92400
10338002 42.5472 43.48373 64.34628 59.75188 57.48852 60.26423 54.81179 53.70885 57.86877 57.02808
GSM1720853 GSM1720854 GSM1720855 GSM1720856
10338001 1822.78600 2014.26000 1737.84200 2001.73400
10338002 59.26121 55.27306 54.36722 49.43959
注释信息的获取可以进行探针和基因的对应,方便后续分析。通过exprs()、pData()和fData()获得的数据都可以利用write.table等进行文件保存。
利用GDS Accession
GDS数据同样可以利用getGEO()函数下载soft文件。例如GDS5881:
> gds = getGEO("GDS5881", GSEMatrix =TRUE, destdir = ".", getGPL = T, AnnotGPL = T)
###destdir设置当前目录
gds可以利用GEOquery包中的Table()获取表达数据,并利用Meta()获得描述信息,其中Meta(gds)$platform可以获得GPL Accession。
> exprSet=Table(gds)
> head(exprSet,1)
ID_REF IDENTIFIER GSM1720845 GSM1720846 GSM1720847 GSM1720848 GSM1720849 GSM1720850 GSM1720851 GSM1720852 GSM1720853
1 10344614 Gm2889 48.4971 47.252 39.3331 49.9048 36.8313 41.9501 37.5569 38.1924 46.0668
GSM1720854 GSM1720855 GSM1720856
1 34.689 38.5762 32.2618
> Meta(gset)$platform
[1] "GPL6246"
针对getGEO返回的gds——GDS数据,可以利用GEOquery包中GDS2Set()和GDS2MA()转变为为ExpressionSets 和limma MALists。
> gds2eSet=GDS2eSet(gds)
> MA=GDS2MA(gds)
再针对返回的gds2eSet,利用exprs()、pData()和fData()同样可以获得表达数据、样品处理分组信息、芯片平台的设计注释信息。返回的MA中涉及大量的描述信息,其中MA$tragets也是样品处理信息。
利用GSM Accession
利用GSM Accession下载的是单样本的表达数据,例如GSM1720833:
> gsm = getGEO("GSM1720833", GSEMatrix =TRUE, destdir = ".", getGPL = T, AnnotGPL = T)
针对gsm,同样是利用GEOquery包中的Table()获取表达数据,并利用Meta()获得描述信息,而获取对应的GSE Accession 和GPL Accsesion利用Meta(gsm)$series_id和Meta(gsm)$platform_id。
利用GPL Accession
针对芯片平台,利用GPL Accession下载得到的数据是芯片的设计和注释信息,可以获得探针组和基因的对应关系,利用Table()函数可以显示annotation,例如GPL6246:
> gpl = getGEO("GPL6246", GSEMatrix =TRUE, destdir = ".", getGPL = T, AnnotGPL = T)
> ann=Table(gpl)
> head(ann,2)
ID Gene title Gene symbol Gene ID UniGene title UniGene symbol UniGene ID
1 10344614 predicted gene 2889 Gm2889 100040658
2 10344616
Nucleotide Title
1 Mus musculus blastocyst blastocyst cDNA, RIKEN full-length enriched library, clone:I1C0009C06 product:hypothetical DeoxyUTP pyrophosphatase/Aspartyl protease, retroviral-type family profile/Retrovirus capsid, C-terminal/Peptidase aspartic/Peptidase aspartic, active site containing protein, full insert sequence///Mus musculus blastocyst blastocyst cDNA, RIKEN full-length enriched library, clone:I1C0042P10 product:hypothetical protein, full insert sequence
2
GI GenBank Accession Platform_CLONEID Platform_ORF Platform_SPOTID Chromosome location
1 74211482///74217103 AK145513///AK145782 chr1:3054233-3054733 18
2 chr1:3102016-3102125
Chromosome annotation GO:Function GO:Process GO:Component GO:Function ID GO:Process ID GO:Component ID
1 Chromosome 18
2
总结
以上数据,最终都可以通过下载整理之后保存成后续分析的文件,相比之手动下载,利用R包会更便捷快速。大家也可以自行整理代码批量下载数据哦!
1. 文章越来越难发?是你没发现新思路,基因家族分析发2-4分文章简单快速,学习链接:基因家族分析实操课程、基因家族文献思路解读
2. 转录组数据理解不深入?图表看不懂?点击链接学习深入解读数据结果文件,学习链接:转录组(有参)结果解读;转录组(无参)结果解读
3. 转录组数据深入挖掘技能-WGCNA,提升你的文章档次,学习链接:WGCNA-加权基因共表达网络分析
4. 转录组数据怎么挖掘?学习链接:转录组标准分析后的数据挖掘、转录组文献解读
5. 微生物16S/ITS/18S分析原理及结果解读、OTU网络图绘制、cytoscape与网络图绘制课程
6. 生物信息入门到精通必修基础课:linux系统使用、perl入门到精通、perl语言高级、R语言入门、R语言画图
7. 医学相关数据挖掘课程,不用做实验也能发文章:TCGA-差异基因分析、GEO芯片数据挖掘、GEO芯片数据标准化、GSEA富集分析课程、TCGA临床数据生存分析、TCGA-转录因子分析、TCGA-ceRNA调控网络分析
8.其他,二代测序转录组数据自主分析、NCBI数据上传、二代测序数据解读
如果觉得我的文章对您有用,请随意打赏。你的支持将鼓励我继续创作!