一、前言

前段時間主管讓我在GEO數(shù)據(jù)庫弄個表達(dá)圖譜毫胜，準(zhǔn)備接下來做各種轉(zhuǎn)錄組分析掖看，現(xiàn)在回過頭來罢荡，覺得這個東西比較簡單灶挟，十分適合上手練習(xí)誉察。思考了一下高诺，決定寫下來分享給大家摘悴。
注意了峭梳，我所說的都是實際工作中會用到的，大家好好看好好學(xué)蹂喻，不過領(lǐng)導(dǎo)有要求不能完全寫工作內(nèi)容葱椭，所以我會有修改。

二口四、什么是GEO數(shù)據(jù)庫

GEO數(shù)據(jù)庫全稱GENE EXPRESSION OMNIBUS孵运，是由美國國立生物技術(shù)信息中心NCBI創(chuàng)建并維護(hù)的基因表達(dá)數(shù)據(jù)庫。它創(chuàng)建于2000年蔓彩，收錄了世界各國研究機(jī)構(gòu)提交的高通量基因表達(dá)數(shù)據(jù)治笨，也就是說只要是目前已經(jīng)發(fā)表的論文，論文中涉及到的基因表達(dá)檢測的數(shù)據(jù)都可以通過這個數(shù)據(jù)庫中找到赤嚼。
我是利用NCBI中GEO DataSets查找GEO數(shù)據(jù)庫資料旷赖，舉個很簡單的例子，我想找胃癌相關(guān)的疾病資料更卒、研究文獻(xiàn)等孵，那么你可以直接搜索gastric carcinoma，當(dāng)然如果這是一種疾病的話逞壁，那你可以直接搜索這種疾病的簡稱流济。

如果你只想關(guān)注人相關(guān)的研究，請在右方選擇類型腌闯。
到這里绳瘟，你必須了解一些比較重要的東西，比如

GEO上有四類數(shù)據(jù)GSM, GSE, GDS, GPL
1.GSM是單個樣本的實驗數(shù)據(jù)
2.GDS是人工整理好的關(guān)于某個話題的GSM的集合姿骏，一個GDS中的GSM的平臺是一樣的
3.GSE是一個實驗項目中的多個芯片實驗糖声，可能使用多個平臺
4.GPL是芯片的平臺，如Affymetrix分瘦， Aglent等

我說簡單點蘸泻，一般我們會利用GSE號在(點擊上圖GES會直接跳轉(zhuǎn)到GEO數(shù)據(jù)庫)GEO數(shù)據(jù)庫中查找資料，那么GSE號就是我們記錄一個實驗項目（就是NCBI每一篇文章啦）的編號嘲玫，如果這個文獻(xiàn)有多個芯片平臺實驗悦施，那么就會產(chǎn)生多個GPL，GPL就是對應(yīng)你所使用的芯片平臺信息去团。GDS用得比較少抡诞，稍微理解一下就行穷蛹。

三、數(shù)據(jù)匯總

這部分跟我的工作有關(guān)昼汗，只是想了解GEO數(shù)據(jù)挖掘使用方法可以跳過不看肴熏。
如果你找的數(shù)據(jù)是一個比較熱門或者比較多人研究的話，你就會得到幾百個或者上千個文獻(xiàn)顷窒，如果我們要對這些文獻(xiàn)進(jìn)行匯總整理蛙吏，這樣做的工程量非常大。你可能需要的東西:GPL鞋吉、GSE鸦做、樣本量、樣本類型坯辩、題目馁龟、樣本處理方法，甚至文章得出什么結(jié)論漆魔。這里我用了我比較擅長的爬蟲來完成這個任務(wù)坷檩，對于大批量資料匯總，能活用自己的技能真的太好不過了改抡，我會另起一篇文章單獨講怎么爬NCBI的數(shù)據(jù)矢炼。

四、GEO數(shù)據(jù)下載

為了完成GEO數(shù)據(jù)挖掘阿纤，做基因表達(dá)譜構(gòu)建句灌，我們必須要下兩個文件（如圖）
1.GPL探針文件
2.表達(dá)矩陣（Series Matrix File）

有萌新可能會說，表達(dá)矩陣人家不是做好了嗎欠拾，我還做啥胰锌？對的，這是別人已經(jīng)處理過的表達(dá)矩陣藐窄，我們只需要再添加上Symbol就可以完成最簡單的GEO數(shù)據(jù)挖掘——基因表達(dá)圖譜了资昧。
是不是很簡單？
對荆忍，就是這么簡單格带，你可以關(guān)掉不看了。

我跟你講刹枉，如果你運氣好叽唱，得到的GPL的文件上就會有對應(yīng)的Symbol，如果沒有微宝，哈棺亭，那就八仙過海各顯神通吧！

五蟋软、數(shù)據(jù)分析（代碼塊）

鴿了這么久居然還有人看這篇文章侦铜，那好专甩，我就把遇到的坑都和大家講一下
我做數(shù)據(jù)分析遇到了很多平臺钟鸵，其他他們的處理方法都不大一樣钉稍，我一一介紹吧
以GPL570為例，這個平臺的數(shù)據(jù)最多棺耍，也是最規(guī)整的贡未。
首先我們讀取GPL文件和GSE文件

setwd("工作路徑")
Sys.setlocale('LC_ALL','C')
#讀取GPL文件
GPL_table = read.table('GPL文件地址',sep = "\t",
                       comment.char = "#", stringsAsFactors = F,
                       header = T, fill = TRUE, quote = "")

