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
2. I got the p.adjust.R from R source. In the method "BH", I am not clear
with the code:
i <- lp:1L
o <- order(p, decreasing = TRUE)
ro <- order(o)
pmin(1, cummin( n / i * p[o] ))[ro]
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 <sm...@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
and provide commented, minimal, self-contained, reproducible code.
