GEO數(shù)據(jù)庫視頻學(xué)習(xí)筆記（ID轉(zhuǎn)換）

上一篇筆記記錄了如何下載芯片數(shù)據(jù)（GEO數(shù)據(jù)庫視頻學(xué)習(xí)筆記（芯片數(shù)據(jù)下載和數(shù)據(jù)讀取））瑟蜈，這篇筆記繼續(xù)跟著視頻學(xué)習(xí)烟逊，視頻地址：here

（一）如何從對象里提取表達(dá)矩陣

# 第二種方法下載數(shù)據(jù)，讀取進(jìn)來
> a = read.table("GSE42872_series_matrix.txt.gz",header = TRUE,sep = "\t", quote= "", fill = T,comment.char = "!")
> class(a) #這里讀取的是一個(gè)data.frame
[1] "data.frame"
> str(a) #這里有6個(gè)樣品铺根，加一個(gè)芯片ID號(hào)宪躯，一共7列
'data.frame':   33297 obs. of  7 variables:
 $ X.ID_REF.    : int  7892501 7892502 7892503 7892504 7892505 7892506 7892507 7892508 7892509 7892510 ...
 $ X.GSM1052615.: num  7.25 6.83 4.4 9.48 4.55 ...
 $ X.GSM1052616.: num  6.81 6.7 4.51 9.68 4.45 ...
 $ X.GSM1052617.: num  7.73 7.02 4.88 9.63 5.12 ...
 $ X.GSM1052618.: num  6.19 6.2 4.36 9.69 4.87 ...
 $ X.GSM1052619.: num  7.05 6.77 4.18 9.91 5.16 ...
 $ X.GSM1052620.: num  7.2 6.24 4.73 9.66 3.99 ...

如果你要是用GEOquery包直接讀取的（上一篇文章里第三種方法），讀取進(jìn)來的是一個(gè)“對象”：

> library(GEOquery)
> gse <- getGEO('gse42872',GSEMatrix = TRUE, AnnotGPL = FALSE) #后面兩個(gè)參數(shù)不是必需的
> gse
$GSE42872_series_matrix.txt.gz
ExpressionSet (storageMode: lockedEnvironment)
assayData: 33297 features, 6 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: GSM1052615 GSM1052616 ... GSM1052620 (6 total)
  varLabels: title geo_accession ... cell type:ch1 (34 total)
  varMetadata: labelDescription
featureData
  featureNames: 7892501 7892502 ... 8180418 (33297 total)
  fvarLabels: ID GB_LIST ... category (12 total)
  fvarMetadata: Column Description labelDescription
experimentData: use 'experimentData(object)'
  pubMedIds: 24469106 
Annotation: GPL6244

“對象”和dataframe是完全不一樣的數(shù)據(jù)結(jié)構(gòu)位迂。dataframe就是行和列組成的访雪。而“對象”像一個(gè)大箱子，里面放著小箱子掂林，每個(gè)小箱子里放著不同信息臣缀，比如表達(dá)矩陣、分組信息等等泻帮。
如果你對一個(gè)“對象”用str：

> str(gse)
List of 1
 $ GSE42872_series_matrix.txt.gz:Formal class 'ExpressionSet' [package "Biobase"] with 7 slots
  .. ..@ experimentData   :Formal class 'MIAME' [package "Biobase"] with 13 slots
  .. .. .. ..@ name             : chr "Tiffany,Joy,Parmenter"
  .. .. .. ..@ lab              : chr ""
  .. .. .. ..@ contact          : chr "tiffanyparmenter@gmail.com"
  .. .. .. ..@ title            : chr "Expression data from BRAFV600E A375 melanoma cells treated with vehicle or vemurafenib"
......#此處省略N行

是非常凌亂的精置。也很復(fù)雜。有多復(fù)雜锣杂。用Jimmy大神的話就是“復(fù)雜到你都沒必要去了解”脂倦。。元莫。
但是你會(huì)發(fā)現(xiàn)這個(gè)對象是list,那么你可以取它的第一個(gè)元素看看：

> gse[[1]]
ExpressionSet (storageMode: lockedEnvironment)
assayData: 33297 features, 6 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: GSM1052615 GSM1052616 ... GSM1052620 (6 total)
  varLabels: title geo_accession ... cell type:ch1 (34 total)
  varMetadata: labelDescription
featureData
  featureNames: 7892501 7892502 ... 8180418 (33297 total)
  fvarLabels: ID GB_LIST ... category (12 total)
  fvarMetadata: Column Description labelDescription
experimentData: use 'experimentData(object)'
  pubMedIds: 24469106 
Annotation: GPL6244

那么怎么從這個(gè)對象里提取表達(dá)矩陣呢赖阻？

> b = gse[[1]]
> View(b)
> a1 = exprs(b)
> View(a1)

這時(shí)a1是從對象里讀取出來的表達(dá)矩陣：

然后和我們前面下載后讀取的表達(dá)矩陣a做個(gè)比較：

> rownames(a) = a[,1]
> a = a[,-1]
> View(a)

現(xiàn)在兩種方法提取的表達(dá)矩陣長得就一樣了。

（二）ID轉(zhuǎn)換

（1）方法一：R包

上面看到行名就是探針的ID踱蠢，那么怎么轉(zhuǎn)成基因ID呢火欧？基因ID轉(zhuǎn)換的步驟需要注意的是，每一種芯片平臺(tái)都對應(yīng)了不同的R包來進(jìn)行ID轉(zhuǎn)換朽基，你需要知道你所用的芯片平臺(tái)對應(yīng)的R包布隔。我們用來練習(xí)的這個(gè)芯片平臺(tái)是Annotation: GPL6244 离陶，那么怎么查呢稼虎？去生信技能樹文章用R獲取芯片探針與基因的對應(yīng)關(guān)系三部曲里查：

