从GEO数据库下载得到表达矩阵一文就够
在第一讲我们详细介绍了GEO数据库的基础知识及规律,也了解了如何利用官方R包GEOquery
来探索GEO数据库,当然,我的生信菜鸟团博客里面也从很多其它角度解析过它,欢迎大家自行搜索学习。总得来说,从GEO数据库里面得到感兴趣数据集的表达矩阵分成两类,最简单的就是直接下载作者归一化好的表达矩阵咯,比较麻烦的就是下载最原始芯片数据,然后根据不同的芯片来一一解读成表达矩阵。
解读GEO数据存放规律及下载,一文就够
解读SRA数据库规律一文就够
通常我们默认作者对其芯片数据处理的步骤是正确的,所以稍微掌握技巧即可下载其归一化的表达矩阵。
而且,我已经把下载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文件其实也是作者处理好的表达矩阵,不过我们需要理解的是如何处理原始芯片数据。
原始的芯片数据处理主要取决于芯片平台,下面我们就一一讲解:
用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")
下面是一些芯片平台的介绍: