GEO数据库视频学习笔记（ID转换）

上一篇笔记记录了如何下载芯片数据（GEO数据库视频学习笔记（芯片数据下载和数据读取）），这篇笔记继续跟着视频学习，视频地址：here

（一）如何从对象里提取表达矩阵

# 第二种方法下载数据，读取进来
> a = read.table("GSE42872_series_matrix.txt.gz",header = TRUE,sep = "\t", quote= "", fill = T,comment.char = "!")
> class(a) #这里读取的是一个data.frame
[1] "data.frame"
> str(a) #这里有6个样品，加一个芯片ID号，一共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包直接读取的（上一篇文章里第三种方法），读取进来的是一个“对象”：

> library(GEOquery)
> gse <- getGEO('gse42872',GSEMatrix = TRUE, AnnotGPL = FALSE) #后面两个参数不是必需的
> 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是完全不一样的数据结构。dataframe就是行和列组成的。而“对象”像一个大箱子，里面放着小箱子，每个小箱子里放着不同信息，比如表达矩阵、分组信息等等。
如果你对一个“对象”用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行

是非常凌乱的。也很复杂。有多复杂。用Jimmy大神的话就是“复杂到你都没必要去了解”。。。
但是你会发现这个对象是list,那么你可以取它的第一个元素看看：

> 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

那么怎么从这个对象里提取表达矩阵呢？

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

这时a1是从对象里读取出来的表达矩阵：

然后和我们前面下载后读取的表达矩阵a做个比较：

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

现在两种方法提取的表达矩阵长得就一样了。

（二）ID转换

（1）方法一：R包

上面看到行名就是探针的ID，那么怎么转成基因ID呢？基因ID转换的步骤需要注意的是，每一种芯片平台都对应了不同的R包来进行ID转换，你需要知道你所用的芯片平台对应的R包。我们用来练习的这个芯片平台是Annotation: GPL6244 ，那么怎么查呢？去生信技能树文章用R获取芯片探针与基因的对应关系三部曲里查：

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

如果你死活装不上这个包，请参考：一个R包怎么也安装不上，那就本地安装呗？

安装后就可以操作了：

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

这里ids这个table里就会有探针的ID和基因的名称：

#看一下基因名这一列里不重复的项（也就是有多少个基因）
> 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)的意思是看每一个基因对应了几个探针
# sort(table(ids$symbol))是从小到大排个序
# tail(sort(table(ids$symbol)))看排序后的表里最后几个基因都对应几个探针
> table(sort(table(ids$symbol)))

    1     2     3     4     5     6     8    10 
18059   599   132    16     5     6     2     2 
#这里看到绝大部分基因只对应一个探针

可以画个图更清晰的看多少基因对应几个探针：

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

> exprSet = a1
#检查a1的探针ID有多少存在在刚才得到的ids表里
> table(rownames(exprSet) %in% ids$probe_id) #%in%的意思是匹配的意思，判断左边的元素是否存在在右边的元素里
FALSE  TRUE 
13483 19814
#只取那些存在于ids表里的那些基因所对应的探针的行
> exprSet = exprSet[rownames(exprSet) %in% ids$probe_id,]
> dim(exprSet)
[1] 19814     6
#同样把上面得到的ids表里也只保留和表达矩阵里存在的基因对应的探针
> 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
#现在表达矩阵里的探针和ids表里的探针就对应起来了

上面的match()函数用法在这里：R实战第五篇：常用函数的用法

下面对表达矩阵分类，这样如果有一个基因名对应多个探针的情况，就会得到一个新的探针列表：

> tmp = by(exprSet,
+          ids$symbol,
+          function(x) rownames(x)[which.max(rowMeans(x))])
# by()分类汇总
# rowMeans(x)行平均值（行是样品，这里实际上求的是某个探针对应的表达量）
# which.max()输出的是括号里最大值对应的序号
> probes = as.character(tmp) #probes就是新的探针列表
> dim(exprSet) #表达矩阵过滤前
[1] 19814     6
> exprSet = exprSet[rownames(exprSet) %in% probes,] #把表达矩阵按照新的探针列表进行过滤
> dim(exprSet) #过滤后
[1] 18821     6

（2）方法二：

你也可以直接下载GPL，GPL是什么东西？参考：解读GEO数据存放规律及下载，一文就够
GPL：每个数据集对应的芯片平台

#考验网速的时候到了：
> gpl <- getGEO('GPL6244',destdir=".")
File stored at: 
./GPL6244.soft

下载完了，这个gpl也是个“对象”：

你可以看这个对象里都有哪些列：

> 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" 
#根据列名你可以提取你想要的列进行提取
#这里Jimmy大神提取的是1,6,7列，但是他用的练习文件不是上面的这个。
#每一个芯片的gpl内容都不一样，对应的列也不一样，这里要注意！！！根据你需要的提取！！！
> head(Table(gpl)[,c(1,6,7)])
#最后保存一下
> write.csv(Table(gpl)[,c(1,6,7)],"GPL****.csv")

每一列都是什么，网站上都告诉你了：