[BioC] SAM analysis : samr package vs siggenes package

[Cécile Laurent]

Hi,

I would like to know why when I run SAM analysis with SamR and Siggenes 
packages, I don't have the same results for the small pvalues.
I test here on golub data, and I don't have the same differential genes 
number when I observe the delta tables (delta.table for samr and 
sam.out at mat.fdr for siggenes) for a FDR 0% (250 are differential with 
samr package and 9 with siggenes package).
On my personnal data, I observe the same difference in differential gene 
number at FDR 5%.
I also observe the same s0, but not the same pi0, in the 2 sam objects.
I think, in my code, the parameters are the same... So, I don't know 
which package used...
Is it something wrong???
Can you explain me these differences in the results..

Thanks,
Cécile


Here is my sessionInfo()

R version 2.6.0 (2007-10-03)
i486-pc-linux-gnu

locale:
LC_CTYPE=fr_FR.UTF-8;LC_NUMERIC=C;LC_TIME=fr_FR.UTF-8;LC_COLLATE=fr_FR.UTF-8;LC_MONETARY=fr_FR.UTF-8;LC_MESSAGES=fr_FR.UTF-8;LC_PAPER=fr_FR.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=fr_FR.UTF-8;LC_IDENTIFICATION=C

attached base packages:
[1] splines   tools     stats     graphics  grDevices utils     datasets
[8] methods   base    

other attached packages:
[1] siggenes_1.12.0 samr_1.25       impute_1.0-5    multtest_1.18.0
[5] survival_2.33   Biobase_1.16.0

loaded via a namespace (and not attached):
[1] rcompgen_0.1-15


### My code :
library(multtest)
library(samr)
library(siggenes)
data(golub)
golubCl=golub.cl
golubCl[which(golub.cl==1)]=2
golubCl[which(golub.cl==0)]=1

## SamR package
samrData=list(x=golub, y=golubCl, geneid=golub.gnames[,3], 
genesnames=golub.gnames[,3], logged2=T)
samr.obj=samr(samrData, resp.type="Two class unpaired", 
testStatistic="standard", nperms=1000, random.seed=123)
delta.table=samr.compute.delta.table(samr.obj, dels=seq(0.1,5,0.05))
pv.samr=samr.pvalues.from.perms(samr.obj$tt, samr.obj$ttstar)

## Siggenes package
sam.out=sam(golub, golub.cl, rand=123,B=1000, 
gene.names=golub.gnames[,3], method="d.stat", var.equal=T, s0=NA, 
include.zero=F, delta=seq(0.1,5,0.05) )

## plot to observe the values (the version with -log() to see the 
difference in the small pvalues)
#plot(sam.out at d[order(sam.out at d)],samr.obj$tt[order(samr.obj$tt)])
#plot(-log(sam.out at d[order(sam.out at d)]),-log(samr.obj$tt[order(samr.obj$tt)]))
plot(sam.out at p.value[order(sam.out at p.value)],pv.samr[order(pv.samr)])
plot(-log(sam.out at p.value[order(sam.out at p.value)]),-log(pv.samr[order(pv.samr)]))