Dear list,

I have been trying to use the SPIA package, but have been encountering some
dificulties in creating my two vectors: one with the log2 fold changes of DE
genes and the other cobntaining the universe of all Entrez gene IDs on the
array. Any help would be greatly apreciated. Please find below my session
info:

R version 2.10.1 (2009-12-14)

Copyright (C) 2009 The R Foundation for Statistical Computing

ISBN 3-900051-07-0



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.



[Previously saved workspace restored]



> library (SPIA)

Loading required package: RCurl

Loading required package: bitops

> data(colorectalcancer)

> DE_Colorectal[1:10]

     3491      2353      1958      1843      3725     23645      9510     84869


 5.960206  5.143502  4.148081  2.429889  1.531126  1.429269  3.937663
-1.147077

     7432      1490

 4.715767  3.447604

> ALL_Colorectal[1:10]

 [1] "3491"  "2353"  "1958"  "1843"  "3725"  "23645" "9510"  "84869" "7432"

[10] "1490"



> head(top)

               ID    logFC  AveExpr        t      P.Value    adj.P.Val
B

10738   201289_at 5.960206 6.226928 23.94388 1.789221e-17 9.782565e-13
25.40124

18604   209189_at 5.143502 7.487305 17.42995 1.560212e-14 2.843486e-10
21.02120

11143 201694_s_at 4.148081 7.038281 16.46040 5.153117e-14 7.043667e-10
20.14963

10490 201041_s_at 2.429889 9.594413 14.06891 1.293706e-12 1.414667e-08
17.66883

10913 201464_x_at 1.531126 8.221044 10.98634 1.685519e-10 1.151947e-06
13.61189

11463   202014_at 1.429269 5.327647 10.45906 4.274251e-10 2.418771e-06
12.80131

> dir()

 [1] "array image K562 Lucena sem VCR.jpeg"

 [2] "Box plot norm data histogram.jpeg"

 [3] "Box plot raw data histogram.jpeg"

 [4] "K562 1.CEL"

 [5] "K562 2_1.CEL"

 [6] "K562 vs Lucena-VCR"

 [7] "limma_completeK562 vs Lucena.xls"

 [8] "Lucena sem VCR 1.CEL"

 [9] "Lucena sem VCR 2.CEL"

[10] "MAplot Norm data.jpeg"

[11] "MAplot Raw data.jpeg"

[12] "QC Stats Graph K562 vc Lucena sem VCR.jpeg"

[13] "RLE NUSE plots.jpeg"

[14] "RNA degradation plot.jpeg"

> library(affy)

Loading required package: Biobase



Welcome to Bioconductor



  Vignettes contain introductory material. To view, type

  'openVignette()'. To cite Bioconductor, see

  'citation("Biobase")' and for packages 'citation(pkgname)'.



> library(tkWidgets)

Loading required package: widgetTools

Loading required package: tcltk

Loading Tcl/Tk interface ... done

Loading required package: DynDoc

Loading required package: tools

> data=ReadAffy(widget=TRUE)

> library(gcrma)

> eset=gcrma(data)

Adjusting for optical effect....Done.

Computing affinitiesLoading required package: AnnotationDbi

.Done.

Adjusting for non-specific binding....Done.

Normalizing

Calculating Expression

> library(genefilter)

> library (hgu133plus2.db)

Loading required package: org.Hs.eg.db

Loading required package: DBI

> esetF = nsFilter (eset, require.entrez=TRUE,remove.dupEntrez=TRUE,
feature.exclude="^AFFX",var.cutof=0.5)$eset

> design = model.matrix (~factor(rep (1:2, each=2)))

> colnames(design)=c("K562", "Lucena")

> design

  K562 Lucena

1    1      0

2    1      0

3    1      1

4    1      1

attr(,"assign")

[1] 0 1

attr(,"contrasts")

attr(,"contrasts")$`factor(rep(1:2, each = 2))`

[1] "contr.treatment"



> library(limma)

> fit =lmFit (esetF, design)

> fit2=eBayes(fit)

> topTable(fit2, coef=2)

              ID     logFC  AveExpr         t      P.Value   adj.P.Val
B

5026   209993_at  8.167898 6.410285  26.34450 7.032016e-07 0.005525598
5.852191

3236  1553436_at  7.512443 6.097913  23.91930 1.175771e-06 0.005525598
5.586671

9476 206488_s_at  7.318355 6.081049  21.61440 2.014549e-06 0.005525598
5.277576

1342   235683_at  6.589716 5.620394  19.19109 3.784627e-06 0.005525598
4.875214

6740 222392_x_at  6.047929 6.736659  18.95527 4.040667e-06 0.005525598
4.830968

8639   210603_at  6.088776 5.773605  18.94361 4.053843e-06 0.005525598
4.828755

3821   202948_at -6.231328 5.580407 -18.55760 4.520484e-06 0.005525598
4.754058

5074   205934_at  5.836512 5.809832  18.26090 4.922783e-06 0.005525598
4.694726

3870 205786_s_at -6.788839 7.957034 -17.55018 6.072079e-06 0.005525598
4.545431

8207   225502_at -6.769984 5.721122 -17.11499 6.933010e-06 0.005525598
4.448722

> library(annotate)

> fit2$genes$EG <- getEG(fit2$genes$ID, "hgu133plus2")

> topTable(fit2, coef=2)

              ID     EG     logFC  AveExpr         t      P.Value
adj.P.Val

5026   209993_at   5243  8.167898 6.410285  26.34450 7.032016e-07
0.005525598

3236  1553436_at 283463  7.512443 6.097913  23.91930 1.175771e-06
0.005525598

9476 206488_s_at    948  7.318355 6.081049  21.61440 2.014549e-06
0.005525598

1342   235683_at 143686  6.589716 5.620394  19.19109 3.784627e-06
0.005525598

6740 222392_x_at  64065  6.047929 6.736659  18.95527 4.040667e-06
0.005525598

8639   210603_at  84779  6.088776 5.773605  18.94361 4.053843e-06
0.005525598

3821   202948_at   3554 -6.231328 5.580407 -18.55760 4.520484e-06
0.005525598

5074   205934_at   5334  5.836512 5.809832  18.26090 4.922783e-06
0.005525598

3870 205786_s_at   3684 -6.788839 7.957034 -17.55018 6.072079e-06
0.005525598

8207   225502_at  81704 -6.769984 5.721122 -17.11499 6.933010e-06
0.005525598

            B

5026 5.852191

3236 5.586671

9476 5.277576

1342 4.875214

6740 4.830968

8639 4.828755

3821 4.754058

5074 4.694726

3870 4.545431

8207 4.448722



> x = hgu133plus2ENTREZID

> fit2$ENTREZ = unlist (as.list (x[fit2$ID]))

> fit2$ENTREZ

NULL

> fit2 = fit2 [! is. na (fit2$ENTREZ),]

Error: unexpected symbol in "fit2 = fit2 [! is. na"

> tg1 = fit2[fit2$adj.P.Val < 0.05,]

> DE_KL = tg1$logFC

> names(DE_KL)= as.vector(tg1$ENTREZ)

> ALL_DE = fit2$ENTREZ

> DE_KL [1:10]

NULL

> DE_KL

NULL

-- 
Marcos B. Pinho
Programa de Engenharia Química - PEQ
Laboratório de Engenharia de Cultivos Celulares- LECC
Universidade Federal do Rio de Janeiro - UFRJ
Instituto Nacional de Câncer - INCA
Rio de Janeiro - Brasil

	[[alternative HTML version deleted]]

