你不得不知道的GEO数据库芯片数据分析方法

大家好，我是阿琛。今天来给大家介绍一个新的专题内容---GEO数据库的使用方法。烹饪需要食材，分析需要数据。数据出发，整个研究的第一步就是数据的下载。对于大部分的研究者而言，拿公开的高通量数据，进行二次分析，是最佳的选择途径。
作为与TCGA数据库齐名的一个大型数据库，GEO数据库包罗万象，对于每个领域的科研工作者很有帮助。GEO数据库是一个储存芯片、二代测序以及其他高通量测序数据的一个数据库。利用这个数据库，我们可以检索到其他一些人上传的一些实验测序数据。
下面，我们来看一下如何使用GEO数据库中的芯片数据进行后续分析。
GEO数据库的简介
GEO数据库，GENE EXPRESSION OMNIBUS，是由美国国立生物技术信息中心NCBI创建并维护的基因表达数据库（官网：https://www.ncbi.nlm.nih.gov/geo/
）。

它最初创建于2000年，主要用于收录各国研究机构提交的高通量基因表达数据，也就是说只要是在目前已经发表的绝大部分论文中，其涉及到的基因表达检测的数据，包括芯片数据，二代测序结果，以及其他形式的高通量检测结果，都可以通过这个数据库中找到。
首先，我们进入GEO数据库中，根据GSE编号，查看一下该数据的一些相关信息。在搜索栏中输入编号，“GSE39582”，然后点击“Search”按钮，进行检索。
当然，熟悉了以后，我们也可以直接输入网址进行快速检索，https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE39582；对于检索页面网址而言，其前面是一样的，唯一的区别是acc=后面的GSE编号，修改编号，可以快速进入对应的结果页面，这样也可以在一定程度上减少由于网速等原因所带来的阻碍。
在结果页面中，我们可以初步看一下该结果对应的发布时间，题目，物种等信息。同时，Experiment type中Expression profiling by array表示该结果是通过芯片获得的表达谱；Overall design中简单介绍了整个研究的设计方案和分组信息。
GEO数据的下载
当获得了芯片的GSE编号后，我们接下来需要将其对应的数据进行下载，从而根据自己的需要进一步分析。关于数据的下载，我们这里主要介绍三种不同的方法。
方法一：下载芯片的原始数据
在检索页面中，一路下拉；在Supplementary file中点击Download中的custom，展开原始数据对应列表；点击“Select all“，然后点击Download，即可将所有样本的原始数据RAW Data文件下载；
虽然这是最直接的方法，但是RAW Data文件相对较大，对下载的网速要求相对较高，而且不同的芯片来源，有不同的处理方法，甚至有些芯片没有处理方法，因为其是对应定制的。所以，一般情况下，不推荐大家下载原始数据。
方法二：下载表达矩阵（series matrix）
在Download family中点击Serier Matrix File(s)，进入下载页面；
待下载完成后，可以直接使用read.table()函数读取进来。
rt <- read.table("GSE39582_series_matrix.txt.gz", sep = '\t', header = T, comment.char = '!')
可以看到，在芯片中，包含了54675个基因探针，586个患者。

方法三：使用R的GEOquery包里面的getGEO()函数直接读取进来（推荐）
library(GEOquery)gset <- getGEO("GSE42872", destdir = ".")
当然，考虑到网速问题，我们可以对参数进行设置，选择不下载平台的注释文件，因为一般来讲注释文件是相对比较大的。

gset <- getGEO("GSE42872", destdir = ".", AnnotGPL = F, getGPL = F)
如果把之前下载的series Matrix文件放在当前目录下，getGEO()函数会直接检测到该文件，并进而直接将其进行读取；

我们可以直接查看一下下载结果gset的变量类型。
class(gset)
可以看到，变量gset是一个列表的形式。

为什么是list格式呢？因为一个GEO芯片项目，是可以对应多个芯片平台的，那么每个平台的数据结果会对应list里面的一个元素。
gset[[1]]
既然是列表，自然可以提取其中的第一个元素出来查看一下。可以看到，其中展示了包含了6个样本，33297个特征，以及相关的临床信息，PMID号，以及注释平台信息。
提取表达和临床信息

3.1 通过pData函数获取分组信息

pdata <- pData(gset[[1]])table(pdata$source_name_ch1)
通过pData()函数，即可提取表达数据中的临床信息；同时，点击Environment中的pdata查看，我们可以查看里面的相关信息。可以看到，其中，非肿瘤的样品19例。

