转录组分析——利用GEO数据挖掘构建基因表达图谱转录组分析——利用GEO数据挖

一、前言前段时间主管让我在GEO数据库弄个表达图谱，准备接下来做各种转录组分析，现在回过头来，觉得这个东西比较简单，十分适合上手练习。思考了一下，决定写下来分享给大家。
注意了，我所说的都是实际工作中会用到的，大家好好看好好学，不过领导有要求不能完全写工作内容，所以我会有修改。
二、什么是GEO数据库 GEO数据库全称GENE EXPRESSION OMNIBUS，是由美国国立生物技术信息中心NCBI创建并维护的基因表达数据库。它创建于2000年，收录了世界各国研究机构提交的高通量基因表达数据，也就是说只要是目前已经发表的论文，论文中涉及到的基因表达检测的数据都可以通过这个数据库中找到。
我是利用NCBI中GEO DataSets查找GEO数据库资料，举个很简单的例子，我想找胃癌相关的疾病资料、研究文献，那么你可以直接搜索gastric carcinoma，当然如果这是一种疾病的话，那你可以直接搜索这种疾病的简称。

文章图片

如果你只想关注人相关的研究，请在右方选择类型。
到这里，你必须了解一些比较重要的东西，比如
GEO上有四类数据GSM, GSE, GDS, GPL
1.GSM是单个样本的实验数据
2.GDS是人工整理好的关于某个话题的GSM的集合，一个GDS中的GSM的平台是一样的
3.GSE是一个实验项目中的多个芯片实验，可能使用多个平台
4.GPL是芯片的平台，如Affymetrix， Aglent等
我说简单点，一般我们会利用GSE号在(点击上图GES会直接跳转到GEO数据库)GEO数据库中查找资料，那么GSE号就是我们记录一个实验项目（就是NCBI每一篇文章啦）的编号，如果这个文献有多个芯片平台实验，那么就会产生多个GPL，GPL就是对应你所使用的芯片平台信息。GDS用得比较少，稍微理解一下就行。
三、数据汇总这部分跟我的工作有关，只是想了解GEO数据挖掘使用方法可以跳过不看。
如果你找的数据是一个比较热门或者比较多人研究的话，你就会得到几百个或者上千个文献，如果我们要对这些文献进行汇总整理，这样做的工程量非常大。你可能需要的东西:GPL、GSE、样本量、样本类型、题目、样本处理方法，甚至文章得出什么结论。这里我用了我比较擅长的爬虫来完成这个任务，对于大批量资料汇总，能活用自己的技能真的太好不过了，我会另起一篇文章单独讲怎么爬NCBI的数据。
四、GEO数据下载为了完成GEO数据挖掘，做基因表达谱构建，我们必须要下两个文件（如图）
1.GPL探针文件
2.表达矩阵（Series Matrix File）

文章图片
有萌新可能会说，表达矩阵人家不是做好了吗，我还做啥？对的，这是别人已经处理过的表达矩阵，我们只需要再添加上Symbol就可以完成最简单的GEO数据挖掘——基因表达图谱了。
是不是很简单？
对，就是这么简单，你可以关掉不看了。
我跟你讲，如果你运气好，得到的GPL的文件上就会有对应的Symbol，如果没有，哈，那就八仙过海各显神通吧！
五、数据分析（代码块）鸽了这么久居然还有人看这篇文章，那好，我就把遇到的坑都和大家讲一下
我做数据分析遇到了很多平台，其他他们的处理方法都不大一样，我一一介绍吧
以GPL570为例，这个平台的数据最多，也是最规整的。
首先我们读取GPL文件和GSE文件

setwd("工作路径") Sys.setlocale('LC_ALL','C') #读取GPL文件 GPL_table = read.table('GPL文件地址',sep = "\t", comment.char = "#", stringsAsFactors = F, header = T, fill = TRUE, quote = "")

第一步，你们先登录到原网页中看看下载的GPL文件大小对不对，因为NCBI外网连接可能存在网络中断的情况，文件具体大小可以再这里看到，比如说54675行、合计7.7582M等

文章图片

