Nothing
#get Jensen-Shannon divergence
JSDt_cluster <- function(phylo,mat,plot=F){
dist.JSD <- function(inMatrix, pseudocount=0.000001, ...) {
KLD <- function(x,y) sum(x*log(x/y))
JSD <- function(x,y) sqrt(0.5*KLD(x,(x+y)/2)+0.5*KLD(y,(x+y)/2))
matrixColSize <- length(colnames(inMatrix))
matrixRowSize <- length(rownames(inMatrix))
colnames <- colnames(inMatrix)
resultsMatrix <- matrix(0, matrixColSize, matrixColSize)
inMatrix = apply(inMatrix,1:2,
function(x) ifelse (x==0,pseudocount,x))
for(i in 1:matrixColSize) {
for(j in 1:matrixColSize) {
resultsMatrix[i,j]=JSD(as.vector(inMatrix[,i]),
as.vector(inMatrix[,j]))
}
}
colnames -> colnames(resultsMatrix) -> rownames(resultsMatrix)
as.dist(resultsMatrix)->resultsMatrix
attr(resultsMatrix, "method") <- "dist"
return(resultsMatrix)
}
#gaussian kernel convolution
dens <- function(x, bw = bw.nrd0, kernel = kernelG, n = 4096,
from = min(x) - 3*sd, to = max(x) + 3*sd, adjust = 1,
...) {
if(has.na <- any(is.na(x))) {
na.omit(x)->x
if(length(x) == 0)
stop('too infinite.')
}
kernelG<-function(x, mean=0, sd=1)
dnorm(x, mean = mean, sd = sd)
x<-log(x)
sd <- (if(is.numeric(bw)) bw[1] else bw(x)) * adjust
X <- seq(from, to, len = n)
M <- outer(X, x, kernel, sd = sd, ...)
structure(list(x = X, y = rowMeans(M), bw = sd,
call = match.call(), n = length(x),
data.name = deparse(substitute(x)),
has.na = has.na), class = "density")
}
#take square-root of Jensen-Shannon divergence
JSDist <- function(x,y) sqrt(dist.JSD(x,y))
x<-lapply(1:dim(mat)[2],function(l){
spectR_t(phylo,mat[,l])$eigenvalues
} )
d<-c()
Ds<-c()
for(i in 1:length(x)){
d[[i]]<-dens(x[[i]])$x
Ds<-as.data.frame(cbind(Ds,d[[i]]))
}
colnames(Ds)<-colnames(mat)
#compute divergence matrix
if(min(Ds)<0){JSD<-as.matrix(JSDist(Ds-min(Ds)+1e-9))}
else{JSD<-as.matrix(JSDist(Ds))}
#cluster on k-medoids
maxCluster=dim(mat)[2]-1
clustersMedoid <- pamk(JSD,krange=1:maxCluster)
#return list of divergence matrix table, clusters with silhouette widths
ret<-list(JSD,clustersMedoid[[1]]$clusinfo,clustersMedoid[[1]]$silinfo$widths[,c(1,3)])
if(plot==T){
#plot heatmap
heatmap(JSD,symm=T)
dev.new()
#plot hierarchical clustering with bootstrap support
clustersHierarchy <- pvclust(JSD,r=seq(0.5,1.5,0.2))
plot(clustersHierarchy,cex=0.3)
pvrect(clustersHierarchy,alpha=0.9)
return(ret)
}
else{return(ret)}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.