Cacher is the code snippet organizer for pro developers

#library(devtools)
#devtools::install_github("PMBio/peer")
library(peer)
library(argparser, quietly=TRUE)

WriteTable <- function(data, filename, index.name) {
  datafile <- file(filename, open = "wt")
  on.exit(close(datafile))
  header <- c(index.name, colnames(data))
  writeLines(paste0(header, collapse="\t"), con=datafile, sep="\n")
  write.table(data, datafile, sep="\t", col.names=FALSE, quote=FALSE)
}




argv = vector(mode='list')
argv$n=250
argv$alphaprior_a = 0.001
  argv$alphaprior_b=0.01
  argv$epsprior_a=0.1
  argv$epsprior_b = 10
  argv$max_iter=1000
  argv$output_dir
  argv$prefix="PEER_gene"

M=t(as.matrix(dge.cpm))
nrows = nrow(M)

# run PEER
n=250
model <- PEER()
invisible(PEER_setNk(model, argv$n))
invisible(PEER_setPhenoMean(model, M))
invisible(PEER_setPriorAlpha(model, argv$alphaprior_a, argv$alphaprior_b))
invisible(PEER_setPriorEps(model, argv$epsprior_a, argv$epsprior_b))
invisible(PEER_setNmax_iterations(model, argv$max_iter))
# if(!is.null(covs)) {
#   invisible(PEER_setCovariates(model, covs))
# }
time <- system.time(PEER_update(model))

X <- PEER_getX(model)  # samples x PEER factors
A <- PEER_getAlpha(model)  # PEER factors x 1
R <- t(PEER_getResiduals(model))  # genes x samples

# add relevant row/column names
c <- paste0("InferredCov",1:ncol(X))
rownames(X) <- rownames(M)
colnames(X) <- c
rownames(A) <- c
colnames(A) <- "Alpha"
A <- as.data.frame(A)
A$Relevance <- 1.0 / A$Alpha
rownames(R) <- colnames(M)
colnames(R) <- rownames(M)

# write results
cat("PEER: writing results ... ")
WriteTable(t(X), file.path(argv$output_dir, paste0(argv$prefix, ".PEER_covariates.txt")), "ID")  # format(X, digits=6)
WriteTable(A, file.path(argv$output_dir, paste0(argv$prefix, ".PEER_alpha.txt")), "ID")
WriteTable(R, file.path(argv$output_dir, paste0(argv$prefix, ".PEER_residuals.txt")), "ID")