从GEO数据库下载得到表达矩阵 一文就够
在第一讲我们详细介绍了GEO数据库的基础知识及规律，也了解了如何利用官方R包GEOquery来探索GEO数据库，当然，我的生信菜鸟团博客里面也从很多其它角度解析过它，欢迎大家自行搜索学习。总得来说，从GEO数据库里面得到感兴趣数据集的表达矩阵分成两类，最简单的就是直接下载作者归一化好的表达矩阵咯，比较麻烦的就是下载最原始芯片数据，然后根据不同的芯片来一一解读成表达矩阵。
首发于生信菜鸟团博客：http://www.bio-info-trainee.com/2087.html
直接下载数据集作者上传的表达矩阵
通常我们默认作者对其芯片数据处理的步骤是正确的，所以稍微掌握技巧即可下载其归一化的表达矩阵。
而且，我已经把下载GEO数据集的表达矩阵这个过程包装成了函数，如下：
downGSE <- function(studyID = "GSE1009", destdir = ".") {
library(GEOquery)
eSet <- getGEO(studyID, destdir = destdir, getGPL = F)
exprSet = exprs(eSet[[1]])
pdata = pData(eSet[[1]])
write.csv(exprSet, paste0(studyID, "_exprSet.csv"))
write.csv(pdata, paste0(studyID, "_metadata.csv"))
return(eSet)
}
参考链接：http://www.bio-info-trainee.com/1085.html
返回的 eSet 是一个对象，需要大家仔细理解，才能进行后续分析。我在生信菜鸟团博客不止一次强调过这个：详细参考https://github.com/bioconductor-china/basic/blob/master/ExpressionSet.md）
先下载原始芯片数据，再自行处理
原始芯片数据作者肯定上传到GEO啦，每个数据集的主页都可以看到各种链接，当然也可以不进入数据集的主页，自己凭空想象出下载链接，因为它们的链接都是有规律的，所以我写了函数来帮助大家获取下载链接，函数如下；
get_GSE_links <-
function(studyID = "GSE1009", down = F, destdir = "./") {
## studyID destdir ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE1nnn/GSE1009/matrix/GSE1009_series_matrix.txt.gz
## ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE1nnn/GSE1009/suppl/GSE1009_RAW.tar
## http://www.ncbi.nlm.nih.gov/geo/browse/?view=samples&mode=csv&series=1009
supp_link = paste0("ftp://ftp.ncbi.nlm.nih.gov/geo/series/", substr(studyID, 1, nchar(studyID) - 3), "nnn/",
studyID, "/suppl/", studyID, "_RAW.tar")
meta_link = paste0("http://www.ncbi.nlm.nih.gov/geo/browse/?view=samples&mode=csv&series=", substr(studyID,
4, nchar(studyID)))
matrix_link = paste0("ftp://ftp.ncbi.nlm.nih.gov/geo/series/", substr(studyID, 1, nchar(studyID) - 3), "nnn/",
studyID, "/matrix/", studyID, "_series_matrix.txt.gz")
print(paste0("The URL for raw data is : ", supp_link))
print(paste0("The URL for metadata is : ", meta_link))
print(paste0("The URL for matrix is : ", matrix_link))
}
根据GEO的ID号，我的函数会自动判断跟它相关的3个链接，大家可以自己复制粘贴到浏览器去下载，也可以用R函数来下载。
> get_GSE_links("GSE42728")
[1] "The URL for raw data is : ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE42nnn/GSE42728/suppl/GSE42728_RAW.tar"
[1] "The URL for metadata is : http://www.ncbi.nlm.nih.gov/geo/browse/?view=samples&mode=csv&series=42728"
[1] "The URL for matrix is : ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE42nnn/GSE42728/matrix/GSE42728_series_matrix.txt.gz"
download.file('ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE42nnn/GSE42872/suppl/GSE42872_RAW.tar',destfile = 'GSE42872_RAW.tar')
可以看到最后那个series_matrix文件其实也是作者处理好的表达矩阵，不过我们需要理解的是如何处理原始芯片数据。
原始的芯片数据处理主要取决于芯片平台，下面我们就一一讲解：
旧版的affymetrix芯片
用affy包，我曾经写过教程：用affy包读取affymetix的基因表达芯片数据-CEL格式数据
library(affy)
#perform mas5 normalization
affy_data = ReadAffy(celfile.path=dir_cels)
eset.mas5 = mas5(affy_data)
exprSet.nologs = exprs(eset.mas5)
exprSet = log(exprSet.nologs, 2) #transform to Log_2 if needed
library(affy)
data <- ReadAffy(celfile.path=dir_cels)
eset <- rma(data)
write.exprs(eset,file="data.txt")
下面是一些芯片平台的介绍：
新版的affymetrix芯片
比如 Human Gene 1.0 St Array ，用oligo，我曾经写过教程：用oligo包来读取affymetix的基因表达芯片数据-CEL格式数据
library(oligo)
celFiles <- list.celfiles()
affyRaw <- read.celfiles(celFiles)
library(pd.mogene.2.0.st)
## 根据芯片平台来载入芯片设计包，没办法自动选择，跟芯片探针包不一样：mogene20sttranscriptcluster.db
eset <- rma(affyRaw)
write.exprs(eset,file="data.txt")
上面是小鼠的芯片，下面是人类的芯片数据处理代码：
studyID='GSE42872';
studyID_probe=paste0(studyID,'_probe')
studyID_gene=paste0(studyID,'_gene')
R_history_data <- paste0(studyID,'.Rdata')
if ( file.exists(R_history_data)){
load( R_history_data )
}else{
library(oligo)
## the directory data should have raw CEL files
geneCELs=list.celfiles('data',listGzipped=T,full.name=T)
affyGeneFS <- read.celfiles(geneCELs)
geneCore <- rma(affyGeneFS, target = "core")
genePS <- rma(affyGeneFS, target = "probeset")
featureData(genePS) <- getNetAffx(genePS, "probeset")
featureData(geneCore) <- getNetAffx(geneCore, "transcript")
}
exprSet=exprs( genePS )
pdata=pData( genePS )
下面是一些芯片平台的介绍：
更新的芯片版本
比如HTA1.0或者HTA2.0，好处比拟转录组测序数据。它们其实有专门的界面版本的软件处理，相关R包目前还没有。
illumina相关芯片
如果illumina的，我也在生信菜鸟团写过博客介绍： 用lumi包来处理illumina的bead系列表达芯片
rm(list=ls())
library(lumi)
# setwd('G:/array/illumina-beadseed-v4/lumi_example')
# fileName <- 'Barnes_gene_profile.txt' # Not Run
## 首先是从illumina的芯片结果文件，自己用R的lumi包来获取表达矩阵。
setwd('G:/array/illumina-beadseed-v4/GSE30669')
fileName <- 'GSE30669_HEK_Sample_Probe_Profile.txt' # Not Run
x.lumi <- lumiR.batch(fileName) ##, sampleInfoFile='sampleInfo.txt')
pData(phenoData(x.lumi))
## Do all the default preprocessing in one step
lumi.N.Q <- lumiExpresso(x.lumi)
### retrieve normalized data
dataMatrix <- exprs(lumi.N.Q)
## 下面是从GEO里面下载表达矩阵
rm(list=ls())
library(GEOquery)
library(limma)
GSE30669 <- getGEO('GSE30669', destdir=".",getGPL = F)
exprSet=exprs(GSE30669[[1]])
GSE30669[[1]]
pdata=pData(GSE30669[[1]])
exprSet=exprs(GSE30669[[1]])
## 很明显可以看到前面得到的dataMatrix 和后面得到的 exprSet 都是我们想要的表达矩阵
重点就是得到表达矩阵，它封装好了一个函数，lumiExpresso可以直接处理LumiBatch对象，这个函数结合了,N,T,B,Q(normalization,transformation,backgroud correction,qulity control)四个步骤，其中Q这个步骤又包括8种统计学图片。在该包的文章有详细说明：http://bioinformatics.oxfordjournals.org/content/24/13/1547.full
而 LumiBatch 对象是通过 lumiR.batch 读取的芯片文件被Illumina Bead Studio toolkit 处理的结果，也就是通常我们从公司或者GEO下载的数据( level 3 的 process data)。
其它数据芯片
再比如agilent公司的， 已经各种各样的实验室定制化芯片，他们的原始芯片数据其实是没有现成的R包的，我们只能相信作者上传的表达矩阵是正确的，直接利用其表达矩阵做下游分析即可。
写在最后
总之，我们要得到所有检测到的基因在所有样本的表达量信息， 就是我们通常说的表达矩阵，才能继续下游分析。
插个题外话，很多朋友会有整合不同的数据集的需求，这个时候需要额外注意批次效应的问题。可以用 RUVseq和svaseq 来找到全新的未知的的因素，再加上我们已知的因素，然后有了这些因素，就可以用 ComBat 和 removeBatchEffect 这样的工具来进行去除批次效应啦。