Bioconductor没想象的那么简单-part8-注释信息必
刘小泽写于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)，此时不需要考虑多重比对的情况，因此可以filter掉
no_na <- mapIds(hugene20sttranscriptcluster.db, ids2, "SYMBOL", "PROBEID",multiVals = "filter") %>% na.omit() %>%length()
no_map <- length(ids2) - no_na #用全部的id数量减去去掉NA之后的数量
欢迎关注我们的公众号～_～
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到jieandze1314@gmail.com
Welcome to our bioinfoplanet!
