118-Bioconductor没想象的那么简单-part8

May 22, 2019 4 min read cnposts

刘小泽写于19.5.22 内容来自Bioc2018的第三章

会接触到的注释资源

芯片注释ChipDb：例如hugene20sttranscriptcluster.db
物种注释OrgDb：如org.Hs.eg.db
位置信息TxDb：它的组成是TxDb.Species.Source.Build.Table 如TxDb.Hsapiens.UCSC.hg19.knownGene ，根据名称就可以推测出，物种是人类，来源是UCSC，版本是hg19，内容是列出已知基因(的位点)；
另外还有一种与之很像，但是是基于Ensembl数据库进行映射注释的EnsDb 包，如：EnsDb.Hsapiens.v79
基因组序列BSgenome，包括物种、提供者、版本、发布日期、发布标准名称、染色体名称等，如BSgenome.Hsapiens.UCSC.hg38
其他：如GO.db、KEGG.db等
在线注释库：AnnotationHub、biomaRt [关于AnnotationHub在part5中介绍过]

与注释包的交互

主要利用select函数，它的使用方法是：select(annopkg, keys, columns, keytype)

annopkg是注释包名称
keys是我们目前有的ID
columns是选择转换输出的列，默认是输出全部的列
keytype指定转换类型，默认使用central keys，也就是核心类型。
那么首先我怎么看有多少种类型呢？使用keytypes(annopkg) 或者columns(annopkg) 另外，我怎么知道其中哪种是核心类型呢？还是有一点规律的：
- 如果是芯片注释包，central keys就是厂商的探针ID；
- 也可以看包名称，比如org.Hs.eg.db中有eg，这就代表了Entrez Gene ID；
- 另外可以加载一个注释包，然后使用head(keys(annopkg)) 看看输出什么内容，帮助猜测；

举个例子

加入现在有一个芯片数据：Affymetrix Human Gene ST 2.0 ，然后想将探针ID转为基因ID

# 比如有5个来自 Affymetrix hugene20 的探针名
ids <- c("16908472","16962185", "16920686","16965513","16819952")
# 转换为一种
library(hugene20sttranscriptcluster.db)
select(hugene20sttranscriptcluster.db, ids, "SYMBOL")
#   PROBEID    SYMBOL
1 16908472 LINC01494
2 16962185      ALG3
3 16920686      <NA>
4 16965513      <NA>
5 16819952      CBFB
#如果要转换成多种，直接可以从keytype参数那里入手
ids <- c('16737401','16657436' ,'16678303') select(hugene20sttranscriptcluster.db, ids, c("SYMBOL","MAP"))
#   PROBEID       SYMBOL     MAP
1  16737401        TRAF6   11p12
2  16657436      DDX11L1 1p36.33
3  16657436 LOC102725121 1p36.33
4  16657436      DDX11L2  2q14.1
5  16657436      DDX11L9 15q26.3
6  16657436     DDX11L10 16p13.3
7  16657436      DDX11L5  9p24.3
8  16657436     DDX11L16    Xq28
9  16657436     DDX11L16    Yq12
10 16678303         ARF1 1q42.13

关于mapIds

从上面例子的结果可以看到，虽然只有三个probe id，但是其中一个probe自己有8个重复项，因为这个探针对应到了多个基因(一般这种情况需要将这个探针去掉)

小tip：如果多个探针对应到一个基因呢？ 可以用最值、平均值、中位值或者MAD值进行过滤

与select函数相似，但mapIds还可以控制重复值这种情况，需要注意的是：

column参数只选择一列
必须指定keytype，而select 对于核心类型可以不指定
增加一个参数：multiVals 用来控制重复值

如果要将probe ID转为Symbol ID，那么就要写成：

mapIds(hugene20sttranscriptcluster.db, ids, "SYMBOL","PROBEID")
 16737401  16657436  16678303 
  "TRAF6" "DDX11L1"    "ARF1"
# 这个结果就没有select返回的许多重复值

上面的例子中没有看到multiVals参数，是因为它有一个默认值first ，表示只保留重复值的第一个，另外还有其他几种：

# 结果返回列表
mapIds(hugene20sttranscriptcluster.db, ids, "SYMBOL", "PROBEID", multiVals = "list")
$`16737401`
[1] "TRAF6"

$`16657436`
[1] "DDX11L1"      "LOC102725121" "DDX11L2"      "DDX11L9"     
[5] "DDX11L10"     "DDX11L5"      "DDX11L16"    

$`16678303`
[1] "ARF1"

# 将重复值作为NA
mapIds(hugene20sttranscriptcluster.db, ids, "SYMBOL", "PROBEID", multiVals = "asNA")
16737401 16657436 16678303 
 "TRAF6"       NA   "ARF1" 

# 过滤掉重复值
mapIds(hugene20sttranscriptcluster.db, ids, "SYMBOL", "PROBEID", multiVals = "filter")
16737401 16678303 
 "TRAF6"   "ARF1"

小练习

Q1: 利用hugene20sttranscriptcluster.db这个包，将Entrez ID为1000的基因对应的基因名输出

>mapIds(hugene20sttranscriptcluster.db, "1000", "SYMBOL", "ENTREZID", multiVals = "asNA")
  1000 
"CDH2"

Q2: PPARG基因对应的Entrez ID是多少

> mapIds(hugene20sttranscriptcluster.db, "PPARG", "ENTREZID", "SYMBOL", multiVals = "asNA")
 PPARG 
"5468"

Q3：GAPDH基因的蛋白数据库ID是多少

# 方法一：mapIds
> mapIds(hugene20sttranscriptcluster.db, "GAPDH", "UNIPROT", "SYMBOL")
'select()' returned 1:many mapping between keys and columns
   GAPDH 
"P04406" 
# 方法二：select
> select(hugene20sttranscriptcluster.db, "GAPDH",  "UNIPROT", "SYMBOL")
'select()' returned 1:many mapping between keys and columns
  SYMBOL UNIPROT
1  GAPDH  P04406
2  GAPDH  V9HVZ4

Q4：给定一个探针列表ids2，如何判断有多少探针只比对到一个基因？如何判断一个基因也没比对上？

# 判断只比对到一个基因数量
tmp <- biomaRt::select(hugene20sttranscriptcluster.db, ids,  "SYMBOL", "PROBEID") %>% na.omit()
sum(table(tmp) == 1 )

# 判断一个基因也没比对上(也就是比对结果为NA)
no_na <- mapIds(hugene20sttranscriptcluster.db, ids2, "SYMBOL", "PROBEID",multiVals = "filter") %>% na.omit() %>%length()
no_map <- length(ids2) - no_na

R Bioconductor