在bioconductor中生成PrimeView Human Gene Expression Array的annotation database

Joe在留言中谈及primeview的基因芯片的处理问题。刚好最近有人来发了一些关于Affymetrix最新的资料,我可以贴一些上来让大家了解一下这种芯片。

PrimeView的最大的特点就是针对成熟标记的人类基因。它全面覆盖人类基因组中已经注释了的基因,包括最新被注释的基因。对于成熟注释的基因,每个探针组都由11个独立的探针组成,而对于其它的基因,每个探针组都由9个探针组成,重复性及可靠性及高。

注释文件来源(版本或者时间)

数据来源 版本号或者时间
UniGene(design time) 219
RefSeq 41
NCBI genome vesion 37
UCSC 19
Ensembl 57
GenBank 177
Entrez 6/2010

它不完全兼容Genome U133以及Human Gene 1.0 ST。想了解更多,请阅读说明文件

对于Affymetrix的PrimeView Human Gene Expression Array,现在在bioconductor.org的注释包当中还只有primeviewcdf及primeviewprobe两个文件,缺少一个primeview.db的注释文件。现在就来使用AnnotationDbi生成这个注释文件。

首先下载注释文件http://www.affymetrix.com/Auth/analysis/downloads/na32/ivt/PrimeView.na32.annot.csv.zip
解压缩之后,利用它生成一个文本文件。这个文本文件的作用就是将primeview中的probeset_id和bioconductor当中的humanchip_db联系起来。中间的联系可以多种多样的,比如说

gb genbank
ug unigene
eg Entrez
refseq refseq

本文使用Entrez为例,来生成这个桥接的文本文件。

 
> probeset<-read.csv("~/Downloads/PrimeView.na32.annot.csv/PrimeView.na32.annot.csv",comment.char="#")
> Entrez.Gene<-probeset[,c("Probe.Set.ID","Entrez.Gene")]
> Entrez.Gene$Entrez.Gene<-gsub("\\s+","",Entrez.Gene$Entrez.Gene)
> Entrez.Gene<-Entrez.Gene[Entrez.Gene$Entrez.Gene!="---",]
> entrez<-strsplit(as.character(Entrez.Gene$Entrez.Gene),"///")
> names(entrez)<-as.character(Entrez.Gene$Probe.Set.ID)
> yapply<-function(X,FUN, ...) { #
+   index <- seq(length.out=length(X)) #
+   namesX <- names(X) #
+   if(is.null(namesX)) namesX <- rep(NA,length(X))#
+   FUN <- match.fun(FUN) #
+   fnames <- names(formals(FUN)) #
+   if( ! "INDEX" %in% fnames ){ #
+     formals(FUN) <- append( formals(FUN), alist(INDEX=) )   } #
+   if( ! "NAMES" %in% fnames ){ #
+     formals(FUN) <- append( formals(FUN), alist(NAMES=) )   } #
+   mapply(FUN,X,INDEX=index, NAMES=namesX,MoreArgs=list(...)) }
> entrez<-yapply(entrez,function(.ele){cbind(rep(NAMES,length(.ele)),gsub(" ","",.ele))})
> entrez<-do.call(rbind,entrez)
> entrez<-unique(entrez)
> write.table(entrez,file="primeviewHumanGeneExprs.txt",sep="\t",col.names=F,row.names=F,quote=F)

如果你希望使用其它的桥接数据库请自己做相应的调整。因为这个注释文件本身就带有我们需要的EntrezID数据,所以不需要自己做太多的工作。如果没有我们需要的ID数据,那可以参与之前的博文

这里有一点需要注意,就是生成的这个桥接文本文件必须是两列的tsv文件,就是tab separated values文件,tab间隔,无引号的那种。第一列是probeset_id,第二列是桥接ID。无行号,无列名。

接下来的工作就是使用AnnotationDbi生成数据库。整个过程会比较缓慢,要有一定的耐心。

> library(AnnotationDbi)
Loading required package: BiocGenerics
 
