R语言处理表达矩阵基因名的几个小知识

测试开头

今天是生信星球陪你的第447天

大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~

就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~

这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！

豆豆写于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