ixxmu
commented
1 year ago 又一个有点难的探针注释 by 生信星球
这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！
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")
如果因为代码看不懂，而跟不上正文的节奏，可以来找我，相当于来一个新手保护期。我的课程都是循环开课。下一期的时间，点进去咨询微信咯 找我们学习生信
https://mp.weixin.qq.com/s/qcdNJ9CHS54zASe78--dmg