Entering edit mode
Marcos Pinho
▴
200
@marcos-pinho-3584
Last seen 9.6 years ago
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]]
