R/nmfconsensus.R
In syspremed/PanNETassigner: A cross-species analysis in pancreatic neuroendocrine tumors reveals molecular subtypes with distinctive clinical, metastatic, developmental, and metabolic characteristics.
nmfconsensus <-function(input.ds, k.init, k.final, num.clusterings, maxniter, error.function, rseed=123456789, stopconvergence = 40, stopfrequency = 10, non.interactive.run = F, doc.string = "",data=FALSE, ...) {
#
# GenePattern Methodology for:
#
# Metagenes and Molecular Pattern Discovery using Matrix Factorization
# Jean-Philippe Brunet, Pablo Tamayo, Todd. R. Golub, and Jill P. Mesirov
#
# Author: Pablo Tamayo (tamayo@genome.wi.mit.edu)
#
# Based on the original matlab version written by Jean-Philippe Brunet (brunet@broad.mit.edu) and
# with additional contributions from: Ted Liefeld (liefeld@broad.mit.edu)
# Date: November 27, 2003
#
# Last change March 3, 2005: modifications to make the output more readable.
#
# Execute from an R console window with this command:
# source("<this file>", echo = TRUE)
# E.g. someoutput <- mynmf2(input.ds="c:\\nmf\\all_aml.res",k.init=2,k.final=5,num.clusterings=20,maxniter=500, error.function="euclidean")
#
# For details on the method see:
#
# Proc. Natl. Acad. Sci. USA 2004 101: 4164-4169
# http://www.broad.mit.edu/cgi-bin/cancer/publications/pub_paper.cgi?mode=view&paper_id=89
#
# Input parameters
#
# input.ds
# input gene expression dataset in GCT or RES format
# k.init
# initial value of k
# k.final
# final value of k
# num.clusterings
# number of NMF clusterings to build consensus matrix
# maxniter
# maximum number of NMF iterations
# error.function
# NMF error function: "divergence" of "euclidean"
# rseed
# random number generator seed
# stopconv
# how many no change checks are needed to stop NMF iterations (convergence)
# stopfreq
# frequency (NMF iterations) of "no change" checks
# non.interactive.run
# flag controlling if the plots are produced interatively (Rgui and saved) or only saved in files
# doc.string
# prefix to be added to the output files
#
# Output files are (prefix with "doc.string")
#
# params.txt
# run parameters and time of execution
# membership.gct
# membership results for samples at all values of K
# cophenetic.txt
# cophenetic values for each K
# cophenetic.plot.jpeg
# plot of cophenetic for each value of K
# consensus.k#.gct (for each value of K)
# consensus matrix for k=#
# consensus.plot.k#.jpeg (for each value of K)
# plot of consensus matrix for k=#
# graphs.k#.jpeg (for each value of K)
# save input parameters
directory = ""
time.string <- as.character(as.POSIXlt(Sys.time(),"GMT"))
if (non.interactive.run == "F"){
filename <- paste(directory, doc.string, ".params.txt", sep="", collapse="")
write(paste("Run of NMF on ", time.string), file=filename)
write(paste("input.ds =", input.ds, sep=" "), file=filename, append=T)
write(paste("k.init = ", k.init, sep=" "), file=filename, append=T)
write(paste("k.final =", k.final, sep=" "), file=filename, append=T)
write(paste("num.clusterings =", num.clusterings, sep=" "), file=filename, append=T)
write(paste("maxniter =", maxniter, sep=" "), file=filename, append=T)
write(paste("error.function =", error.function, sep=" "), file=filename, append=T)
write(paste("rseed =", rseed, sep=" "), file=filename, append=T)
write(paste("directory =", directory, sep=" "), file=filename, append=T)
write(paste("stopconv =", stopconvergence, sep=" "), file=filename, append=T)
write(paste("stopfreq =", stopfrequency, sep=" "), file=filename, append=T)
write(paste("non.interctive.run =", non.interactive.run, sep=" "), file=filename, append=T)
write(paste("doc.string =", doc.string, sep=" "), file=filename, append=T)
}
######################################################################################
######### Include all functions required by nmfconcensus function here
######################################################################################
read.dataset <- function(file) {
result <- regexpr(paste(".txt","$",sep=""), tolower(file))
if(result[[1]] != -1)
return(read.gct(file))
result <- regexpr(paste(".res","$",sep=""), tolower(file))
if(result[[1]] != -1)
return(read.res(file))
stop("Input is not a res or gct file.")
}
matrix.abs.plot <- function(V, axes = F, log = F, norm = T, transpose = T, matrix.order = T, max.v = 1, min.v = 0, main = " ", sub = " ", xlab = " ", ylab = " ") {
rows <- length(V[,1])
cols <- length(V[1,])
if (log == T) {
V <- log(V)
}
B <- matrix(0, nrow=rows, ncol=cols)
for (i in 1:rows) {
for (j in 1:cols) {
if (matrix.order == T) {
k <- rows - i + 1
} else {
k <- i
}
if (norm == T) {
if ((max.v == 1) && (min.v == 0)) {
max.val <- max(V)
min.val <- min(V)
} else {
max.val = max.v
min.val = min.v
}
}
B[k, j] <- max.val - V[i, j] + min.val
}
}
if (transpose == T) {
B <- t(B)
}
if (norm == T) {
#removed gamma = 1.5
image(z = B, zlim = c(min.val, max.val), axes = axes, col = rainbow(100, s = 1.0, v = 0.75, start = 0.0, end = 0.75), main = main, sub = sub, xlab = xlab, ylab = ylab)
} else {
image(z = B, axes = axes, col = rainbow(100, s = 1, v = 0.6, start = 0.1, end = 0.9, gamma = 1), main = main, sub = sub, xlab = xlab, ylab = ylab)
}
return(list(B, max.val, min.val))
}
metagene.plot <- function(H, main = " ", sub = " ", xlab = "samples ", ylab = "amplitude") {
k <- length(H[,1])
S <- length(H[1,])
index <- 1:S
maxval <- max(H)
minval <- min(H)
plot(index, H[1,], xlim=c(1, S), ylim=c(minval, maxval), main = main, sub = sub, ylab = ylab, xlab = xlab, type="n")
for (i in 1:k) {
lines(index, H[i,], type="l", col = i, lwd=2)
}
}
NMF <- function(V, k, maxniter = 2000, seed = 123456, stopconv = 40, stopfreq = 10) {
N <- length(V[,1])
M <- length(V[1,])
set.seed(seed)
W <- matrix(runif(N*k), nrow = N, ncol = k) # Initialize W and H with random numbers
H <- matrix(runif(k*M), nrow = k, ncol = M)
VP <- matrix(nrow = N, ncol = M)
error.v <- vector(mode = "numeric", length = maxniter)
new.membership <- vector(mode = "numeric", length = M)
old.membership <- vector(mode = "numeric", length = M)
eps <- .Machine$double.eps
for (t in 1:maxniter) {
VP = W %*% H
H <- H * (crossprod(W, V)/crossprod(W, VP)) + eps
VP = W %*% H
H.t <- t(H)
W <- W * (V %*% H.t)/(VP %*% H.t) + eps
error.v[t] <- sqrt(sum((V - VP)^2))/(N * M)
if (t %% stopfreq == 0) {
for (j in 1:M) {
class <- order(H[,j], decreasing=T)
new.membership[j] <- class[1]
}
if (sum(new.membership == old.membership) == M) {
no.change.count <- no.change.count + 1
} else {
no.change.count <- 0
}
if (no.change.count == stopconv) break
old.membership <- new.membership
}
}
return(list(W = W, H = H, t = t, error.v = error.v))
}
NMF.div <- function(V, k, maxniter = 2000, seed = 123456, stopconv = 40, stopfreq = 10) {
N <- length(V[,1])
M <- length(V[1,])
set.seed(seed)
W <- matrix(runif(N*k), nrow = N, ncol = k) # Initialize W and H with random numbers
H <- matrix(runif(k*M), nrow = k, ncol = M)
VP <- matrix(nrow = N, ncol = M)
error.v <- vector(mode = "numeric", length = maxniter)
new.membership <- vector(mode = "numeric", length = M)
old.membership <- vector(mode = "numeric", length = M)
no.change.count <- 0
eps <- .Machine$double.eps
for (t in 1:maxniter) {
VP = W %*% H
W.t <- t(W)
H <- H * (W.t %*% (V/VP)) + eps
norm <- apply(W, MARGIN=2, FUN=sum)
for (i in 1:k) {
H[i,] <- H[i,]/norm[i]
}
VP = W %*% H
H.t <- t(H)
W <- W * ((V/VP) %*% H.t) + eps
norm <- apply(H, MARGIN=1, FUN=sum)
for (i in 1:k) {
W[,i] <- W[,i]/norm[i]
}
error.v[t] <- sum(V * log((V + eps)/(VP + eps)) - V + VP)/(M * N)
if (t %% stopfreq == 0) {
for (j in 1:M) {
class <- order(H[,j], decreasing=T)
new.membership[j] <- class[1]
}
if (sum(new.membership == old.membership) == M) {
no.change.count <- no.change.count + 1
} else {
no.change.count <- 0
}
if (no.change.count == stopconv) break
old.membership <- new.membership
}
}
return(list(W = W, H = H, t = t, error.v = error.v))
}
ConsPlot <- function(V, col.labels, col.names, main = " ", sub = " ", xlab=" ", ylab=" ") {
# Plots a heatmap plot of a consensus matrix
cols <- length(V[1,])
B <- matrix(0, nrow=cols, ncol=cols)
max.val <- max(V)
min.val <- min(V)
for (i in 1:cols) {
for (j in 1:cols) {
k <- cols - i + 1
B[k, j] <- max.val - V[i, j] + min.val
}
}
col.names2 <- rev(col.names)
col.labels2 <- rev(col.labels)
D <- matrix(0, nrow=(cols + 1), ncol=(cols + 1))
col.tag <- vector(length=cols, mode="numeric")
current.tag <- 0
col.tag[1] <- current.tag
for (i in 2:cols) {
if (col.labels[i] != col.labels[i - 1]) {
current.tag <- 1 - current.tag
}
col.tag[i] <- current.tag
}
col.tag2 <- rev(col.tag)
D[(cols + 1), 2:(cols + 1)] <- ifelse(col.tag %% 2 == 0, 1.02, 1.01)
D[1:cols, 1] <- ifelse(col.tag2 %% 2 == 0, 1.02, 1.01)
D[(cols + 1), 1] <- 1.03
D[1:cols, 2:(cols + 1)] <- B[1:cols, 1:cols]
#gamma = 1.5 removed
col.map <- c(rainbow(100, s = 1.0, v = 0.75, start = 0.0, end = 0.75), "#BBBBBB", "#333333", "#FFFFFF")
image(1:(cols + 1), 1:(cols + 1), t(D), col = col.map, axes=FALSE, main=main, sub=sub, xlab= xlab, ylab=ylab)
for (i in 1:cols) {
col.names[i] <- paste(" ", substr(col.names[i], 1, 12), sep="")
col.names2[i] <- paste(substr(col.names2[i], 1, 12), " ", sep="")
}
axis(2, at=1:cols, labels=col.names2, adj= 0.5, tick=FALSE, las = 1, cex.axis=0.50, font.axis=1, line=-1)
axis(2, at=1:cols, labels=col.labels2, adj= 0.5, tick=FALSE, las = 1, cex.axis=0.65, font.axis=1, line=-1)
axis(3, at=2:(cols + 1), labels=col.names, adj= 1, tick=FALSE, las = 3, cex.axis=0.50, font.axis=1, line=-1)
axis(3, at=2:(cols + 1), labels=as.character(col.labels), adj = 1, tick=FALSE, las = 1, cex.axis=0.65, font.axis=1, line=-1)
return()
}
read.res <- function(filename = "NULL") {
#
# Reads a gene expression dataset in RES format and converts it into an R data frame
#
header.cont <- readLines(filename, n = 1)
temp <- unlist(strsplit(header.cont, "\t"))
colst <- length(temp)
header.labels <- temp[seq(3, colst, 2)]
ds <- read.delim(filename, header=F, row.names = 2, sep="\t", quote="", skip=3, blank.lines.skip=T, comment.char="", as.is=T)
colst <- length(ds[1,])
cols <- (colst - 1)/2
rows <- length(ds[,1])
A <- matrix(nrow=rows - 1, ncol=cols)
A <- ds[1:rows, seq(2, colst, 2)]
table1 <- data.frame(A)
names(table1) <- header.labels
return(table1)
}
read.gct <- function(filename = "NULL") {
#
# Reads a gene expression dataset in GCT format and converts it into an R data frame
#
ds <- read.delim(filename, header=T, quote="", row.names=1, blank.lines.skip=T, comment.char="", as.is=T,na.strings="NA")
# ds <- read.delim(filename, header=T, quote="", skip=2, row.names=1, blank.lines.skip=T, comment.char="", as.is=T)
# ds <- ds[-1]
return(ds)
}
write.gct <- function (gct, filename)
{
f <- file(filename, "w")
gct_1<-cbind(rownames(gct),gct)
colnames(gct_1)[1]<-"Name"
write.table(gct_1,file = f, row.names=FALSE, sep="\t")
# cat("#1.2", "\n", file = f, append = TRUE, sep = "")
# cat(dim(gct)[1], "\t", dim(gct)[2], "\n", file = f, append = TRUE, sep = "")
# cat("Name", file = f, append = TRUE, sep = "")
# cat("Description", file = f, append = TRUE, sep = "")
# names <- names(gct)
# for (j in 1:length(names)) {
# cat("\t", names[j], file = f, append = TRUE, sep = "")
# }
# cat("\n", file = f, append = TRUE, sep = "")
# oldWarn <- options(warn = -1)
# m <- matrix(nrow = dim(gct)[1], ncol = dim(gct)[2] + 1)
# m[, 1] <- row.names(gct)
# m[, 2] <- row.names(gct)
# index <- 2
# for (i in 1:dim(gct)[2]) {
# m[, index] <- gct[, i]
# index <- index + 1
# }
# write.table(m, file = f, append = TRUE, quote = FALSE, sep = "\t", eol = "\n", col.names = FALSE, row.names = FALSE)
# close(f)
# options(warn = 0)
return(gct)
}
########################################
k.init<-as.integer(k.init)
k.final<-as.integer(k.final)
num.clusterings<-as.integer(num.clusterings)
n.iter<-as.integer(maxniter)
if (!is.na(rseed)){
seed <- as.integer(rseed)
}
stopfreq <- as.integer(stopfrequency)
stopconv <- as.integer(stopconvergence)
# library(mva)
# library(MASS)
# library(GenePattern)
D <- data
A <- data.matrix(D)
# Threshold negative values to small quantity
eps <- .Machine$double.eps
A[A < 0] <- eps
cols <- length(A[1,])
rows <- length(A[,1])
col.names <- colnames(D) # edited to colnames instead of names(D)
num.k <- k.final - k.init + 1
rho <- vector(mode = "numeric", length = num.k)
k.vector <- vector(mode = "numeric", length = num.k)
k.index <- 1
connect.matrix.ordered <- array(0, c(num.k, cols, cols))
for (k in k.init:k.final) {
if (non.interactive.run == F) {
if (.Platform$OS.type == "windows") {
filename <- paste(directory, doc.string, "graphs.k", k, sep="", collapse="")
windows(width = 9, height = 11)
} else if (.Platform$OS.type == "unix") {
filename <- paste(directory, doc.string, "graphs.k", k, ".pdf", sep="", collapse="")
pdf(file=filename, width = 9, height = 11)
}
} else {
if (.Platform$OS.type == "unix") {
filename <- paste(directory, doc.string, "graphs.k", k, ".pdf", sep="", collapse="")
pdf(file=filename, width = 9, height = 11)
} else if (.Platform$OS.type == "windows") {
filename <- paste(directory, doc.string, "graphs.k", k, ".pdf", sep="", collapse="")
pdf(file=filename, width = 9, height = 11)
}
}
nf <- layout(matrix(c(1,2,3,4,5,6,7,8), 4, 2, byrow=T), c(1, 1, 1, 1), c(1, 1), TRUE)
assign <- matrix(0, nrow = num.clusterings, ncol = cols)
for (i in 1:num.clusterings) {
if (non.interactive.run == "F"){
print(paste("Computing clustering number=", i, " for k=", k, sep="", collapse=" "))
}
if (error.function == "divergence"){
NMF.out <- NMF.div(V = A, k = k, maxniter = n.iter, seed = seed + i, stopconv = stopconv, stopfreq = stopfreq)
} else if (error.function == "euclidean"){
NMF.out <- NMF(V = A, k = k, maxniter = n.iter, seed = seed + i, stopconv = stopconv, stopfreq = stopfreq)
} else {
stop(paste("Un-supported error function=", error.function, sep=""))
}
if (non.interactive.run == "F"){
print(paste(NMF.out$t, " NMF iterations performed", sep="", collapse=" "))
}
for (j in 1:cols) { # Find membership
class <- order(NMF.out$H[,j], decreasing=T)
assign[i, j] <- class[1]
}
if (i == 3) { # '2''20' example for first clustering iteration
H.saved <- NMF.out$H
sub.string <- paste(doc.string, " k=", k, sep="")
plot(1:NMF.out$t, NMF.out$error.v[1:NMF.out$t], pch = 20, cex = 1.5, col = 1, xlab="time", ylab="NMF error", sub=sub.string, main=paste("Example of NMF convergence plot k=", k, sep=""))
if (rows < 1000) {
W <- NMF.out$W
} else {
W <- NMF.out$W[sample(x = 1:rows, size = 1000),]
}
sub.string <- paste(doc.string, " k=", k, sep="")
matrix.abs.plot(W, sub = sub.string, log = F, main = "Example W matrix (orig. ordering)", ylab = "genes", xlab ="metasamples")
matrix.abs.plot(H.saved, sub = sub.string, log = F, main = "Example H matrix (orig. ordering)", ylab = "metagenes", xlab ="samples")
metagene.plot(H = H.saved, main = "Example metagenes (orig. ordering)", sub = sub.string, xlab = "samples", ylab = "metagenes")
}
# Anguraj added to find the metagenes
filename <- paste(directory, doc.string, "metagenes.k.",k,".",i, ".txt", sep="", collapse="")
write.table(NMF.out$W, filename,sep="\t")
filename1 <- paste(directory, doc.string, "metasamples.k.",k,".",i, ".txt", sep="", collapse="")
write.table(NMF.out$H, filename1,sep="\t")
rm(NMF.out)
} ## end for (i in 1:num.clusterings)
# compute consensus matrix
connect.matrix <- matrix(0, nrow = cols, ncol = cols)
for (i in 1:num.clusterings) {
for (j in 1:cols) {
for (p in 1:cols) {
if (j != p) {
if (assign[i, j] == assign[i, p]) {
connect.matrix[j, p] <- connect.matrix[j, p] + 1
}
} else {
connect.matrix[j, p] <- connect.matrix[j, p] + 1
}
}
}
}
connect.matrix <- connect.matrix / num.clusterings
dist.matrix <- 1 - connect.matrix
dist.matrix <- as.dist(dist.matrix)
HC <- hclust(dist.matrix, method="average")
dist.coph <- cophenetic(HC)
k.vector[k.index] <- k
rho[k.index] <- cor(dist.matrix, dist.coph)
rho[k.index] <- signif(rho[k.index], digits = 4)
# connect.matrix.ordered <- matrix(0, nrow=cols, ncol = cols)
for (i in 1:cols) {
for (j in 1:cols) {
connect.matrix.ordered[k.index, i, j] <- connect.matrix[HC$order[i], HC$order[j]]
}
}
# compute consensus clustering membership
membership <- cutree(HC, k = k)
max.k <- max(membership)
items.names.ordered <- col.names[HC$order]
membership.ordered <- membership[HC$order]
results <- data.frame(cbind(membership.ordered, items.names.ordered))
if (k > k.init){
all.membership <- cbind(all.membership, membership);
} else {
all.membership <- cbind(membership);
}
sub.string <- paste(doc.string, " k=", k, sep="")
matrix.abs.plot(connect.matrix.ordered[k.index,,], sub=sub.string, log = F, main = "Ordered Consensus Matrix", ylab = "samples", xlab ="samples")
plot(HC, xlab="samples", cex = 0.75, labels = col.names, sub = sub.string, col = "blue", main = paste("Ordered Linkage Tree. Coph=", rho[k.index]))
resultsGct <- data.frame(membership.ordered)
row.names(resultsGct) <- items.names.ordered
filename <- paste(directory, doc.string, "consensus.k.",k, ".gct", sep="", collapse="")
write.gct(resultsGct, filename)
H.sorted <- H.saved[,HC$order]
sub.string <- paste(doc.string, " k=", k, sep="")
matrix.abs.plot(H.sorted, sub = sub.string, log = F, main = "Example H matrix (ordered)", ylab = "metagenes", xlab ="samples")
metagene.plot(H = H.sorted, sub = sub.string, main = "Example metagenes (ordered)", xlab = "samples", ylab = "metagenes")
if (non.interactive.run == F) {
if (.Platform$OS.type == "windows") {
savePlot(filename = filename, type ="jpeg", device = dev.cur())
} else if (.Platform$OS.type == "unix") {
dev.off()
}
} else {
dev.off()
}
if (non.interactive.run == F) {
if (.Platform$OS.type == "windows") {
filename <- paste(directory, doc.string, "consensus.plot.k", k, sep="", collapse="")
windows(width = 8.5, height = 11)
} else if (.Platform$OS.type == "unix") {
filename <- paste(directory, doc.string, "consensus.plot.k", k, ".pdf", sep="", collapse="")
pdf(file=filename, width = 8.5, height = 11)
}
} else {
if (.Platform$OS.type == "unix") {
filename <- paste(directory, doc.string, "consensus.plot.k", k, ".pdf", sep="", collapse="")
pdf(file=filename, width = 8.5, height = 11)
} else if (.Platform$OS.type == "windows") {
filename <- paste(directory, doc.string, "consensus.plot.k", k, ".pdf", sep="", collapse="")
pdf(file=filename, width = 8.5, height = 11)
}
}
nf <- layout(matrix(c(1), 1, 1, byrow=T), c(1, 1), c(1, 1), TRUE)
conlabel <- paste("Consensus k =", k, sep=" ", collapse="")
sub.string <- paste("Consensus matrix k=", k, "; dataset= ", input.ds, sep="")
ConsPlot(connect.matrix.ordered[k.index,,], col.labels = membership.ordered, col.names = items.names.ordered, main = " ", sub=sub.string, xlab=" ", ylab=" ")
if (non.interactive.run == F) {
if (.Platform$OS.type == "windows") {
savePlot(filename = filename, type ="jpeg", device = dev.cur())
} else if (.Platform$OS.type == "unix") {
dev.off()
}
} else {
dev.off()
}
k.index <- k.index + 1
} # end of loop over k
# Save consensus matrices in one file
if (non.interactive.run == F) {
if (.Platform$OS.type == "windows") {
filename <- paste(directory, doc.string, "consensus.all.k.plot", sep="")
windows(width = 8.5, height = 11)
} else if (.Platform$OS.type == "unix") {
filename <- paste(directory, doc.string, "consensus.all.k.plot.pdf", sep="")
pdf(file=filename, width = 8.5, height = 11)
}
} else {
if (.Platform$OS.type == "unix") {
filename <- paste(directory, doc.string, "consensus.all.k.plot.pdf", sep="")
pdf(file=filename, width = 8.5, height = 11)
} else if (.Platform$OS.type == "windows") {
filename <- paste(directory, doc.string, "consensus.all.k.plot.pdf", sep="")
pdf(file=filename, width = 8.5, height = 11)
}
}
nf <- layout(matrix(c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16), 4, 4, byrow=T), c(1, 1, 1, 1), c(1, 1, 1, 1), TRUE)
for (k in 1:num.k) {
matrix.abs.plot(connect.matrix.ordered[k,,], log = F, main = paste("k=", k.vector[k]),
sub = paste("Cophenetic coef.=", rho[k]), ylab = "samples", xlab ="samples")
}
y.range <- c(1 - 2*(1 - min(rho)), 1)
plot(k.vector, rho, main ="Cophenetic Coefficient", xlim=c(k.init, k.final), ylim=y.range,
xlab = "k", ylab="Cophenetic correlation", type = "n")
lines(k.vector, rho, type = "l", col = "black")
points(k.vector, rho, pch=22, type = "p", cex = 1.25, bg = "black", col = "black")
if (non.interactive.run == F) {
if (.Platform$OS.type == "windows") {
savePlot(filename = filename, type ="jpeg", device = dev.cur())
} else if (.Platform$OS.type == "unix") {
dev.off()
}
} else {
dev.off()
}
if (non.interactive.run == F) {
if (.Platform$OS.type == "windows") {
filename <- paste(directory, doc.string, "cophenetic.plot", sep="")
windows(width = 8.5, height = 11)
} else if (.Platform$OS.type == "unix") {
filename <- paste(directory, doc.string, "cophenetic.plot.pdf", sep="")
pdf(file=filename, width = 8.5, height = 11)
}
} else {
if (.Platform$OS.type == "unix") {
filename <- paste(directory, doc.string, "cophenetic.plot.pdf", sep="")
pdf(file=filename, width = 8.5, height = 11)
} else if (.Platform$OS.type == "windows") {
filename <- paste(directory, doc.string, "cophenetic.plot.pdf", sep="")
pdf(file=filename, width = 8.5, height = 11)
}
}
# Write the membership matrix
resultsmembership <- data.frame(all.membership)
row.names(resultsmembership) <- col.names
filename <- paste(directory, doc.string, "membership", ".gct", sep="", collapse="")
write.gct(resultsmembership , filename)
y.range <- c(1 - 2*(1 - min(rho)), 1)
plot(k.vector, rho, main ="Cophenetic Coefficient", xlim=c(k.init, k.final), ylim=y.range, xlab = "k", ylab="Cophenetic correlation", type = "n")
lines(k.vector, rho, type = "l", col = "black")
points(k.vector, rho, pch=22, type = "p", cex = 1.25, bg = "black", col = "black")
if (non.interactive.run == F) {
if (.Platform$OS.type == "windows") {
savePlot(filename = filename, type ="jpeg", device = dev.cur())
} else if (.Platform$OS.type == "unix") {
dev.off()
}
} else {
dev.off()
}
xx <- cbind(k.vector, rho)
write(xx, file= paste(directory, doc.string, ".", "cophenetic.txt", sep=""))
}
syspremed/PanNETassigner documentation built on May 31, 2019, 12:48 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
nmfconsensus <-function(input.ds, k.init, k.final, num.clusterings, maxniter, error.function, rseed=123456789, stopconvergence = 40, stopfrequency = 10, non.interactive.run = F, doc.string = "",data=FALSE, ...) {
#
# GenePattern Methodology for:
#
# Metagenes and Molecular Pattern Discovery using Matrix Factorization
# Jean-Philippe Brunet, Pablo Tamayo, Todd. R. Golub, and Jill P. Mesirov
#
# Author: Pablo Tamayo (tamayo@genome.wi.mit.edu)
#
# Based on the original matlab version written by Jean-Philippe Brunet (brunet@broad.mit.edu) and
# with additional contributions from: Ted Liefeld (liefeld@broad.mit.edu)
# Date: November 27, 2003
#
# Last change March 3, 2005: modifications to make the output more readable.
#
# Execute from an R console window with this command:
# source("<this file>", echo = TRUE)
# E.g. someoutput <- mynmf2(input.ds="c:\\nmf\\all_aml.res",k.init=2,k.final=5,num.clusterings=20,maxniter=500, error.function="euclidean")
#
# For details on the method see:
#
# Proc. Natl. Acad. Sci. USA 2004 101: 4164-4169
# http://www.broad.mit.edu/cgi-bin/cancer/publications/pub_paper.cgi?mode=view&paper_id=89
#
# Input parameters
#
# input.ds
# input gene expression dataset in GCT or RES format
# k.init
# initial value of k
# k.final
# final value of k
# num.clusterings
# number of NMF clusterings to build consensus matrix
# maxniter
# maximum number of NMF iterations
# error.function
# NMF error function: "divergence" of "euclidean"
# rseed
# random number generator seed
# stopconv
# how many no change checks are needed to stop NMF iterations (convergence)
# stopfreq
# frequency (NMF iterations) of "no change" checks
# non.interactive.run
# flag controlling if the plots are produced interatively (Rgui and saved) or only saved in files
# doc.string
# prefix to be added to the output files
#
# Output files are (prefix with "doc.string")
#
# params.txt
# run parameters and time of execution
# membership.gct
# membership results for samples at all values of K
# cophenetic.txt
# cophenetic values for each K
# cophenetic.plot.jpeg
# plot of cophenetic for each value of K
# consensus.k#.gct (for each value of K)
# consensus matrix for k=#
# consensus.plot.k#.jpeg (for each value of K)
# plot of consensus matrix for k=#
# graphs.k#.jpeg (for each value of K)
# save input parameters
directory = ""
time.string <- as.character(as.POSIXlt(Sys.time(),"GMT"))
if (non.interactive.run == "F"){
filename <- paste(directory, doc.string, ".params.txt", sep="", collapse="")
write(paste("Run of NMF on ", time.string), file=filename)
write(paste("input.ds =", input.ds, sep=" "), file=filename, append=T)
write(paste("k.init = ", k.init, sep=" "), file=filename, append=T)
write(paste("k.final =", k.final, sep=" "), file=filename, append=T)
write(paste("num.clusterings =", num.clusterings, sep=" "), file=filename, append=T)
write(paste("maxniter =", maxniter, sep=" "), file=filename, append=T)
write(paste("error.function =", error.function, sep=" "), file=filename, append=T)
write(paste("rseed =", rseed, sep=" "), file=filename, append=T)
write(paste("directory =", directory, sep=" "), file=filename, append=T)
write(paste("stopconv =", stopconvergence, sep=" "), file=filename, append=T)
write(paste("stopfreq =", stopfrequency, sep=" "), file=filename, append=T)
write(paste("non.interctive.run =", non.interactive.run, sep=" "), file=filename, append=T)
write(paste("doc.string =", doc.string, sep=" "), file=filename, append=T)
}
######################################################################################
######### Include all functions required by nmfconcensus function here
######################################################################################
read.dataset <- function(file) {
result <- regexpr(paste(".txt","$",sep=""), tolower(file))
if(result[[1]] != -1)
return(read.gct(file))
result <- regexpr(paste(".res","$",sep=""), tolower(file))
if(result[[1]] != -1)
return(read.res(file))
stop("Input is not a res or gct file.")
}
matrix.abs.plot <- function(V, axes = F, log = F, norm = T, transpose = T, matrix.order = T, max.v = 1, min.v = 0, main = " ", sub = " ", xlab = " ", ylab = " ") {
rows <- length(V[,1])
cols <- length(V[1,])
if (log == T) {
V <- log(V)
}
B <- matrix(0, nrow=rows, ncol=cols)
for (i in 1:rows) {
for (j in 1:cols) {
if (matrix.order == T) {
k <- rows - i + 1
} else {
k <- i
}
if (norm == T) {
if ((max.v == 1) && (min.v == 0)) {
max.val <- max(V)
min.val <- min(V)
} else {
max.val = max.v
min.val = min.v
}
}
B[k, j] <- max.val - V[i, j] + min.val
}
}
if (transpose == T) {
B <- t(B)
}
if (norm == T) {
#removed gamma = 1.5
image(z = B, zlim = c(min.val, max.val), axes = axes, col = rainbow(100, s = 1.0, v = 0.75, start = 0.0, end = 0.75), main = main, sub = sub, xlab = xlab, ylab = ylab)
} else {
image(z = B, axes = axes, col = rainbow(100, s = 1, v = 0.6, start = 0.1, end = 0.9, gamma = 1), main = main, sub = sub, xlab = xlab, ylab = ylab)
}
return(list(B, max.val, min.val))
}
metagene.plot <- function(H, main = " ", sub = " ", xlab = "samples ", ylab = "amplitude") {
k <- length(H[,1])
S <- length(H[1,])
index <- 1:S
maxval <- max(H)
minval <- min(H)
plot(index, H[1,], xlim=c(1, S), ylim=c(minval, maxval), main = main, sub = sub, ylab = ylab, xlab = xlab, type="n")
for (i in 1:k) {
lines(index, H[i,], type="l", col = i, lwd=2)
}
}
NMF <- function(V, k, maxniter = 2000, seed = 123456, stopconv = 40, stopfreq = 10) {
N <- length(V[,1])
M <- length(V[1,])
set.seed(seed)
W <- matrix(runif(N*k), nrow = N, ncol = k) # Initialize W and H with random numbers
H <- matrix(runif(k*M), nrow = k, ncol = M)
VP <- matrix(nrow = N, ncol = M)
error.v <- vector(mode = "numeric", length = maxniter)
new.membership <- vector(mode = "numeric", length = M)
old.membership <- vector(mode = "numeric", length = M)
eps <- .Machine$double.eps
for (t in 1:maxniter) {
VP = W %*% H
H <- H * (crossprod(W, V)/crossprod(W, VP)) + eps
VP = W %*% H
H.t <- t(H)
W <- W * (V %*% H.t)/(VP %*% H.t) + eps
error.v[t] <- sqrt(sum((V - VP)^2))/(N * M)
if (t %% stopfreq == 0) {
for (j in 1:M) {
class <- order(H[,j], decreasing=T)
new.membership[j] <- class[1]
}
if (sum(new.membership == old.membership) == M) {
no.change.count <- no.change.count + 1
} else {
no.change.count <- 0
}
if (no.change.count == stopconv) break
old.membership <- new.membership
}
}
return(list(W = W, H = H, t = t, error.v = error.v))
}
NMF.div <- function(V, k, maxniter = 2000, seed = 123456, stopconv = 40, stopfreq = 10) {
N <- length(V[,1])
M <- length(V[1,])
set.seed(seed)
W <- matrix(runif(N*k), nrow = N, ncol = k) # Initialize W and H with random numbers
H <- matrix(runif(k*M), nrow = k, ncol = M)
VP <- matrix(nrow = N, ncol = M)
error.v <- vector(mode = "numeric", length = maxniter)
new.membership <- vector(mode = "numeric", length = M)
old.membership <- vector(mode = "numeric", length = M)
no.change.count <- 0
eps <- .Machine$double.eps
for (t in 1:maxniter) {
VP = W %*% H
W.t <- t(W)
H <- H * (W.t %*% (V/VP)) + eps
norm <- apply(W, MARGIN=2, FUN=sum)
for (i in 1:k) {
H[i,] <- H[i,]/norm[i]
}
VP = W %*% H
H.t <- t(H)
W <- W * ((V/VP) %*% H.t) + eps
norm <- apply(H, MARGIN=1, FUN=sum)
for (i in 1:k) {
W[,i] <- W[,i]/norm[i]
}
error.v[t] <- sum(V * log((V + eps)/(VP + eps)) - V + VP)/(M * N)
if (t %% stopfreq == 0) {
for (j in 1:M) {
class <- order(H[,j], decreasing=T)
new.membership[j] <- class[1]
}
if (sum(new.membership == old.membership) == M) {
no.change.count <- no.change.count + 1
} else {
no.change.count <- 0
}
if (no.change.count == stopconv) break
old.membership <- new.membership
}
}
return(list(W = W, H = H, t = t, error.v = error.v))
}
ConsPlot <- function(V, col.labels, col.names, main = " ", sub = " ", xlab=" ", ylab=" ") {
# Plots a heatmap plot of a consensus matrix
cols <- length(V[1,])
B <- matrix(0, nrow=cols, ncol=cols)
max.val <- max(V)
min.val <- min(V)
for (i in 1:cols) {
for (j in 1:cols) {
k <- cols - i + 1
B[k, j] <- max.val - V[i, j] + min.val
}
}
col.names2 <- rev(col.names)
col.labels2 <- rev(col.labels)
D <- matrix(0, nrow=(cols + 1), ncol=(cols + 1))
col.tag <- vector(length=cols, mode="numeric")
current.tag <- 0
col.tag[1] <- current.tag
for (i in 2:cols) {
if (col.labels[i] != col.labels[i - 1]) {
current.tag <- 1 - current.tag
}
col.tag[i] <- current.tag
}
col.tag2 <- rev(col.tag)
D[(cols + 1), 2:(cols + 1)] <- ifelse(col.tag %% 2 == 0, 1.02, 1.01)
D[1:cols, 1] <- ifelse(col.tag2 %% 2 == 0, 1.02, 1.01)
D[(cols + 1), 1] <- 1.03
D[1:cols, 2:(cols + 1)] <- B[1:cols, 1:cols]
#gamma = 1.5 removed
col.map <- c(rainbow(100, s = 1.0, v = 0.75, start = 0.0, end = 0.75), "#BBBBBB", "#333333", "#FFFFFF")
image(1:(cols + 1), 1:(cols + 1), t(D), col = col.map, axes=FALSE, main=main, sub=sub, xlab= xlab, ylab=ylab)
for (i in 1:cols) {
col.names[i] <- paste(" ", substr(col.names[i], 1, 12), sep="")
col.names2[i] <- paste(substr(col.names2[i], 1, 12), " ", sep="")
}
axis(2, at=1:cols, labels=col.names2, adj= 0.5, tick=FALSE, las = 1, cex.axis=0.50, font.axis=1, line=-1)
axis(2, at=1:cols, labels=col.labels2, adj= 0.5, tick=FALSE, las = 1, cex.axis=0.65, font.axis=1, line=-1)
axis(3, at=2:(cols + 1), labels=col.names, adj= 1, tick=FALSE, las = 3, cex.axis=0.50, font.axis=1, line=-1)
axis(3, at=2:(cols + 1), labels=as.character(col.labels), adj = 1, tick=FALSE, las = 1, cex.axis=0.65, font.axis=1, line=-1)
return()
}
read.res <- function(filename = "NULL") {
#
# Reads a gene expression dataset in RES format and converts it into an R data frame
#
header.cont <- readLines(filename, n = 1)
temp <- unlist(strsplit(header.cont, "\t"))
colst <- length(temp)
header.labels <- temp[seq(3, colst, 2)]
ds <- read.delim(filename, header=F, row.names = 2, sep="\t", quote="", skip=3, blank.lines.skip=T, comment.char="", as.is=T)
colst <- length(ds[1,])
cols <- (colst - 1)/2
rows <- length(ds[,1])
A <- matrix(nrow=rows - 1, ncol=cols)
A <- ds[1:rows, seq(2, colst, 2)]
table1 <- data.frame(A)
names(table1) <- header.labels
return(table1)
}
read.gct <- function(filename = "NULL") {
#
# Reads a gene expression dataset in GCT format and converts it into an R data frame
#
ds <- read.delim(filename, header=T, quote="", row.names=1, blank.lines.skip=T, comment.char="", as.is=T,na.strings="NA")
# ds <- read.delim(filename, header=T, quote="", skip=2, row.names=1, blank.lines.skip=T, comment.char="", as.is=T)
# ds <- ds[-1]
return(ds)
}
write.gct <- function (gct, filename)
{
f <- file(filename, "w")
gct_1<-cbind(rownames(gct),gct)
colnames(gct_1)[1]<-"Name"
write.table(gct_1,file = f, row.names=FALSE, sep="\t")
# cat("#1.2", "\n", file = f, append = TRUE, sep = "")
# cat(dim(gct)[1], "\t", dim(gct)[2], "\n", file = f, append = TRUE, sep = "")
# cat("Name", file = f, append = TRUE, sep = "")
# cat("Description", file = f, append = TRUE, sep = "")
# names <- names(gct)
# for (j in 1:length(names)) {
# cat("\t", names[j], file = f, append = TRUE, sep = "")
# }
# cat("\n", file = f, append = TRUE, sep = "")
# oldWarn <- options(warn = -1)
# m <- matrix(nrow = dim(gct)[1], ncol = dim(gct)[2] + 1)
# m[, 1] <- row.names(gct)
# m[, 2] <- row.names(gct)
# index <- 2
# for (i in 1:dim(gct)[2]) {
# m[, index] <- gct[, i]
# index <- index + 1
# }
# write.table(m, file = f, append = TRUE, quote = FALSE, sep = "\t", eol = "\n", col.names = FALSE, row.names = FALSE)
# close(f)
# options(warn = 0)
return(gct)
}
########################################
k.init<-as.integer(k.init)
k.final<-as.integer(k.final)
num.clusterings<-as.integer(num.clusterings)
n.iter<-as.integer(maxniter)
if (!is.na(rseed)){
seed <- as.integer(rseed)
}
stopfreq <- as.integer(stopfrequency)
stopconv <- as.integer(stopconvergence)
# library(mva)
# library(MASS)
# library(GenePattern)
D <- data
A <- data.matrix(D)
# Threshold negative values to small quantity
eps <- .Machine$double.eps
A[A < 0] <- eps
cols <- length(A[1,])
rows <- length(A[,1])
col.names <- colnames(D) # edited to colnames instead of names(D)
num.k <- k.final - k.init + 1
rho <- vector(mode = "numeric", length = num.k)
k.vector <- vector(mode = "numeric", length = num.k)
k.index <- 1
connect.matrix.ordered <- array(0, c(num.k, cols, cols))
for (k in k.init:k.final) {
if (non.interactive.run == F) {
if (.Platform$OS.type == "windows") {
filename <- paste(directory, doc.string, "graphs.k", k, sep="", collapse="")
windows(width = 9, height = 11)
} else if (.Platform$OS.type == "unix") {
filename <- paste(directory, doc.string, "graphs.k", k, ".pdf", sep="", collapse="")
pdf(file=filename, width = 9, height = 11)
}
} else {
if (.Platform$OS.type == "unix") {
filename <- paste(directory, doc.string, "graphs.k", k, ".pdf", sep="", collapse="")
pdf(file=filename, width = 9, height = 11)
} else if (.Platform$OS.type == "windows") {
filename <- paste(directory, doc.string, "graphs.k", k, ".pdf", sep="", collapse="")
pdf(file=filename, width = 9, height = 11)
}
}
nf <- layout(matrix(c(1,2,3,4,5,6,7,8), 4, 2, byrow=T), c(1, 1, 1, 1), c(1, 1), TRUE)
assign <- matrix(0, nrow = num.clusterings, ncol = cols)
for (i in 1:num.clusterings) {
if (non.interactive.run == "F"){
print(paste("Computing clustering number=", i, " for k=", k, sep="", collapse=" "))
}
if (error.function == "divergence"){
NMF.out <- NMF.div(V = A, k = k, maxniter = n.iter, seed = seed + i, stopconv = stopconv, stopfreq = stopfreq)
} else if (error.function == "euclidean"){
NMF.out <- NMF(V = A, k = k, maxniter = n.iter, seed = seed + i, stopconv = stopconv, stopfreq = stopfreq)
} else {
stop(paste("Un-supported error function=", error.function, sep=""))
}
if (non.interactive.run == "F"){
print(paste(NMF.out$t, " NMF iterations performed", sep="", collapse=" "))
}
for (j in 1:cols) { # Find membership
class <- order(NMF.out$H[,j], decreasing=T)
assign[i, j] <- class[1]
}
if (i == 3) { # '2''20' example for first clustering iteration
H.saved <- NMF.out$H
sub.string <- paste(doc.string, " k=", k, sep="")
plot(1:NMF.out$t, NMF.out$error.v[1:NMF.out$t], pch = 20, cex = 1.5, col = 1, xlab="time", ylab="NMF error", sub=sub.string, main=paste("Example of NMF convergence plot k=", k, sep=""))
if (rows < 1000) {
W <- NMF.out$W
} else {
W <- NMF.out$W[sample(x = 1:rows, size = 1000),]
}
sub.string <- paste(doc.string, " k=", k, sep="")
matrix.abs.plot(W, sub = sub.string, log = F, main = "Example W matrix (orig. ordering)", ylab = "genes", xlab ="metasamples")
matrix.abs.plot(H.saved, sub = sub.string, log = F, main = "Example H matrix (orig. ordering)", ylab = "metagenes", xlab ="samples")
metagene.plot(H = H.saved, main = "Example metagenes (orig. ordering)", sub = sub.string, xlab = "samples", ylab = "metagenes")
}
# Anguraj added to find the metagenes
filename <- paste(directory, doc.string, "metagenes.k.",k,".",i, ".txt", sep="", collapse="")
write.table(NMF.out$W, filename,sep="\t")
filename1 <- paste(directory, doc.string, "metasamples.k.",k,".",i, ".txt", sep="", collapse="")
write.table(NMF.out$H, filename1,sep="\t")
rm(NMF.out)
} ## end for (i in 1:num.clusterings)
# compute consensus matrix
connect.matrix <- matrix(0, nrow = cols, ncol = cols)
for (i in 1:num.clusterings) {
for (j in 1:cols) {
for (p in 1:cols) {
if (j != p) {
if (assign[i, j] == assign[i, p]) {
connect.matrix[j, p] <- connect.matrix[j, p] + 1
}
} else {
connect.matrix[j, p] <- connect.matrix[j, p] + 1
}
}
}
}
connect.matrix <- connect.matrix / num.clusterings
dist.matrix <- 1 - connect.matrix
dist.matrix <- as.dist(dist.matrix)
HC <- hclust(dist.matrix, method="average")
dist.coph <- cophenetic(HC)
k.vector[k.index] <- k
rho[k.index] <- cor(dist.matrix, dist.coph)
rho[k.index] <- signif(rho[k.index], digits = 4)
# connect.matrix.ordered <- matrix(0, nrow=cols, ncol = cols)
for (i in 1:cols) {
for (j in 1:cols) {
connect.matrix.ordered[k.index, i, j] <- connect.matrix[HC$order[i], HC$order[j]]
}
}
# compute consensus clustering membership
membership <- cutree(HC, k = k)
max.k <- max(membership)
items.names.ordered <- col.names[HC$order]
membership.ordered <- membership[HC$order]
results <- data.frame(cbind(membership.ordered, items.names.ordered))
if (k > k.init){
all.membership <- cbind(all.membership, membership);
} else {
all.membership <- cbind(membership);
}
sub.string <- paste(doc.string, " k=", k, sep="")
matrix.abs.plot(connect.matrix.ordered[k.index,,], sub=sub.string, log = F, main = "Ordered Consensus Matrix", ylab = "samples", xlab ="samples")
plot(HC, xlab="samples", cex = 0.75, labels = col.names, sub = sub.string, col = "blue", main = paste("Ordered Linkage Tree. Coph=", rho[k.index]))
resultsGct <- data.frame(membership.ordered)
row.names(resultsGct) <- items.names.ordered
filename <- paste(directory, doc.string, "consensus.k.",k, ".gct", sep="", collapse="")
write.gct(resultsGct, filename)
H.sorted <- H.saved[,HC$order]
sub.string <- paste(doc.string, " k=", k, sep="")
matrix.abs.plot(H.sorted, sub = sub.string, log = F, main = "Example H matrix (ordered)", ylab = "metagenes", xlab ="samples")
metagene.plot(H = H.sorted, sub = sub.string, main = "Example metagenes (ordered)", xlab = "samples", ylab = "metagenes")
if (non.interactive.run == F) {
if (.Platform$OS.type == "windows") {
savePlot(filename = filename, type ="jpeg", device = dev.cur())
} else if (.Platform$OS.type == "unix") {
dev.off()
}
} else {
dev.off()
}
if (non.interactive.run == F) {
if (.Platform$OS.type == "windows") {
filename <- paste(directory, doc.string, "consensus.plot.k", k, sep="", collapse="")
windows(width = 8.5, height = 11)
} else if (.Platform$OS.type == "unix") {
filename <- paste(directory, doc.string, "consensus.plot.k", k, ".pdf", sep="", collapse="")
pdf(file=filename, width = 8.5, height = 11)
}
} else {
if (.Platform$OS.type == "unix") {
filename <- paste(directory, doc.string, "consensus.plot.k", k, ".pdf", sep="", collapse="")
pdf(file=filename, width = 8.5, height = 11)
} else if (.Platform$OS.type == "windows") {
filename <- paste(directory, doc.string, "consensus.plot.k", k, ".pdf", sep="", collapse="")
pdf(file=filename, width = 8.5, height = 11)
}
}
nf <- layout(matrix(c(1), 1, 1, byrow=T), c(1, 1), c(1, 1), TRUE)
conlabel <- paste("Consensus k =", k, sep=" ", collapse="")
sub.string <- paste("Consensus matrix k=", k, "; dataset= ", input.ds, sep="")
ConsPlot(connect.matrix.ordered[k.index,,], col.labels = membership.ordered, col.names = items.names.ordered, main = " ", sub=sub.string, xlab=" ", ylab=" ")
if (non.interactive.run == F) {
if (.Platform$OS.type == "windows") {
savePlot(filename = filename, type ="jpeg", device = dev.cur())
} else if (.Platform$OS.type == "unix") {
dev.off()
}
} else {
dev.off()
}
k.index <- k.index + 1
} # end of loop over k
# Save consensus matrices in one file
if (non.interactive.run == F) {
if (.Platform$OS.type == "windows") {
filename <- paste(directory, doc.string, "consensus.all.k.plot", sep="")
windows(width = 8.5, height = 11)
} else if (.Platform$OS.type == "unix") {
filename <- paste(directory, doc.string, "consensus.all.k.plot.pdf", sep="")
pdf(file=filename, width = 8.5, height = 11)
}
} else {
if (.Platform$OS.type == "unix") {
filename <- paste(directory, doc.string, "consensus.all.k.plot.pdf", sep="")
pdf(file=filename, width = 8.5, height = 11)
} else if (.Platform$OS.type == "windows") {
filename <- paste(directory, doc.string, "consensus.all.k.plot.pdf", sep="")
pdf(file=filename, width = 8.5, height = 11)
}
}
nf <- layout(matrix(c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16), 4, 4, byrow=T), c(1, 1, 1, 1), c(1, 1, 1, 1), TRUE)
for (k in 1:num.k) {
matrix.abs.plot(connect.matrix.ordered[k,,], log = F, main = paste("k=", k.vector[k]),
sub = paste("Cophenetic coef.=", rho[k]), ylab = "samples", xlab ="samples")
}
y.range <- c(1 - 2*(1 - min(rho)), 1)
plot(k.vector, rho, main ="Cophenetic Coefficient", xlim=c(k.init, k.final), ylim=y.range,
xlab = "k", ylab="Cophenetic correlation", type = "n")
lines(k.vector, rho, type = "l", col = "black")
points(k.vector, rho, pch=22, type = "p", cex = 1.25, bg = "black", col = "black")
if (non.interactive.run == F) {
if (.Platform$OS.type == "windows") {
savePlot(filename = filename, type ="jpeg", device = dev.cur())
} else if (.Platform$OS.type == "unix") {
dev.off()
}
} else {
dev.off()
}
if (non.interactive.run == F) {
if (.Platform$OS.type == "windows") {
filename <- paste(directory, doc.string, "cophenetic.plot", sep="")
windows(width = 8.5, height = 11)
} else if (.Platform$OS.type == "unix") {
filename <- paste(directory, doc.string, "cophenetic.plot.pdf", sep="")
pdf(file=filename, width = 8.5, height = 11)
}
} else {
if (.Platform$OS.type == "unix") {
filename <- paste(directory, doc.string, "cophenetic.plot.pdf", sep="")
pdf(file=filename, width = 8.5, height = 11)
} else if (.Platform$OS.type == "windows") {
filename <- paste(directory, doc.string, "cophenetic.plot.pdf", sep="")
pdf(file=filename, width = 8.5, height = 11)
}
}
# Write the membership matrix
resultsmembership <- data.frame(all.membership)
row.names(resultsmembership) <- col.names
filename <- paste(directory, doc.string, "membership", ".gct", sep="", collapse="")
write.gct(resultsmembership , filename)
y.range <- c(1 - 2*(1 - min(rho)), 1)
plot(k.vector, rho, main ="Cophenetic Coefficient", xlim=c(k.init, k.final), ylim=y.range, xlab = "k", ylab="Cophenetic correlation", type = "n")
lines(k.vector, rho, type = "l", col = "black")
points(k.vector, rho, pch=22, type = "p", cex = 1.25, bg = "black", col = "black")
if (non.interactive.run == F) {
if (.Platform$OS.type == "windows") {
savePlot(filename = filename, type ="jpeg", device = dev.cur())
} else if (.Platform$OS.type == "unix") {
dev.off()
}
} else {
dev.off()
}
xx <- cbind(k.vector, rho)
write(xx, file= paste(directory, doc.string, ".", "cophenetic.txt", sep=""))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.