又一个有点难的探针注释

作者: 小洁忘了怎么分身 | 来源:发表于2023-08-11 12:52 被阅读0次

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")

又一个有点难的探针注释

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

2.读取它

3.解析

4.逐个击破

4.1. Ace View

4.2.ENSEMBL

4.2.GenBank

4.4. RefSeq

4.5 Havana transcript

4.6 UCSC Genes

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读