多种方法实现基因ID SYMBOL 与 ENTREZID转换
2022-08-16 本文已影响0人
小白要变大神
方法一
ID转换，ENTREZID是进行GO分析最好的ID类型（针对clusterProfiler）
ID转换用到的是bitr()函数，但我的SYMBOL 转 ENTREZID 接近40%失败？不知道什么原因？
bitr()的使用方法：
head(deg1) # 差异基因
# 我们的deg1，都没有重复的 SYMBOL，总共 42175个SYMBOL （ENSEMBL）基因
table(duplicated(deg1$SYMBOL))；table(duplicated(deg1$ENSEMBL))
# FALSE
# 42175
# deg1 <- deg1[!duplicated(deg1$SYMBOL),] # 如有重复的SYMBOL 去重
df <- bitr(deg1$SYMBOL,
fromType = "SYMBOL",
toType = c("ENSEMBL","ENTREZID"),
OrgDb = org.Hs.eg.db)
# 合并deg1 df 为d1
d1 <- inner_join(deg1,df); head(d1) #合并
dim(d1)
# [1] 24525 12
table(duplicated(d1$SYMBOL))
# FALSE
# 24525
table(duplicated(d1$ENSEMBL))
# FALSE
# 24525
table(duplicated(d1$ENTREZID))
# FALSE
# 24525
方法二
预先安装AnnotationDbi 和 org.Hs.eg.db，加载org.Hs.eg.db
library(org.Hs.eg.db)
keytypes(org.Hs.eg.db) #查看基因编号系统名称
image.png
k = keys(org.Hs.eg.db, keytype = "ENSEMBL")
head(k,5)
[1] "ENSG00000121410" "ENSG00000175899" "ENSG00000256069" "ENSG00000171428" "ENSG00000156006"
基于提取的ENSEMBL ID，提取对应的所有Gene ID(ENTREZID)，(以及Symbol），并查看一下提取的内容。
list=select(org.Hs.eg.db, keys=k, columns = c("ENTREZID","SYMBOL"), keytype="ENSEMBL")
#'select()' returned 1:many mapping between keys and columns
dim(list)
[1] 29140 3
head(list,5)
ENSEMBL ENTREZID SYMBOL
1 ENSG00000121410 1 A1BG
2 ENSG00000175899 2 A2M
3 ENSG00000256069 3 A2MP1
4 ENSG00000171428 9 NAT1
5 ENSG00000156006 10 NAT2
预先准备的需要转的 ENSEMBL ID，如何找到他们对应的Gene ID(ENTREZID)和Symbol，例如ID 中的，获得的对应关系：ID_list
head( deg1$ENSEMBL) # 我们自己的准备要转的
# [1] "ENSG00000256069" "ENSG00000127837" "ENSG00000129673" "ENSG00000276016" "ENSG00000075624" "ENSG00000204262"
ID_list=list[match(deg1$ENSEMBL，list[,"ENSEMBL"]),]
table(duplicated(ID_list$ENSEMBL))
# FALSE TRUE
# 26852 15323
table(duplicated(ID_list$SYMBOL))
# FALSE TRUE
# 26794 15381
table(duplicated(ID_list$ENTREZID))
#FALSE TRUE
# 26794 15381
ID_list
ENSEMBL ENTREZID SYMBOL
3 ENSG00000256069 3 A2MP1
8 ENSG00000127837 14 AAMP
9 ENSG00000129673 15 AANAT
30 ENSG00000276016 29 ABR
59 ENSG00000075624 60 ACTB
1017 ENSG00000204262 1290 COL5A2
3856 ENSG00000149294 4684 NCAM1
7605 ENSG00000069943 9488 PIGB
8058 ENSG00000173992 9973 CCS
10155 ENSG00000166171 25911 DPCD
17531 ENSG00000177201 127064 OR2T12
# 再合并deg1, ID_list, 为 d2
d2 <- inner_join(deg1, ID_list, by="SYMBOL" ) #合并
head(d2)
dim(d2)
[1] 24560 13
table(duplicated(d2$SYMBOL)) #有重复
#FALSE TRUE
# 24507 53
table(duplicated(d2$ENSEMBL.x))
# FALSE TRUE
# 24507 53
table(duplicated(d2$ENTREZID))
# FALSE TRUE
# 24507 53
方法三
自己下载：
NCBI的基因entrez ID相关文件介绍 , https://www.plob.org/article/9798.html
访问NCBI ftp 下载数据
https://ftp.ncbi.nlm.nih.gov/
https://ftp.ncbi.nlm.nih.gov/gene/DATA/
NCBI下载数据
gene2ensembl，gene2accession， gene2pubmed，gene2go, 以及 gene_info信息文件，它们的核心连接是gene的entrez ID号
人类的：
https://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO/Mammalia/Homo_sapiens.gene_info.gz
最主要的前几列
# 下载 Homo_sapiens.gene_info
# [https://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO/Mammalia/Homo_sapiens.gene_info.gz](https://links.jianshu.com/go?to=https%3A%2F%2Fftp.ncbi.nlm.nih.gov%2Fgene%2FDATA%2FGENE_INFO%2FMammalia%2FHomo_sapiens.gene_info.gz)
less -SN Homo_sapiens.gene_info
#tax_id #物种编号9606 是人类
#GeneID 基因ID 最新的；所以用旧的ID 无法转换，可以尝试参考中方法或者爬虫
#Symbol 基因名
#LocusTag 别名
#tax_id GeneID Symbol LocusTag
9606 1 A1BG - A1B|ABG|GAB|HYST2477
image.png
参考 ：基因 ID 和 Symbol 转换 https://www.jianshu.com/p/7fde2fc6efa4
# 导入Homo_sapiens.gene_info 到 Rstudio，完成id转换
gene_info= fread("Homo_sapiens.gene_info.cp", sep="\t")
head(gene_info)；dim(gene_info)
# 75533 16
names(gene_info)
[1] "#tax_id"
[2] "GeneID"
[3] "Symbol"
[4] "LocusTag"
[5] "Synonyms"
[6] "dbXrefs"
[7] "chromosome"
[8] "map_location"
[9] "description"
[10] "type_of_gene"
[11] "Symbol_from_nomenclature_authority"
[12] "Full_name_from_nomenclature_authority"
[13] "Nomenclature_status"
[14] "Other_designations"
[15] "Modification_date"
[16] "Feature_type"
allID_NEW= data.frame(gene_info[, c("GeneID","Symbol")])
colnames(allID_NEW)=c("ENTREZID","SYMBOL")
table(duplicated(allID_NEW$SYMBOL))
#FALSE TRUE
# 75379 154
table(duplicated(allID_NEW$ENTREZID))
# FALSE
# 75533
head(allID_NEW)
ENTREZID SYMBOL
1 1 A1BG
2 2 A2M
3 3 A2MP1
4 9 NAT1
5 10 NAT2
6 11 NATP
# 我们的deg1，都没有重复的 SYMBOL
table(duplicated(deg1$SYMBOL))；table(duplicated(deg1$ENSEMBL))
# FALSE
# 42175
table(deg1$SYMBOL %in% allID_NEW$SYMBOL)
# FALSE TRUE
# 14399 27776
# 把我们的准备要转的deg1, 和下载的gene_info 合并转ID
d3 = inner_join(deg1, allID_NEW, by="SYMBOL")
dim(d3)
# 27780 12
table(duplicated(d3$SYMBOL)) # SYMBOL有4个重复
# FALSE TRUE
# 27776 4
table(duplicated(d3$ENSEMBL))
# FALSE TRUE
# 27776 4
table(duplicated(d3$ENTREZID)) #没有重复
# FALSE
# 27780
d3$SYMBOL[duplicated(d3$ENSEMBL)]
[1] "TEC" "HBD" "MEMO1" "MMD2"
d3 %>% filter(., SYMBOL == "HBD")
同一个SYMBOL、ENSEMBL 对应有 不同的 ENTREZID？
SYMBOL有重复，ENTREZID不重复
方法四：成功率最高！
head(deg1) # 差异基因
dim(deg1)
head(deg1) # 差异基因
需要提取gtf文件中ensembl_id(gene_id)和gene symbol(gene_name)的对应关系，此gtf文件来源于前面featurecounts 计数时用的gtf，长这样：
gtf文件
从gtf 中提取 gene_id、gene_name 、 gene_type。
此步在linux或者R中操作都可以，我个人比较喜欢用linux命令，因此示范一下在linux中的操作，最后会得到g2s_vm25_gencode.txt文件
cat gtf_geneid2symbol_gencode.sh
# 以下为bash脚本内容，在linux下运行：bash gtf_geneid2symbol_gencode.sh
gtf="gencode.v32.annotation_nochr.gtf"
### 从gtf 中提取 gene_id、gene_name 、 gene_type
grep 'gene_id' $gtf | awk -F 'gene_id \"' '{print $2}' |awk -F '\"' '{print $1}' >gene_id_tmp
grep 'gene_id' $gtf | awk -F 'gene_name \"' '{print $2}' |awk -F '\"' '{print $1}' >gene_name_tmp
grep 'gene_id' $gtf | awk -F 'gene_type \"' '{print $2}' |awk print $1}' >gene_type_tmp
paste gene_id_tmp gene_name_tmp gene_type_tmp >last_tmp
## 生成 gencode.v32.annotation.txt， 包含3列： gene_id gene_name gene_type
uniq last_tmp > gencode.v32.annotation.txt
rm *_tmp
head gencode.v32.annotation.txt
head gencode.v32.annotation.txt
## 多个ensembl id 可 对应与一个gene symbol
gs=read.table("gencode.v32.annotation.txt", header=T)
table(gs$gene_type) # gene_type
table(duplicated(gs$ENSEMBL)); table(duplicated(gs$SYMBOL)) # gene symbol有重复
dim(gs)
# 60609 3
dim(allID_NEW) # 方法三得到的，从genecode下载、提取的
# 75533 2
head(gs)
ENSEMBL SYMBOL gene_type
1 ENSG00000223972 DDX11L1 transcribed_unprocessed_pseudogene
2 ENSG00000227232 WASH7P unprocessed_pseudogene
3 ENSG00000278267 MIR6859-1 miRNA
4 ENSG00000243485 MIR1302-2HG lncRNA
5 ENSG00000284332 MIR1302-2 miRNA
6 ENSG00000237613 FAM138A lncRNA
head(allID_NEW) # 方法三得到的，从genecode下载、提取的
ENTREZID SYMBOL
1 1 A1BG
2 2 A2M
3 3 A2MP1
4 9 NAT1
5 10 NAT2
6 11 NATP
# 合并从gtf 提取出的 gs 和 genecode 的 gene_info 提取的 allID_NEW
table( gs$SYMBOL %in% allID_NEW$SYMBOL )
# FALSE TRUE
# 23350 37259
ID_conversion = inner_join(gs, allID_NEW, by="SYMBOL")
dim(f)
# 37263 4
table(duplicated(ID_conversion$SYMBOL)) # 136个SYMBOL重复
# FALSE TRUE
# 37127 136
table(duplicated(ID_conversion$ENSEMBL)) # 40个ENSEMBL重复
# FALSE TRUE
# 37223 40
table(duplicated(ID_conversion$ENTREZID)) # 132个ENTREZID重复
# FALSE TRUE
# 37131 132
# 保存，方便以后使用
write.table(ID_conversion, "gene_ID_conversion.txt", sep="\t", col.names=T,row.names=T , quote=F)
head(ID_conversion)
ENSEMBL SYMBOL gene_type ENTREZID
1 ENSG00000223972 DDX11L1 transcribed_unprocessed_pseudogene 100287102
2 ENSG00000227232 WASH7P unprocessed_pseudogene 653635
3 ENSG00000278267 MIR6859-1 miRNA 102466751
4 ENSG00000243485 MIR1302-2HG lncRNA 107985730
5 ENSG00000284332 MIR1302-2 miRNA 100302278
6 ENSG00000237613 FAM138A lncRNA 645520
dim(deg1)
# 36769 11
# 根据制作出的 ID_conversion 注释合并我们自己的deg1
d4= inner_join(deg1, ID_conversion )
dim(d4)
大功告成！自己下载的代码转换，比R包转换注释的基因 ENTREZID更多， 转化效率提高！
# 方法一
table(duplicated(d1$ENTREZID))
# FALSE
# 24525
# 方法二
table(duplicated(d2$ENTREZID))
# FALSE TRUE
# 24507 53
#方法三
table(duplicated(d3$ENTREZID))
#FALSE
#27780
# 方法三比方法一多注释 3506个 ID
length(setdiff(d3$ENTREZID, d1$ENTREZID))
# 3506
dif= setdiff(d3$ENTREZID, d1$ENTREZID) #多注释出来的ID
setdiff= d3 %>% filter(., ENTREZID %in% dif)
dim(setdiff)
# 3506 12
head(setdiff)
# 方法4
dim(d4)
多注释出来的ID head(setdiff)
# 看看这些多注释出来的ID 有多少是 Up Down 基因
setdiff %>% filter(., change != "Stable")%>% dim()
# [1] 410 12 # 410个是差异基因，赚了
DEG_setdiff= setdiff %>% filter(., change != "Stable")
# 看看这些多注释出来的ID 都是什么，主要是lncRNA
table(DEG_setdiff$gene_type)
多注释出来的ID 都是什么
附：表达矩阵换行名
参考：gene ID / Gene Symbol / Ensembl ID https://www.jianshu.com/p/3b27c32fa392
###合并
# 去掉ENSEMBL的小数点后的版本号
deg1 = deg1 %>% separate(ENSEMBL, into= c("ENSEMBL", 'foo')) %>% dplyr::select(.,-foo) ; head(deg1)
# 多个ensembl_id 可对应与一个gene symbol
table(duplicated(deg1$ENSEMBL)); table(duplicated(deg1$SYMBOL)) # 我的gene symbol有1000+多重复的，可以去重！
head(deg1);dim(deg1)
head(deg1)
head(gs) # 前面根据gtf提出来的 ，包含3列： gene_id gene_name gene_type
# lncRNA = c("3prime_overlapping_ncRNA", "antisense", "bidirectional_promoter_lncRNA", "lincRNA", "macro_lncRNA", "non_coding", "processed_transcript", "sense_intronic" , "sense_overlapping")
#gtf中共14826个lnc
k1 = gs$gene_type %in% "lncRNA" ; table(k1)
# FALSE TRUE
# 43760 16849
#gtf中共19814个mRNA
k2 = gs$gene_type == "protein_coding";table(k2)
# k2
# FALSE TRUE
# 40644 19965
# deg中有多少lncRNA?
k3 = deg1$gene_type %in% "lncRNA" ;table(k3) #deg中共7501个lnc
# k3
# FALSE TRUE
# 22035 2490
#deg中有多少个mRNA？
k4 = deg1$gene_type =="protein_coding";table(k4)
# k4
# FALSE TRUE
# 6203 18322
# 差异的mRNA和lncRNA 各有多少
table(deg1$change)
# Down Stable Up
# 1477 20505 2543
k5 = deg1$change !="Stable"
#有差异的lncRNA
table(k3&k5)
# FALSE TRUE
# 24079 446
# 有差异的mRNA
table(k4&k5)
# FALSE TRUE
# FALSE TRUE
21254 3271
表达矩阵的行名id转换
head(gs) #包含3列： gene_id 、 gene_name 、gene_type
exp[1:6,1:6] # 表达矩阵
exp = exp[rownames(exp) %in% gs$gene_id,] #从表达矩阵中的行名=gtf中的gene_id
gs = gs[match( rownames(exp), gs$gene_id) , ] #将gs的顺序调整成和表达矩阵行名一致
identical(gs$gene_id, rownames(exp)) #检查是否一致
# [1] TRUE
k = !duplicated(gs$gene_name);table(k) #gene_symbol做行名不能重复，!duplicated() 去重
# k
# FALSE TRUE
# 193 30152
gs = gs[k,]
exp = exp[k,]
rownames(exp) = gs$gene_name
参考： 基因ID转换方法总结 https://www.jianshu.com/p/6b0b5c55058f
