又一个有点难的探针注释

1.下载GPL页面的表格文件

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL23159

2.读取它

b = read.delim("GPL23159-69552.txt",
               check.names = F,
               comment.char = "#")
colnames(b)

##  [1] "ID"           "probeset_id"  "seqname"      "strand"      
##  [5] "start"        "stop"         "total_probes" "category"    
##  [9] "SPOT_ID"      "SPOT_ID"

b[1,10]

## [1] "NM_001005484 // RefSeq // Homo sapiens olfactory receptor, family 4, subfamily F, member 5 (OR4F5), mRNA. // chr1 // 100 // 100 // 0 // --- // 0 /// ENST00000335137 // ENSEMBL // olfactory receptor, family 4, subfamily F, member 5 [gene_biotype:protein_coding transcript_biotype:protein_coding] // chr1 // 100 // 100 // 0 // --- // 0 /// OTTHUMT00000003223 // Havana transcript // olfactory receptor, family 4, subfamily F, member 5[gene_biotype:protein_coding transcript_biotype:protein_coding] // chr1 // 100 // 100 // 0 // --- // 0 /// uc001aal.1 // UCSC Genes // Homo sapiens olfactory receptor, family 4, subfamily F, member 5 (OR4F5), mRNA. // chr1 // 100 // 100 // 0 // --- // 0 /// CCDS30547.1 // ccdsGene // olfactory receptor, family 4, subfamily F, member 5 [Source:HGNC Symbol;Acc:HGNC:14825] // chr1 // 100 // 100 // 0 // --- // 0"

3.解析

有用信息在第10列,集成了各种基因的信息。总觉得应该有啥现成工具可以解析吧。没查到,硬上好了。。

library(tidyverse)
tmp = str_split(b[,10]," // ",simplify = T)
dim(tmp)

## [1] 27189  2043

tmp[1,1:3]

## [1] "NM_001005484"                                                                   
## [2] "RefSeq"                                                                         
## [3] "Homo sapiens olfactory receptor, family 4, subfamily F, member 5 (OR4F5), mRNA."

有两个办法,第一反应是提取第三列,括号里面的内容即可。但是两万多行看过去,哟很多特例。放弃治疗。

第二个办法是用第一列内容id转换为symbol。

table(tmp[,2])
## 
##                            Ace View          circbase 
##              5741               402               435 
##           ENSEMBL           GenBank Havana transcript 
##              1128                42               124 
##        lncRNAWiki            RefSeq        UCSC Genes 
##                42             19265                10

发现。。。基因ID的类型还挺多的啊,其中大部队是 RefSeq。

a = data.frame(b[,1],tmp[,1:2])
colnames(a) = c("probe_id","gene_id","bank")
table(a$bank)
## 
##                            Ace View          circbase 
##              5741               402               435 
##           ENSEMBL           GenBank Havana transcript 
##              1128                42               124 
##        lncRNAWiki            RefSeq        UCSC Genes 
##                42             19265                10

a = filter(a,bank!="") %>%
  arrange(bank)
m = unique(a$bank);m
## [1] "Ace View"          "ENSEMBL"           "GenBank"          
## [4] "Havana transcript" "RefSeq"            "UCSC Genes"       
## [7] "circbase"          "lncRNAWiki"
a2 = split(a,a$bank) #把这个数据拆成列表。

一个一个类型的处理,挺麻烦的。

4.逐个击破

4.1. Ace View

数据库介绍直接偷懒GPT

Ace View是一个基因注释和基因结构数据库,提供了人类和其他一些物种的全面基因注释信息。它由NCBI(美国国家生物技术信息中心)开发和维护。

以下是Ace View数据库的主要特点和功能:

1.基因注释信息:Ace View提供了广泛的基因注释信息,包括基因结构、表达模式、转录变体、蛋白质编码序列等。这些注释数据来自多个来源,包括NCBI RefSeq、ENSEMBL、UniProt等。