因此，根据临床信息，我们可以对样品进行分组，分为肿瘤组和正常组；
library(stringr)group_list <- ifelse(str_detect(pdata$source_name_ch1, "Adenocarcinoma"), "tumor", "normal")#设置参考水平group_list = factor(group_list, levels = c("normal","tumor"))table(group_list)
最终，得到正常组19例样品，肿瘤组566例样品。

3.2 通过exprs()函数获取表达矩阵并校正
整理完临床信息后，我们需要提取对应的表达数据。对于表达数据，除了下载Series Matrix后直接使用read.table()函数进行读取外，我们也可以直接从GEOquery下载得到的变量gset中进行提取。
exp <- exprs(gset[[1]])boxplot(exp,outline=FALSE, notch=T,col=group_list, las=2)
使用exprs()函数可以从gset[[1]]提取表达信息；同时，我们可以使用boxplot()函数先简单看一下整体样本的表达情况。

由于每一次技术重复的时候，都会有误差，芯片的原始数据是由仪器读取的，不同的读取时间，或者扫描仪光线的强弱都会导致同一类型的样本出现误差。正式分析前，我们需要对其进行人工校正一下。这里我们用limma包内置的一个函数，
normalizeBetweenArrays()函数。
library(limma) exp=normalizeBetweenArrays(exp)boxplot(exp,outline=FALSE, notch=T,col=group_list, las=2)
可以看到，经过校正，整个表达水平基本趋于一致。

range(exp)
此外，使用range()函数查看一下表达数据exp的取值范围；一般而言，范围在20以内的表达值基本已经经过了log对数转换。

ID变换
整理好了表达矩阵以后，我们需要将探针的id转换成为基因的Gene symbol。对于探针id的转换过程，目前主要是通过R包来进行转换。接下来，我们来看一下如何进行芯片探针id的转换过程。

方法一：使用R包转换

随着芯片平台的普遍使用，其基因的注释信息也被整理成了不同的R包；因此，通常情况下我们使用R包来注释。不同的平台，对应着不同的R包。首先，我们来看一下这个数据集使用的平台类型。
index = gset[[1]]@annotation
通过提取列表gset[[1]]中的注释信息，可以看到，该芯片使用的是我们最常见的平台，GPL570。

对于GPL570，其对应的R包是hgu133plus2.db包；查找显示，其储存在Bioconductor中，下载并进行安装R包。
if(!require("hgu133plus2.db")) BiocManager::install("hgu133plus2.db")library(hgu133plus2.db)
首先，我们来看看，在hgu133plus2.db包中，包含了哪些信息；

ls("package:hgu133plus2.db")
可以看到，除了Symbol信息外，在其中还包含了Ensemble id，Entrez id等信息，可以需要进行提取。

ids <- toTable(hgu133plus2SYMBOL)head(ids)
提取其中的Symbol信息，可以看到，最终获得了probe id和Gene symbol的对应信息。

length(unique(ids$symbol))table(sort(table(ids$symbol)))
其中，经过去重复，一共存在20174个不同的Gene symbol，且部分基因存在多条探针的对应关系。

接下来，我们需要将其进行一一对应匹配。
library(tidyverse)exp <- as.data.frame(exp)exp <- exp %>% mutate(probe_id=rownames(exp)) %>% inner_join(ids,by="probe_id") %>% select(probe_id, symbol, everything())exp <- exp[!duplicated(exp$symbol),]rownames(exp) <- exp$symbolexp <- exp[,-(1:2)]
经过id匹配，并去重复，最终得到了20174个基因的表达结果；

同时，我们可以查看一下前3行前3列的表达情况。
exp[1:3,1:3]

当然，除了R包注释外，还有其他的注释方法，比如使用网页下载的soft文件进行注释，或者有些特殊的芯片内容，需要自己手工比对注释。

方法二：使用soft文件注释

方法三：手工注释

PCA分析
表达矩阵到此基本整理完成。接下来，在正式的差异分析之前，我们首先可以做一个PCA分析，整体水平查看正常和肿瘤两组样品直接是否存在显著的差异。
library(FactoMineR)library(factoextra)dat=as.data.frame(t(exp))dat.pca <- PCA(dat, graph = FALSE)pca_plot <- fviz_pca_ind(dat.pca, geom.ind = "point", col.ind = group_list, palette = c("#00AFBB", "#E7B800"), addEllipses = TRUE, legend.title = "Groups")
结果显示：
在该芯片中，癌和癌旁组织的表达水平存在一定的差异。

差异分析
对于芯片数据的差异分析，我们一般使用limma包来进行。关于差异分析的输入文件，主要是两个，第一是整理好的表达矩阵，其中行名为基因名，列名为样本名；第二是分组信息（group list）。
library(limma)design=model.matrix(~group_list)fit=lmFit(exp,design)fit=eBayes(fit)deg=topTable(fit,coef=2,number = Inf)
最终，使用topTable()函数提取所有基因的差异分析。


colnames(deg)
可以看到，在结果表格中，包含了6块的内容，包括我们常见的logFC值，以及P.Value，adj.P.Val等等。

logFC=1.5P.Value = 0.05type1 = (deg$P.Value < P.Value)&(deg$logFC < -logFC)type2 = (deg$P.Value < P.Value)&(deg$logFC > logFC)deg$change = ifelse(type1,"down",ifelse(type2,"up","stable"))table(deg$change)
接下来，我们需要根据设定的阈值，｜logFC｜>1.5和P值<0.05，将显著差异表达的基因进行分组。结果显示：显著下调的基因有427个，显著上调的基因有242个。

可视化分析
1、火山图
对于差异分析结果，火山图和热图是两种最常见的展示方式。首先，我们来看一下火山图的绘制方法。对于火山图的绘制，大家可以参考之前的推文（生信最重要的图之一，十分钟帮你搞定！建议收藏！！）。

方法一：基于ggpubr包绘制火山图

#install.packages("ggpubr")#install.packages("ggthemes")library(ggpubr)library(ggthemes)deg$logP <- -log10(deg$P.Value)ggscatter(deg, x = "logFC", y = "logP", color = "change", palette = c("blue", "black", "red"), size = 1) + theme_base() + geom_hline(yintercept = -log10(P.Value), linetype = "dashed") + geom_vline(xintercept = c(-logFC, logFC), linetype = "dashed")
结果显示：
接下来，我们可以进一步给火山图添加标签，把显著上调和显著下调基因中前5名的基因名进行标注；
###添加标签deg$Label = ""#新加一列labeldeg <- deg[order(deg$P.Value), ]
up.genes <- head(rownames(deg)[which(deg$change == "up")], 5)down.genes <- head(rownames(deg)[which(deg$change == "down")], 5)#将up.genes和down.genes合并，并加入到Label中deg.top5.genes <- c(as.character(up.genes), as.character(down.genes))deg$Label[match(deg.top5.genes, rownames(deg))] <- deg.top5.genes
ggscatter(deg, x = "logFC", y = "logP", color = "change", palette = c("blue", "black", "red"), size = 1, label = deg$Label, font.label = 8, repel = T, xlab = "log2FoldChange", ylab = "-log10(P-value)") + theme_base() + geom_hline(yintercept = -log10(P.Value), linetype = "dashed") + geom_vline(xintercept = c(-logFC, logFC), linetype = "dashed")
结果显示：

方法二：基于ggplot2包绘制火山图

library(ggplot2)p <-ggplot(data = deg,aes(x = logFC,y = -log10(P.Value))) +geom_point(alpha=0.4,size=3.5,aes(color=change)) +ylab("-log10(Pvalue)")+scale_color_manual(values=c("blue", "grey","red"))+geom_vline(xintercept=c(-logFC,logFC),lty=4,col="black",lwd=0.8) +geom_hline(yintercept = -log10(P.Value),lty=4,col="black",lwd=0.8) +theme_bw()p
结果显示：
当然，我们也可以对其进行添加标签的操作；
x1 = deg %>% filter(change == "up") %>% head(5)x2 = deg %>% filter(change == "down") %>% head(5)label = rbind(x1,x2)volcano_plot <- p + geom_point(size = 3, shape = 1, data = label) + ggrepel::geom_label_repel(data = label, aes(label = rownames(label)), color="black")
结果显示：
2、热图
cg = rownames(deg)[deg$change !="stable"]diff=exp[cg,]
提取差异表达基因的表达情况；

library(pheatmap)annotation_col=data.frame(group=group_list)rownames(annotation_col)=colnames(diff) pheatmap(diff,annotation_col=annotation_col,scale = "row",show_rownames = F,show_colnames =F,color = colorRampPalette(c("navy", "white", "red"))(50),fontsize = 10,fontsize_row=3,fontsize_col=3)
结果显示：
到此，GEO数据库芯片数据的下载，probe id的转换，差异分析已经基本完成了，整个文章中最难的80%内容也已经基本解决。接下来，就是针对这些差异基因的常规分析，包括GO分析，KEGG分析，GSEA分析，蛋白-蛋白互作网络（Protein-protein interaction，PPI）。
—END—
继续阅读
阅读原文