第一步，你們先登錄到原網(wǎng)頁中看看下載的GPL文件大小對不對蒙袍，因為NCBI外網(wǎng)連接可能存在網(wǎng)絡(luò)中斷的情況俊卤，文件具體大小可以再這里看到，比如說54675行害幅、合計7.7582M等

順便提一下消恍，下載GPL是下面那個Download full table...按鈕，如果不提供直接下載以现，而是加載到網(wǎng)頁上的平臺狠怨，你們直接CTRL+A復(fù)制到txt里，一樣的邑遏，不過注意一定要等NCBI加載完佣赖。GSE同理，下載完數(shù)據(jù)記得都要核查數(shù)據(jù)大小记盒。

第二步憎蛤，核查GPL探針文件信息。為了完成基因表達(dá)圖譜纪吮，就一定要看探針文件是否有Gene Symbol的信息俩檬，我特地拿570來說，不僅是這個平臺的資料最多碾盟，而且信息也很完整棚辽，因為有一些平臺只有GB_ACC，網(wǎng)上有方法回對GB_ACC到基因名稱巷疼，這里就點一下晚胡。

Gene Symbol

第三步，準(zhǔn)備合并表達(dá)矩陣嚼沿。細(xì)分下來的步驟就是在GPL中找出基因列和ID列估盘。merge到GSE中，替換掉原來的探針列骡尽。這個是我做好的代碼遣妥，大家改一改自己的路徑就可以用了。

setwd("路徑")
Sys.setlocale('LC_ALL','C')
#讀取GPL文件
GPL_table = read.table('GPL570.txt',sep = "\t",
                       comment.char = "#", stringsAsFactors = F,
                       header = T, fill = TRUE, quote = "")  
#讀取GSE文件
GSE4100 <- read.table('GSE4100_series_matrix.txt',sep = "\t",
                        comment.char = "!", stringsAsFactors = F,header = T, fill=TRUE)#43931

#表達(dá)矩陣制作（已有完整基因名稱）
ID_Sybmol = GPL_table[,c(1,11)] #GPL對應(yīng)ID列
colnames(ID_Sybmol)[2]="Symbol" #更改名稱為Symbol攀细，主要是為了對其下面的求平均函數(shù)

#合并ID與基因列
Exp = merge(ID_Sybmol,GSE4100,by.x = "ID",by.y = "ID_REF",all=T)  #45782個
Exp = Exp[,-1]

到這里箫踩，就到表達(dá)矩陣最核心的地方爱态，數(shù)據(jù)居然有4萬多行，但基因?qū)嶋H上只有2萬多個境钟，這是咋回事锦担？
其實你看一下數(shù)據(jù)就會發(fā)現(xiàn)，很多基因是重復(fù)的慨削，而且有些是一對多的探針洞渔。
1.所謂一對多，就是一個探針文件對應(yīng)N個基因缚态，這種探針我們的方法是直接去除磁椒。
2.基因重復(fù)，我們會以求平均的方式去除過多的行玫芦。
3.空值浆熔，0值，我都會去除桥帆。
4.數(shù)據(jù)歸一化（Normalization）医增，目的是為了后去的去批次差異（batch effect），如果你不懂环葵，那請你百度调窍。
因此第四步，數(shù)據(jù)清洗张遭，我上完整的數(shù)據(jù)清洗代碼

#數(shù)據(jù)過濾
Exp = Exp[Exp$Symbol != "",] #45782
Exp = na.omit(Exp)   #45782

#去除
Exp1 = data.frame(Exp[-grep("/",Exp$"Symbol"),]) #去一對多邓萨，grep是包含的意思，-就是不包含菊卷，symbol列不包含 #42984

#求平均值
meanfun <- function(x) {
  x1 <- data.frame(unique(x[,1]))
  colnames(x1) <- c("Symbol")
  for (i in 2:ncol(x)){
    x2 <- data.frame(tapply(x[,i],x[,1],mean))
    x2[,2] <- rownames(x2)
    colnames(x2) <- c(colnames(x)[i], "Symbol")
    x1 <- merge(x1,x2,by.x = "Symbol",by.y = "Symbol",all=FALSE)
  }
  return(x1)
}


Exp2<- meanfun(Exp1)#23520
# 查看數(shù)據(jù)
par(cex = 0.7)
n.sample=ncol(Exp2[,-1])
if(n.sample>40) par(cex = 0.5)
cols <- rainbow(n.sample*1.2)
boxplot(Exp2[,-1], col = cols,main="expression value",las=2)

write.table(Exp2,"Exp_Original.txt",row.names = F,quote = F,sep="\t")

# 校正
row.names(Exp2) = Exp2[,1]
Exp2 = log(Exp2[,-1])  # 我使用的是Log缔恳、這里面大家用的會是log2，目的一樣洁闰，但是最后出的數(shù)據(jù)集中度不一樣歉甚，大家一定要注意

# 測試校正數(shù)據(jù)
par(cex = 0.7)
n.sample=ncol(Exp2)
if(n.sample>40) par(cex = 0.5)
cols <- rainbow(n.sample*1.2)
boxplot(Exp2, col = cols,main="expression value",las=2)

Symbol = row.names(Exp2)
Exp_test = cbind(Symbol,Exp2)

write.table(Exp_test,"Exp.txt",row.names = F,quote = F,sep="\t")

完了我們的表達(dá)矩陣就大功告成啦一共23520個。

轉(zhuǎn)錄組分析——利用GEO數(shù)據(jù)挖掘構(gòu)建基因表達(dá)圖譜