任务:有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/
博主您好!我想做一些GEO上面的芯片的差异分析并且做出热图。我想请教以下几个问题:
1.是需要下载CEL的原始数据文件呢,还是下载txt文件呢。如果txt文件是标准化了的,是不是可以直接用了。有必要做cel的分析吗。
2.如果用CEL文件中的数据,每次都要进行质量控制分析吗。如果质量不好的,就不可用是吗。
谢谢博主。
您好。
1. 需要分析情况,如果你只是在一个GEO的数据中进行分析,可以只下载标准化了的txt表达文件。此时可以不从头分析。但是你必须100%相信作者上传了正确的文件。
2. 质量控制这件事情很难说。有很多时候,我面对的研究者都不希望听到我说他们的芯片质量不好。当我们说质量不好时,必须有足够的证据。因为分析芯片所需要的时间比较短,所以我一般都会做两次,一次去除质量不好的,一次不去除,比较一下两次的结果,哪次更符合预期。这种做法是很不提倡的,因为它带来了倾向性的结果。
博主您好,为什么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连接打不开,请问下去哪里具体下载MoEx-1_0-st-v1.na31.mm9.probeset.csv呢?谢谢
去这里找找试试运气吧?http://www.affymetrix.com/support/technical/byproduct.affx?product=moexon-st