139-R语言处理基因注释信息的几个小知识
豆豆写于19.9.14 使用的示例数据是一个单细胞的sce（SingleCellExperiment）对象。如果对这个对象不了解也没有太大关系，在这次的内容中，我们只需要关心它的基因名，使用一般的表达矩阵也是一样的。后台回复“sce”获得测试数据
我们上游定量得到的基因一般都是Ensembl或Entrez的ID，这样保证了ID号的准确性和稳定性。但是，只看这样的ID我们还是不能知道这个是什么基因，因为常用的还是Symbol ID，像“TP53“、”BRCA1“等，因此就涉及到了一个常用知识点：
基因ID转换
这里是小鼠数据，因此使用小鼠的注释包
**第一种方法：**mapIds(得到的是一个字符串)，并且得到的基因名和原来矩阵的基因名是一一对应的
# 加载测试数据和物种注释库
load('bioinfoplanet_sce.Rdata')
library(org.Mm.eg.db)
symb <- mapIds(org.Mm.eg.db, keys=rownames(sce), keytype="ENSEMBL", column="SYMBOL")
# 看看数据类型
> head(symb)
ENSMUSG00000102693 ENSMUSG00000064842 ENSMUSG00000051951 ENSMUSG00000102851
NA NA "Xkr4" NA
ENSMUSG00000103377 ENSMUSG00000104017
NA NA
# 检测一下它们的Ensembl ID是否一样
> identical(names(symb),rownames(sce))
[1] TRUE
rowData(sce)$SYMBOL <- symb
rowData(sce)$ENSEMBL <- rownames(sce)
**第二种方法：**select(得到的是一个数据框）,它得到的基因名顺序和原矩阵基因名顺序不同，需要再match
library(org.Mm.eg.db)
tmp <- select(org.Mm.eg.db, keys=rownames(sce), keytype="ENSEMBL", column="SYMBOL")
# 看看数据类型
> head(tmp)
ENSEMBL SYMBOL
1 ENSMUSG00000102693 <NA>
2 ENSMUSG00000064842 <NA>
3 ENSMUSG00000051951 Xkr4
4 ENSMUSG00000102851 <NA>
5 ENSMUSG00000103377 <NA>
6 ENSMUSG00000104017 <NA>
# 检测一下它们的Ensembl ID是否一样
> identical(tmp$ENSEMBL,rownames(sce))
[1] FALSE
# 再match一下，按谁匹配谁放第一个
> identical(tmp$ENSEMBL[match(rownames(sce),tmp$ENSEMBL)],rownames(sce))
[1] TRUE
# 和第一种方法的比价，结果就一样啦
> identical(tmp$SYMBOL[match(rownames(sce),tmp$ENSEMBL)],as.character(symb))
[1] TRUE
# 于是第二种方法的结果就是
rowData(sce)$SYMBOL <- tmp$SYMBOL[match(rownames(sce),tmp$ENSEMBL)]
rowData(sce)$ENSEMBL <- rownames(sce)
综上两点，看来使用mapIds更加方便
好，最后看一下添加的基因信息：
> head(rowData(sce))
DataFrame with 6 rows and 3 columns
GeneLength ENSEMBL SYMBOL
<integer> <character> <character>
ENSMUSG00000102693 1070 ENSMUSG00000102693 NA
ENSMUSG00000064842 110 ENSMUSG00000064842 NA
ENSMUSG00000051951 6094 ENSMUSG00000051951 Xkr4
ENSMUSG00000102851 480 ENSMUSG00000102851 NA
ENSMUSG00000103377 2819 ENSMUSG00000103377 NA
ENSMUSG00000104017 2233 ENSMUSG00000104017 NA
上面👆的转Symbol name的操作确实很直观，但是毕竟不是所有的Ensembl ID都有Symbol name对应的。
- 从上面的几行结果看到，没有Symbol 对应的就会被填成”NA“；
- 另外如果出现重复的symbol，我们还要去重。
因此，可以了解下scater包带的一个操作：
整合基因名（ID+Name）
uniquifyFeatureNames() make gene name unique and valid：
- 如果有symbol name的，它会将Ensembl ID替换为symbol name，没有name的仍然使用ID；
- 如果name相同、ID不同（即两个Ensembl ID对应同一个Symbol name），它会将ID与name组合保证特异性
library(scater)
# 参数很简单，第一个是ID，第二是names
# 看这个函数的帮助文档，大体就是能做这么件事
uniquifyFeatureNames(
ID=paste0("ENSG0000000", 1:5),
names=c("A", NA, "B", "C", "A"))
[1] "A_ENSG00000001" "ENSG00000002" "B" "C"
[5] "A_ENSG00000005"
# 对sce的基因名操作
rownames(sce) <- uniquifyFeatureNames(rowData(sce)$ENSEMBL, rowData(sce)$SYMBOL)
head(rownames(sce))
## [1] "ENSMUSG00000102693" "ENSMUSG00000064842" "Xkr4"
## [4] "ENSMUSG00000102851" "ENSMUSG00000103377" "ENSMUSG00000104017"
添加基因所在染色体位置信息
利用这个包： TxDb.Mmusculus.UCSC.mm10.ensGene
可以参考：http://bioinfoblog.it/tag/txdb/
BiocManager::install('TxDb.Mmusculus.UCSC.mm10.ensGene')
library(TxDb.Mmusculus.UCSC.mm10.ensGene)
# TxDb storing transcript annotations（包括CDS、Exon、Transcript的start、end、chr-location、strand）
columns(TxDb.Mmusculus.UCSC.mm10.ensGene)
# 返回一个GRanges对象
head(transcripts(TxDb.Mmusculus.UCSC.mm10.ensGene,columns=c('CDSCHROM')))
# 补充：如果是想看org.db其中的内容
if(F){
keytypes(org.Hs.eg.db)
head(mappedkeys(org.Hs.egENSEMBL))
}
location <- mapIds(TxDb.Mmusculus.UCSC.mm10.ensGene,
keys=rowData(sce)$ENSEMBL,
column="CDSCHROM", keytype="GENEID")
rowData(sce)$CHR <- location
summary(location=="chrM")
## Mode FALSE TRUE NA's
## logical 22428 13 24255
结果会发现有大量的NA，然后比对到线粒体的有13个，在单细胞数据中，太多基因比对到线粒体是不好的现象。具体原因会在后续的细胞质控介绍。
