Some useful R codes written by myself:

Network#

diffusion kernel#

diffusion.kernel <- function(adj, beta = 1)
{
 lap <- adj - diag(rowSums(adj))
 eig <- eigen(lap)
 dif <- eig$vectors %*% diag(exp(beta * eig$values)) %*% t(eig$vectors)
 dif
}

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))) 
}