上次更新了一个GEO数据框的帖子，包含了GEO分析的主要组成部分，其中有些区域为了代码复用，不是很容易懂。现在对其进行注释，并增加KEGG分析的内容，力求能让他变成最有诚意的GEO数据分析教程，文末还录制了一个友好的导学视频，希望物尽其用。
GEO数据库包罗万象，对于每个领域的科研工作者很有帮助。
我最开始分析GEO芯片主要使用别人的代码，跑流程。但是绝对不能出现报错，一报错就束手无策。
遇到比较困难的地方就使用别人的小工具来代替。比如，ID转换，GO分析等。直到去年我才能够完全使用R语言实现。
这两年，网上的教程帖子井喷一样涌现，大有教师超过学员的趋势。
但是，GEO的分析跟TCGA不一样，里面要考虑的细节比较多。
比如我们今天要解决以下几个问题：
GEO数据如何方便地下载？
数据什么时候需要标准化？
什么时候做log转换？
芯片探针如何转换？
差异基因如何获得？
配对样本的对比矩阵如何构建？
配对样本如何作图展示？
热图火山图究竟在什么？
GO分析，KEGG分析怎么做，如何解释？
看文献得到GSE号，这是数据的名片，GSE32575，
在这个网址输入这个号我们可以知道这个芯片的信息
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi
点击进去后，会有这个芯片的介绍
比如，这个芯片一开始出现的时候是干什么的？
用别人的数据支持自己的科研假设。
在往下看还有平台信息GPL6102，这个信息可以帮助我们来注释文件。
Analyze with GEO2R按钮，就可以按照以下帖子介绍的方法，无代码分析这个芯片。黑色区域是代码块，可以一整块一整块复制到Rstudio的脚本中，然后一行行运行，代码框中的#号表示注释，也就是这一行的内容不会被当成代码执行。
关于如何配置自己的R语言环境，可以看这篇帖子。
学习R语言，从这一课开始
捣鼓好了之后，我们就可以往下操作了，而当你完成这个帖子的那一刻，你也会一脚踏进数据挖掘的大门，至少说，见了世面，以后不会被别人那些云山雾罩的话给忽悠。
## 加载R包
library(GEOquery)
## 下载数据，如果文件夹中有会直接读入
gset = getGEO('GSE32575', destdir='.',getGPL = F)
## 获取ExpressionSet对象，包括的表达矩阵和分组信息
gset=gset[[1]]
以后只要更改GSE号，就可以直接下载别的GEO数据，很方便。
## 获取分组信息,点开查阅信息
pdata=pData(gset)
## 只要后36个,本次选择的这36个是配对的。
## 所以，别人的芯片我们也不是全要，我们只要适合自己的数据
group_list=c(rep('before',18),rep('after',18))
group_list=factor(group_list)
## 强制限定顺序
group_list <>'before')
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)
我们通过分组信息已经知道，前12个样本本次不需要，所以先去掉
exprSet = as.data.frame(exprSet)[,-seq(1,12)]
打开现在的表达矩阵，是这个样子的
ex <>
qx <>
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)] <>
exprSet <>
print('log2 transform finished')}else{print('log2 transform not needed')}
这个语句会自动判断，如果已经log话就跳过，并提示
'log2 transform not needed'
如果没有log话，他自动log2，并且提示：
'log2 transform finished'
这时候再来看这个数据就规则多了
这步会很顺利，也会很难，是技术上的限速环节。
通常情况下我们使用R包来注释，R包和平台的注释信息对应关系我整理了一个表格 platformMap，是一个文本文件，在文末有获取方法。
platformMap <>'platformMap.txt')
我们根据上述帖子里面提到的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='https://mirrors.ustc.edu.cn/bioc/')
## 没有这个包就给你装，有就不装
if(!require('illuminaHumanv2.db')) BiocManager::install('illuminaHumanv2.db',update = F,ask = F)
## 加载R包
library(illuminaHumanv2.db)
获取探针对应的symbol信息
## 获取表达矩阵的行名，就是探针名称
probeset <>
## 使用lookup函数，找到探针在illuminaHumanv2.db中的对应基因名称
## 如果分析别的芯片数据，把illuminaHumanv2.db更换即可
SYMBOL <>'illuminaHumanv2.db', 'SYMBOL')
## 转换为向量
symbol = as.vector(unlist(SYMBOL))
制作probe2symbol转换文件
probe2symbol = data.frame(probeset,symbol,stringsAsFactors = F)
一个基因会对应对个探针，有些基因名称会是重复的，这些都需要处理。
对于，多个探针，我们选取在样本中平均表达量最高的探针作为对应基因的表达量。一下代码完成所有事情，而且可以复用。
library(dplyr)
library(tibble)
exprSet <- exprset %>%
rownames_to_column(var='probeset') %>%
#合并探针的信息
inner_join(probe2symbol,by='probeset') %>%
#去掉多余信息
select(-probeset) %>%
#重新排列
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')
现在数据变成这个样子
这里是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个。
这里配对和非配对的差异分析，只相差一步，也就是是否有配对信息
差异分析的结果，配对和非配对理解起来比较抽象，我们用图片直观地感受一下。
首先我们需要得到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)
作图展示:
挑选了一个基因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')
我们还可以批量地作图，这段不必要，可以自己一张张画，我只是展示一下如何批量画出差异最大的8个基因：
library(dplyr)
library(tibble)
allDiff_arrange <- alldiff_pair %>%
rownames_to_column(var='genesymbol') %>%
arrange(desc(abs(logFC)))
genes <>1:8]
plotlist <>function(x){
data =data.frame(data_plot[,c('pairinfo','group')],gene=data_plot[,x])
ggplot(data, aes(group,gene,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 ',x))+
theme_classic()+
theme(legend.position = 'none')
})
library(cowplot)
plot_grid(plotlist=plotlist, ncol=4,labels = LETTERS[1:8])
而配对样本的展示，我们用TCGA的数据演示过很多次
找出TCGA中的配对样本并正确展示数据
配对样本数据如何计算及展示？
差异分析后还有很多操作，
画一张热图
library(pheatmap)
## 设定差异基因阈值，减少差异基因用于提取表达矩阵
allDiff_pair=topTable(fit,adjust='BH',coef='group_listafter',number=Inf,p.value=0.05,lfc =0.5)
##提前部分数据用作热图绘制
heatdata <>
##制作一个分组信息用于注释
annotation_col <>
rownames(annotation_col) <>
#如果注释出界，可以通过调整格子比例和字体修正
pheatmap(heatdata, #热图的数据
cluster_rows = TRUE,#行聚类
cluster_cols = TRUE,#列聚类，可以看出样本之间的区分度
annotation_col =annotation_col, #标注样本分类
annotation_legend=TRUE, # 显示注释
show_rownames = F,# 显示行名
show_colnames = F,# 显示列名
scale = 'row', #以行来标准化，这个功能很不错
color =colorRampPalette(c('blue', 'white','red'))(100))
画热图的意义在哪？
第一看样本质量：本来before和after两组应该完全分开的，但是热图里面after有两个样本跟bfefore分不开，要考虑是不是测量失误，还是本身样本就特殊
我们还可以画一个火山图
library(ggplot2)
library(ggrepel)
library(dplyr)
libr
data - topTable(fit,adjust='BH',coef='group_listafter',number=Inf)
data$significant - as.factor(data$adj.P.Val0.05 & abs(data$logFC) > 0.5)
data$gene - rownames(data)
ggplot(data=data, aes(x=logFC, y =-log10(adj.P.Val),color=significant)) +
geom_point(alpha=0.8, size=1.2,col='black')+
geom_point(data=subset(data, logFC > 0.5),alpha=0.8, size=1.2,col='red')+
geom_point(data=subset(data, logFC <>-0.5),alpha=0.6, size=1.2,col='blue')+
labs(x='log2 (fold change)',y='-log10 (adj.P.Val)')+
theme(plot.title = element_text(hjust = 0.4))+
geom_hline(yintercept = -log10(0.05),lty=4,lwd=0.6,alpha=0.8)+
geom_vline(xintercept = c(0.5,-0.5),lty=4,lwd=0.6,alpha=0.8)+
theme_bw()+
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = 'black')) +
geom_point(data=subset(data, abs(logFC) > 1),alpha=0.8, size=3,col='green')+
geom_text_repel(data=subset(data, abs(logFC) > 1),
aes(label=gene),col='black',alpha = 0.8)
火山图画起来简单，难的是如何定制化展示。比如在图上有不同的颜色，不同的点来表示数据。有什么意义呢？我其实理解的也不是很透彻，但是考虑到他的横坐标是变化倍数，纵坐标是p值的负数，那么p值越小，倍数越大的基因就会出现在左右上方。
数据是死的，而图形化展示是活的，火山图可能是大家觉得比较好的展示差异基因的方法。从这个图中我们还可以看出，两边是否对称，比如有的药物就是抑制基因转录的，那么加了药物的两组，可能会出现，上调和下调基因不对称的情形，这个从数据中可以分析出来，但是用图就一目了然了。
此外还有耳熟能详的GO分析，KEGG分析，我觉得都是名气很大用的很少的聚类方法。
做这个的意义在于，我们获得了几千个差异基因，但是我们确定最重要的基因呢？因人而异，不同的课题组，想的不一样。但是大体的思路是聚类，假如P53通路有100个基因，对于人类20000个基因而言，随便抓一把基因，出现P53通路基因的概率是100/20000,是0.005，现在3000个差异基因中就出现了50个P53通路的基因，这个概率是0.016，明显超过了他随机出现的概率，我们就可以说，p53通路被富集了，这个通路可能在你研究的课题中有重要作用。
Y叔的神包clusterprofiler可以做出特别好看的图。
比如GO分析，有三个组成部分，分别是CC，cellular compartment 细胞组分，BP，biological process，生物进程，MF，molecular function，分子功能。
## 加载R包
suppressMessages(library(clusterProfiler))
#获得基因列表
gene <>
#基因名称转换，返回的是数据框
gene = bitr(gene, fromType='SYMBOL', toType='ENTREZID', OrgDb='org.Hs.eg.db')
de <>
## GO分析
go <>'org.Hs.eg.db', ont='all')
library(ggplot2)
p <>'ONTOLOGY') +facet_grid(ONTOLOGY~., scale='free')
p
如何看这个图，如何解释？看图容易解释难。简单说来，一看颜色，二看大小。颜色越红p值越小越可信，点越打表示富集的基因比例越多，越有可能被验证。
举个例子，BP那一项中，有个叫neutrophil mediate immunity的GO通路被富集了，看字面意思应该是中性粒细胞介导的免疫反应，联系到这到实验设计，就是一开始的summary和文章部分，我觉得这个是合理的，因为他研究的就是长期慢性炎症对肥胖的影响，这里我们发现，炎症导致的免疫反应可能可以解释这个事情。还有很多其他的GO通路，需要自己根据背景去解释。去发现线索。
还有一个是广为人知的KEGG通路富集分析，Y叔的神包依然可以胜任。
EGG - enrichKEGG(gene= gene$ENTREZID,
organism = 'hsa',
pvalueCutoff = 0.05)
dotplot(EGG)
这一看，全是大咖类的通路，每一个都是身价百万，说明啥问题，这篇文章的实验设计能够从多个方面解释。细胞周期，免疫反应，凋亡，等等，假如是我们自己的结果，我们就需要用自己的背景去思考，这么多通路是并联的还是串联的，是以此发生的级联反应，还是单兵作战呢？我觉得大部分情况下，我们会尝试用一个通路去说明。这是需要背景知识的。
假如我们对凋亡感兴趣，我们还可以看看我们有多少基因被富集在了凋亡通路，Y叔的神包clusterprofiler也是可以做的。而且做的很云淡风轻。
首先我们把富集的结果变成数据框，便于自己查看。
test <>
browseKEGG(EGG, 'hsa04210')
那么这时候我们就会想，差异基因是有高有低的，而通路富集的时候缺失了这一个信息。一种解决方案是，分别用高表达和低表达基因去做这个分析。我觉得信息是冗余和偏误的，因为任何一个通路的激活，都是伴随着基因的高表达和低表达，本身这就是一个整体，没必要分开分析。
另外一种方法是，先总体分析，然后再标注出上调和下调的基因。那么Y叔的神包也可以做的。
那么我们这个GEO的分析就讲完了。来看看我们获得了哪些图片呢？
也是可以的，我自己觉得这个分析不需要人为设定差异基因的阈值，比如我们认为大于2倍就叫差异基因，这是不科学的，GSEA分析更加靠谱，可以代替GO和KEGG分析。
假如我们通过以上的方法最终锁定了一两个基因，那么我们还可以对这个基因进行批量的相关性分析。
单基因批量相关性分析的妙用
R语言充其量只是个工具，神器不神，主要看人。假如你现在已经读过很多文献，做过很多实验，R语言一定有用，因为你首先想的会是用R语言解决自己的科学问题。
假如你现在还没有怎么读过文献，先不要急着要学习R语言。经常有学生加我微信，上来就说我要学习生信，我一听心里都很慌张，因为我知道，就我掌握的东西离生信还离得远，我就是从科研人员的角度功利地去运用生信知识解决问题。所以学习生信之前，先冷静几天，不要拿自己的爱好跟别人的专业去比较。
不分类别地去多看10分左右的文献可能是更好的策略，再高分的就高屋建瓴谈哲学对初学者不友好，只能学其行不能学其神，低分的容易把初学者带入科研皮条客的大坑，10分左右的文献，起手不高，论证严谨，大部分情况下比许多导师还可靠(这里用词很谨慎)。
有时候我会偏执地认为，读文献是性价比最高的提升科学思维的方法，一个科研人员的水平高低可以简单地用文献阅读量来比较。然后在回过头来看看这个帖子凝练自己的科学问题。有了科学问题我们就有了灵魂，才知道如何在海量资源中筛选有用的数据。
课题设计：收不完的病人查不完的房，临床医生如何快速地设计一个靠谱的课题？
下周，我会更新一个长长的TCGA的帖子，这样我的2018年就这么过去了。
因为考虑到有些朋友是新手，所以我还录制了一个友好的导学视频，教大家如何把帖子里面的内容化为己有，连同R语言安装教程还有platformMap文件一起打包。有需要的朋友，加我微信guotosky，简单介绍一下自己，即可获取。
不需要任何赞赏，如果觉得这个帖子好看，就在右下角点击好看，这是微信的新功能。
祝大家新年快乐。
联系客服