顺便提一下，下载GPL是下面那个Download full table...按钮，如果不提供直接下载，而是加载到网页上的平台，你们直接CTRL+A复制到txt里，一样的，不过注意一定要等NCBI加载完。GSE同理，下载完数据记得都要核查数据大小。
第二步，核查GPL探针文件信息。为了完成基因表达图谱，就一定要看探针文件是否有Gene Symbol的信息，我特地拿570来说，不仅是这个平台的资料最多，而且信息也很完整，因为有一些平台只有GB_ACC，网上有方法回对GB_ACC到基因名称，这里就点一下。

文章图片
Gene Symbol 【转录组分析——利用GEO数据挖掘构建基因表达图谱】第三步，准备合并表达矩阵。细分下来的步骤就是在GPL中找出基因列和ID列。merge到GSE中，替换掉原来的探针列。这个是我做好的代码，大家改一改自己的路径就可以用了。

setwd("路径") Sys.setlocale('LC_ALL','C') #读取GPL文件 GPL_table = read.table('GPL570.txt',sep = "\t", comment.char = "#", stringsAsFactors = F, header = T, fill = TRUE, quote = "") #读取GSE文件 GSE4100 <- read.table('GSE4100_series_matrix.txt',sep = "\t", comment.char = "!", stringsAsFactors = F,header = T, fill=TRUE)#43931#表达矩阵制作（已有完整基因名称） ID_Sybmol = GPL_table[,c(1,11)] #GPL对应ID列 colnames(ID_Sybmol)[2]="Symbol" #更改名称为Symbol，主要是为了对其下面的求平均函数#合并ID与基因列 Exp = merge(ID_Sybmol,GSE4100,by.x = "ID",by.y = "ID_REF",all=T)#45782个 Exp = Exp[,-1]

到这里，就到表达矩阵最核心的地方，数据居然有4万多行，但基因实际上只有2万多个，这是咋回事？
其实你看一下数据就会发现，很多基因是重复的，而且有些是一对多的探针。
1.所谓一对多，就是一个探针文件对应N个基因，这种探针我们的方法是直接去除。
2.基因重复，我们会以求平均的方式去除过多的行。
3.空值，0值，我都会去除。
4.数据归一化（Normalization），目的是为了后去的去批次差异（batch effect），如果你不懂，那请你百度。
因此第四步，数据清洗，我上完整的数据清洗代码

#数据过滤 Exp = Exp[Exp$Symbol != "",] #45782 Exp = na.omit(Exp)#45782#去除 Exp1 = data.frame(Exp[-grep("/",Exp$"Symbol"),]) #去一对多，grep是包含的意思，-就是不包含，symbol列不包含 #42984#求平均值 meanfun <- function(x) { x1 <- data.frame(unique(x[,1])) colnames(x1) <- c("Symbol") for (i in 2:ncol(x)){ x2 <- data.frame(tapply(x[,i],x[,1],mean)) x2[,2] <- rownames(x2) colnames(x2) <- c(colnames(x)[i], "Symbol") x1 <- merge(x1,x2,by.x = "Symbol",by.y = "Symbol",all=FALSE) } return(x1) }Exp2<- meanfun(Exp1)#23520 # 查看数据 par(cex = 0.7) n.sample=ncol(Exp2[,-1]) if(n.sample>40) par(cex = 0.5) cols <- rainbow(n.sample*1.2) boxplot(Exp2[,-1], col = cols,main="expression value",las=2)write.table(Exp2,"Exp_Original.txt",row.names = F,quote = F,sep="\t")# 校正 row.names(Exp2) = Exp2[,1] Exp2 = log(Exp2[,-1])# 我使用的是Log、这里面大家用的会是log2，目的一样，但是最后出的数据集中度不一样，大家一定要注意# 测试校正数据 par(cex = 0.7) n.sample=ncol(Exp2) if(n.sample>40) par(cex = 0.5) cols <- rainbow(n.sample*1.2) boxplot(Exp2, col = cols,main="expression value",las=2)Symbol = row.names(Exp2) Exp_test = cbind(Symbol,Exp2)write.table(Exp_test,"Exp.txt",row.names = F,quote = F,sep="\t")

完了我们的表达矩阵就大功告成啦一共23520个。