作者：jzhang Update：2020-07-12 Emali：jzhang910@qq.com
前言：我们都知道很多人在进行GEO数据库挖掘的时候，首先遇到的第一个痛点就是探针ID转换成gene symbol的问题，gene symbol因为被大家熟悉所以得到广泛采用。针对此问题，我进行了一个总结贴。
一般ID转换分为以下几个办法：
- 使用GEO数据库提供的GPL文件
- 金标准当然是去基因芯片的厂商的官网直接去下载啦
- 使用bioconductor里面的芯片探针注释包
- 使用探针序列进行序列比对注释（以上办法都不能解决的时候）
- 生信技能树jimmy针对以上问题开发的四个包：idmap1，idmap2，idmap3，AnnoProbe
1.GEO官网
使用ftp地址下载
# 以下网址存在一些特点
https://ftp.ncbi.nlm.nih.gov/geo/platforms/GPL15nnn/GPL15207
https://ftp.ncbi.nlm.nih.gov/geo/platforms/GPL14nnn/GPL14951
https://ftp.ncbi.nlm.nih.gov/geo/platforms/GPL15nnn/GPL15314
R代码版本
options('download.file.method.GEOquery'='auto')
options('GEOquery.inmemory.gpl'=FALSE)
library(GEOquery)
library(Biobase)
# Download GPL file, put it in the current directory, and load it:
gpl <- getGEO("GPL16699", destdir=".")
colnames(Table(gpl))
head(Table(gpl)[,c(1,10,13)])
probe2symbol <- Table(gpl)[,c(1,13)]
2.芯片官网
根据芯片名字搜索得到注释文件。
3.注释包
3.1 bioconductor
library(hgu133plus2.db)
# 查看有哪些对应关系
ls("package:hgu133plus2.db")
# 得到对应关系
probe2symbol2 <- toTable(hgu133plus2SYMBOL)
probe2symbol2[1:6,]
probe2entrezid <- toTable(hgu133plus2ENTREZID)
probe2entrezid[1:6,]
probe2id <- merge(probe2symbol2,probe2entrezid,by="probe_id")
3.2 idmap1
idmap1解决了bioconductor包下载困难的问题，是整合全部的bioconductor里面的芯片探针注释包。
library(devtools)
devtools::install_github("jmzeng1314/idmap1")
library(idmap1)
ids1 <- idmap1::getIDs("gpl570")
head(ids1)
# 你想知道我们支持哪些平台吗，当然是可以看的
# 查看一个数据集合中有哪些数据集
data(package="idmap1")
data(p2s_df)
gpl <- unique(p2s_df[,3])
3.3 idmap2
idmap2解决了GPL平台的soft文件下载困难，把全部表达芯片的soft文件下载好后，提前它们的芯片和基因symbol对应关系并存储好。不同的gpl平台的soft信息不一样，一般来说，大家关心的其实就是探针的ID，以及基因的symbol列。
library(devtools)
devtools::install_github("jmzeng1314/idmap2")
library(idmap2)
ids2 <- idmap2::get_soft_IDs("gpl570")
head(ids2)
# 你想知道我们支持哪些平台吗，当然是可以看的
data(gpl_list)
gpl <- unique(gpl_list$gpl)
3.4 idmap3
idmap3解决了那些并不提供探针的注释信息的平台，大量的GPL平台的soft文件并不提供基因注释信息。是下载全部的GPL的soft文件里面的探针碱基序列比对后重新注释的包。
比如：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL21827
能拿到的信息就探针序列，所以是idmap1和idmap2都是无法注释的，这个时候就需要我们的idmap3啦！
library(devtools)
devtools::install_github("jmzeng1314/idmap3")
library(idmap3)
ids3 <- idmap3::get_pipe_IDs("GPL21827")
head(ids3)
# 你想知道我们支持哪些平台吗，当然是可以看的
data(gpl_list)
gpl_list[,2:1]
3.5 AnnoProbe
整合了前面三个包的内容，根据type参数分为bioc对应bioconductor的注释包，soft对应GEO的soft格式的芯片注释文件，pipe对应重新注释的文件。
library(devtools)
devtools::install_github("jmzeng1314/AnnoProbe")
# 仅仅是一句话，就拿到了这个平台的探针的注释信息。需要注意的是，这个函数的type参数，其实是有3个选择，这里我演示的是选择soft这个来源的基因注释信息。
library(AnnoProbe)
ids <- AnnoProbe::idmap("GPL570",type = "soft")
head(ids)
# 你可以指定ID_type，目前只能是选择 "ENSEMBL" or "SYMBOL"，然后这个函数就会为你进行ID转换及坐标，还有基因类型的注释。
library(AnnoProbe)
IDs <- c("DDX11L1", "MIR6859-1", "OR4G4P", "OR4F5")
ID_type = "SYMBOL"
annoGene(IDs, ID_type)
annoGene(IDs, ID_type,out_file ='tmp.html')
annoGene(IDs, ID_type,out_file ='tmp.csv')
4.实战-重注释
针对只有探针序列的GPL平台芯片，进行重注释 https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA
芯片注释
比如芯片GPL21827，Agilent-079487 Arraystar Human LncRNA microarray V4 在GEO里面仅仅是提供了探针序列
GPL21827
4.1 得到探针文件的fasta序列格式文件
## part1：得到探针序列的fa格式文件
# 这个包需要注意两个配置，一般来说自动化的配置是足够的。
options('download.file.method.GEOquery'='auto')
options('GEOquery.inmemory.gpl'=FALSE)
library(GEOquery)
# 首先下载序列
gset <- getGEO("GPL21827", destdir="../Analysis/data-GPL-anno" )
# 保存到本地
save(gset,file="../Analysis/data-GPL-anno/GPL21827_eSet.Rdata")
# 载入数据
lname <- load("../Analysis/data-GPL-anno/GPL21827_eSet.Rdata")
lname
colnames(Table(gset))
# 得到探针和序列
probe2seq <- Table(gset)[,c(1,4)]
# 保存成fasta格式
all_recs <- paste(apply(probe2seq,1,function(x) paste0('>',x[1],'\n',x[2])),collapse = '\n')
write.table(all_recs,file = "../Analysis/data-GPL-anno/GPL21827.fasta",row.names = F,col.names = F,quote = F)
4.2 比对得到bam文件并提取比对好的坐标
# 比对并排序
index=/teach/database/genome/Ensembl/Homo_sapiens/GRCh38_release95/Homo_sapiens.GRCh38_release95.genome
fa=/teach/project/GPL-anno/GPL21827/fasta/GPL21827.fasta
# 比对
hisat2 -p 4 -x $index -f $fa -S GPL21827.Hisat_aln.sam
# 排序
samtools sort -@ 3 -o GPL21827.Hisat_aln.sorted.bam GPL21827.Hisat_aln.sam
# 比对好的bam文件映射到对应参考基因组的注释信息
bam=GPL21827.Hisat_aln.sorted.bam
# -F 4 过滤掉没有mapping上的reads
samtools view -F 4 $bam | awk 'BEGIN{OFS="\t";}{print $3,$4-60,$4+60,$1,$5,$6}' |perl -alne '{print unless $F[0] =~ /v1/}' >GPL21827_human_align.bed
4.3 根据参考基因组的gtf文件提取已知注释坐标
# 格式化gff文件
wget -c ftp://ftp.ensembl.org/pub/release-95/gff3/homo_sapiens/Homo_sapiens.GRCh38.95.gff3.gz
zcat Homo_sapiens.GRCh38.95.gff3.gz |perl -alne '{next unless $F[2] eq "gene";/gene_id \"(.*?)\";.*?gene_name \"(.*?)\";.*?gene_biotype \"(.*?)\";/;print "$2\t$3\t$1\t$F[0]\t$F[3]\t$F[4]"}' >Homo_sapiens.GRCh38.95.gff3.txt
## 制作bed文件
# 编码蛋白基因的bed
grep -w protein_coding Homo_sapiens.GRCh38.95.gff3.txt |awk 'BEGIN{OFS="\t";}{print $4,$5,$6,$1,$2,$3}' >human_pro.bed
# 非编码蛋白的基因bed
grep -v -w protein_coding Homo_sapiens.GRCh38.95.gff3.txt |awk 'BEGIN{OFS="\t";}{print $4,$5,$6,$1,$2,$3}' >human_non.bed
# 全部的bed
cat Homo_sapiens.GRCh38.95.gff3.txt |awk 'BEGIN{OFS="\t";}{print $4,$5,$6,$1,$2,$3}' >human_all.bed
4.4 提取交集得到探针注释信息
bedtools intersect -a GPL21827_human_align.bed -b human_all.bed -f 1 -wa -wb >GPL21827_human_anno.bed
最后注释结果为：
GPL21827芯片注释结果