2.多样化的物种覆盖:除了人类基因组的注释外,Ace View还提供了其他一些物种的基因注释,如小鼠、斑马鱼、果蝇等。这使得研究人员可以在不同物种之间进行比较和分析。

3.基因结构可视化:Ace View通过可视化工具,如基因结构图,帮助用户直观地理解基因的结构和转录变异。这些图形显示了外显子、内含子、UTR(5' UTR和3' UTR)以及其他关键结构元素。

4.集成的功能注释:除了基本的基因结构信息,Ace View还提供了额外的功能注释,如基因功能、GO(Gene Ontology)注释、亚细胞定位等。这些数据有助于研究人员理解基因的功能和调控。

5.查询和下载工具:Ace View提供了灵活的查询和下载工具,使用户可以根据不同的准则搜索和访问所需的注释信息。用户可以通过基因名、序列、类别等方式进行查询,并选择将结果以文本或图形形式下载。

总而言之,Ace View是一个全面的基因注释和基因结构数据库,为研究人员提供了丰富的基因注释信息,并帮助他们理解基因的结构和功能。它是一个有用的工具,可以支持生物学研究、基因组分析和功能注释等领域的工作。

拿出一个ID看看

a2$`Ace View`$gene_id[1]
## [1] "MIIP.lAug10-unspliced"

"MICAL3.vdAug10-unspliced" 是一个命名标识符,可能与基因、转录本或蛋白质相关。让我们逐个解释这个名字的各个组成部分:

MICAL3:这是一个基因的名称,代表某个具体的基因。

vdAug10:这部分通常表示数据版本或更新的时间戳。在这种情况下,"vdAug10" 可能表示数据集的更新日期为 "August 2010"(即2010年8月)。

unspliced:这个术语通常指未剪接(unspliced)的转录本。转录是基因的DNA序列在转录过程中生成RNA的过程。在此过程中,原始RNA转录后需要进行剪接,去除非编码区域(intron)并连接编码区域(exon)。然而,有时会存在未经剪接的转录本,即保留了一些或全部内含子的转录本。

所以以点号为分隔符拆分,保留他的第一部分就行了。

str_split_i(a2$`Ace View`$gene_id[1],"\\.",1)
## [1] "MIIP"
#全部搞掉。
a2$`Ace View`$gene_id = str_split_i(a2$`Ace View`$gene_id,"\\.",1)

4.2.ENSEMBL

a2$ENSEMBL$gene_id[1] 
## [1] "ENST00000415919"

也是转录本ID。可以用bitr转换为symbol

library(clusterProfiler)
library(org.Hs.eg.db)
library(tidyverse)
g = bitr(a2$ENSEMBL$gene_id,fromType = "ENSEMBLTRANS",
         toType = "SYMBOL",OrgDb = "org.Hs.eg.db")
a2$ENSEMBL = inner_join(a2$ENSEMBL,g,by = c("gene_id"="ENSEMBLTRANS"))%>%
  mutate(gene_id = SYMBOL)%>%
  dplyr::select(-4)

4.2.GenBank

只有42个。。。忽略也行的。

搜了一下如何转换,发现一个R包rentrez可以搞定。

又遇到一对多匹配的情况,就直接只保留第一个了。