Attaching package: ‘BiocGenerics’
 
The following object(s) are masked from ‘package:stats’:
 
    xtabs
 
The following object(s) are masked from ‘package:base’:
 
    anyDuplicated, cbind, colnames, duplicated, eval, Filter, Find, get, intersect, lapply, Map, mapply, mget, order, paste, pmax, pmax.int, pmin, pmin.int, Position, rbind, Reduce,
    rep.int, rownames, sapply, setdiff, table, tapply, union, unique
 
Loading required package: Biobase
Welcome to Bioconductor
 
    Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see 'citation("Biobase")', and for packages 'citation("pkgname")'.
 
> available.chipdbschemas()
 [1] "ARABIDOPSISCHIP_DB" "BOVINECHIP_DB"      "CANINECHIP_DB"      "CHICKENCHIP_DB"    
 [5] "ECOLICHIP_DB"       "FLYCHIP_DB"         "HUMANCHIP_DB"       "HUMANCROSSCHIP_DB" 
 [9] "MOUSECHIP_DB"       "PIGCHIP_DB"         "RATCHIP_DB"         "WORMCHIP_DB"       
[13] "YEASTCHIP_DB"       "ZEBRAFISHCHIP_DB"  
> source("http://bioconductor.org/biocLite.R")
> biocLite("human.db0")
> dir.create("primeview")
> populateDB("HUMANCHIP_DB",affy=F,prefix="primeview",fileName="primeviewHumanGeneExprs.txt",metaDataSrc=c("DBSCHEMA"="HUMANCHIP_DB","ORGANISM"="Homo sapiens","SPECIES"="Human","MANUFACTURER"="Affymetrix","CHIPNAME"="PrimeView Human Gene Expression Array","MANUFACTURERURL"="http://www.affymetrix.com"),baseMapType="eg",outputDir="primeview")
baseMapType is egPrepending MetadataCreating Genes tableAppending ProbesFound 0 Probe AccessionsAppending Gene InfoFound 47378 Gene NamesFound 47378 Gene SymbolsAppending ChromosomesAppending Cytogenetic LocationsAppending OmimAppending RefSeqAppending PubmedAppending UnigeneAppending ChrLengthsAppending 3 GO tablesAppending 3 GO ALL tablesAppending KEGGAppending ECAppending Chromosome LocationsAppending PfamAppending PrositeAppending AliasAppending EnsemblAppending UniprotAppending Metadata 
 
simplifying probes tabledropping redundant data[1] TRUE
> seed<-new("AnnDbPkgSeed",Package = "primeview.db",Version="1.0.0",PkgTemplate="HUMANCHIP.DB",AnnObjPrefix="primeview")
> makeAnnDbPkg(seed,file.path("primeview","primeview.sqlite"),dest_dir="primeview")
Creating package in primeview/primeview.db

我们在当前文件夹下新建了一个primeview的文件夹,其下会有一个primeview.db的文件夹,这就是生成的注释文件。
最后的工作就是安装这个文件。

> install.packages("primeview/primeview.db",repos=NULL,type="source")
* installing *source* package ‘primeview.db’ ...
** R
** inst
** preparing package for lazy loading
** help
*** installing help indices
** building package indices
** testing if installed package can be loaded
*** arch - i386
*** arch - x86_64
 
* DONE (primeview.db)
> library(primeview.db)
Loading required package: org.Hs.eg.db

