刘小泽写于19.5.22
内容来自Bioc2018的第三章~注释信息必知必会
Bioconductor没想象的那么简单(part1)
Bioconductor没想象的那么简单(part2)
Bioconductor没想象的那么简单(part3)
Bioconductor没想象的那么简单-(part4)
Bioconductor没想象的那么简单-(part5)
Bioconductor没想象的那么简单(part6)
Bioconductor没想象的那么简单-part7
会接触到的注释资源
芯片注释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还可以控制重复值这种情况，需要注意的是：
如果要将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
