欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

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

程序员文章站 2024-03-18 21:28:46
...

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.

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


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

                 https://www.jianshu.com/p/4d0812195b65