a2$GenBank
##               probe_id     gene_id    bank
## 1531 TC0100007910.hg.1    BC072388 GenBank
## 1532 TC0100009833.hg.1    BC012069 GenBank
## 1533 TC0100013466.hg.1    BC008063 GenBank
## 1534 TC0100015017.hg.1    BC006095 GenBank
## 1535 TC0100016321.hg.1    BC000790 GenBank
## 1536 TC0100018481.hg.1    BC090943 GenBank
## 1537 TC0300007739.hg.1  BC000512_2 GenBank
## 1538 TC0400009894.hg.1    BC140710 GenBank
## 1539 TC0400011637.hg.1  BC015832_2 GenBank
## 1540 TC0500007434.hg.1  BC157112_9 GenBank
## 1541 TC0500007552.hg.1    BC157872 GenBank
## 1542 TC0600010770.hg.1    BC113541 GenBank
## 1543 TC0600013898.hg.1  BC003627_2 GenBank
## 1544 TC0600014094.hg.1    BC007661 GenBank
## 1545 TC0700006836.hg.1  BC001603_3 GenBank
## 1546 TC0700009292.hg.1 BC157112_10 GenBank
## 1547 TC0800010103.hg.1  BC044587_2 GenBank
## 1548 TC0Y00006573.hg.1  BC075016_4 GenBank
## 1549 TC0Y00007325.hg.1  BC121113_2 GenBank
## 1550 TC1100008612.hg.1    BC000354 GenBank
## 1551 TC1100010222.hg.1    BC001781 GenBank
## 1552 TC1100010987.hg.1    BC019385 GenBank
## 1553 TC1200008677.hg.1    BC007512 GenBank
## 1554 TC1200010857.hg.1  BC073964_2 GenBank
## 1555 TC1200011543.hg.1  BC000455_2 GenBank
## 1556 TC1300008149.hg.1    BC140945 GenBank
## 1557 TC1300009638.hg.1  BC105276_2 GenBank
## 1558 TC1500008658.hg.1    BC066981 GenBank
## 1559 TC1500010237.hg.1    BC000483 GenBank
## 1560 TC1500010688.hg.1    BC136597 GenBank
## 1561 TC1500010707.hg.1    BC105788 GenBank
## 1562 TC1600006459.hg.1    BC146939 GenBank
## 1563 TC1600011373.hg.1    BC014471 GenBank
## 1564 TC1700006471.hg.1    BC064534 GenBank
## 1565 TC1700009870.hg.1    BC066948 GenBank
## 1566 TC1700012423.hg.1    BC107850 GenBank
## 1567 TC1800008333.hg.1    BC070280 GenBank
## 1568 TC1900009454.hg.1    BC172325 GenBank
## 1569 TC2100007508.hg.1    BC036452 GenBank
## 1570 TC2200006432.hg.1  BC156846_2 GenBank
## 1571 TC2200008578.hg.1  BC157112_7 GenBank
## 1572 TC2200009350.hg.1    BC006490 GenBank
library(rentrez)
gene_id <- a2$GenBank$gene_id  
re = list()
for(i in 1:length(gene_id)){
  search_result<- entrez_search(db = "gene", term = gene_id[[i]])
  if(length(search_result$ids)==1){
    summary_result <- entrez_summary(db = "gene", id = search_result$ids)
    re[[i]] = summary_result$name
  }else if (length(search_result$ids)==0){
    re[[i]] <- NA
  }else if (length(search_result$ids)>1){
    summary_result <- entrez_summary(db = "gene", id = search_result$ids)
    re[[i]] = summary_result[[1]]$name
  }
}
a2$GenBank$gene_id = unlist(re)
a2$GenBank = na.omit(a2$GenBank)

4.4. RefSeq

RefSeq是一个由美国国家生物技术信息中心(NCBI)维护的基因组序列和注释数据库。它提供了一种标准化的命名和编号系统,用于唯一标识基因、转录本和蛋白质序列。

在RefSeq中,存在不同类型的RefSeq ID,其中最常见的是以下两种:

1.NM_XXXXXX:这是RefSeq mRNA(messenger RNA)记录的标识符,用于表示基因的转录本。每个转录本都被分配一个唯一的NM ID。

2.NP_XXXXXX:这是RefSeq蛋白质(protein)记录的标识符,用于表示基因的编码蛋白质序列。每个蛋白质序列都被分配一个唯一的NP ID。

这些RefSeq ID通过数字(XXXXXX)来唯一标识每个具体的转录本或蛋白质。NM和NP ID之间存在着相应的对应关系,可以通过匹配前缀和数字部分进行链接。

