GEO數(shù)據(jù)挖掘

第二遍看GEO數(shù)據(jù)挖掘的視頻，比之前更容易看懂很多，也有很多東西是需要重新掌握的，發(fā)現(xiàn)了很多第一遍完全沒有聽到的東西。

P1 0-項(xiàng)目總覽及Github介紹

 # 以文件夾以及Project的形式來管理代碼及文件
 # 讀文獻(xiàn)理解文獻(xiàn)的實(shí)驗(yàn)設(shè)計(jì)
 # 在github中有課程的代碼：https://github.com/jmzeng1314/GEO/tree/master/task1-check-specific-genes

P2 1-通用文獻(xiàn)閱讀及規(guī)律

 # 一個(gè)基因?qū)?yīng)多個(gè)轉(zhuǎn)錄本，即一個(gè)基因?qū)?yīng)多個(gè)探針，之后操作的時(shí)候要進(jìn)行篩選
 # 2015年，芯片→測序儀
 # 不同的芯片的數(shù)據(jù)可以組成一個(gè)大的數(shù)據(jù)，不是重點(diǎn)

P3 2-了解GEO數(shù)據(jù)庫

 # GEO數(shù)據(jù)庫
 # GEO數(shù)據(jù)庫主頁：樣品、平臺(tái)(定制的平臺(tái)可能會(huì)沒有探針以及基因?qū)?yīng)的關(guān)系)
 # ref seq ID: NM開頭，可以跟基因?qū)?yīng)
 # 芯片的基礎(chǔ)知識(shí)
 # illumina與affymatrix的芯片是不一樣的
 # GPL570: hgu133plus2是使用最多的一款芯片
 # 看GSE相關(guān)的文章

P4 3-數(shù)據(jù)下載的3種方式

# 1. Raw data下載 # 不同的芯片的處理方法不一樣，不推薦，分析不好
    ## 用Read.Affy來讀取，舊版的芯片
    ## 不同的版本用不同的函數(shù)來提取
# 2. Series Matrix files # 要注意是否下載完整 Linux 可以用md5檢驗(yàn)
a=read.table('GSE42972_series_matrix.txt.gz',sep='\t',comment.char='!',fill=T,header=T,quote="")
    ## 思路：先用編輯器打開看一下，找到相應(yīng)的函數(shù)來讀取數(shù)據(jù)。
# 3. Geoquerry包下載：下載過來的是一個(gè)對(duì)象，可以把表達(dá)矩陣、分組信息提取出來
library(GEOquerry)
gset<-getGEO('GSE42872',destdir='.',
             AnnotGPL=F,
             getGPL=F)  #避免下載平臺(tái)，因?yàn)楸容^大
gset

P5 4-ID轉(zhuǎn)換大全

library(GEOquerry)
gset<-getGEO('GSE42872',destdir='.',
             AnnotGPL=F,
             getGPL=F)  #避免下載平臺(tái)，因?yàn)楸容^大
gset
class(gset)
b= gset[[1]] 
a1= exprs(b) ##這個(gè)取出來的跟前面的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這個(gè)包就是做了上面的這句話。

