[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