来完成你的⽣信作业，这是最有诚意的GEO数据库教程
上次更新了⼀个GEO数据框的帖⼦，包含了GEO分析的主要组成部分，其中有些区域为了代码复⽤，不是很容易懂。现在对其进⾏注释，并增加KEGG分析的内容，⼒求能让他变成最有诚意的GEO数据分析教程，⽂以下是正⽂： GEO数据库包罗万象，对于每个领域的科研⼯作者很有帮助。
我最开始分析GEO芯⽚主要使⽤别⼈的代码，跑流程。但是绝对不能出现报错，⼀报错就束⼿⽆策。
遇到⽐较困难的地⽅就使⽤别⼈的⼩⼯具来代替。⽐如，ID转换，GO分析等。直到去年我才能够完全使⽤R语⾔实现。
这两年，⽹上的教程帖⼦井喷⼀样涌现，⼤有教师超过学员的趋势。
但是，GEO的分析跟TCGA不⼀样，⾥⾯要考虑的细节⽐较多。
⽐如我们今天要解决以下⼏个问题：
1. GEO数据如何⽅便地下载？
3. 什么时候做log转换？
4. 芯⽚探针如何转换？
6. 配对样本的对⽐矩阵如何构建？
7. 配对样本如何作图展⽰？
8. 热图⽕⼭图究竟在什么？
9. GO分析，KEGG分析怎么做，如何解释？
看⽂献得到GSE号，这是数据的名⽚，GSE32575，
⽂献中提及的GSE号
GEO官⽹
点击进去后，会有这个芯⽚的介绍
⽐如，这个芯⽚⼀开始出现的时候是⼲什么的？
summary
认真读⼀读就知道作者为什么要做这个芯⽚。
如果还不是很清楚，可以往下看，可以看到实验设计，甚⾄还有这个芯⽚当时发表的⽂章。下载下来读⼀读，可以看看当时这个芯⽚的哪些信息被使⽤了，还有哪些没被使⽤，我们可以提出新的科研假设，⽤别⼈的信⽤别⼈的数据⽀持⾃⼰的科研假设。 在往下看还有平台信息GPL6102，这个信息可以帮助我们来注释⽂件。
图中的48表⽰这个芯⽚做了48个样本，每个样本都写明是什么样本。
假如点击Analyze with GEO2R按钮，就可以按照以下帖⼦介绍的⽅法，⽆代码分析这个芯⽚。
⽆代码GEO芯⽚分析图⽂教程
但是今天，我们要⽤更加⽅便快捷的⽅法，假如你有⼀点耐⼼，完全可以复制这篇帖⼦的每⼀步操作，
最后得到能发表的图⽚。
⿊⾊区域是代码块，可以⼀整块⼀整块复制到Rstudio的脚本中，然后⼀⾏⾏运⾏，代码框中的#号表⽰注释，也就是这⼀⾏的内容不会被当成代码执⾏。
关于如何配置⾃⼰的R语⾔环境，可以看这篇帖⼦。
学习R语⾔，从这⼀课开始
捣⿎好了之后，我们就可以往下操作了，⽽当你完成这个帖⼦的那⼀刻，你也会⼀脚踏进数据挖掘的⼤门，⾄少说，见了世⾯，以后不会被别⼈那些云⼭雾罩的话给忽悠。
1.加载R包获取数据。
## 加载R包
library(GEOquery)
## 下载数据，如果⽂件夹中有会直接读⼊
gset = getGEO('GSE32575', destdir='.',getGPL = F)
## 获取ExpressionSet对象，包括的表达矩阵和分组信息
gset=gset[[1]]
以后只要更改GSE号，就可以直接下载别的GEO数据，很⽅便。
2.通过pData函数获取分组信息
## 获取分组信息，点开查阅信息
pdata=pData(gset)
## 只要后36个，本次选择的这36个是配对的。
## 所以，别⼈的芯⽚我们也不是全要，我们只要适合⾃⼰的数据
group_list=c(rep('before',18),rep('after',18))
group_list=factor(group_list)
## 强制限定顺序
group_list <- relevel(group_list, ref='before')
3.通过exprs函数获取表达矩阵并校正
exprSet=exprs(gset)
可以先简单看⼀下整体样本的表达情况，其中exprSet[,-c(1:12)]的意思是去掉这个数据的前12个，因为我们需要的后⾯36个
boxplot(exprSet[,-c(1:12)],outline=FALSE, notch=T,col=group_list, las=2)
未校正前
需要⼈⼯校正⼀下，⽤的⽅法类似于Quntile Normalization
这⾥我们⽤limma包内置的⼀个函数
library(limma)
exprSet=normalizeBetweenArrays(exprSet)
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)
校正后
4.判断是否需要进⾏数据转换
我们通过分组信息已经知道，前12个样本本次不需要，所以先去掉
exprSet = as.data.frame(exprSet)[,-seq(1,12)]
打开现在的表达矩阵，是这个样⼦的
表达矩阵局部
⾁眼看到数字很⼤，这时候需要log转换(选log2)。
此外还有⼀个脚本可以⾃动判断是否需要log转换，这个脚本来⾃于GEO2R，我改写了⼀点点。
ex <- exprSet
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))LogC <- (qx[5] > 100) ||
(qx[6]-qx[1] > 50 && qx[2] > 0) ||
(qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
if (LogC) { ex[which(ex <= 0)] <- NaN
exprSet <- log2(ex)
print('log2 transform finished')}else{print('log2 transform not needed')}
这个语句会⾃动判断，如果已经log话就跳过，并提⽰
'log2 transform not needed'
如果没有log话，他⾃动log2，并且提⽰：
'log2 transform finished'
这时候再来看这个数据就规则多了
log2转换后
但是我们发现⼀个问题，这个表达数据⾏名是探针，不是我们熟悉的基因名称如果TP53，BRCA1这样的，所以我们需要注释他。
5.获取注释信息
这步会很顺利，也会很难，是技术上的限速环节。
通常情况下我们使⽤R包来注释，R包和平台的注释信息对应关系我整理了⼀个表格 platformMap，是⼀个⽂本⽂件，在⽂末有获取⽅法。
platformMap <- data.table::fread('')
我们根据上述帖⼦⾥⾯提到的platformMap，到对应的R包是：
'illuminaHumanv2.db'
这⾥⾯的GPL6102是平台名称，我们在⼀开始已经在⽹页中知道了，如果忘记了，不要紧，这⾏代码可以⾃动获取
index = gset@annotation
我们也发现，表格中对应的 'illuminaHumanv2'，跟我们需要的R包 'illuminaHumanv2.db'，中间相差个“.db”，如果是单个，就⾃⼰加上，我们也可以使⽤代码的⽅法，⾃动获取：index = gset@annotation
platformDB = paste0(platformMap$bioc_package[grep(index,platformMap$gpl)],'.db')
安装R包并加载
## 第⼀句是为了增加镜像
if(length(getOption('BioC_mirror'))==0) options(BioC_mirror='mirrors./bioc/')
## 没有这个包就给你装，有就不装
if(!require('illuminaHumanv2.db')) BiocManager::install('illuminaHumanv2.db',update = F,ask = F)
## 加载R包
library(illuminaHumanv2.db)
获取探针对应的symbol信息
## 获取表达矩阵的⾏名，就是探针名称
probeset <- rownames(exprSet)
## 使⽤lookup函数，到探针在illuminaHumanv2.db中的对应基因名称
## 如果分析别的芯⽚数据，把illuminaHumanv2.db更换即可
SYMBOL <- annotate::lookUp(probeset,'illuminaHumanv2.db', 'SYMBOL')
## 转换为向量
symbol = as.vector(unlist(SYMBOL))
制作probe2symbol转换⽂件
probe2symbol = data.frame(probeset,symbol,stringsAsFactors = F)
probe2symbol部分结果
有了这个⽂件，我们就可以进⾏探针转换了
有些平台是没有对应的R包的，可以⾃⼰下载平台信息来做
skr！GEO芯⽚数据的探针ID转换
有些GEO平台的探针转换⽐较⿇烦
正则表达式是我们认识世界的哲学
6.探针转换与基因去重
⼀个基因会对应对个探针，有些基因名称会是重复的，这些都需要处理。
对于，多个探针，我们选取在样本中平均表达量最⾼的探针作为对应基因的表达量。⼀下代码完成所有事情，⽽且可以复⽤。
library(dplyr)
library(tibble)
exprSet <- exprSet %>%
rownames_to_column(var='probeset') %>%
#合并探针的信息
inner_join(probe2symbol,by='probeset') %>%
#去掉多余信息
select(-probeset) %>%
#重新排列
select(symbol,everything()) %>%
select(symbol,everything()) %>%
#求出平均数(这边的点号代表上⼀步产出的数据)
mutate(rowMean =rowMeans(.[grep('GSM', names(.))])) %>%
#去除symbol中的NA
filter(symbol != 'NA') %>%
#把表达量的平均值按从⼤到⼩排序
arrange(desc(rowMean)) %>%
# symbol留下第⼀个
distinct(symbol,.keep_all = T) %>%
#反向选择去除rowMean这⼀列
select(-rowMean) %>%
# 列名变成⾏名
column_to_rownames(var = 'symbol')
现在数据变成这个样⼦
探针去重与转换
⽽且，⾏数从原来的48701变成18962⾏。
7.差异分析
这⾥是limma包的核⼼环节，也是理解以后其他差异分析⼯具如Deseq2的基础。
差异分析的矩阵构建有两种⽅式
我们选取格式⽐较简单的。如果没有配对信息，差异分析这样做：
design=model.matrix(~ group_list)
fit=lmFit(exprSet,design)
fit=eBayes(fit)
allDiff=topTable(fit,adjust='fdr',coef='group_listafter',number=Inf,p.value=0.05)
其中allDiff得到6799⾏。
因为我们这⾥实际上是有配对信息的，差异分析应该这样做：
pairinfo = factor(rep(1:18,2))
design=model.matrix(~ pairinfo+group_list)
fit=lmFit(exprSet,design)
fit=eBayes(fit)
allDiff_pair=topTable(fit,adjust='BH',coef='group_listafter',number=Inf,p.value=0.05)
得到的差异分析结果allDiff_pair有7650个，⽐直接做差异分析要多接近1000个。
这⾥配对和⾮配对的差异分析，只相差⼀步，也就是是否有配对信息
配对和⾮配对⽐较
8.作图验证(⾮必要)
差异分析的结果，配对和⾮配对理解起来⽐较抽象，我们⽤图⽚直观地感受⼀下。
⾸先我们需要得到ggplot2喜欢的数据格式，⾏是观测，列是变量，这也叫clean data，tide data, 清洁数据。⽽我们学习R语⾔的过程中，⼤部分时间都是在调整格式，这也是线下课上特别强调的⼀点。⾸先基因名称需要在列，所以需要⽤t函数，⾏列转置。
data_plot = as.data.frame(t(exprSet))
data_plot = data.frame(pairinfo=rep(1:18,2),
group=group_list,
data_plot,stringsAsFactors = F)
data_plot局部
作图展⽰:
挑选了⼀个基因CAMKK2
library(ggplot2)
ggplot(data_plot, aes(group,CAMKK2,fill=group)) +
geom_jitter(aes(colour=group), size=2, alpha=0.5)+
xlab('') +
ylab(paste('Expression of ','CAMKK2'))+
theme_classic()+
theme(legend.position = 'none')
⽆配对信息
看起来杂乱⼀⽚，根本得不到有⽤信息，我们加⼊他的配对信息，再来作图：
library(ggplot2)
ggplot(data_plot, aes(group,CAMKK2,fill=group)) +
geom_boxplot() +
geom_point(size=2, alpha=0.5) +
geom_line(aes(group=pairinfo), colour='black', linetype='11') +
xlab('') +
ylab(paste('Expression of ','CAMKK2'))+
theme_classic()+
theme(legend.position = 'none')
加⼊配对信息
突然间飞到枝头变凤凰。⽽这些基因就是普通差异分析会遗漏的部分，配对差异分析的时候他们⼜变得有意义。
再来看看真正有差异的基因吧，⽐如：
library(ggplot2)
ggplot(data_plot, aes(group,CH25H,fill=group)) +
geom_boxplot() +
geom_point(size=2, alpha=0.5) +
geom_line(aes(group=pairinfo), colour='black', linetype='11') +
xlab('') +
ylab(paste('Expression of ','CH25H'))+
theme_classic()+
theme(legend.position = 'none')
CH25H
我们还可以批量地作图，这段不必要，可以⾃⼰⼀张张画，我只是展⽰⼀下如何批量画出差异最⼤的8个基因：
library(dplyr)
library(tibble)
allDiff_arrange <- allDiff_pair %>%
rownames_to_column(var='genesymbol') %>%
值越⼩，倍
明。这是需要背景知识的。
假如我们对凋亡感兴趣，我们还可以看看我们有多少基因被富集在了凋亡通路，Y叔的神包clusterprofiler也是可以做的。⽽且做的很云淡风轻。
⾸先我们把富集的结果变成数据框，便于⾃⼰查看。
test <- data.frame(EGG)
我们发现想看的凋亡通路牌号是，hsa04210，那我们就⽤下⾯的代码浏览⼀下，会跳出⼀个⽹页。
browseKEGG(EGG, 'hsa04210')
apoptosis通路
其中标红的就是富集到的通路。
那么这时候我们就会想，差异基因是有⾼有低的，⽽通路富集的时候缺失了这⼀个信息。⼀种解决⽅案是，分别⽤⾼表达和低表达基因去做这个分析。我觉得信息是冗余和偏误的，因为任何⼀个通路的
激活，都是伴随另外⼀种⽅法是，先总体分析，然后再标注出上调和下调的基因。那么Y叔的神包也可以做的。
好像上了⼀点档次，留给⼤家⾃⼰去实现。
那么我们这个GEO的分析就讲完了。来看看我们获得了哪些图⽚呢？
只要是参加过线下课的同学，或者只要R语⾔基本能⼒过关的同学，轻⽽易举就能重复这个帖⼦的所有内容。
即使是零基础，安装这个帖⼦准备好⼯作环境，也能完成。
学习R语⾔，从这⼀课开始
如果继续做GESA，基因富集分析
也是可以的，我⾃⼰觉得这个分析不需要⼈为设定差异基因的阈值，⽐如我们认为⼤于2倍就叫差异基因，这是不科学的，GSEA分析更加靠谱，可以代替GO和KEGG分析。
GSEA汇总图
正向调控
负向调控
假如我们通过以上的⽅法最终锁定了⼀两个基因，那么我们还可以对这个基因进⾏批量的相关性分析。
单基因批量相关性分析的妙⽤
那么现在的问题是，技术上实现了，理论上该如何整合。
R语⾔充其量只是个⼯具，神器不神，主要看⼈。假如你现在已经读过很多⽂献，做过很多实验，R语⾔⼀定有⽤，因为你⾸先想的会是⽤R语⾔解决⾃⼰的科学问题。
假如你现在还没有怎么读过⽂献，先不要急着要学习R语⾔。经常有学⽣加我，上来就说我要学习⽣信，我⼀听⼼⾥都很慌张，因为我知道，就我掌握的东西离⽣信还离得远，我就是从科研⼈员的⾓度功利地去运不分类别地去多看10分左右的⽂献可能是更好的策略，再⾼分的就⾼屋建瓴谈哲学对初学者不友好，只能学其⾏不能学其神，低分的容易把初学者带⼊科研⽪条客的⼤坑，10分左右的⽂献，起⼿不⾼，论证严谨，⼤
有时候我会偏执地认为，读⽂献是性价⽐最⾼的提升科学思维的⽅法，⼀个科研⼈员的⽔平⾼低可以简单地⽤⽂献阅读量来⽐较。然后在回过头来看看这个帖⼦凝练⾃⼰的科学问题。有了科学问题我们就有了灵魂，才课题设计：收不完的病⼈查不完的房，临床医⽣如何快速地设计⼀个靠谱的课题？
下周，我会更新⼀个长长的TCGA的帖⼦，这样我的2018年就这么过去了。
因为考虑到有些朋友是新⼿，所以我还录制了⼀个友好的导学视频，教⼤家如何把帖⼦⾥⾯的内容化为⼰有，连同R语⾔安装教程还有platformMap⽂件⼀起打包。有需要的朋友，加我guotosky，简单介绍⼀
不需要任何赞赏，如果觉得这个帖⼦好看，就在右下⾓点击好看，这是的新功能。
祝⼤家新年快乐。
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议