如何为affymetrix exon array chip创建一个R下的annotation package 3

任务:有MoEx-1_0-st-v1芯片的注释文件,MoEx-1_0-st-v1.na31.mm9.probeset.csv。现在需要为其创建一个注释包,以方便注释及使用。

MoEx-1_0-st-v1.na31.mm9.probeset.csv可以从http://www.affymetrix.com/estore/browse/products.jsp?productId=131474&categoryId=35765&productName=GeneChip-Mouse-Exon-1.0-ST-Array#1_1下载到。

第一步,获取所有probeset的genbank号。

probeset<-read.csv("MoEx-1_0-st-v1/MoEx-1_0-st-v1.na31.mm9.probeset.csv",comment.char="#")
#所有的genbank号都可以在gene_assignment中找到。在gene_assignment中,有两种分隔符,“///”以及"//"
#///用于分隔不同的注释行 //用于分隔一行内不同的注释内容。其中genbank号在注释行的第一格。
probeset<-probeset[,c("probeset_id","gene_assignment")]
gc()
genbank<-as.character(probeset$gene_assignment)
genbank<-strsplit(genbank,"///")
names(genbank)<-probeset$probeset_id
genbank<-sapply(genbank,function(x){x<-gsub(" ","",x); x<-do.call(rbind,sapply(x,strsplit,"//")); x[,1]})
#下面将genbank变成一对一的列表,需要使用yapply
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(...)) }
 
genbank<-yapply(genbank,function(.ele){cbind(rep(NAMES,length(.ele)),.ele)})
genbank<-do.call(rbind,genbank)
colnames(genbank)<-c("probeset_id","genbank")
rownames(genbank)<-1:nrow(genbank)
#因为其中夹杂着有ensembl id,所以需要把这些号转换成genbank id,这里使用biomaRt来完成转换
ensembl_t_ids<-genbank[grepl("^ENS",genbank[,2]),2]
library(biomaRt)
ensembl = useMart("ensembl")
ensembl = useDataset("mmusculus_gene_ensembl",mart=ensembl)
access<-getBM(attributes = c("ensembl_transcript_id","ox_refseq_mrna__dm_dbprimary_acc_1074"), filters = "ensembl_transcript_id",values=ensembl_t_ids,mart=ensembl)
access<-merge(genbank,access,by.x="genbank",by.y="ensembl_transcript_id")
access<-access[,c(2,3)]
colnames(access)<-colnames(genbank)
genbank<-rbind(genbank[!grepl("^ENS",genbank[,2]),],access)

第二步,将所有probeset的genbank号排除重复之后写入一个文本文件

genbank<-genbank[genbank[,2]!="---",]
genbank<-genbank[genbank[,2]!="",]
genbank<-unique(genbank)
write.table(genbank, file = "MoEx-1_0-st-v1_probeset.txt", sep = "\t", col.names = F, row.names=F, qmethod = "double")

第三步,使用AnnotationDBI来建立注释包

library(AnnotationDbi)
source("http://bioconductor.org/biocLite.R")
biocLite("mouse.db0")
Infile <- "MoEx-1_0-st-v1_probeset.txt"
OutDir <- "where-you-want-to-save-your-package"
library(org.Mm.eg.db)
makeDBPackage("MOUSECHIP_DB",
              affy = FALSE,
              prefix = "moex10stv1probeset",
              fileName = Infile,
              outputDir=OutDir,
              baseMapType = "gbNRef",
              version = "1.0.0",
              manufacturer = "Affymetrix",
              chipName = "moex10stv1",
              manufacturerUrl = "http://www.affymetrix.com",
              author = "Jianhong Ou",
              maintainer = "qiuworld.com")

第四步,试验安装包。在OutDir找到新生成的moex10stv1probeset.db
R CMD INSTALL –build moex10stv1probeset.db

之后你就可以用library(moex10stv1probeset.db)来调用它了.

全文参考了:http://bioinformaticsagatha.wordpress.com/2011/03/14/how-to-create-the-annotation-package-for-mouse-gene-1-1-st-array/

3 thoughts on “如何为affymetrix exon array chip创建一个R下的annotation package

  1. Pingback: 在bioconductor中生成PrimeView Human Gene Expression Array的annotation database | 糗世界

  2. Reply 波波 12月 23,2014 11:54 下午

    博主您好!我想做一些GEO上面的芯片的差异分析并且做出热图。我想请教以下几个问题:
    1.是需要下载CEL的原始数据文件呢,还是下载txt文件呢。如果txt文件是标准化了的,是不是可以直接用了。有必要做cel的分析吗。
    2.如果用CEL文件中的数据,每次都要进行质量控制分析吗。如果质量不好的,就不可用是吗。

    谢谢博主。

    • Reply admin 12月 26,2014 10:15 上午

      您好。
      1. 需要分析情况,如果你只是在一个GEO的数据中进行分析,可以只下载标准化了的txt表达文件。此时可以不从头分析。但是你必须100%相信作者上传了正确的文件。
      2. 质量控制这件事情很难说。有很多时候,我面对的研究者都不希望听到我说他们的芯片质量不好。当我们说质量不好时,必须有足够的证据。因为分析芯片所需要的时间比较短,所以我一般都会做两次,一次去除质量不好的,一次不去除,比较一下两次的结果,哪次更符合预期。这种做法是很不提倡的,因为它带来了倾向性的结果。

Leave a Reply

  

  

  

%d 博主赞过: