生信学习——GEO数据挖掘

写在前面——按照生信技能树的学习路线，学完R语言就该学习GEO数据挖掘了。有人说GEO数据挖掘可以快速发文（https://zhuanlan.zhihu.com/p/36303146），不知道靠不靠谱。反正学一学总没有坏处。看完Jimmy老师的视频，写一篇总结方便日后复习。这里有很多操作在 《生信人的20个R语言习题》都可以见到，那里写的更加详细。

视频教程：https://www.bilibili.com/video/BV1is411H7Hq?p=1
代码地址：https://github.com/jmzeng1314/GEO

STEP1：表达矩阵ID转换

首先理解下面的4个概念：
GEO Platform (GPL)
GEO Sample (GSM)
GEO Series (GSE)
GEO Dataset (GDS)
理解起来也很容易。一篇文章可以有一个或者多个GSE数据集，一个GSE里面可以有一个或者多个GSM样本。多个研究的GSM样本可以根据研究目的整合为一个GDS，不过GDS本身用的很少。而每个数据集都有着自己对应的芯片平台，就是GPL。
​

用R获取芯片探针与基因的对应关系三部曲-bioconductor
http://www.bio-info-trainee.com/1399.html

# setwd(dir = "geo_learn/")

### step 1 ###
# 获得GSE数据集的表达矩阵
if(F){
   
  suppressPackageStartupMessages(library(GEOquery))
  gset <- getGEO('GSE42872', destdir=".",
                 AnnotGPL = F,
                 getGPL = F)
  save(gset,file='GSE42872_gset.Rdata')
}
load('GSE42872_gset.Rdata')

exprSet <- exprs(gset[[1]])
pdata <- pData(gset[[1]])
group_list <- c(rep('control', 3), rep('case', 3))

# 以下操作等同于exprs(gset[[1]])
# a <- read.table("GSE42872_series_matrix.txt.gz", 
#                 sep = "\t", quote = "", fill = T,
#                 comment.char = "!", header = T)
# rownames(a) <- a[,1]
# a <- a[,-1]

### step 2 ###
# 根据gset中的Annotation: GPL6244找到对应的R包，安装并使用
# BiocManager::install("hugene10sttranscriptcluster.db")
suppressPackageStartupMessages(library(hugene10sttranscriptcluster.db))

# 找不到对应R包的可以使用下面这种方法
# gpl <- getGE0('GPL6480', destdir=".")
# colnames(Table(gpl))##[1] 41108 17
# head(Table(gpL)[,c(1,6,7)])
# ## you need to check this,which column do you need
# write.csv(Table(gpl)[,c(1,6,7)],"GPL6400.csv")

### step 3 ###
# 获得探针与基因的对应关系，对表达矩阵进行ID转换
ls("package:hugene10sttranscriptcluster.db")
ids <- toTable(hugene10sttranscriptclusterSYMBOL)

# 将表达矩阵中没有对应基因名字的探针去除
table(rownames(exprSet) %in% ids$probe_id)
dim(exprSet)
exprSet <- exprSet[(rownames(exprSet) %in% ids$probe_id),]
dim(exprSet)

# 将exprSet与ids的数据顺序一一对应
ids <- ids[match(rownames(exprSet),ids$probe_id),]
dim(ids)

# 整合表达矩阵
# 多个探针对应一个基因的情况下，只保留在所有样本里面平均表达量最大的那个探针。
tmp <- by(exprSet,ids$symbol,
          function(x) rownames(x)[which.max(rowMeans(x))])
tmp[1:20]
probes <- as.character(tmp) 
exprSet <- exprSet[rownames(exprSet) %in% probes, ]
dim(exprSet)
dim(ids)

rownames(exprSet) <- ids[match(rownames(exprSet),ids$probe_id),2]

save(exprSet, group_list, file = 'GSE42872_new_exprSet.Rdata')

转换前的exprSet

转换后的exprSet

STEP2：差异分析

load('GSE42872_new_exprSet.Rdata')

# 绘制boxplot，看数据分布是否整齐
library(reshape2)
m_exprSet <- melt(exprSet)
head(m_exprSet)
colnames(m_exprSet) <- c("symbol", "sample", "value")
head(m_exprSet)
m_exprSet$group <- rep(group_list, each = nrow(exprSet))
head(m_exprSet)

library(ggplot2)
ggplot(m_exprSet, aes(x = sample, y = value, fill = group)) + geom_boxplot()

# clustering
# 看聚类效果，效果好则说明数据可用
colnames(exprSet) <- paste(group_list,1:6,sep='')
hc <- hclust(dist(t(exprSet)))
nodePar <- list(lab.cex = 0.6<