方法一
ID轉換，ENTREZID是進行GO分析最好的ID類型（針對clusterProfiler）
ID轉換用到的是bitr()函數(shù)，但我的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) #查看基因編號系統(tǒng)名稱
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），并查看一下提取的內(nèi)容。
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 下載數(shù)據(jù)
https://ftp.ncbi.nlm.nih.gov/
https://ftp.ncbi.nlm.nih.gov/gene/DATA/
NCBI下載數(shù)據(jù)
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 轉換 http://www.itdecent.cn/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 計數(shù)時用的gtf，長這樣：
gtf文件
從gtf 中提取 gene_id、gene_name 、 gene_type。
此步在linux或者R中操作都可以，我個人比較喜歡用linux命令，因此示范一下在linux中的操作，最后會得到g2s_vm25_gencode.txt文件
cat gtf_geneid2symbol_gencode.sh
# 以下為bash腳本內(nèi)容，在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
# 根據(jù)制作出的 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 http://www.itdecent.cn/p/3b27c32fa392
###合并
# 去掉ENSEMBL的小數(shù)點后的版本號
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) # 前面根據(jù)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的順序調(diào)整成和表達矩陣行名一致
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轉換方法總結 http://www.itdecent.cn/p/6b0b5c55058f
