经过比对发现hgu133plus2.db的注释数据与GEO上的注释数据比,多出了一些,想知道这些多出来的注释数据的正确性。
总结来说,就是hgu133plus2.db的注释数据比GEO和affy的注释数据都要多,多出来的那是怎么回事?
还有就是 hgu133plus2注释数据中没有1个探针对应多个entrez id的情况,但是其他两个注释文件都有,这种情况应该如何取舍?选用哪个注释数据比较好?
具体代码如下:
library("hgu133plus2.db")
library(GEOquery)
hug133plus2=toTable(hgu133plus2ENTREZID)
grep(pattern = '/',x=hug133plus2$gene_id,fixed = T)
which(hug133plus2$gene_id=='')
nrow(hug133plus2)
#没有1对多 行数为42307
##取GEO注释文件
GEOanno=getGEO(file ='Mycloud/课题/辐射_GSEA/DataSource/Annotation_file/GPL570.annot.gz')
GEOanno=Table(GEOanno)
colnames(GEOanno)
GEOanno=GEOanno[,c(1,4)]
##看一下有多少空的注释
length(which(GEOanno$`Gene ID`==''))
#9557行空注释 删除!
rindex=which(GEOanno$`Gene ID`=='')
GEOanno=GEOanno[-rindex,]
geo_probe_multigene=length(grep(pattern = '/',x=GEOanno$`Gene ID`,fixed = T))
nrow(GEOanno)
geo_probe_multigene/nrow(GEOanno)
# 1对多的探针占0.049比例 所占比例不多
#affy官方的
affyanno=read.csv(file = 'Mycloud/课题/辐射_GSEA/DataSource/Annotation_file/HG-U133_Plus_2-na36-annot-csv/HG-U133_Plus_2.na36.annot.csv',header = T,as.is = T,comment.char = '#')
colnames(affyanno)
affyanno_entrez=affyanno[,c(1,19)]
length(which(affyanno_entrez$Entrez.Gene=='---'))
rindex=which(affyanno_entrez$Entrez.Gene=='---')
affyanno_entrez=affyanno_entrez[-rindex,]
nrow(affyanno_entrez)
affy_probe_multigene=length(grep(pattern = '/',x=affyanno_entrez$Entrez.Gene,fixed = T))
affy_probe_multigene/nrow(affyanno_entrez)
## 0.03的比例 更低
##注释文件比较
##hug133与GEO文件比较
hugmerge=paste0(hug133plus2$probe_id,hug133plus2$gene_id)
geomerge=paste0(GEOanno$ID,GEOanno$`Gene ID`)
nrow(hug133plus2)
nrow(GEOanno)
## GEO行数比hug多
table(geomerge %in% hugmerge)
dif_geo_hug=GEOanno[!(geomerge %in% hugmerge),]
nrow(dif_geo_hug)/nrow(hug133plus2)
##8%的探针不同。。
##查看后大部分是1对多的探针的 删掉1对多再看看
GEOanno=GEOanno[-(grep(pattern = '/',x=GEOanno$`Gene ID`,fixed = T)),]
nrow(GEOanno)-nrow(hug133plus2)
# 还是多接近600行 再进行比对
geomerge=paste0(GEOanno$ID,GEOanno$`Gene ID`)
table(geomerge %in% hugmerge)
dif_geo_hug=GEOanno[!(geomerge %in% hugmerge),]
nrow(dif_geo_hug)/nrow(hug133plus2)
#这次只有2%的探针不同了
dif_hug_geo=hug133plus2[!(hugmerge %in% geomerge),]
dif_hug_geo[,1]
##为什么这些探针注释不同?
dif_hug_geo[1,1]
hug133plus2[which(hug133plus2$probe_id=='1552258_at'),]
# 结果为 112597
GEOanno[which(GEOanno$ID=='1552258_at'),]
which(GEOanno$ID=='1552258_at')
## 没有这个探针 说明一开始这个探针没有对应的enrezid 被当成空探针删除了
##那么hug133的注释数据对不对?
## 测试一下affy官方的数据
affyanno_entrez[which(affyanno_entrez$Probe.Set.ID=='1552258_at'),]
## 结果也没有 why?
阅读全文
收起全文