[BioC] SPIA package

[Seth Falcon]

See below...

On 3/17/10 10:52 AM, Marcos Pinho wrote:
> 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]))

At this point, you should inspect fit2$ID and determine its type using 
class.

And you might inspect what you get from:

   as.list(x[head(fit2$ID)])

to try and understand why the result of the above in your session is 
returning NULL.

Also, with all of your output, I couldn't find sessionInfo() which might 
be useful in helping you any further.

+ seth

>
>> fit2$ENTREZ
>
> NULL

...

-- 
Seth Falcon
Bioconductor Core Team | FHCRC