[BioC] Questions about multtest

[michael watson (IAH-C)]

OK, I using the multtest package to analyse my data, following the instructions in multtest.pdf.

I run:

> t <- mt.teststat(data[,6:12], c(0,0,0,1,1,1,1), test="t")

which calculates the t statistic for my data.  The t statistic for my first gene comes up as:

> t[1]
[1] 40.60158

Presumably, this is equivalent to me running t.test:

> t.test(data[1,9:12], data[1,6:8], var.equal=FALSE, alternative="two.sided")

        Welch Two Sample t-test

data:  data[1, 9:12] and data[1, 6:8]
t = 40.6016, df = 2, p-value = 0.0006061
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 1.713804 2.120092
sample estimates:
    mean of x     mean of y
-1.596190e-15 -1.916948e+00

So I want to know how I can get p-values for the t statistics I have just calculated using mt.teststat.

This is where I get confused - multtest.pdf says I should "compute raw nominal two-sided p-values for the 3,051 test statistics using the standard Gaussian distribution":

> rawp0 <- 2 * (1 - pnorm(abs(t)))

However, rawp0:

> rawp0[1]
[1] 0

So according to this procedure, my raw (unadjusted) p-value is zero, whereas according to t.test(), my p-value is 0.0006061.

Soldiering on, I want to calculate adjusted p-values accoridng to Benjamini and Hochberg:

>res <- mt.rawp2adjp(rawp0, "BH")
>adjp <- res$adjp[order(res$index), ]
>adjp[1]
[1] 0

Does this mean that my adjusted p-value from the Benjamini and Hochberg procedure for my first gene is zero?  Haven't I lost some information along the way when 0.0006061 was converted to just plain old zero?  How do I get accurate raw p-values from the multtest package?

Or am I doing something wrong?

Thanks
Mick