> BiocManager::install("hugene10sttranscriptcluster.db")
> library("hugene10sttranscriptcluster.db")

如果你死活裝不上這個(gè)包，請參考：一個(gè)R包怎么也安裝不上招刨，那就本地安裝唄霎俩？

安裝后就可以操作了：

> ids = toTable(hugene10sttranscriptclusterSYMBOL)
> View(ids)

這里ids這個(gè)table里就會(huì)有探針的ID和基因的名稱：

#看一下基因名這一列里不重復(fù)的項(xiàng)（也就是有多少個(gè)基因）
> length(unique(ids$symbol))
[1] 18821
> tail(sort(table(ids$symbol)))

  RPL41  UBTFL1  CDK11B  UBE2D3    IGKC LRRFIP1 
      6       6       8       8      10      10 
#這里table(ids$symbol)的意思是看每一個(gè)基因?qū)?yīng)了幾個(gè)探針
# sort(table(ids$symbol))是從小到大排個(gè)序
# tail(sort(table(ids$symbol)))看排序后的表里最后幾個(gè)基因都對應(yīng)幾個(gè)探針
> table(sort(table(ids$symbol)))

    1     2     3     4     5     6     8    10 
18059   599   132    16     5     6     2     2 
#這里看到絕大部分基因只對應(yīng)一個(gè)探針

可以畫個(gè)圖更清晰的看多少基因?qū)?yīng)幾個(gè)探針：

> plot(table(sort(table(ids$symbol))))

> exprSet = a1
#檢查a1的探針I(yè)D有多少存在在剛才得到的ids表里
> table(rownames(exprSet) %in% ids$probe_id) #%in%的意思是匹配的意思，判斷左邊的元素是否存在在右邊的元素里
FALSE  TRUE 
13483 19814
#只取那些存在于ids表里的那些基因所對應(yīng)的探針的行
> exprSet = exprSet[rownames(exprSet) %in% ids$probe_id,]
> dim(exprSet)
[1] 19814     6
#同樣把上面得到的ids表里也只保留和表達(dá)矩陣?yán)锎嬖诘幕驅(qū)?yīng)的探針
> ids = ids[match(rownames(exprSet),ids$probe_id),]
> head(ids)
  probe_id    symbol
1  7896759 LINC01128
2  7896761    SAMD11
3  7896779    KLHL17
4  7896798   PLEKHN1
5  7896817     ISG15
6  7896822      AGRN
> exprSet[1:5,1:5]
        GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619
7896759    8.75126    8.61650    8.81149    8.32067    8.41445
7896761    8.39069    8.52617    8.43338    9.17284    9.10216
7896779    8.20228    8.30886    8.18518    8.13322    8.06453
7896798    8.41004    8.37679    8.27521    8.34524    8.35557
7896817    7.72204    7.74572    7.78022    7.72308    7.53797
#現(xiàn)在表達(dá)矩陣?yán)锏奶结樅蚷ds表里的探針就對應(yīng)起來了

上面的match()函數(shù)用法在這里：R實(shí)戰(zhàn) 第五篇：常用函數(shù)的用法

下面對表達(dá)矩陣分類，這樣如果有一個(gè)基因名對應(yīng)多個(gè)探針的情況打却，就會(huì)得到一個(gè)新的探針列表：

> tmp = by(exprSet,
+          ids$symbol,
+          function(x) rownames(x)[which.max(rowMeans(x))])
# by()分類匯總
# rowMeans(x)行平均值（行是樣品杉适，這里實(shí)際上求的是某個(gè)探針對應(yīng)的表達(dá)量）
# which.max()輸出的是括號(hào)里最大值對應(yīng)的序號(hào)
> probes = as.character(tmp) #probes就是新的探針列表
> dim(exprSet) #表達(dá)矩陣過濾前
[1] 19814     6
> exprSet = exprSet[rownames(exprSet) %in% probes,] #把表達(dá)矩陣按照新的探針列表進(jìn)行過濾
> dim(exprSet) #過濾后
[1] 18821     6

（2）方法二：

你也可以直接下載GPL，GPL是什么東西柳击？參考：解讀GEO數(shù)據(jù)存放規(guī)律及下載猿推，一文就夠
GPL：每個(gè)數(shù)據(jù)集對應(yīng)的芯片平臺(tái)

#考驗(yàn)網(wǎng)速的時(shí)候到了：
> gpl <- getGEO('GPL6244',destdir=".")
File stored at: 
./GPL6244.soft

下載完了，這個(gè)gpl也是個(gè)“對象”：

你可以看這個(gè)對象里都有哪些列：

> colnames(Table(gpl))
 [1] "ID"              "GB_LIST"         "SPOT_ID"        
 [4] "seqname"         "RANGE_GB"        "RANGE_STRAND"   
 [7] "RANGE_START"     "RANGE_STOP"      "total_probes"   
[10] "gene_assignment" "mrna_assignment" "category" 
#根據(jù)列名你可以提取你想要的列進(jìn)行提取
#這里Jimmy大神提取的是1,6,7列捌肴，但是他用的練習(xí)文件不是上面的這個(gè)蹬叭。
#每一個(gè)芯片的gpl內(nèi)容都不一樣，對應(yīng)的列也不一樣状知，這里要注意；辔濉！饥悴！根據(jù)你需要的提忍勾！Ｎ魃琛瓣铣！
> head(Table(gpl)[,c(1,6,7)])
#最后保存一下
> write.csv(Table(gpl)[,c(1,6,7)],"GPL****.csv")

每一列都是什么，網(wǎng)站上都告訴你了：