11 thoughts on “在bioconductor中生成PrimeView Human Gene Expression Array的annotation database

  1. 非常感谢您的这篇博文,我已经安装上述的生成了db。
    qc()时还是出现报错,显示:
    > rawData.qc<-qc(rawData)
    错误于.qc.file.setQCArray(l) : Array name should be a single string.

    我也试图将文件名改成简单的字符串 但是依旧不行,不知是何原因,向您请教了。

  2. hi 又过来请教您问题了。
    想比较一下不同normalization方法,查到affycomp这个包
    看了说明 google的 但是还是没有搞懂怎么用这个包

    难道这个包真的就是这样“This package isn’t designed to work with your data, it is designed to work with the Affy and GeneLogic spike-in datasets. The goal of the package is to make it easy for people to compare their preferred summarization methods easily.”的目的?

    那如何对自己的数据 进行有针对性的方法选择 比较呢?

  3. 你好!我在下载最新的注释文件后,生成db文件时出现了错误:
    >makeAnnDbPkg(seed,file.path(“primeview”,”primeview.sqlite”),dest_dir=”primeview”)
    Error in match(x, table, nomatch = 0L) :
    object ‘doc_template_names’ not found

    请问是错误在哪儿呢?不太明白。谢谢

    1. 因为这里的代码是早期的代码,现在最新的包已经变更了很多。您需要查看AnnotationForge
      说明中的最新的代码。你会很快发现,全新的代码下,如果是affymetrix的话,桥接文件都不需要自己准备。祝你好运。

  4. 问一个比较弱弱的问题,在做这个时间想的。
    affy给出的annot注释文件了,里面已经有很详细的注释信息了。那为啥还要在做一个db呢?直接从注释的文件里将ID和gene对应上不行吗?
    还有个问题,由此衍生出来的。不同的CDF文件可能对找差异基因的影响较大,尤其是在brainarray下载的cdf文件,很多都是处理后的一个probe对应一个gene,这样岂不是很容易miss一些基因?但是如果用affy的多个probe,出现一个基因对应的probe变化趋势(一个上、一个下)如何处理呢?或者比如10 probe中,只有2-3个probe有差异,如何判定是否是差异基因呢?
    有些琐碎,谢谢!

    1. 1. 为什么不直接使用affymetrix公司提供的注释信息。
      答:可以直接使用,但是因为affymetrix公司为了将注释信息放在一个xls表单里,所以有很多项都整合在同一个单元格内了,所以不方便直接使用。需要经过转换。bioconductor的db库就是方便了人们,让大家不需要自己再做转换。

      2. 为什么原本的多个probe对应一个gene,转换至db库中就变成1对1的关系了。
      答:这里有个概念需要区别:什么是probe,什么是probeset。前者指的是一个探针,而后者指的是一组探针。一个探针可能会非特异性结合多个基因,也有可能多个基因共享同一段DNA片段。这在设计的时候不可能完全避免。但是对于测定一个基因表达量的一组探针而言,在经过统计运算之后,比如说线性拟合,或者使用贝叶斯判断等,都可以排除不准确的探针而确定基因的表达量。

  5. 如果我理解没错的话,用affy的注释xls自己制作的db,其仍然保留了所有的probe,在挑选差异基因时,是将单个probe进行处理,其实得到的是差异probe,转化成基因后很可能是某个基因有重复。
    而在brainarray下载的cdf,其已经合并了,但是问题就是其选择的某个基因的probe或probeset是否准确,我总是担心他会miss一些信息。
    我现在不成熟的想法是,根据probe来挑差异基因,得到的差异probe对应到gene上,但在合并时需要考虑probe的权重还是什么?(没想好),不仅仅是均值。不知是否合理。

  6. 呵呵,还有个问题讨扰啊。用新的AnnotationForge包时,的确可以用makeAnnDbPkg直接将affy给的注释包做成db文件,但是有个问题,在一个probe对应下的gene或者eg有的会有xxxx///xxx///xxx///的形式,使得最终得到的结果就是这些probe没有被标注gene信息,可能是未识别。如果我要更改这个格式,是将该probe下对应的symbol取出一个就可以了呢?还是将所有的symbol都取出,然后对应到同一个probe下面(和你之前做的桥联文件相似)?若是后者,新的包似乎又不认,比较困惑。。。

    1. 你只需要取出一个就可以了。关键是取出的都必须是一个数据库的注释,否则没办法做出效果比较好的桥联文件。

发表评论

电子邮件地址不会被公开。 必填项已用*标注