Dear Michael,

Thanks very much.

Jiang
On Tue, Oct 4, 2011 at 2:39 PM, R. Michael Weylandt <
michael.weylandt@gmail.com> wrote:

> On Mon, Oct 3, 2011 at 10:08 PM, chunjiang he <camelbbs@gmail.com> wrote:
> > Hi,
> >
> > I have two questions want to ask.
> >
> > 1. If I have a matrix like this, and I want to figure out the rows whose
> > value in the 3rd column are less than 0.05. How can I do it with R.
> > hsa-let-7a--MBTD1    0.528239197    2.41E-05
> > hsa-let-7a--APOBEC1    0.507869409    5.51E-05
> > hsa-let-7a--PAPOLA    0.470451884    0.000221774
> > hsa-let-7a--NF2    0.469280186    0.000231065
> > hsa-let-7a--SLC17A5    0.454597978    0.000381713
> > hsa-let-7a--THOC2    0.447714054    0.000479322
> > hsa-let-7a--SMG7    0.444972282    0.000524129
> >
>
> Suppose your data is "d": then try which(d[,3] < 0.05)
>
> > 2. I got the p.adjust.R from R source. In the method "BH", I am not clear
> > with the code:
> >           i <- lp:1L
>
>
> # Just the same as seq(lp, 1 , by = -1)
>
>
> >           o <- order(p, decreasing = TRUE)
> >           ro <- order(o)
> >           pmin(1, cummin( n / i * p[o] ))[ro]
>
> # pmin does parallel minimums, p[o] is the same as sort(p) and
> ordering by [ro] puts the outputted values in reverse order than the
> went in.
>
> As an exercise, I'd suggest you get the original paper, see how the
> calculation is done there, implement it in R as best you can, even if
> it seems loop-y, and refine it down to R Core's implementation. One of
> the best ways I know to learn to think vectorwise.
>
> Sorry I can't help more, but I don't know the method so I dont want to
> read too much into the code and say something that I havent thought
> through (Lord knows I do that enough on this list!!)
>
> Michael
>
> >
>
>
> > How to explain the first and the fourth row.
> > ====================p.adjust.R=======================================
> > p.adjust.methods <-
> >    c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr",
> "none")
> > p.adjust <- function(p, method = p.adjust.methods, n = length(p))
> > {
> >    ## Methods 'Hommel', 'BH', 'BY' and speed improvements contributed by
> >    ## Gordon Smyth <smyth@wehi.edu.au>.
> >    method <- match.arg(method)
> >    if(method == "fdr") method <- "BH" # back compatibility
> >    nm <- names(p)
> >    p <- as.numeric(p); names(p) <- nm
> >    p0 <- p
> >    if(all(nna <- !is.na(p))) nna <- TRUE
> >    p <- p[nna]
> >    lp <- length(p)
> >    stopifnot(n >= lp)
> >    if (n <= 1) return(p0)
> >    if (n == 2 && method == "hommel") method <- "hochberg"
> >    p0[nna] <-
> >  switch(method,
> >        bonferroni = pmin(1, n * p),
> >        holm = {
> >     i <- seq_len(lp)
> >     o <- order(p)
> >     ro <- order(o)
> >     pmin(1, cummax( (n - i + 1L) * p[o] ))[ro]
> >        },
> >        hommel = { ## needs n-1 >= 2 in for() below
> >     if(n > lp) p <- c(p, rep.int(1, n-lp))
> >     i <- seq_len(n)
> >     o <- order(p)
> >     p <- p[o]
> >     ro <- order(o)
> >     q <- pa <- rep.int( min(n*p/i), n)
> >     for (j in (n-1):2) {
> >         ij <- seq_len(n-j+1)
> >         i2 <- (n-j+2):n
> >         q1 <- min(j*p[i2]/(2:j))
> >         q[ij] <- pmin(j*p[ij], q1)
> >         q[i2] <- q[n-j+1]
> >         pa <- pmax(pa,q)
> >     }
> >     pmax(pa,p)[if(lp < n) ro[1:lp] else ro]
> >        },
> >        hochberg = {
> >     i <- lp:1L
> >     o <- order(p, decreasing = TRUE)
> >     ro <- order(o)
> >     pmin(1, cummin( (n - i + 1L) * p[o] ))[ro]
> >        },
> >        BH = {
> >     i <- lp:1L
> >     o <- order(p, decreasing = TRUE)
> >     ro <- order(o)
> >     pmin(1, cummin( n / i * p[o] ))[ro]
> >        },
> >        BY = {
> >     i <- lp:1L
> >     o <- order(p, decreasing = TRUE)
> >     ro <- order(o)
> >     q <- sum(1L/(1L:n))
> >     pmin(1, cummin(q * n / i * p[o]))[ro]
> >        },
> >        none = p)
> >    p0
> > }
> > ============================================================
> >
> >
> > I wrote a code to do my work in BH correction like the following:
> >
> > rm(list=ls())
> > a<-read.csv("test.txt",sep="\t",header=F,quote="")
> > b<-a[order(a[,3],decreasing=TRUE),]
> > c<-p.adjust(b[,3],method="BH")
> > b[,4]<-c
> > write.table(b,"zz.txt",sep="\t")
> >
> > Is that right? Thanks for all.
> >
> > Jiang
> >
> >        [[alternative HTML version deleted]]
> >
> > ______________________________________________
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html<http://www.r-project.org/posting-guide.html>
> > and provide commented, minimal, self-contained, reproducible code.
> >
>

	[[alternative HTML version deleted]]