## 包裝成函數(shù)  
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包來獲取對(duì)應(yīng)GEO數(shù)據(jù)的實(shí)驗(yàn)信息，參考
#  http://www.bio-info-trainee.com/1085.html
downGSE('GSE42872') #對(duì)多兩個(gè)csv表格
## ID轉(zhuǎn)換的背景知識(shí)
# 探針需要轉(zhuǎn)換成基因名，多個(gè)探針可能對(duì)應(yīng)同一個(gè)基因，需要策略篩選
# 1. 找到GPL平臺(tái)對(duì)應(yīng)的bioconductor的包
# 2. 找到GPL自己的注釋
library(hugene10sttranscriptcluster.db)
?hugene10sttranscriptcluster.db
ls('package:hugene10sttranscriptcluster.db')
ids=toTable(hugene10sttranscriptclusterSYMBOL) #蛋白編碼的基因就2萬個(gè)左右
length(unique(ids$symbol)) #看一下有多少個(gè)基因
tail(sort(table(ids$symbol))) #看一下基因?qū)?yīng)多個(gè)探針的后面幾位
table(sort(table(ids$symbol))) #看一下基因?qū)?yīng)的探針情況 #作者的數(shù)據(jù)可能已經(jīng)把數(shù)據(jù)過濾了，所以要默認(rèn)它是對(duì)的，不然要自己處理。
plot(table(sort(table(ids$symbol))))
exprSet=a1
table(rownames(exprSet))%in%ids$probe_id)
dim(exprSet) #只有經(jīng)過注釋的探針我們才能進(jìn)行分析，這個(gè)代碼用來看過濾之前的探針數(shù)目
exprSet=exprSet[rownames(exprSet)%in%ids$probe_id,]
dim(exprSet) # 看一下過濾完之后的探針數(shù)目
ids=ids[match(rownames(exprSet),ids$probe_id),]# 將ids順序調(diào)整為與表達(dá)矩陣一樣
# 為什么要順序調(diào)整成一樣，就是為了通過邏輯值來判斷，一個(gè)的正確與否跟另外一個(gè)一樣。
head(ids)
exprSet[1:5,1:5]
tmp=by(exprSet,ids$symbol,function(x) rownames(x)[which.max(rowMeans(x)])
#需要對(duì)表達(dá)矩陣進(jìn)行分割，達(dá)到一個(gè)新的探針列表，這里取的是最大值
probes=as.character(tmp)
# 根據(jù)probe來進(jìn)行過濾
dim(exprSet)
exprSet=exprSet[rownames(exprSet)%in%probes,]                                                  
dim(exprSet)                                                                

GPL找不到相對(duì)應(yīng)的包的解決方法

# 參考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-了解你的表達(dá)矩陣

# 使用R語言20題來進(jìn)行演示
## 檢驗(yàn)常見的基因，如內(nèi)參、自己敲低或者過表處理的基因，GAPDH，ACTB的表達(dá)
## 查看分布圖，各個(gè)樣本的表達(dá)的boxplot
## PCA圖與Hclust圖
plot(hclust(dist(t(exprSet))))
# 代碼在R語言的20題，這邊要對(duì)代碼進(jìn)行修改以及合并。

P7 6-差異分析

# 參考微信文章《用Limma對(duì)芯片數(shù)據(jù)做差異分析》
library(limma)
design<- model.matrix(~0+factor(group_list))#構(gòu)建design矩陣，分組矩陣
# 即tmp=data.frame(case=c(0,0,0,1,1,1),control=c(1,1,1,0,0,0)),樣本多就需要用之前的來進(jìn)行操作。
colnames(design)=levels(factor(group_list))#列名是分組信息，并且賦值
rownames(design)=colnames(exprSet) #樣品名是design的行名
## 構(gòu)建比較矩陣
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
## 進(jìn)行分析
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好看的火山圖代碼
# 富集分析，對(duì)上調(diào)下調(diào)基因進(jìn)行注釋，知道這些基因是干什么用的。
# 參考微信文章 《結(jié)果注釋一文就夠》
# 你的通路有多少基因，你的差異基因中有多少被抽中了
# 1000/20000,基因被抽中的概率是1/20,cell cycle的基因假設(shè)有100個(gè)，那么被抽中的個(gè)數(shù)應(yīng)該是5個(gè)，但是如果抽到了20個(gè)，說明cell cycle這個(gè)通路發(fā)生了變化。
## KEGG 分析
library(clusterProfiler)
## 需要將基因轉(zhuǎn)變?yōu)镋NTREZID
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這個(gè)通路距離，它的BgRatio為36/7431,而例子中18/425，差異基因中富集了18個(gè)，說明這個(gè)通路被明顯影響了
vignette('clusterProfiler') #查看包的說明書，需要從頭看到尾
## GO analysis

P9 8-KEGG-GO等數(shù)據(jù)庫的注釋及GSEA分析

# KEGG,GO的缺點(diǎn)
# 超幾何分布檢驗(yàn)的弊端，無法考量每一個(gè)基因的本身作用，事實(shí)上基因是有重要與不重要的區(qū)分的
# GSEA分析
data(geneList,package='DOSE')
geneList
boxplot(geneList) #差異基因的LogFC,要看包的說明書，知道是哪里來的geneList
nrDEG$logFC
geneList=nrDEG$logFC
names(geneList)=rownames(nrDEG)

## 會(huì)報(bào)錯(cuò)，需要轉(zhuǎn)換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)
## 需要自己調(diào)整，將自己的數(shù)據(jù)調(diào)整別人輪子需要的數(shù)據(jù)。
## GSEA軟件需要制作練個(gè)輸入文件

P10 9- 收尾的幾點(diǎn)建議

# 看github的代碼

P11 10- 批量生存分析代碼大放送

# 人群太大了，異質(zhì)性太高，沒有一個(gè)決定性的分子能夠決定死與活，病人情況太復(fù)雜，生存分析只是其中的一種。
# 參考公眾號(hào)《生信技能樹 生存分析》
# 跑R語言的習(xí)題