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

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. Reply Joe 5月 4,2012 3:15 上午

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

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

  2. Reply Joe 6月 4,2012 1:14 上午

    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. Reply Joe 7月 25,2013 9:20 下午

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

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

    • Reply admin 7月 29,2013 8:21 上午

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

  4. Reply Joe 7月 31,2013 2:08 上午

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

    • Reply admin 8月 1,2013 7:59 上午

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

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

  5. Reply Joe 8月 1,2013 8:03 下午

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

  6. Reply Joe 8月 2,2013 7:46 下午

    非常感谢!多有讨扰!

  7. Reply Joe 8月 2,2013 7:51 下午

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

    • Reply admin 8月 7,2013 9:07 上午

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

Leave a Reply

  

  

  

%d 博主赞过: