[BioC] a question on SPIA package

[Jing Huang]

Dear Adi and other members,

I am trying to apply my microarray data to SPIA R package that is created by Adi Tarca. I only got one pathway. I am wondering if I did anything wrong. Could you help me on this?

I attached my topTable of fit data that contains lists of genes that is differently expressed, LFC , p.value, adj.p.value and so on.

Here is my R:

> library(lattice)
> m <- exprs(eset)
> A <- unname(rowMeans(m))
> M <- as.vector(m) -A
> Sample <- colnames(m)[col(m)]
> df <- data.frame(A=A, M=M, Sample=Sample)
> treatments=factor(c(1,1,1,2,2,2,3,3,3,4,4,4), labels=c("CTRL","HIF1a","HIF2a","HIF1a2a"))
> design=model.matrix(~treatments)
> colnames(design)=c("CTRL","HIF1a","HIF2a","HIF1a2a")
> design
   CTRL HIF1a HIF2a HIF1a2a
1     1     0     0       0
2     1     0     0       0
3     1     0     0       0
4     1     1     0       0
5     1     1     0       0
6     1     1     0       0
7     1     0     1       0
8     1     0     1       0
9     1     0     1       0
10    1     0     0       1
11    1     0     0       1
12    1     0     0       1
attr(,"assign")
[1] 0 1 1 1
attr(,"contrasts")
attr(,"contrasts")$treatments
[1] "contr.treatment"

> fit=lmFit(eset,design)
> fit=eBayes(fit)
> results=classifyTestsF(fit, p.value=0.001)
> summary(results)
    CTRL HIF1a HIF2a HIF1a2a
-1     0   440    84     577
0      0 50437 53683   50649
1  54675  3798   908    3449
> options(digits=3)
> colnames(topTable(fit,coef="HIF1a",n=20))
[1] "ID"                    "Gene.title"            "Gene.symbol"
 [4] "Gene.ID"               "UniGene.title"         "UniGene.symbol"
 [7] "UniGene.ID"            "Nucleotide.Title"      "GI"
[10] "GenBank.Accession"     "Platform_CLONEID"      "Platform_ORF"
[13] "Platform_SPOTID"       "Chromosome.location"   "Chromosome.annotation"
[16] "GO.Function"           "GO.Process"            "GO.Component"
[19] "GO.Function.1"         "GO.Process.1"          "GO.Component.1"
[22] "logFC"                 "AveExpr"               "t"
[25] "P.Value"               "adj.P.Val"             "B"
> library(SPIA)
> x=topTable(fit,coef="HIF1a")
> library(hgu133plus2.db)
Loading required package: AnnotationDbi
Loading required package: org.Hs.eg.db
Loading required package: DBI


> y=hgu133plus2ENTREZID
> x$ENTREZ=unlist(as.list(y[x$ID]))
> x=x[!is.na(x$ENTREZ),]
> x=x[!duplicated(x$ENTREZ),]
> tg1=x[x$adj.P.Val<0.1,]
> HIF1a=tg1$logFC
> names(HIF1a)=as.vector(tg1$ENTREZ)
> HIF1a1=x$ENTREZ
> res=spia(de=HIF1a,all=HIF1a1,organism="hsa",nB=2000,plots=F,beta=NULL,combine="fisher",verbose=F)
> res
           Name    ID pSize NDE pNDE tA pPERT pG pGFdr pGFWER    Status
1 RNA transport 03013     1   1    1  0    NA  1     1      1 Inhibited
                                                    KEGGLINK
1 http://www.genome.jp/dbget-bin/show_pathway?hsa03013+54913


Many Many Thanks

Jing Huang PhD

Oregon Health $ Sciences University

Portland Oregon