Search
Question: limma2annaffy
0
gravatar for Guest User
4.6 years ago by
Guest User12k
Guest User12k wrote:
Dear Jim, I am performing differential expression analysis using affymetrix human gene 2.0 ST arrays. I get an error message from limma2annaffy() function. Error in .checkKeysAreWellFormed(keys) : keys must be supplied in a character vector with no NAs I get the annotated HTML and txt tables for each desired contrast by using topTable and annaffy tools but I would like to know what is going wrong with limma2annaffy since I find this function truly useful and convenient. Thank you in advance for your time and your help, Jose -- output of sessionInfo(): R version 3.0.2 (2013-09-25) -- "Frisbee Sailing" Copyright (C) 2013 The R Foundation for Statistical Computing Platform: x86_64-apple-darwin10.8.0 (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. [R.app GUI 1.62 (6558) x86_64-apple-darwin10.8.0] [History restored from /Users/Jose/.Rapp.history] > setwd("/Users/Jose/Documents/Oscar/MDA_261213/MDA_271213") > library(limma) > library(oligo) Loading required package: BiocGenerics Loading required package: parallel Attaching package: ???????BiocGenerics??????? The following objects are masked from ???????package:parallel???????: clusterApply, clusterApplyLB, clusterCall, clusterEvalQ, clusterExport, clusterMap, parApply, parCapply, parLapply, parLapplyLB, parRapply, parSapply, parSapplyLB The following object is masked from ???????package:limma???????: plotMA The following object is masked from ???????package:stats???????: xtabs The following objects are masked from ???????package:base???????: anyDuplicated, append, as.data.frame, as.vector, cbind, colnames, duplicated, eval, evalq, Filter, Find, get, intersect, is.unsorted, lapply, Map, mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce, rep.int, rownames, sapply, setdiff, sort, table, tapply, union, unique, unlist Loading required package: oligoClasses Welcome to oligoClasses version 1.24.0 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")'. ====================================================================== ====================================== Welcome to oligo version 1.26.0 ====================================================================== ====================================== Attaching package: ???????oligo??????? The following object is masked from ???????package:BiocGenerics???????: normalize The following object is masked from ???????package:limma???????: backgroundCorrect > Data<-read.celfiles(list.celfiles()) Loading required package: pd.hugene.2.0.st Loading required package: RSQLite Loading required package: DBI Platform design info loaded. Reading in : OC_MDAC1_(HuGene-2_0-st).CEL Reading in : OC_MDAC2_(HuGene-2_0-st).CEL Reading in : OC_MDAC3_(HuGene-2_0-st).CEL Reading in : OC_MDAP+T1_(HuGene-2_0-st).CEL Reading in : OC_MDAP+T2_(HuGene-2_0-st).CEL Reading in : OC_MDAP+T3_(HuGene-2_0-st).CEL Reading in : OC_MDAP1_(HuGene-2_0-st).CEL Reading in : OC_MDAP2_(HuGene-2_0-st).CEL Reading in : OC_MDAP3_(HuGene-2_0-st).CEL Reading in : OC_MDAT1_(HuGene-2_0-st).CEL Reading in : OC_MDAT2_(HuGene-2_0-st).CEL Reading in : OC_MDAT3_(HuGene-2_0-st).CEL > Data GeneFeatureSet (storageMode: lockedEnvironment) assayData: 2598544 features, 12 samples element names: exprs protocolData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: exprs dates varMetadata: labelDescription channel phenoData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: index varMetadata: labelDescription channel featureData: none experimentData: use 'experimentData(object)' Annotation: pd.hugene.2.0.st > eset<-rma(Data) Background correcting Normalizing Calculating Expression > eset ExpressionSet (storageMode: lockedEnvironment) assayData: 53617 features, 12 samples element names: exprs protocolData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: exprs dates varMetadata: labelDescription channel phenoData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: index varMetadata: labelDescription channel featureData: none experimentData: use 'experimentData(object)' Annotation: pd.hugene.2.0.st > library(affycoretools) Loading required package: affy Attaching package: ???????affy??????? The following objects are masked from ???????package:oligo???????: intensity, MAplot, mm, mm<-, mmindex, pm, pm<-, pmindex, probeNames, rma The following object is masked from ???????package:oligoClasses???????: list.celfiles Loading required package: GO.db Loading required package: AnnotationDbi Loading required package: KEGG.db KEGG.db contains mappings based on older data because the original resource was removed from the the public domain before the most recent update was produced. This package should now be considered deprecated and future versions of Bioconductor may not have it available. Users who want more current data are encouraged to look at the KEGGREST or reactome.db packages KernSmooth 2.23 loaded Copyright M. P. Wand 1997-2009 > eset_m <- getMainProbes(eset) > eset_m ExpressionSet (storageMode: lockedEnvironment) assayData: 46255 features, 12 samples element names: exprs protocolData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: exprs dates varMetadata: labelDescription channel phenoData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: index varMetadata: labelDescription channel featureData: none experimentData: use 'experimentData(object)' Annotation: pd.hugene.2.0.st ## getMainProbes subsets main (44629) but also normgene->exon (1626) probes from eset??? > library("hugene20sttranscriptcluster.db") Loading required package: org.Hs.eg.db > annotation(eset_m) <- "hugene20sttranscriptcluster.db" > eset_m ExpressionSet (storageMode: lockedEnvironment) assayData: 46255 features, 12 samples element names: exprs protocolData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: exprs dates varMetadata: labelDescription channel phenoData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: index varMetadata: labelDescription channel featureData: none experimentData: use 'experimentData(object)' Annotation: hugene20sttranscriptcluster.db > library(arrayQualityMetrics) > arrayQualityMetrics(expressionset=eset, outdir="QC normalized", force=TRUE) The directory 'QC normalized' has been created. Error in tmp[i] : invalid subscript type 'list' > hist(eset_m) > dev.copy2eps(file="eset_m.eps") quartz 2 > library(genefilter) > f1 = kOverA(3,2) > filt = filterfun(f1) > index = genefilter(eset_m,filt) > eset_filt = eset_m[index,] > dim(eset_filt) Features Samples 45290 12 > dim(eset_m) Features Samples 46255 12 > eset_filt ExpressionSet (storageMode: lockedEnvironment) assayData: 45290 features, 12 samples element names: exprs protocolData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: exprs dates varMetadata: labelDescription channel phenoData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: index varMetadata: labelDescription channel featureData: none experimentData: use 'experimentData(object)' Annotation: hugene20sttranscriptcluster.db > hist(eset_filt) > dev.copy2eps(file="eset_filt.eps") quartz 2 > pData(eset_filt) index OC_MDAC1_(HuGene-2_0-st).CEL 1 OC_MDAC2_(HuGene-2_0-st).CEL 2 OC_MDAC3_(HuGene-2_0-st).CEL 3 OC_MDAP+T1_(HuGene-2_0-st).CEL 4 OC_MDAP+T2_(HuGene-2_0-st).CEL 5 OC_MDAP+T3_(HuGene-2_0-st).CEL 6 OC_MDAP1_(HuGene-2_0-st).CEL 7 OC_MDAP2_(HuGene-2_0-st).CEL 8 OC_MDAP3_(HuGene-2_0-st).CEL 9 OC_MDAT1_(HuGene-2_0-st).CEL 10 OC_MDAT2_(HuGene-2_0-st).CEL 11 OC_MDAT3_(HuGene-2_0-st).CEL 12 > label <- read.delim("label.txt") > label array index label 1 OC_MDAC1_(HuGene-2_0-st).CEL 1 MDA_C 2 OC_MDAC2_(HuGene-2_0-st).CEL 2 MDA_C 3 OC_MDAC3_(HuGene-2_0-st).CEL 3 MDA_C 4 OC_MDAP+T1_(HuGene-2_0-st).CEL 4 MDA_PT 5 OC_MDAP+T2_(HuGene-2_0-st).CEL 5 MDA_PT 6 OC_MDAP+T3_(HuGene-2_0-st).CEL 6 MDA_PT 7 OC_MDAP1_(HuGene-2_0-st).CEL 7 MDA_P 8 OC_MDAP2_(HuGene-2_0-st).CEL 8 MDA_P 9 OC_MDAP3_(HuGene-2_0-st).CEL 9 MDA_P 10 OC_MDAT1_(HuGene-2_0-st).CEL 10 MDA_T 11 OC_MDAT2_(HuGene-2_0-st).CEL 11 MDA_T 12 OC_MDAT3_(HuGene-2_0-st).CEL 12 MDA_T > plotPCA(eset_filt, groups=rep(1:4, each=3), groupnames=unique(paste(label$label))) > dev.copy2eps(file="pca.eps") quartz 2 > plotPCA(eset_filt, groups=rep(1:4, each=3), groupnames=unique(paste(label$label)), screeplot=TRUE) > dev.copy2eps(file="pca_scree.eps") quartz 2 > plotPCA(eset_filt, groups=rep(1:4, each=3), groupnames=unique(paste(label$label)), squarepca=TRUE) > dev.copy2eps(file="pca_sq.eps") quartz 2 > plotPCA(eset_filt, groups=rep(1:4, each=3), groupnames=unique(paste(label$label)), squarepca=TRUE, plot3d=TRUE, pcs=c(1,2,3)) Sometimes rgl doesn't plot the first time. If there isn't anything in the plotting window, close it and re-run plotPCA(). > rgl.postscript("pca_3d.eps","eps",drawText=TRUE) > sampleNames(eset_filt) [1] "OC_MDAC1_(HuGene-2_0-st).CEL" "OC_MDAC2_(HuGene-2_0-st).CEL" "OC_MDAC3_(HuGene-2_0-st).CEL" [4] "OC_MDAP+T1_(HuGene-2_0-st).CEL" "OC_MDAP+T2_(HuGene-2_0-st).CEL" "OC_MDAP+T3_(HuGene-2_0-st).CEL" [7] "OC_MDAP1_(HuGene-2_0-st).CEL" "OC_MDAP2_(HuGene-2_0-st).CEL" "OC_MDAP3_(HuGene-2_0-st).CEL" [10] "OC_MDAT1_(HuGene-2_0-st).CEL" "OC_MDAT2_(HuGene-2_0-st).CEL" "OC_MDAT3_(HuGene-2_0-st).CEL" > cond.factor <- factor(c("MDAC","MDAC","MDAC","MDA_PT","MDA_PT","MDA_ PT","MDA_P","MDA_P","MDA_P","MDA_T","MDA_T","MDA_T"), levels = c("MDAC","MDA_PT","MDA_P","MDA_T")) > design <- model.matrix(~0+cond.factor) > colnames(design)<-c("MDAC","MDA_PT","MDA_P","MDA_T") > design MDAC MDA_PT MDA_P MDA_T 1 1 0 0 0 2 1 0 0 0 3 1 0 0 0 4 0 1 0 0 5 0 1 0 0 6 0 1 0 0 7 0 0 1 0 8 0 0 1 0 9 0 0 1 0 10 0 0 0 1 11 0 0 0 1 12 0 0 0 1 attr(,"assign") [1] 1 1 1 1 attr(,"contrasts") attr(,"contrasts")$cond.factor [1] "contr.treatment" > lm <- lmFit(eset_filt, design) > contrast.matrix <- makeContrasts(MDA_PT-MDAC, MDA_P-MDAC, MDA_T- MDAC, levels=design) > contrast.matrix Contrasts Levels MDA_PT - MDAC MDA_P - MDAC MDA_T - MDAC MDAC -1 -1 -1 MDA_PT 1 0 0 MDA_P 0 1 0 MDA_T 0 0 1 > fit2 <- contrasts.fit(lm, contrast.matrix) > eb <- eBayes(fit2) > results <- decideTests(eb) > vc <- vennCounts2(results, method ="same") > vennDiagram(vc, cex=0.8) > dev.copy2eps(file="venn.eps") quartz 2 > vennSelect(eset_filt,design,results,contrast.matrix,eb) > limma2annaffy(eset_filt,eb,design,contrast.matrix,annotation(eset_fi lt),pfilt=0.05,fldfilt=0.585,save=TRUE) You are going to output 3 tables, With this many genes in each: 2134 3415 275 Do you want to accept or change these values? [ a/c ] a Error in .checkKeysAreWellFormed(keys) : keys must be supplied in a character vector with no NAs ## gene numbers are ok but function does not work out further??? > limma2annaffy(eset_filt,eb,design,contrast.matrix,annotation(eset_fi lt),number=10) You are going to output 3 tables, With this many genes in each: 10 10 10 Do you want to accept or change these values? [ a/c ] a Error in .checkKeysAreWellFormed(keys) : keys must be supplied in a character vector with no NAs > limma2annaffy(eset_filt,eb,design,contrast.matrix,annotation(eset_fi lt),number=10) You are going to output 3 tables, With this many genes in each: 10 10 10 Do you want to accept or change these values? [ a/c ] c Do you want to change the p-value? [ y/n ] n Do you want to change the fold filter? [ y/n ] n You are going to output 3 tables, With this many genes in each: 30 30 30 Do you want to accept or change these values? [ a/c ] a Error in .checkKeysAreWellFormed(keys) : keys must be supplied in a character vector with no NAs ## error with gene number modification (10 -> 30 (default setting)) and function??? > fStat <- topTableF(eb, number=Inf) > head(fStat) MDA_PT...MDAC MDA_P...MDAC MDA_T...MDAC AveExpr F P.Value adj.P.Val 16859205 -1.286312 -4.2941076 0.01920662 9.660757 810.8735 1.102616e-14 4.993749e-10 16672452 -3.395777 -3.6027376 -3.30340265 5.170242 719.3912 2.336551e-14 5.291120e-10 16853028 5.559876 0.7256685 0.12395099 4.229877 662.0674 3.932584e-14 5.936891e-10 16833204 4.791469 0.4702840 0.61199749 6.136889 567.6096 1.031520e-13 1.167939e-09 17052685 -4.434038 -1.9167135 -1.44742803 8.153717 511.1014 1.988363e-13 1.801059e-09 17067314 -3.330336 -3.0812685 -0.35449045 5.620444 390.7517 1.064119e-12 7.861126e-09 > write.table(fStat,file="Fstat.txt",sep = "\t", row.names = TRUE) > dim(fStat) [1] 45290 7 > dim(eset_filt) Features Samples 45290 12 > MDA_PT_sig <- topTable(eb, coef = 1, number = Inf, p.value = 0.05,lfc=0.585) > dim(MDA_PT_sig) #2134 [1] 2134 6 > head(MDA_PT_sig) logFC AveExpr t P.Value adj.P.Val B 17052685 -4.434038 8.153717 -38.39125 1.882904e-14 3.930583e-10 20.70057 16853028 5.559876 4.229877 38.06801 2.093170e-14 3.930583e-10 20.64453 16672452 -3.395777 5.170242 -37.41013 2.603610e-14 3.930583e-10 20.52750 16833204 4.791469 6.136889 36.18503 3.949201e-14 4.471483e-10 20.29845 16743091 2.715139 7.201918 27.99236 9.712218e-13 8.797327e-09 18.29467 17095194 2.547995 5.724030 27.07882 1.466936e-12 1.002695e-08 18.00646 > MDA_PT_sig$tcid <- rownames(MDA_PT_sig) > head(MDA_PT_sig) logFC AveExpr t P.Value adj.P.Val B tcid 17052685 -4.434038 8.153717 -38.39125 1.882904e-14 3.930583e-10 20.70057 17052685 16853028 5.559876 4.229877 38.06801 2.093170e-14 3.930583e-10 20.64453 16853028 16672452 -3.395777 5.170242 -37.41013 2.603610e-14 3.930583e-10 20.52750 16672452 16833204 4.791469 6.136889 36.18503 3.949201e-14 4.471483e-10 20.29845 16833204 16743091 2.715139 7.201918 27.99236 9.712218e-13 8.797327e-09 18.29467 16743091 17095194 2.547995 5.724030 27.07882 1.466936e-12 1.002695e-08 18.00646 17095194 > library(annaffy) > aaf.handler() [1] "Probe" "Symbol" "Description" "Chromosome" [5] "Chromosome Location" "GenBank" "Gene" "Cytoband" [9] "UniGene" "PubMed" "Gene Ontology" "Pathway" > anncols <- aaf.handler()[c(1:3,6:7,9:12)] > anntable <- aafTableAnn(MDA_PT_sig$tcid, "hugene20sttranscriptcluster.db", anncols) > head(MDA_PT_sig) logFC AveExpr t P.Value adj.P.Val B tcid 17052685 -4.434038 8.153717 -38.39125 1.882904e-14 3.930583e-10 20.70057 17052685 16853028 5.559876 4.229877 38.06801 2.093170e-14 3.930583e-10 20.64453 16853028 16672452 -3.395777 5.170242 -37.41013 2.603610e-14 3.930583e-10 20.52750 16672452 16833204 4.791469 6.136889 36.18503 3.949201e-14 4.471483e-10 20.29845 16833204 16743091 2.715139 7.201918 27.99236 9.712218e-13 8.797327e-09 18.29467 16743091 17095194 2.547995 5.724030 27.07882 1.466936e-12 1.002695e-08 18.00646 17095194 > class(MDA_PT_sig) [1] "data.frame" > testtable <- aafTable(logFC = round(MDA_PT_sig$logFC,2), tstat=round(MDA_PT_sig$t,1),adjPval =signif(MDA_PT_sig$adj.P.Val,1), probeids= MDA_PT_sig$tcid) > table <- merge(anntable,testtable) > exprtable <- aafTableInt(eset_filt,probeids=MDA_PT_sig$tcid) > tableF <- merge(table,exprtable) > head(tableF) An object of class "aafTable" Slot "probeids": [1] "17052685" Slot "table": $Probe An object of class "aafList" [[1]] An object of class "aafProbe" [1] "17052685" $Symbol An object of class "aafList" [[1]] character(0) attr(,"class") [1] "aafSymbol" $Description An object of class "aafList" [[1]] character(0) attr(,"class") [1] "aafDescription" $GenBank An object of class "aafList" [[1]] character(0) attr(,"class") [1] "aafGenBank" $Gene An object of class "aafList" [[1]] An object of class "aafLocusLink" integer(0) $UniGene An object of class "aafList" [[1]] character(0) attr(,"class") [1] "aafUniGene" $PubMed An object of class "aafList" [[1]] An object of class "aafPubMed" integer(0) $`Gene Ontology` An object of class "aafList" [[1]] An object of class "aafGO" list() $Pathway An object of class "aafList" [[1]] An object of class "aafPathway" list() $logFC An object of class "aafList" [[1]] [1] -4.43 $tstat An object of class "aafList" [[1]] [1] -38.4 $adjPval An object of class "aafList" [[1]] [1] 4e-10 $`OC_MDAC1_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 9.961219 $`OC_MDAC2_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 10.22801 $`OC_MDAC3_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 10.12056 $`OC_MDAP+T1_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 5.455098 $`OC_MDAP+T2_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 5.849167 $`OC_MDAP+T3_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 5.703406 $`OC_MDAP1_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 8.109914 $`OC_MDAP2_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 8.325627 $`OC_MDAP3_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 8.124103 $`OC_MDAT1_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 8.66443 $`OC_MDAT2_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 8.719271 $`OC_MDAT3_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 8.5838 > saveHTML(tableF,"tableF.html",title="MDA_PT pval<0.05 fc>1.5", open=TRUE) > saveText(tableF,"table_MDA_PT.txt") ## idem for the other two contrast... ##volcanos and MA... ##heatmap... > sessionInfo() R version 3.0.2 (2013-09-25) Platform: x86_64-apple-darwin10.8.0 (64-bit) locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] parallel stats graphics grDevices utils datasets methods base other attached packages: [1] gplots_2.12.1 annaffy_1.34.0 [3] rgl_0.93.963 genefilter_1.44.0 [5] arrayQualityMetrics_3.18.0 hugene20sttranscriptcluster.db_2.13.0 [7] org.Hs.eg.db_2.10.1 affycoretools_1.34.0 [9] KEGG.db_2.10.1 GO.db_2.10.1 [11] AnnotationDbi_1.24.0 affy_1.40.0 [13] pd.hugene.2.0.st_3.8.0 RSQLite_0.11.4 [15] DBI_0.2-7 oligo_1.26.0 [17] Biobase_2.22.0 oligoClasses_1.24.0 [19] BiocGenerics_0.8.0 limma_3.18.7 loaded via a namespace (and not attached): [1] affxparser_1.34.0 affyio_1.30.0 affyPLM_1.38.0 annotate_1.40.0 [5] AnnotationForge_1.4.4 beadarray_2.12.0 BeadDataPackR_1.14.0 BiocInstaller_1.12.0 [9] biomaRt_2.18.0 Biostrings_2.30.1 biovizBase_1.10.7 bit_1.1-11 [13] bitops_1.0-6 BSgenome_1.30.0 Cairo_1.5-5 Category_2.28.0 [17] caTools_1.16 cluster_1.14.4 codetools_0.2-8 colorspace_1.2-4 [21] DESeq2_1.2.8 dichromat_2.0-0 digest_0.6.4 edgeR_3.4.2 [25] ff_2.2-12 foreach_1.4.1 Formula_1.1-1 gcrma_2.34.0 [29] gdata_2.13.2 GenomicFeatures_1.14.2 GenomicRanges_1.14.4 ggbio_1.10.10 [33] ggplot2_0.9.3.1 GOstats_2.28.0 graph_1.40.1 grid_3.0.2 [37] gridExtra_0.9.1 GSEABase_1.24.0 gtable_0.1.2 gtools_3.1.1 [41] Hmisc_3.13-0 hwriter_1.3 IRanges_1.20.6 iterators_1.0.6 [45] KernSmooth_2.23-10 labeling_0.2 lattice_0.20-24 latticeExtra_0.6-26 [49] locfit_1.5-9.1 MASS_7.3-29 Matrix_1.1-0 munsell_0.4.2 [53] PFAM.db_2.10.1 plyr_1.8 preprocessCore_1.24.0 proto_0.3-10 [57] R.methodsS3_1.5.2 R.oo_1.15.8 R.utils_1.28.4 R2HTML_2.2.1 [61] RBGL_1.38.0 RColorBrewer_1.0-5 Rcpp_0.10.6 RcppArmadillo_0.3.930.1 [65] RCurl_1.95-4.1 ReportingTools_2.2.0 reshape2_1.2.2 Rsamtools_1.14.2 [69] rtracklayer_1.22.0 scales_0.2.3 setRNG_2011.11-2 splines_3.0.2 [73] stats4_3.0.2 stringr_0.6.2 survival_2.37-4 SVGAnnotation_0.93-1 [77] tools_3.0.2 VariantAnnotation_1.8.8 vsn_3.30.0 XML_3.95-0.2 [81] xtable_1.7-1 XVector_0.2.0 zlibbioc_1.8.0 -- Sent via the guest posting facility at bioconductor.org.
ADD COMMENTlink modified 4.6 years ago by James W. MacDonald46k • written 4.6 years ago by Guest User12k
0
gravatar for James W. MacDonald
4.6 years ago by
United States
James W. MacDonald46k wrote:
Hi Jose, I am not in the office this week, so I don't have easy access to a computer running R. I don't see anything obviously wrong with your code, and as you note you can get results using annaffy directly. Internally limma2annaffy is just getting an index of significant probesets and then getting the IDs by probeids <- featureNames(eset_filt)[index] So you might check to see that you don't have any NA values for the probeset IDs. Not sure if that is the problem, but worth checking I suppose. Lately I have been working away from annaffy and towards ReportingTools, which is still actively maintained and a bit nicer for the end user. It's not that fully integrated, but see for example ?vennPage for html Venn diagrams with clickable links. To replace limma2annaffy I use the internal function tableFilt. Something like library(ReportingTools) eb$genes = select(hugene20sttranscriptcluster.db, feature names(eset_filt), c("SYMBOL","ENTREZID","GENENAME")) #add more if you like tab = lapply(1:ncol(contrast), function(x) affycoretools:::tableFilt(eb, x, pfilt=0.05, fldfilt=0.585)) out =lapply(seq(along=tab), function(x) { htab = HTMLReport(gsub(" ", "_", colnames(contrast)[x]), colnames(contrast)[x], reportDirectory = "./reports") publish(tab[x], htab, .modifyDF = list(affyLinks, entrezLinks)) finish(htab) #you can use write.table here to output text version return(htab) }) Note this is not tested and was typed on a tablet, so there may well be some autocorrect errors in there. Then indx = HTMLReport("index", "Put a header here") publish(Link(out, indx), indx) finish(indx) Which will give you a nice HTML page with links to each table. The vennPage function is also intended to be passed in via the publish() step in that last code block, so you can have a single index page that acts as a point of entry. Also note that ReportingTools will also accept an MarrayLM object as input and will create an HTML page. See the vignettes for that package. I just don't really like the stock output you get, so have rolled my own. Best, Jim Dear Jim, I am performing differential expression analysis using affymetrix human gene 2.0 ST arrays. I get an error message from limma2annaffy() function. Error in .checkKeysAreWellFormed(keys) : keys must be supplied in a character vector with no NAs I get the annotated HTML and txt tables for each desired contrast by using topTable and annaffy tools but I would like to know what is going wrong with limma2annaffy since I find this function truly useful and convenient. Thank you in advance for your time and your help, Jose -- output of sessionInfo(): R version 3.0.2 (2013-09-25) -- "Frisbee Sailing" Copyright (C) 2013 The R Foundation for Statistical Computing Platform: x86_64-apple-darwin10.8.0 (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. [R.app GUI 1.62 (6558) x86_64-apple-darwin10.8.0] [History restored from /Users/Jose/.Rapp.history] > setwd("/Users/Jose/Documents/Oscar/MDA_261213/MDA_271213") > library(limma) > library(oligo) Loading required package: BiocGenerics Loading required package: parallel Attaching package: ‘BiocGenerics’ The following objects are masked from ‘package:parallel’: clusterApply, clusterApplyLB, clusterCall, clusterEvalQ, clusterExport, clusterMap, parApply, parCapply, parLapply, parLapplyLB, parRapply, parSapply, parSapplyLB The following object is masked from ‘package:limma’: plotMA The following object is masked from ‘package:stats’: xtabs The following objects are masked from ‘package:base’: anyDuplicated, append, as.data.frame, as.vector, cbind, colnames, duplicated, eval, evalq, Filter, Find, get, intersect, is.unsorted, lapply, Map, mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce, rep.int, rownames, sapply, setdiff, sort, table, tapply, union, unique, unlist Loading required package: oligoClasses Welcome to oligoClasses version 1.24.0 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")'. ====================================================================== ====================================== Welcome to oligo version 1.26.0 ====================================================================== ====================================== Attaching package: ‘oligo’ The following object is masked from ‘package:BiocGenerics’: normalize The following object is masked from ‘package:limma’: backgroundCorrect > Data<-read.celfiles(list.celfiles()) Loading required package: pd.hugene.2.0.st Loading required package: RSQLite Loading required package: DBI Platform design info loaded. Reading in : OC_MDAC1_(HuGene-2_0-st).CEL Reading in : OC_MDAC2_(HuGene-2_0-st).CEL Reading in : OC_MDAC3_(HuGene-2_0-st).CEL Reading in : OC_MDAP+T1_(HuGene-2_0-st).CEL Reading in : OC_MDAP+T2_(HuGene-2_0-st).CEL Reading in : OC_MDAP+T3_(HuGene-2_0-st).CEL Reading in : OC_MDAP1_(HuGene-2_0-st).CEL Reading in : OC_MDAP2_(HuGene-2_0-st).CEL Reading in : OC_MDAP3_(HuGene-2_0-st).CEL Reading in : OC_MDAT1_(HuGene-2_0-st).CEL Reading in : OC_MDAT2_(HuGene-2_0-st).CEL Reading in : OC_MDAT3_(HuGene-2_0-st).CEL > Data GeneFeatureSet (storageMode: lockedEnvironment) assayData: 2598544 features, 12 samples element names: exprs protocolData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: exprs dates varMetadata: labelDescription channel phenoData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: index varMetadata: labelDescription channel featureData: none experimentData: use 'experimentData(object)' Annotation: pd.hugene.2.0.st > eset<-rma(Data) Background correcting Normalizing Calculating Expression > eset ExpressionSet (storageMode: lockedEnvironment) assayData: 53617 features, 12 samples element names: exprs protocolData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: exprs dates varMetadata: labelDescription channel phenoData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: index varMetadata: labelDescription channel featureData: none experimentData: use 'experimentData(object)' Annotation: pd.hugene.2.0.st > library(affycoretools) Loading required package: affy Attaching package: ‘affy’ The following objects are masked from ‘package:oligo’: intensity, MAplot, mm, mm<-, mmindex, pm, pm<-, pmindex, probeNames, rma The following object is masked from ‘package:oligoClasses’: list.celfiles Loading required package: GO.db Loading required package: AnnotationDbi Loading required package: KEGG.db KEGG.db contains mappings based on older data because the original resource was removed from the the public domain before the most recent update was produced. This package should now be considered deprecated and future versions of Bioconductor may not have it available. Users who want more current data are encouraged to look at the KEGGREST or reactome.db packages KernSmooth 2.23 loaded Copyright M. P. Wand 1997-2009 > eset_m <- getMainProbes(eset) > eset_m ExpressionSet (storageMode: lockedEnvironment) assayData: 46255 features, 12 samples element names: exprs protocolData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: exprs dates varMetadata: labelDescription channel phenoData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: index varMetadata: labelDescription channel featureData: none experimentData: use 'experimentData(object)' Annotation: pd.hugene.2.0.st ## getMainProbes subsets main (44629) but also normgene->exon (1626) probes from eset > library("hugene20sttranscriptcluster.db") Loading required package: org.Hs.eg.db > annotation(eset_m) <- "hugene20sttranscriptcluster.db" > eset_m ExpressionSet (storageMode: lockedEnvironment) assayData: 46255 features, 12 samples element names: exprs protocolData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: exprs dates varMetadata: labelDescription channel phenoData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: index varMetadata: labelDescription channel featureData: none experimentData: use 'experimentData(object)' Annotation: hugene20sttranscriptcluster.db > library(arrayQualityMetrics) > arrayQualityMetrics(expressionset=eset, outdir="QC normalized", force=TRUE) The directory 'QC normalized' has been created. Error in tmp[i] : invalid subscript type 'list' > hist(eset_m) > dev.copy2eps(file="eset_m.eps") quartz 2 > library(genefilter) > f1 = kOverA(3,2) > filt = filterfun(f1) > index = genefilter(eset_m,filt) > eset_filt = eset_m[index,] > dim(eset_filt) Features Samples 45290 12 > dim(eset_m) Features Samples 46255 12 > eset_filt ExpressionSet (storageMode: lockedEnvironment) assayData: 45290 features, 12 samples element names: exprs protocolData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: exprs dates varMetadata: labelDescription channel phenoData rowNames: OC_MDAC1_(HuGene-2_0-st).CEL OC_MDAC2_(HuGene-2_0-st).CEL ... OC_MDAT3_(HuGene-2_0-st).CEL (12 total) varLabels: index varMetadata: labelDescription channel featureData: none experimentData: use 'experimentData(object)' Annotation: hugene20sttranscriptcluster.db > hist(eset_filt) > dev.copy2eps(file="eset_filt.eps") quartz 2 > pData(eset_filt) index OC_MDAC1_(HuGene-2_0-st).CEL 1 OC_MDAC2_(HuGene-2_0-st).CEL 2 OC_MDAC3_(HuGene-2_0-st).CEL 3 OC_MDAP+T1_(HuGene-2_0-st).CEL 4 OC_MDAP+T2_(HuGene-2_0-st).CEL 5 OC_MDAP+T3_(HuGene-2_0-st).CEL 6 OC_MDAP1_(HuGene-2_0-st).CEL 7 OC_MDAP2_(HuGene-2_0-st).CEL 8 OC_MDAP3_(HuGene-2_0-st).CEL 9 OC_MDAT1_(HuGene-2_0-st).CEL 10 OC_MDAT2_(HuGene-2_0-st).CEL 11 OC_MDAT3_(HuGene-2_0-st).CEL 12 > label <- read.delim("label.txt") > label array index label 1 OC_MDAC1_(HuGene-2_0-st).CEL 1 MDA_C 2 OC_MDAC2_(HuGene-2_0-st).CEL 2 MDA_C 3 OC_MDAC3_(HuGene-2_0-st).CEL 3 MDA_C 4 OC_MDAP+T1_(HuGene-2_0-st).CEL 4 MDA_PT 5 OC_MDAP+T2_(HuGene-2_0-st).CEL 5 MDA_PT 6 OC_MDAP+T3_(HuGene-2_0-st).CEL 6 MDA_PT 7 OC_MDAP1_(HuGene-2_0-st).CEL 7 MDA_P 8 OC_MDAP2_(HuGene-2_0-st).CEL 8 MDA_P 9 OC_MDAP3_(HuGene-2_0-st).CEL 9 MDA_P 10 OC_MDAT1_(HuGene-2_0-st).CEL 10 MDA_T 11 OC_MDAT2_(HuGene-2_0-st).CEL 11 MDA_T 12 OC_MDAT3_(HuGene-2_0-st).CEL 12 MDA_T > plotPCA(eset_filt, groups=rep(1:4, each=3), groupnames=unique(paste(label$label))) > dev.copy2eps(file="pca.eps") quartz 2 > plotPCA(eset_filt, groups=rep(1:4, each=3), groupnames=unique(paste(label$label)), screeplot=TRUE) > dev.copy2eps(file="pca_scree.eps") quartz 2 > plotPCA(eset_filt, groups=rep(1:4, each=3), groupnames=unique(paste(label$label)), squarepca=TRUE) > dev.copy2eps(file="pca_sq.eps") quartz 2 > plotPCA(eset_filt, groups=rep(1:4, each=3), groupnames=unique(paste(label$label)), squarepca=TRUE, plot3d=TRUE, pcs=c(1,2,3)) Sometimes rgl doesn't plot the first time. If there isn't anything in the plotting window, close it and re-run plotPCA(). > rgl.postscript("pca_3d.eps","eps",drawText=TRUE) > sampleNames(eset_filt) [1] "OC_MDAC1_(HuGene-2_0-st).CEL" "OC_MDAC2_(HuGene-2_0-st).CEL" "OC_MDAC3_(HuGene-2_0-st).CEL" [4] "OC_MDAP+T1_(HuGene-2_0-st).CEL" "OC_MDAP+T2_(HuGene-2_0-st).CEL" "OC_MDAP+T3_(HuGene-2_0-st).CEL" [7] "OC_MDAP1_(HuGene-2_0-st).CEL" "OC_MDAP2_(HuGene-2_0-st).CEL" "OC_MDAP3_(HuGene-2_0-st).CEL" [10] "OC_MDAT1_(HuGene-2_0-st).CEL" "OC_MDAT2_(HuGene-2_0-st).CEL" "OC_MDAT3_(HuGene-2_0-st).CEL" > cond.factor <- factor(c("MDAC","MDAC","MDAC","MDA_PT","MDA_PT","MDA_PT","MDA_P","MDA_ P","MDA_P","MDA_T","MDA_T","MDA_T"), levels = c("MDAC","MDA_PT","MDA_P","MDA_T")) > design <- model.matrix(~0+cond.factor) > colnames(design)<-c("MDAC","MDA_PT","MDA_P","MDA_T") > design MDAC MDA_PT MDA_P MDA_T 1 1 0 0 0 2 1 0 0 0 3 1 0 0 0 4 0 1 0 0 5 0 1 0 0 6 0 1 0 0 7 0 0 1 0 8 0 0 1 0 9 0 0 1 0 10 0 0 0 1 11 0 0 0 1 12 0 0 0 1 attr(,"assign") [1] 1 1 1 1 attr(,"contrasts") attr(,"contrasts")$cond.factor [1] "contr.treatment" > lm <- lmFit(eset_filt, design) > contrast.matrix <- makeContrasts(MDA_PT-MDAC, MDA_P-MDAC, MDA_T- MDAC, levels=design) > contrast.matrix Contrasts Levels MDA_PT - MDAC MDA_P - MDAC MDA_T - MDAC MDAC -1 -1 -1 MDA_PT 1 0 0 MDA_P 0 1 0 MDA_T 0 0 1 > fit2 <- contrasts.fit(lm, contrast.matrix) > eb <- eBayes(fit2) > results <- decideTests(eb) > vc <- vennCounts2(results, method ="same") > vennDiagram(vc, cex=0.8) > dev.copy2eps(file="venn.eps") quartz 2 > vennSelect(eset_filt,design,results,contrast.matrix,eb) > limma2annaffy(eset_filt,eb,design,contrast.matrix,annotation(eset_filt ),pfilt=0.05,fldfilt=0.585,save=TRUE) You are going to output 3 tables, With this many genes in each: 2134 3415 275 Do you want to accept or change these values? [ a/c ] a Error in .checkKeysAreWellFormed(keys) : keys must be supplied in a character vector with no NAs ## gene numbers are ok but function does not work out further > limma2annaffy(eset_filt,eb,design,contrast.matrix,annotation(eset_filt ),number=10) You are going to output 3 tables, With this many genes in each: 10 10 10 Do you want to accept or change these values? [ a/c ] a Error in .checkKeysAreWellFormed(keys) : keys must be supplied in a character vector with no NAs > limma2annaffy(eset_filt,eb,design,contrast.matrix,annotation(eset_filt ),number=10) You are going to output 3 tables, With this many genes in each: 10 10 10 Do you want to accept or change these values? [ a/c ] c Do you want to change the p-value? [ y/n ] n Do you want to change the fold filter? [ y/n ] n You are going to output 3 tables, With this many genes in each: 30 30 30 Do you want to accept or change these values? [ a/c ] a Error in .checkKeysAreWellFormed(keys) : keys must be supplied in a character vector with no NAs ## error with gene number modification (10 -> 30 (default setting)) and function > fStat <- topTableF(eb, number=Inf) > head(fStat) MDA_PT...MDAC MDA_P...MDAC MDA_T...MDAC AveExpr F P.Value adj.P.Val 16859205 -1.286312 -4.2941076 0.01920662 9.660757 810.8735 1.102616e-14 4.993749e-10 16672452 -3.395777 -3.6027376 -3.30340265 5.170242 719.3912 2.336551e-14 5.291120e-10 16853028 5.559876 0.7256685 0.12395099 4.229877 662.0674 3.932584e-14 5.936891e-10 16833204 4.791469 0.4702840 0.61199749 6.136889 567.6096 1.031520e-13 1.167939e-09 17052685 -4.434038 -1.9167135 -1.44742803 8.153717 511.1014 1.988363e-13 1.801059e-09 17067314 -3.330336 -3.0812685 -0.35449045 5.620444 390.7517 1.064119e-12 7.861126e-09 > write.table(fStat,file="Fstat.txt",sep = "\t", row.names = TRUE) > dim(fStat) [1] 45290 7 > dim(eset_filt) Features Samples 45290 12 > MDA_PT_sig <- topTable(eb, coef = 1, number = Inf, p.value = 0.05,lfc=0.585) > dim(MDA_PT_sig) #2134 [1] 2134 6 > head(MDA_PT_sig) logFC AveExpr t P.Value adj.P.Val B 17052685 -4.434038 8.153717 -38.39125 1.882904e-14 3.930583e-10 20.70057 16853028 5.559876 4.229877 38.06801 2.093170e-14 3.930583e-10 20.64453 16672452 -3.395777 5.170242 -37.41013 2.603610e-14 3.930583e-10 20.52750 16833204 4.791469 6.136889 36.18503 3.949201e-14 4.471483e-10 20.29845 16743091 2.715139 7.201918 27.99236 9.712218e-13 8.797327e-09 18.29467 17095194 2.547995 5.724030 27.07882 1.466936e-12 1.002695e-08 18.00646 > MDA_PT_sig$tcid <- rownames(MDA_PT_sig) > head(MDA_PT_sig) logFC AveExpr t P.Value adj.P.Val B tcid 17052685 -4.434038 8.153717 -38.39125 1.882904e-14 3.930583e-10 20.70057 17052685 16853028 5.559876 4.229877 38.06801 2.093170e-14 3.930583e-10 20.64453 16853028 16672452 -3.395777 5.170242 -37.41013 2.603610e-14 3.930583e-10 20.52750 16672452 16833204 4.791469 6.136889 36.18503 3.949201e-14 4.471483e-10 20.29845 16833204 16743091 2.715139 7.201918 27.99236 9.712218e-13 8.797327e-09 18.29467 16743091 17095194 2.547995 5.724030 27.07882 1.466936e-12 1.002695e-08 18.00646 17095194 > library(annaffy) > aaf.handler() [1] "Probe" "Symbol" "Description" "Chromosome" [5] "Chromosome Location" "GenBank" "Gene" "Cytoband" [9] "UniGene" "PubMed" "Gene Ontology" "Pathway" > anncols <- aaf.handler()[c(1:3,6:7,9:12)] > anntable <- aafTableAnn(MDA_PT_sig$tcid, "hugene20sttranscriptcluster.db", anncols) > head(MDA_PT_sig) logFC AveExpr t P.Value adj.P.Val B tcid 17052685 -4.434038 8.153717 -38.39125 1.882904e-14 3.930583e-10 20.70057 17052685 16853028 5.559876 4.229877 38.06801 2.093170e-14 3.930583e-10 20.64453 16853028 16672452 -3.395777 5.170242 -37.41013 2.603610e-14 3.930583e-10 20.52750 16672452 16833204 4.791469 6.136889 36.18503 3.949201e-14 4.471483e-10 20.29845 16833204 16743091 2.715139 7.201918 27.99236 9.712218e-13 8.797327e-09 18.29467 16743091 17095194 2.547995 5.724030 27.07882 1.466936e-12 1.002695e-08 18.00646 17095194 > class(MDA_PT_sig) [1] "data.frame" > testtable <- aafTable(logFC = round(MDA_PT_sig$logFC,2), tstat=round(MDA_PT_sig$t,1),adjPval =signif(MDA_PT_sig$adj.P.Val,1), probeids= MDA_PT_sig$tcid) > table <- merge(anntable,testtable) > exprtable <- aafTableInt(eset_filt,probeids=MDA_PT_sig$tcid) > tableF <- merge(table,exprtable) > head(tableF) An object of class "aafTable" Slot "probeids": [1] "17052685" Slot "table": $Probe An object of class "aafList" [[1]] An object of class "aafProbe" [1] "17052685" $Symbol An object of class "aafList" [[1]] character(0) attr(,"class") [1] "aafSymbol" $Description An object of class "aafList" [[1]] character(0) attr(,"class") [1] "aafDescription" $GenBank An object of class "aafList" [[1]] character(0) attr(,"class") [1] "aafGenBank" $Gene An object of class "aafList" [[1]] An object of class "aafLocusLink" integer(0) $UniGene An object of class "aafList" [[1]] character(0) attr(,"class") [1] "aafUniGene" $PubMed An object of class "aafList" [[1]] An object of class "aafPubMed" integer(0) $`Gene Ontology` An object of class "aafList" [[1]] An object of class "aafGO" list() $Pathway An object of class "aafList" [[1]] An object of class "aafPathway" list() $logFC An object of class "aafList" [[1]] [1] -4.43 $tstat An object of class "aafList" [[1]] [1] -38.4 $adjPval An object of class "aafList" [[1]] [1] 4e-10 $`OC_MDAC1_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 9.961219 $`OC_MDAC2_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 10.22801 $`OC_MDAC3_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 10.12056 $`OC_MDAP+T1_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 5.455098 $`OC_MDAP+T2_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 5.849167 $`OC_MDAP+T3_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 5.703406 $`OC_MDAP1_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 8.109914 $`OC_MDAP2_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 8.325627 $`OC_MDAP3_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 8.124103 $`OC_MDAT1_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 8.66443 $`OC_MDAT2_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 8.719271 $`OC_MDAT3_(HuGene-2_0-st).CEL` An object of class "aafList" [[1]] An object of class "aafIntensity" [1] 8.5838 > saveHTML(tableF,"tableF.html",title="MDA_PT pval<0.05 fc>1.5", open=TRUE) > saveText(tableF,"table_MDA_PT.txt") ## idem for the other two contrast... ##volcanos and MA... ##heatmap... > sessionInfo() R version 3.0.2 (2013-09-25) Platform: x86_64-apple-darwin10.8.0 (64-bit) locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] parallel stats graphics grDevices utils datasets methods base other attached packages: [1] gplots_2.12.1 annaffy_1.34.0 [3] rgl_0.93.963 genefilter_1.44.0 [5] arrayQualityMetrics_3.18.0 hugene20sttranscriptcluster.db_2.13.0 [7] org.Hs.eg.db_2.10.1 affycoretools_1.34.0 [9] KEGG.db_2.10.1 GO.db_2.10.1 [11] AnnotationDbi_1.24.0 affy_1.40.0 [13] pd.hugene.2.0.st_3.8.0 RSQLite_0.11.4 [15] DBI_0.2-7 oligo_1.26.0 [17] Biobase_2.22.0 oligoClasses_1.24.0 [19] BiocGenerics_0.8.0 limma_3.18.7 loaded via a namespace (and not attached): [1] affxparser_1.34.0 affyio_1.30.0 affyPLM_1.38.0 annotate_1.40.0 [5] AnnotationForge_1.4.4 beadarray_2.12.0 BeadDataPackR_1.14.0 BiocInstaller_1.12.0 [9] biomaRt_2.18.0 Biostrings_2.30.1 biovizBase_1.10.7 bit_1.1-11 [13] bitops_1.0-6 BSgenome_1.30.0 Cairo_1.5-5 Category_2.28.0 [17] caTools_1.16 cluster_1.14.4 codetools_0.2-8 colorspace_1.2-4 [21] DESeq2_1.2.8 dichromat_2.0-0 digest_0.6.4 edgeR_3.4.2 [25] ff_2.2-12 foreach_1.4.1 Formula_1.1-1 gcrma_2.34.0 [29] gdata_2.13.2 GenomicFeatures_1.14.2 GenomicRanges_1.14.4 ggbio_1.10.10 [33] ggplot2_0.9.3.1 GOstats_2.28.0 graph_1.40.1 grid_3.0.2 [37] gridExtra_0.9.1 GSEABase_1.24.0 gtable_0.1.2 gtools_3.1.1 [41] Hmisc_3.13-0 hwriter_1.3 IRanges_1.20.6 iterators_1.0.6 [45] KernSmooth_2.23-10 labeling_0.2 lattice_0.20-24 latticeExtra_0.6-26 [49] locfit_1.5-9.1 MASS_7.3-29 Matrix_1.1-0 munsell_0.4.2 [53] PFAM.db_2.10.1 plyr_1.8 preprocessCore_1.24.0 proto_0.3-10 [57] R.methodsS3_1.5.2 R.oo_1.15.8 R.utils_1.28.4 R2HTML_2.2.1 [61] RBGL_1.38.0 RColorBrewer_1.0-5 Rcpp_0.10.6 RcppArmadillo_0.3.930.1 [65] RCurl_1.95-4.1 ReportingTools_2.2.0 reshape2_1.2.2 Rsamtools_1.14.2 [69] rtracklayer_1.22.0 scales_0.2.3 setRNG_2011.11-2 splines_3.0.2 [73] stats4_3.0.2 stringr_0.6.2 survival_2.37-4 SVGAnnotation_0.93-1 [77] tools_3.0.2 VariantAnnotation_1.8.8 vsn_3.30.0 XML_3.95-0.2 [81] xtable_1.7-1 XVector_0.2.0 zlibbioc_1.8.0 -- Sent via the guest posting facility at bioconductor.org. [[alternative HTML version deleted]]
ADD COMMENTlink written 4.6 years ago by James W. MacDonald46k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 249 users visited in the last hour