学会了GEO的数据处理，又能怎样？

学会了GEO的数据处理，又能怎样？
GEO数据库包罗万象，对于每个领域的科研工作者很有帮助。

我最开始分析GEO芯片主要使用别人的代码，跑流程。

但是绝对不能报错，一报错就束手无策。

遇到比较困难的地方就使用别人的小工具来代替。

比如，ID转换，聚类分析等。

直到去年我才能够完全使用R语言实现。

这两年，网上的教程帖子井喷一样涌现，大有教师超过学员的趋势。

但是，GEO的分析跟TCGA不一样，里面要考虑的细节比较多。

比如我们今天要解决以下几个问题：
1.GEO数据如何方便地下载？
2.数据什么时候需要标准化？
3.什么时候做log转换？
4.芯片探针如何转换？
5.差异基因如何获得？
6.配对样本的对比矩阵如何构建？
7.配对样本如何作图展示？
8.热图火山图究竟在什么？
9.GO分析，KEGG分析能干什么事情？
看文献得到GSE号，这是数据的名片，GSE32575，
文献中提及的GSE号
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)
可以先简单看一下整体样本的表达情况
boxplot(exprSet,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)
校正后
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 no t needed')}
这个语句会自动判断，如果已经log话就跳过，并提示
'log2 transform not needed'
如果没有log话，他自动log2，并且提示：
'log2 transform finished'
这时候再来看这个数据就规则多了
log2转换后
但是我们发现一个问题，这个表达数据行名是探针，不是我们熟悉的基因名称如果TP53，BRCA1这样的，所以我们需要注释他。

5.获取注释信息
这步会很顺利，也会很难，是技术上的限速环节。

skr！GEO芯片数据的探针ID转换
有些GEO平台的探针转换比较麻烦
正则表达式是我们认识世界的哲学
我们根据上述帖子里面提到的platformMap，找到对应的R包是：'illuminaHumanv2.db'
安装R包并加载
if(length(getOption('BioC_mirror'))==0) options(BioC_mirror ='https:///bioc/')
if(!require('illuminaHumanv2.db')) BiocManager::install('illu minaHumanv2.db',update = F,ask = F)
library(illuminaHumanv2.db, character.only=TRUE)
获取探针对应的symbol信息
probeset <- rownames(exprSet)
SYMBOL <- annotate::lookUp(probeset, platformDB, 'SYMB OL')
symbol = as.vector(unlist(SYMBOL))
制作probe2symbol转换文件
probe2symbol = data.frame(probeset,symbol,stringsAsFact ors = F)
probe2symbol部分结果
有了这个文件，我们就可以进行探针转换了
6.探针转换与基因去重
一个基因会对应对个探针，有些基因名称会是重复的，这些都需要处理。

对于，多个探针，我们选取在样本中平均表达量最高的探针作为对应基因的表达量。

一下代码完成所有事情，而且可以复用。

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')
现在数据变成这个样子
探针去重与转换
而且，行数从原来的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',numb
er=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',n umber=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
而配对样本的展示，我们用TCGA的数据演示过很多次
找出TCGA中的配对样本并正确展示数据
配对样本数据如何计算及展示？
9.后续分析
差异分析后还有很多操作，
画一张热图
library(pheatmap)
## 设定差异基因阈值，减少差异基因用于提取表达矩阵
allDiff_pair=topTable(fit,adjust='BH',coef='group_listafter',n umber=Inf,p.value=0.05,lfc =0.5)
##提前部分数据用作热图绘制
heatdata <- exprSet[rownames(allDiff_pair),]
##制作一个分组信息用于注释
annotation_col <- data.frame(group_list)
rownames(annotation_col) <- colnames(heatdata)
#如果注释出界，可以通过调整格子比例和字体修正
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)
data <- topTable(fit,adjust='BH',coef='group_listafter',numb er=Inf)
data$significant <- as.factor(data$adj.P.Val<0.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, si ze=3,col='green')+
geom_text_repel(data=subset(data, abs(logFC) > 1),
aes(label=gene),col='black',alpha = 0.8)
定制火山图
火山图画起来简单，难的是如何定制化展示。

比如在图上有不同的颜色，不同的点来表示数据。

此外还有耳熟能详的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 <- rownames(allDiff)
#基因名称转换，返回的是数据框
gene = bitr(gene, fromType='SYMBOL', toType='ENTREZID', OrgDb='org.Hs.eg.db')
de <- gene$ENTREZID
## GO分析
go <- enrichGO(gene = de, OrgDb = 'org.Hs.eg.db', ont='all ')
library(ggplot2)
p <- dotplot(go, split='ONTOLOGY') +facet_grid(ONTOLOG
Y~., scale='free')
p
GO分析
KEGG分析
EGG <- enrichKEGG(gene= gene$ENTREZID, organism = 'hsa',
pvalueCutoff = 0.05)
dotplot(EGG)
KEGG分析
只要是参加过线下课的同学，或者只要R语言基本能力过关的同学，轻而易举就能重复这个帖子的所有内容。

即使是零基础，安装这个帖子准备好工作环境，也能完成。

学习R语言，从这一课开始
但是，又能怎么样？Y叔的这个clusterprofiler包还可以作出更多好看的图片，现在关键的问题是，我们如何去是解释, 如何提取有用信息为自己所用。

我们临床医生学习生信很功利，是为了解决自己的问题，首先，你自己得有需要解决的科学问题。

而我自己产生科学问题的方法就三个：
•Search，
•看文献，
•预实验。

这可能是为什么我跟好很多人交流后发现
这样一种现象，上了年纪的老师，博一博二的
师兄师姐，更倾向于使用生信解决自己的问题。

而年级相对偏小的研一研二的朋友(实际上他
们年龄不小)，更容易沉迷技术，喜欢画出云
山雾罩的毛线图。

所以如果你现在是研一，研二，多看文献，培养科学思维才是最重要的。

说到底，生信并不是刚需啊。