Matrix#

merge rows with same name#

merge.rows <- function(x, method=mean, na.rm=FALSE) {
	x <- as.matrix(x[!is.na(rownames(x)), ])
	r <- unique(rownames(x))
	m <- length(r)
	if (dim(x)[1] != m) {
		y <- matrix(nrow=m, ncol=dim(x)[2])
		rownames(y) <- r
		colnames(y) <- colnames(x)
		for (i in 1:m) {
			n <- sum(rownames(x) == r[i])
			if (n > 1) {
				s <- as.matrix(x[rownames(x) == r[i], ])
				y[i, ] <- apply(s, 2, method, na.rm)
			} else {
				y[i, ] <- x[r[i], ]
			}
		}
	} else {
		y <- x
	}
	y
}

Statistics#

robust fold change#

min.abs <- function(x) {
	x[which.min(abs(x))]
}

get.robust.fc <- function(case, control) {
	m <- dim(case)[2]
	n <- dim(control)[2]
	fc.raw <- matrix(nrow=dim(case)[1], ncol=m*n)
	for (i in 1:m) {
		for (j in 1:n) {
			fc.raw[, (i-1)*m+j] <- case[,i] - control[,j]
		}
	}
	fc <- apply(fc.raw, 1, min.abs)
}

fold change with confidence interval#

fc <- function(x, y, confidence.level=90, var.equal=F) {
	fc.interval(length(x), mean(x), var(x), length(y), mean(y), var(y), confidence.level, var.equal)
}

fc.interval <- function(x.n, x.mu, x.var, y.n, y.mu, y.var, confidence.level=90, var.equal=F) {
	mu <- x.mu - y.mu
	if (var.equal) {
		nu <- x.n + y.n - 2
		se <- sqrt(((x.n-1)*x.var + (y.n-1)*y.var) / nu) * sqrt(1/x.n + 1/y.n)
	} else {
		nu <- (x.var/x.n + y.var/y.n)^2 / (x.var^2/x.n^2/(x.n-1) + y.var^2/y.n^2/(y.n-1))
		se <- sqrt(x.var/x.n + y.var/y.n)
	}
	t <- -qt((1-confidence.level/100)/2, df=nu)
	if (mu >= 0) {
		mu.lower <- max(0, mu - t*se)
		mu.upper <- mu + t*se
	} else {
		mu.lower <- min(0, mu + t*se)
		mu.upper <- mu - t*se
	}
	data.frame(fc=mu, df=nu, lower=mu.lower, upper=mu.upper)
}

p-value#

p.value <- function(cdf, z, params=numeric(0), side=0) {
	n <- length(params)
	p <- switch(n+1,
			cdf(z),
			cdf(z, params),
			cdf(z, params[1], params[2]),
			cdf(z, params[1], params[2], params[3])
		)
	if (side < 0) p
	else if (side > 0) 1-p
	else if (p < 1/2) 2*p
	else 2*(1-p)
}

t-test#

one sample#

my.test1.t <- function(x.n, x.mu, x.var, mu, side=0) {
	nu <- x.n - 1
	t <- (x.mu - mu) / sqrt(x.var / x.n)
	p <- p.value(pt, t, params=nu, side=side)
	data.frame(mean=x.mu, df=nu, statistic=t, p.value=p)
}

two samples#

my.test2.t <- function(x.n, x.mu, x.var, y.n, y.mu, y.var, var.equal=F, side=0) {
	mu <- x.mu - y.mu
	if (var.equal) {
		nu <- x.n + y.n - 2
		s <- sqrt(((x.n-1)*x.var + (y.n-1)*y.var) / nu)
		t <- mu / (s * sqrt(1/x.n + 1/y.n))
	} else {
		nu <- (x.var/x.n + y.var/y.n)^2 / (x.var^2/x.n^2/(x.n-1) + y.var^2/y.n^2/(y.n-1))
		t <- mu / sqrt(x.var/x.n + y.var/y.n)
	}
	p <- p.value(pt, t, params=nu, side=side)
	data.frame(mean=mu, df=nu, statistic=t, p.value=p)
}

normal test#

one sample#

my.test1.norm <- function(x.n, x.mu, x.var, mu, side=0) {
	nu <- x.n - 1
	z <- (x.mu - mu) / sqrt(x.var / x.n)
	p <- p.value(pnorm, z, side=side)
	data.frame(mean=x.mu, df=nu, statistic=z, p.value=p)
}

two samples#

my.test2.norm <- function(x.n, x.mu, x.var, y.n, y.mu, y.var, side=0) {
	mu <- x.mu - y.mu
	nu <- x.n + y.n
	z <- mu / sqrt(x.var/x.n + y.var/y.n)
	p <- p.value(pnorm, z, side=side)
	data.frame(mean=mu, df=nu, statistic=z, p.value=p)
}

Graphics#

Color function to generate green-red heat maps#

updown.colors <- function(n = 50, down.col = 1/3, up.col=1, saturation = 1) { 
	c(hsv(down.col, saturation, seq(1,0,length=n)), hsv(up.col, saturation, seq(0,1,length=n))) 
}