Dear R-Mac-List,
I find that a script I worked did not
work on the command line R but worked when I
pasted the same script into the RGui or sourced it
from the GUI. Although this is a Bioconductor script,
I suspect that it is a R-for-MAC problem, either in how
I am using R for Mac or if there is a problem with the
this implementation of R.
I would greatly appreciate any help it running scripts
like the following from the command line.
Here is the script:
library(affy)
library(annaffy)
library(gcrma)
library(limma)
library(hgu133plus2.db)
raw<-ReadAffy()
gcrmanm<-gcrma(raw,fast=FALSE)
#write.exprs(gcrmanm,"thirdbatchexprs.txt")
#allready done
targets<-readTargets("Targets.txt")
targets
pressure<-factor(targets$pressure,levels=c("ctrl","tr"))
pressure
patient<-factor(targets$patient)
patient
design<-model.matrix(~pressure+patient)
design
arrayw<-arrayWeights(gcrmanm, design, method = "genebygene", maxiter = 50, tol
= 1e-10)
fit<-lmFit(gcrmanm,design,weights=arrayw)
fit<-eBayes(fit)
comp1<-topTable(fit,coef="pressuretr",number=54675,adjust.method="BH",sort.by="B")
names(comp1)
names(comp1)[1]="Probe"
probeids<-comp1[,1]
aaf.handler()
anncols<-aaf.handler()[c(1:5,7:8,11:12)]
anncols
anntable<-aafTableAnn(probeids, "hgu133plus2.db", colnames=anncols)
comp1table<-aafTableFrame(comp1,names(comp1),probeids=probeids)
comp1ann<-merge(comp1table,anntable)
saveText(comp1ann,"trVSctrl.txt")
Here is my terminal command:
ccrfml1:cels friedman$ R CMD BATCH limpresallSC.R
The file did not run to completion.
Here is my error message:
ccrfml1:cels friedman$ cat *.Rout
###########################################################
R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-apple-darwin9.8.0/x86_64 (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.
Error: object geneNames is not exported by 'namespace:Biobase'
Execution halted
###########################################################
ccrfml1:cels friedman$
Here is a record of my session pasted into the GUI:
R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i386-apple-darwin9.8.0/i386 (32-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.53 (6335) i386-apple-darwin9.8.0]
[Workspace restored from /Users/friedman/.RData]
[History restored from /Users/friedman/.Rapp.history]
> sessionInfo()
R version 2.15.2 (2012-10-26)
Platform: i386-apple-darwin9.8.0/i386 (32-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] stats graphics grDevices utils datasets methods base
> library(affy)
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")'.
> library(annaffy)
Loading required package: GO.db
Loading required package: AnnotationDbi
Loading required package: DBI
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. One possible alternative to consider is to look at the
reactome.db package
> library(gcrma)
> library(limma)
> library(hgu133plus2.db)
Loading required package: org.Hs.eg.db
> sessionInfo()
R version 2.15.2 (2012-10-26)
Platform: i386-apple-darwin9.8.0/i386 (32-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] stats graphics grDevices utils datasets methods base
other attached packages:
[1] hgu133plus2.db_2.8.0 org.Hs.eg.db_2.8.0 limma_3.12.1
gcrma_2.30.0
[5] annaffy_1.28.0 KEGG.db_2.7.1 GO.db_2.7.1
RSQLite_0.11.1
[9] DBI_0.2-5 AnnotationDbi_1.20.2 affy_1.36.0
Biobase_2.18.0
[13] BiocGenerics_0.4.0
loaded via a namespace (and not attached):
[1] affyio_1.24.0 BiocInstaller_1.8.3 Biostrings_2.26.2
IRanges_1.16.4
[5] parallel_2.15.2 preprocessCore_1.18.0 splines_2.15.2
stats4_2.15.2
[9] zlibbioc_1.2.0
> raw<-ReadAffy()
Error in AllButCelsForReadAffy(..., filenames = filenames, widget = widget, :
No cel filennames specified and no cel files in specified
directory:/Users/friedman
> raw<-ReadAffy()
> gcrmanm<-gcrma(raw,fast=FALSE)
Adjusting for optical effect........................Done.
Computing affinities.Done.
Adjusting for non-specific binding........................Done.
Normalizing
Calculating Expression
>
>
> targets<-readTargets("Targets.txt")
> targets
FileName Name patient pressure
1 01_1.CEL 01_1 1 ctrl
2 01_2.CEL 01_2 1 tr
3 02_1.CEL 02_1 2 ctrl
4 02_1.CEL 02_1 2 tr
5 03_1.CEL 03_1 3 ctrl
6 03_2.CEL 03_2 3 tr
7 04_1.CEL 04_1 4 ctrl
8 04_2.CEL 04_2 4 tr
9 05_1.CEL 05_1 5 ctrl
10 05_2.CEL 05_2 5 tr
11 06_1.CEL 06_1 6 ctrl
12 06_2.CEL 06_2 6 tr
13 07_1.CEL 07_1 7 ctrl
14 07_2.CEL 07_2 7 tr
15 08_1.CEL 08_1 8 ctrl
16 08_2.CEL 08_2 8 tr
17 09_1.CEL 09_1 9 ctrl
18 09_2.CEL 09_2 9 tr
19 10_1.CEL 10_1 10 ctrl
20 10_2.CEL 10_2 10 tr
21 11_1.CEL 11_1 11 ctrl
22 11_2.CEL 11_2 11 tr
23 12_1.CEL 12_1 12 ctrl
24 12_2.CEL 12_2 12 tr
> pressure<-factor(targets$pressure,levels=c("ctrl","tr"))
> pressure
[1] ctrl tr ctrl tr ctrl tr ctrl tr ctrl tr ctrl tr ctrl tr ctrl
tr ctrl tr ctrl
[20] tr ctrl tr ctrl tr
Levels: ctrl tr
> patient<-factor(targets$patient)
> patient
[1] 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10 10 11 11 12 12
Levels: 1 2 3 4 5 6 7 8 9 10 11 12
>
> design<-model.matrix(~pressure+patient)
> design
(Intercept) pressuretr patient2 patient3 patient4 patient5 patient6 patient7
patient8 patient9
1 1 0 0 0 0 0 0 0
0 0
2 1 1 0 0 0 0 0 0
0 0
3 1 0 1 0 0 0 0 0
0 0
4 1 1 1 0 0 0 0 0
0 0
5 1 0 0 1 0 0 0 0
0 0
6 1 1 0 1 0 0 0 0
0 0
7 1 0 0 0 1 0 0 0
0 0
8 1 1 0 0 1 0 0 0
0 0
9 1 0 0 0 0 1 0 0
0 0
10 1 1 0 0 0 1 0 0
0 0
11 1 0 0 0 0 0 1 0
0 0
12 1 1 0 0 0 0 1 0
0 0
13 1 0 0 0 0 0 0 1
0 0
14 1 1 0 0 0 0 0 1
0 0
15 1 0 0 0 0 0 0 0
1 0
16 1 1 0 0 0 0 0 0
1 0
17 1 0 0 0 0 0 0 0
0 1
18 1 1 0 0 0 0 0 0
0 1
19 1 0 0 0 0 0 0 0
0 0
20 1 1 0 0 0 0 0 0
0 0
21 1 0 0 0 0 0 0 0
0 0
22 1 1 0 0 0 0 0 0
0 0
23 1 0 0 0 0 0 0 0
0 0
24 1 1 0 0 0 0 0 0
0 0
patient10 patient11 patient12
1 0 0 0
2 0 0 0
3 0 0 0
4 0 0 0
5 0 0 0
6 0 0 0
7 0 0 0
8 0 0 0
9 0 0 0
10 0 0 0
11 0 0 0
12 0 0 0
13 0 0 0
14 0 0 0
15 0 0 0
16 0 0 0
17 0 0 0
18 0 0 0
19 1 0 0
20 1 0 0
21 0 1 0
22 0 1 0
23 0 0 1
24 0 0 1
attr(,"assign")
[1] 0 1 2 2 2 2 2 2 2 2 2 2 2
attr(,"contrasts")
attr(,"contrasts")$pressure
[1] "contr.treatment"
attr(,"contrasts")$patient
[1] "contr.treatment"
> arrayw<-arrayWeights(gcrmanm, design, method = "genebygene", maxiter = 50,
> tol = 1e-10)
> fit<-lmFit(gcrmanm,design,weights=arrayw)
> fit<-eBayes(fit)
>
> comp1<-topTable(fit,coef="pressuretr",number=54675,adjust.method="BH",sort.by="B")
> names(comp1)
[1] "ID" "logFC" "AveExpr" "t" "P.Value" "adj.P.Val" "B"
> names(comp1)[1]="Probe"
>
> probeids<-comp1[,1]
>
> aaf.handler()
[1] "Probe" "Symbol" "Description"
"Chromosome"
[5] "Chromosome Location" "GenBank" "Gene"
"Cytoband"
[9] "UniGene" "PubMed" "Gene Ontology"
"Pathway"
>
> anncols<-aaf.handler()[c(1:5,7:8,11:12)]
> anncols
[1] "Probe" "Symbol" "Description"
"Chromosome"
[5] "Chromosome Location" "Gene" "Cytoband" "Gene
Ontology"
[9] "Pathway"
>
> anntable<-aafTableAnn(probeids, "hgu133plus2.db", colnames=anncols)
> comp1table<-aafTableFrame(comp1,names(comp1),probeids=probeids)
> comp1ann<-merge(comp1table,anntable)
> saveText(comp1ann,"trVSctrl.txt")
There were 50 or more warnings (use warnings() to see the first 50)
###########################################################
It ran to completion very nicely.
The script also ran when I sourced it from the GUI.
I would appreciate any suggestions that could help me run the script
from the command line,
Richard A. Friedman, PhD
Associate Research Scientist,
Biomedical Informatics Shared Resource
Herbert Irving Comprehensive Cancer Center (HICCC)
Lecturer,
Department of Biomedical Informatics (DBMI)
Educational Coordinator,
Center for Computational Biology and Bioinformatics (C2B2)/
National Center for Multiscale Analysis of Genomic Networks (MAGNet)/
Columbia Initiative in Systems Biology
Room 824
Irving Cancer Research Center
Columbia University
1130 St. Nicholas Ave
New York, NY 10032
(212)851-4765 (voice)
[email protected]
http://cancercenter.columbia.edu/~friedman/
In memoriam, Ray Bradbury
[[alternative HTML version deleted]]
_______________________________________________
R-SIG-Mac mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-mac
