GEO数据挖掘
第二遍看GEO数据挖掘的视频，比之前更容易看懂很多，也有很多东西是需要重新掌握的，发现了很多第一遍完全没有听到的东西。
P1 0-项目总览及Github介绍
# 以文件夹以及Project的形式来管理代码及文件
# 读文献理解文献的实验设计
# 在github中有课程的代码：https://github.com/jmzeng1314/GEO/tree/master/task1-check-specific-genes
P2 1-通用文献阅读及规律
# 一个基因对应多个转录本，即一个基因对应多个探针，之后操作的时候要进行筛选
# 2015年，芯片→测序仪
# 不同的芯片的数据可以组成一个大的数据，不是重点
P3 2-了解GEO数据库
# GEO数据库
# GEO数据库主页：样品、平台(定制的平台可能会没有探针以及基因对应的关系)
# ref seq ID: NM开头，可以跟基因对应
# 芯片的基础知识
# illumina与affymatrix的芯片是不一样的
# GPL570: hgu133plus2是使用最多的一款芯片
# 看GSE相关的文章
P4 3-数据下载的3种方式
# 1. Raw data下载 # 不同的芯片的处理方法不一样，不推荐，分析不好
## 用Read.Affy来读取，旧版的芯片
## 不同的版本用不同的函数来提取
# 2. Series Matrix files # 要注意是否下载完整 Linux 可以用md5检验
a=read.table('GSE42972_series_matrix.txt.gz',sep='\t',comment.char='!',fill=T,header=T,quote="")
## 思路：先用编辑器打开看一下，找到相应的函数来读取数据。
# 3. Geoquerry包下载：下载过来的是一个对象，可以把表达矩阵、分组信息提取出来
library(GEOquerry)
gset<-getGEO('GSE42872',destdir='.',
AnnotGPL=F,
getGPL=F) #避免下载平台，因为比较大
gset
P5 4-ID转换大全
library(GEOquerry)
gset<-getGEO('GSE42872',destdir='.',
AnnotGPL=F,
getGPL=F) #避免下载平台，因为比较大
gset
class(gset)
b= gset[[1]]
a1= exprs(b) ##这个取出来的跟前面的a相同
#a=read.table('GSE42972_series_matrix.txt.gz',sep='\t',comment.char='!',fill=T,header=T,quote="")
# rownames(a)=a[,1]
# a=a[,-1]
## GEOquerry这个包就是做了上面的这句话。
## 包装成函数
downGSE <- function(studyID = "GSE1009", destdir = ".") {
library(GEOquery)
eSet <- getGEO(studyID, destdir = destdir, getGPL = F)
exprSet = exprs(eSet[[1]])
pdata = pData(eSet[[1]])
write.csv(exprSet, paste0(studyID, "_exprSet.csv"))
write.csv(pdata, paste0(studyID, "_metadata.csv"))
return(eSet)
}
# 使用GEOmetadb包来获取对应GEO数据的实验信息，参考
# http://www.bio-info-trainee.com/1085.html
downGSE('GSE42872') #对多两个csv表格
## ID转换的背景知识
# 探针需要转换成基因名，多个探针可能对应同一个基因，需要策略筛选
# 1. 找到GPL平台对应的bioconductor的包
# 2. 找到GPL自己的注释
library(hugene10sttranscriptcluster.db)
?hugene10sttranscriptcluster.db
ls('package:hugene10sttranscriptcluster.db')
ids=toTable(hugene10sttranscriptclusterSYMBOL) #蛋白编码的基因就2万个左右
length(unique(ids$symbol)) #看一下有多少个基因
tail(sort(table(ids$symbol))) #看一下基因对应多个探针的后面几位
table(sort(table(ids$symbol))) #看一下基因对应的探针情况 #作者的数据可能已经把数据过滤了，所以要默认它是对的，不然要自己处理。
plot(table(sort(table(ids$symbol))))
exprSet=a1
table(rownames(exprSet))%in%ids$probe_id)
dim(exprSet) #只有经过注释的探针我们才能进行分析，这个代码用来看过滤之前的探针数目
exprSet=exprSet[rownames(exprSet)%in%ids$probe_id,]
dim(exprSet) # 看一下过滤完之后的探针数目
ids=ids[match(rownames(exprSet),ids$probe_id),]# 将ids顺序调整为与表达矩阵一样
# 为什么要顺序调整成一样，就是为了通过逻辑值来判断，一个的正确与否跟另外一个一样。
head(ids)
exprSet[1:5,1:5]
tmp=by(exprSet,ids$symbol,function(x) rownames(x)[which.max(rowMeans(x)])
#需要对表达矩阵进行分割，达到一个新的探针列表，这里取的是最大值
probes=as.character(tmp)
# 根据probe来进行过滤
dim(exprSet)
exprSet=exprSet[rownames(exprSet)%in%probes,]
dim(exprSet)
GPL找不到相对应的包的解决方法
# 参考github：https://github.com/jmzeng1314/my-R/tree/master/9-microarray-examples
gpl<- getGEO('GPL6480',destdir='.')
colnames(Table(gpl)) ## 41108 17
head(Table(gpl))[,c(1,6,7)] # 需要自己检查，需要哪几列
write.csv(Table(gpl)[,c(1,6,7)],"GPL6480.csv")
P6 5-了解你的表达矩阵
# 使用R语言20题来进行演示
## 检验常见的基因，如内参、自己敲低或者过表处理的基因，GAPDH，ACTB的表达
## 查看分布图，各个样本的表达的boxplot
## PCA图与Hclust图
plot(hclust(dist(t(exprSet))))
# 代码在R语言的20题，这边要对代码进行修改以及合并。
P7 6-差异分析
# 参考微信文章《用Limma对芯片数据做差异分析》
library(limma)
design<- model.matrix(~0+factor(group_list))#构建design矩阵，分组矩阵
# 即tmp=data.frame(case=c(0,0,0,1,1,1),control=c(1,1,1,0,0,0)),样本多就需要用之前的来进行操作。
colnames(design)=levels(factor(group_list))#列名是分组信息，并且赋值
rownames(design)=colnames(exprSet) #样品名是design的行名
## 构建比较矩阵
contrast.matrix<-makeContrasts(paste0(unique(group_list),collpase='-'),levels=design)
contrast.matrix
# paste0(unique(group_list),collpase='-'),知道是谁跟谁比，case比control，还是control比case，a/b → a=1,b=-1
## 进行分析
fit<- lmFit(exprSet,design)
fit2<- contrasts.fit(fit,contrast.matrix)
fit2<- eBayes(fit2)
tempOutpu=topTable(fit2,coef=1,n=Inf)
nrDEG=na.omit(tempOutput)
head(nrDEG)
## 画热图
library(pheatmap)
choose_gene=head(rownames(nrDEG),25)
choose_matrix=exprSet[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))
pheatmap(choose_matrix)
# 参考 《差异分析 一文不够》
P8 7-火山图及热图制作及美化
# volcano plot
colnames(nrDEG)
## plot(nrDEG$logFC,-log10(nrDEG$P.value))
## 可以google好看的火山图代码
# 富集分析，对上调下调基因进行注释，知道这些基因是干什么用的。
# 参考微信文章 《结果注释一文就够》
# 你的通路有多少基因，你的差异基因中有多少被抽中了
# 1000/20000,基因被抽中的概率是1/20,cell cycle的基因假设有100个，那么被抽中的个数应该是5个，但是如果抽到了20个，说明cell cycle这个通路发生了变化。
## KEGG 分析
library(clusterProfiler)
## 需要将基因转变为ENTREZID
gene.df<- bitr(gene,fromType='SYMBOL',
toType=c('ENSEMBL',"ENTREZID"),
OrgDb=org.Hs.eg.db)
head(gene.df)
## ENTREZID
kk<-enrichKEGG(gene=gene.df$ENTREZID,
organism='hsa',
pvalueCutoff=0.05)
head(kk)[,1:6]
## 拿DNA replication这个通路距离，它的BgRatio为36/7431,而例子中18/425，差异基因中富集了18个，说明这个通路被明显影响了
vignette('clusterProfiler') #查看包的说明书，需要从头看到尾
## GO analysis
P9 8-KEGG-GO等数据库的注释及GSEA分析
# KEGG,GO的缺点
# 超几何分布检验的弊端，无法考量每一个基因的本身作用，事实上基因是有重要与不重要的区分的
# GSEA分析
data(geneList,package='DOSE')
geneList
boxplot(geneList) #差异基因的LogFC,要看包的说明书，知道是哪里来的geneList
nrDEG$logFC
geneList=nrDEG$logFC
names(geneList)=rownames(nrDEG)
## 会报错，需要转换ID
gene.df<-bitr(names(geneList),fromType='SYMBOL',
toType=c('ENSEMBL','ENTREZID'),
OrgDb=org.Hs.eg.db)
head(gene.df)
tmp=data.frame(SYMBOL=names(geneList),
logFc=as.numeric(geneList))
tmp=merge(tmp,gene.df,by='SYMBOL')
geneList=tmp$logFc
names(geneList)=gene.df$ENREZID
geneList=sort(geneList,decreasing=T)
## 主要是为了满足kk2的输入文件的条件
kk2<- gseKEGG(geneList =geneList,
organism ='hsa',
nPerm =1000,
minGSSize=120,
pvalueCutoff=0.05,
verbose= FALSE)
head(kk2)
## 需要自己调整，将自己的数据调整别人轮子需要的数据。
## GSEA软件需要制作练个输入文件
P10 9- 收尾的几点建议
# 看github的代码
P11 10- 批量生存分析代码大放送
# 人群太大了，异质性太高，没有一个决定性的分子能够决定死与活，病人情况太复杂，生存分析只是其中的一种。
# 参考公众号《生信技能树 生存分析》
# 跑R语言的习题
