Hi Holger,

Thanks for the response (I am sorry to answer so late as I was on a business
trip).

I will look into the description to try and understand it better.

I already thought that the number is too low. I am still not sure though
about the first two runs.

What about the R.fold parameter?
Does it make sense to set it before I run the SAM analysis, or is it better
to first run it with all the genes and than exclude those genes with a fold
induction under a certain threshold?

In the second run I get almost the same amount of called genes with a
"better" FDR value.
Is there a way to say which FDR value can be accepted?

Thanks

Assa


On Tue, Jun 28, 2011 at 13:11, Holger Schwender <holger.schw@gmx.de> wrote:

> Hi Assa,
>
> siggenes is described in
>
> http://cran.r-project.org/doc/Rnews/Rnews_2006-5.pdf
>
> on pages 45ff, and with a bit more technical details in
>
> https://eldorado.tu-dortmund.de/bitstream/2003/23306/1/diss_schwender.pdf
>
> Section 6.3.2. Another possibility is to look into the code for sam, in
> particular into the functions d.stat and siggenes:::stats.cal.
>
> The reason why the FDR is so low in the third example is the very, very low
> estimate for p0. Take a look at the number of falsely called genes (False)
> for about the same number of called genes. These numbers are close to each
> other. The FDR is estimated by
>
> p0 * False/Called
>
> so the FDR gets very low if p0 is very small.
>
> I would thus not trust the results of the third analysis, as the p0
> estimation has not worked properly here. Not really sure why, but the reason
> might be that too many null genes are removed when use.dm=FALSE.
>
> Best,
> Holger
>
> -------- Original-Nachricht --------
> > Datum: Tue, 28 Jun 2011 09:41:53 +0200
> > Von: Assa Yeroslaviz <frymor@gmail.com>
> > An: bioconductor <bioconductor@stat.math.ethz.ch>, holger.schw@gmx.de
> > Betreff: siggenes parameters
>
> > Hi everybody,
> >
> > I have a question about the behavior of siggenes in R (2.12) with
> > different
> > parameters.
> > Maybe I need to understand better how siggenes calculated the
> > differentially
> > regulated genes, but I would like to ask you for your help in
> > understanding
> > it.
> >
> > Maybe you can also direct me to where I can find the answers for this
> > problem(s)
> >
> > I am running a sam analysis of wild type vs. mutant in drosophila miRNA
> > microarrays.
> >
> > I first read the files into an affybatch and than normalized them (RMA).
> > After normalizing the arrays I extracted all non-drosophila miRNA from my
> > expressionSet and ran the sam analysis only with the 186 dme-miRNA
> probes.
> >
> > cl <- c(0,0,0,1,1,1) # "wt1", "wt2", "wt3", "mt1", "mt2", "mt3"
> > sam.out <- sam(dme_rma,cl, var.equal=FALSE, B=100, include.zero=FALSE,
> > gene.names=Probe_names, na.replace=FALSE, rand=123)
> >
> > I run the sam analysis with different parameters and got *very* different
> > results.
> >
> > the first run was with the default parameters. Here I got  these results:
> > SAM Analysis for the Two-Class Unpaired Case Assuming Unequal Variances
> >
> >  s0 = 0.3176  (The 55 % quantile of the s values.)
> >
> >    Delta    p0 False Called   FDR cutlow cutup j2  j1
> > 1   0.10 0.882 74.85     92 0.717 -0.130 2.135 91 186
> > 2   0.11 0.882 74.85     92 0.717 -0.130 2.135 91 186
> > 3   0.12 0.882 28.25     39 0.639 -0.575 2.135 38 186
> > 4   0.13 0.882 24.35     38 0.565 -0.633 2.135 37 186
> > ...
> > 10  0.19 0.882 21.65     37 0.516 -0.672 2.135 36 186
> > 11  0.20 0.882  12.4     24 0.456 -0.860 2.135 23 186
> > 12  0.21 0.882  12.4     24 0.456 -0.860 2.135 23 186
> > 13  0.22 0.882  12.4     24 0.456 -0.860 2.135 23 186
> > 14  0.23 0.882  10.9     22 0.437 -0.905 2.135 21 186
> > 15 0.24 0.882 0.4      1 0.353   -Inf 2.135  0 186
> > ...
> > 20  0.29 0.882   0.4      1 0.353   -Inf 2.135  0 186
> > 21  0.30 0.882     0      0     0   -Inf   Inf  0 187
> > 22  0.31 0.882     0      0     0   -Inf   Inf  0 187
> >
> > I get here a very high FDR value for very few miRNAs.
> >
> > For my second run I added the parameter R.fold = 1.2, to exclude all
> > miRNAs
> > under this value from the analysis.
> > Here I get much better results with lower FDR values:
> >
> >  Number of variables having a fold change >= 1.2 or <= 0.8333 : 90
> >
> >  s0 = 0.0477  (The 0 % quantile of the s values.)
> >
> >    Delta    p0 False Called    FDR cutlow cutup j2 j1
> > 1    0.1 0.378  55.5     82 0.2557 -0.708 0.425 45 54
> > 2    0.2 0.378 52.25     80 0.2467 -0.708 0.562 45 56
> > 3    0.3 0.378  48.4     79 0.2314 -0.708 0.676 45 57
> > 4    0.5 0.378  24.9     48 0.1960 -0.708 2.504 45 88
> > 5    0.6 0.378  23.7     45 0.1990 -0.708   Inf 45 91
> > 6    0.7 0.378 10.65     31 0.1298 -1.170   Inf 31 91
> > 7    0.8 0.378  6.75     25 0.1020 -1.513   Inf 25 91
> > 8    1.0 0.378  0.85      4 0.0803 -2.712   Inf  4 91
> > 9    1.1 0.378  0.25      1 0.0944 -4.810   Inf  1 91
> > 10   1.2 0.378     0      0      0   -Inf   Inf  0 91
> >
> > I than ran a third sam with the parameter use.dm=FALSE. This gave me
> again
> > different results with much lower FDR values:
> >
> >  Number of variables having a fold change >= 1.2 or <= 0.8333 : 97
> >
> >  s0 = 0.1488  (The 5 % quantile of the s values.)
> >
> >    Delta    p0 False Called     FDR cutlow cutup j2 j1
> > 1    0.1 0.009  53.5     81 0.00612 -0.541 0.545 48 65
> > 2    0.2 0.009 50.65     80 0.00587 -0.541 0.616 48 66
> > 3    0.3 0.009 26.85     48 0.00519 -0.541   Inf 48 98
> > 4    0.4 0.009 26.85     48 0.00519 -0.541   Inf 48 98
> > 5    0.5 0.009 26.85     48 0.00519 -0.541   Inf 48 98
> > 6    0.6 0.009     0      0       0   -Inf   Inf  0 98
> > 7    0.7 0.009     0      0       0   -Inf   Inf  0 98
> > 8    0.8 0.009     0      0       0   -Inf   Inf  0 98
> > 9    0.9 0.009     0      0       0   -Inf   Inf  0 98
> > 10   1.0 0.009     0      0       0   -Inf   Inf  0 98
> >
> >
> > As you can see my obvious problem.  Which is the better (right???) way of
> > running this analysis.
> > I think, to understand that, one needs to know how SAM works. I don't
> > understand why I get such different results in the FDR between the three
> > possibilities.
> >
> > As far as I understand it, the lower the FDR value, the better. But how
> > can
> > it be, that I have in the third run over 80 miRNAs with such a low FDR,
> > while I have only few DE miRNAs in the first run with a much higher FDR
> > value?
> >
> > I would appreciate any help,
> >
> > thanks
> > Assa
>
> --
> Empfehlen Sie GMX DSL Ihren Freunden und Bekannten und wir
> belohnen Sie mit bis zu 50,- Euro! https://freundschaftswerbung.gmx.de
>

	[[alternative HTML version deleted]]

