TCGA中RNA-counts数据中ensemble_ID转换成gene_ID的方法

1.测试数据的准备和相关包的安装。library(stringr)> d1 <- read.table('test.txt', sep = '\t', header = TRUE)> d1tagtcga1 ENSG00000000003.13 2969 4725 1350 16672ENSG00...

随风而逝*

20864人浏览 · 2020-01-16 17:26:30

随风而逝* · 2020-01-16 17:26:30 发布

1.测试数据的准备和相关包的安装。

library(stringr)
> d1 <- read.table('test.txt', sep = '\t', header = TRUE)
> d1
                 tag    t    c    g    a
1 ENSG00000000003.13 2969 4725 1350 1667
2  ENSG00000000005.5    5   14    2    0
3 ENSG00000000419.11 1608 1588  749  888

安装biomaRt（因为我的R版本是3.5以上的。因此安装方式如下）
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(version = "3.10")
安装biomaRt包
BiocManager::install(c("biomaRt"))

2.去除ensemble_ID的版本号。version其实就是小数点后面的部分，ensembl_ID与其他基因ID进行转换时是不带有小数点的。如果你拿到的数据有小数点，那就没法顺利的merge了，所以就要把它去掉。现在有两种方法可以直接使用。

方法一：

> d1$ensemble_id=unlist(str_split(d1$tag,"[.]",simplify=T))[,1]
> d1
                 tag    t    c    g    a     ensemble_id
1 ENSG00000000003.13 2969 4725 1350 1667 ENSG00000000003
2  ENSG00000000005.5    5   14    2    0 ENSG00000000005
3 ENSG00000000419.11 1608 1588  749  888 ENSG00000000419

方法二：

#将_后面的数字替换为空赋值给a
a<- gsub("\\_\\d*", "", d1$tag)
#将.后面的数字替换为空赋值给ENSEMBL
ENSEMBL <- gsub("\\.\\d*", "", a)
# 将ENSEMBL重新添加到raw_count_filt1矩阵
row.names(d1) <- ENSEMBL
d1 <- cbind(ENSEMBL,d1)
colnames(d1)[1] <- c("ensembl_gene_id")


d1
                ensembl_gene_id                tag    t    c    g    a
ENSG00000000003 ENSG00000000003 ENSG00000000003.13 2969 4725 1350 1667
ENSG00000000005 ENSG00000000005  ENSG00000000005.5    5   14    2    0
ENSG00000000419 ENSG00000000419 ENSG00000000419.11 1608 1588  749  888

3.对基因进行注释-获取gene_symbol，用bioMart对ensembl_id转换成gene_symbol

library(biomaRt)
1.显示一下能连接的数据库
listMarts()

               biomart               version
1 ENSEMBL_MART_ENSEMBL      Ensembl Genes 99
2   ENSEMBL_MART_MOUSE      Mouse strains 99
3     ENSEMBL_MART_SNP  Ensembl Variation 99
4 ENSEMBL_MART_FUNCGEN Ensembl Regulation 99



这里我们选择ensembl数据库
2.用useMart函数选定数据库

plant<-useMart("ensembl")

3.用listDatasets()函数显示当前数据库所含的基因组注释

                           dataset
1         acalliptera_gene_ensembl
2       acarolinensis_gene_ensembl
3        acchrysaetos_gene_ensembl
4        acitrinellus_gene_ensembl
5        amelanoleuca_gene_ensembl
6          amexicanus_gene_ensembl
7            ampachon_gene_ensembl
8          anancymaae_gene_ensembl
9     applatyrhynchos_gene_ensembl
10       atestudineus_gene_ensembl
11            bbbison_gene_ensembl
12         bgrunniens_gene_ensembl
13           bihybrid_gene_ensembl
14             bmutus_gene_ensembl
15         bsplendens_gene_ensembl
16            btaurus_gene_ensembl
17           bthybrid_gene_ensembl
18        cabingdonii_gene_ensembl
19        capalliatus_gene_ensembl
20            caperea_gene_ensembl
21              catys_gene_ensembl
... 有很多

这里我们要获取的基因注释的基因是人类基因，所以选择hsapiens_gene_ensembl

4.用useDataseq()函数选定数据库中的基因组

>mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))  

##这条语句的意思是：选定ensembl数据库中的hsapiens_gene_ensembl基因组


5.选定我们需要获得的注释类型


用lsitFilters()函数查看可选择的类型，选定要获取的注释类型，以及已知注释的类型

listFilters(mart)

                                     name
1                         chromosome_name
2                                   start
3                                     end
4                              band_start
5                                band_end
6                            marker_start
7                              marker_end
8                           encode_region
9                                  strand
10                     chromosomal_region
11                              with_ccds
12                            with_chembl
13          with_clone_based_ensembl_gene
14    with_clone_based_ensembl_transcript
15                            with_dbass3
16                            with_dbass5
17                 with_ens_hs_transcript
18                with_ens_hs_translation
19             with_entrezgene_trans_name
20                              with_embl
21                      with_arrayexpress
22                            with_genedb
...  有很多

选择好数据库，基因组，要获得的注释类型，和已知的注释类型，就可以开始获取注释了

6.用getBM()函数获取注释

hg_symbols<- getBM(attributes=c('ensembl_gene_id','hgnc_symbol',"chromosome_name", "start_position","end_position", "band"), filters= 'ensembl_gene_id', values = d1$ensembl_gene_id, mart = mart)

结果如下：
  ensembl_gene_id hgnc_symbol chromosome_name start_position end_position   band
1 ENSG00000000003      TSPAN6               X      100627108    100639991  q22.1
2 ENSG00000000005        TNMD               X      100584936    100599885  q22.1
3 ENSG00000000419        DPM1              20       50934867     50958555 q13.13

这个函数有
4个参数
attributers()里面的值为我们要获取的注释类型
filters（）里面的值为我们已知的注释类型
values= 这个值就是我们已知的注释类型的数据，把上面我们通过数据处理得到的ensembl基因序号作为ensembl_gene_id 的值
mart= 这个值是我们所选定的数据库的基因组
mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
获取完注释就可以把注释文件和基因表达量文件合并起来了
注释就完成了！

本博主新开公众号，希望大家能扫码关注一下，十分感谢大家。

本文参考：https://www.jianshu.com/p/72ee40bed5b4