[BioC] help with RankProd plaase
Kimpel, Mark William
mkimpel at iupui.edu
Thu Nov 16 17:47:46 CET 2006
I am a new user of the RankProd package and am having difficulty getting it to work with my data, whereas I have no problem with the golub data used in the vignette. The RankProd functions run fine and do not produce errors with my data, but I get no significant genes, whereas the same data with Limma and SAM produce over a thousand sig. genes at FDR<0.1. I can't help but think I am making a silly mistake in setting up the cl or origen vectors for RankProd or perhaps in something else.
Below is my sessionInfo() and script for both Limma and RankProd. Any help would be appreciated.
Thanks,
Mark
sessionInfo()
R version 2.4.0 (2006-10-03)
i386-pc-mingw32
locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
attached base packages:
[1] "splines" "stats" "graphics" "grDevices" "datasets" "utils" "methods" "tools" "base"
other attached packages:
rgu34acdf svMisc multtest affycoretools biomaRt RCurl XML GOstats
"1.14.0" "0.9-5" "1.12.0" "1.6.0" "1.8.0" "0.7-0" "0.99-93" "2.0.3"
Category genefilter survival KEGG RBGL annotate GO graph
"2.0.3" "1.12.0" "2.29" "1.14.1" "1.10.0" "1.12.0" "1.14.1" "1.12.0"
RankProd RWinEdt limma affy affyio Biobase
"2.6.0" "1.7-5" "2.9.1" "1.12.1" "1.2.0" "1.12.2"
#####################################
#The Limma method
require(limma)
factor.mat<-factor(pData(affy.object.preprocessed.BB01$eSet)$Treatment) #affy.object.preprocessed.BB01$eSet is the
# expressionSet created with justRMA
treatments <- factor(factor.mat) #Treatment has two levels, P and NP
design <- model.matrix(~ -1+treatments)
colNames <- unique(factor.mat); o<-order(colNames); colNames<-colNames[o]; colnames(design)<-colNames
fit <- lmFit(exprs(affy.object.preprocessed.BB01$eSet),des=design)
########################################################################
#Make contrast matrix
contrast <-makeContrasts(
(P - NP),
levels=design)
###########################################################################
fit2 <- contrasts.fit(fit,contrast)
fit2<-eBayes(fit2)
rslt <- decideTests(fit2, method="separate", adjust.method="BH",p.value=0.1)
topTable(fit2,coef=1,number=10,genelist=fit$genes,adjust.method="BH",sort.by="B",resort.by=NULL)
#########################################################################
#The RankProd method
RP.out<-RPadvance(data=exprs(affy.object.preprocessed.BB01$eSet),
cl=c(as.numeric(unclass(factor(pData(affy.object.preprocessed.BB01$eSet)$Treatment))) - 1),
origin=rep(1, ncol(exprs(affy.object.preprocessed.BB01$eSet))),
num.perm=100,
logged=TRUE,
na.rm=FALSE,
gene.names=NULL,
plot=TRUE,
rand=NULL)
plotRP(RP.out, cutoff=0.1)
topGene(RP.out,cutoff = 0.1, logged = TRUE, logbase = 2, gene.names=affy.object.preprocessed.BB01$annotations[,4])
Mark W. Kimpel MD
Official Business Address:
Department of Psychiatry
Indiana University School of Medicine
PR M116
Institute of Psychiatric Research
791 Union Drive
Indianapolis, IN 46202
Preferred Mailing Address:
15032 Hunter Court
Westfield, IN 46074
(317) 490-5129 Work, & Mobile
(317) 663-0513 Home (no voice mail please)
1-(317)-536-2730 FAX
More information about the Bioconductor
mailing list