RefSeq ID是广泛使用的标识符,用于引用和访问NCBI数据库中的基因和蛋白质注释信息。通过这些ID,研究人员可以快速准确地定位特定的基因转录本和蛋白质序列,并进行相应的功能注释和研究。

也是可以用bitr转换的。

head(a2$RefSeq$gene_id)
## [1] "NM_001005484" "NM_152486"    "NM_198317"    "NM_001160184"
## [5] "NM_005101"    "NM_001305275"
g = bitr(a2$RefSeq$gene_id,fromType = "REFSEQ",
     toType = "SYMBOL",OrgDb = "org.Hs.eg.db")
a2$RefSeq = inner_join(a2$RefSeq,g,by = c("gene_id"="REFSEQ"))%>%
  mutate(gene_id = SYMBOL)%>%
  dplyr::select(-4)

4.5 Havana transcript

试了一下木有成功。。。只当把方法记下来。

library(biomaRt)

# 连接到Ensembl数据库
ensembl <- useMart("ensembl")
# 选择人类物种和数据集
ensembl <- useDataset(dataset = "hsapiens_gene_ensembl", mart = ensembl)
# 定义Havana转录本ID列表
havana_transcripts <- a2$`Havana transcript`$gene_id

# 查询Havana转录本对应的基因符号
results <- getBM(attributes = c("ensembl_transcript_id", "external_gene_name"),
                 filters = "ensembl_transcript_id",
                 values = havana_transcripts,
                 mart = ensembl)

# 查看结果
print(results)
## [1] ensembl_transcript_id external_gene_name   
## <0 rows> (or 0-length row.names)

4.6 UCSC Genes

同样搞不定..只能说ID类型真是多

# 加载所需的包
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(GenomicFeatures)
library(org.Hs.eg.db)
# 创建TxDb对象
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
# 定义UCSC转录本ID
ucsc_transcript_id <- a2$`UCSC Genes`$gene_id

# 获取 UCSC 转录本的基因信息
gene_info <- select(txdb, keys = ucsc_transcript_id,
                    keytype = "TXID",
                    columns = c("TXID", "GENEID"),
                    multiVals = "first")
## Error in .testForValidKeys(x, keys, keytype, fks): None of the keys entered are valid keys for 'TXID'. Please use the keys method to see a listing of valid arguments.

"uc011nam.2" 是一个转录本标识符,通常用于描述人类的 alternative splicing isoform(转录本的剪接异构体)。这个标识符是由UCSC基因组浏览器 (UCSC Genome Browser) 分配的。

UCSC基因组浏览器是一个在线的基因组注释和可视化工具,提供了大量的基因、转录本和变体信息。"uc011nam.2" 代表某个特定基因的一个转录本剪接异构体。然而,仅通过这个标识符无法确定所属的基因或其他相关信息。

要获取关于 "uc011nam.2" 转录本的更多详细信息,建议使用UCSC基因组浏览器进行查询,并检索与该转录本相关的注释、序列、表达模式等数据。通过浏览器的搜索功能,您可以输入 "uc011nam.2" 来获取特定转录本的相关信息。

最终搞定的是这几个,也已经对应到两万多基因。

ids = rbind(a2$`Ace View`,
            a2$GenBank,
            a2$ENSEMBL,
            a2$RefSeq)
ids = ids[,-3]
colnames(ids)[2] = "symbol"
head(ids)
##            probe_id   symbol
## 1 TC0100006877.hg.1     MIIP
## 2 TC0100006906.hg.1 C1orf158
## 3 TC0100007305.hg.1    KDM1A
## 4 TC0100007922.hg.1     PPIE
## 5 TC0100008151.hg.1   RAD54L
## 6 TC0100008342.hg.1     PODN
save(ids,file = "ids.Rdata")
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 203,456评论 5 477
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,370评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 150,337评论 0 337
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,583评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,596评论 5 365
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,572评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,936评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,595评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,850评论 1 297
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,601评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,685评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,371评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,951评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,934评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,167评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 43,636评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,411评论 2 342

推荐阅读更多精彩内容