###
library(aricode)
library(SNFtool)
library(cluster)
library(iClusterPlus)
library(PINSPlus)
library(HCfused)
source("~/GitHub/HC-fused/NEMO.R")
source("~/GitHub/HC-fused/sim.R")
mode="HC"
this_method="ward.D"
n.iter <- 100
my_var <- 5
VVV <- c(0.1,0.5,1,5,10,20) #c(10,15,20,30,50)
OURS <- vector("list",length(VVV))
SNF <- vector("list",length(VVV))
CON <- vector("list",length(VVV))
BAY <- vector("list",length(VVV))
PINS <- vector("list",length(VVV))
X_NEMO <- vector("list",length(VVV))
OURS_sil <- vector("list",length(VVV))
SNF_sil <- vector("list",length(VVV))
CON_sil <- vector("list",length(VVV))
BAY_sil <- vector("list",length(VVV))
PINS_sil <- vector("list",length(VVV))
NEMO_sil <- vector("list",length(VVV))
truelab <- c(1,1,2,2,2,2,2,2,3,3)
count <- 0
for(my_var in VVV){
count <- count + 1
for(xx in 1:n.iter){
print(my_var)
print(xx)
#HC-fused
res <- sim2(FALSE, my_var, mode=mode)
mat1 <- res$mat1
mat2 <- res$mat2
#obj <- as.list(1:dim(mat2)[1])
#res <- HC_fused(obj, mat1, mat2, n.iter=10)
#m <- (res$NETWORK/max(res$NETWORK))
HCres <- HC_fused_subtyping(omics = list(mat1 , mat2),
max.k = 9, this_method = "ward.D", HC.iter = 50)
P <- HCres$P
hc_final <- hclust(as.dist(P), method=this_method)
sil <- calc.SIL(as.dist(P),9, method=this_method)
k <- as.numeric(names(which.max(sil)))
cl <- cutree(hc_final,k)
#print(cl)
OURS[[count]] <- c(OURS[[count]],clustComp(cl,truelab)$ARI) #rand.index(cl,truelab)
OURS_sil[[count]] <- c(OURS_sil[[count]],sil[2])
#SNF
K=9
alpha=0.5
TTT=20
C=3
res <- sim2(FALSE, my_var)
Data1 <- res$mat1
Data2 <- res$mat2
Data1 = standardNormalization(Data1)
Data2 = standardNormalization(Data2)
Dist1 = (dist2(as.matrix(Data1),as.matrix(Data1)))^(1/2)
Dist2 = (dist2(as.matrix(Data2),as.matrix(Data2)))^(1/2)
W1 = affinityMatrix(Dist1, K, alpha)
W2 = affinityMatrix(Dist2, K, alpha)
W = SNF(list(W1,W2), K, TTT)
sil <- calc.SIL(as.dist(1-W),9)
#k <- as.numeric(names(which.max(sil)))
estimationResult <- estimateNumberOfClustersGivenGraph(W, 2:9);
k <- estimationResult[[1]]
group = spectralClustering(W,k)
SNF[[count]] <- c(SNF[[count]],clustComp(group,truelab)$ARI)
SNF_sil[[count]] <- c(SNF_sil[[count]],sil[2])
#rand.index(group,truelab)
#concatenate
res <- sim2(FALSE, my_var)
Data1 <- res$mat1
Data2 <- res$mat2
Data <- cbind(Data1,Data2)
hc_con <- hclust(dist(Data))
sil <- calc.SIL(dist(Data),9)
k <- as.numeric(names(which.max(sil)))
cl_con <- cutree(hc_con,k)
#print(cl)
CON[[count]] <- c(CON[[count]],clustComp(cl_con,truelab)$ARI) #rand.index(cl,truelab)
CON_sil[[count]] <- c(CON_sil[[count]],sil[2])
#iClusterPlus
#BIC <- rep(NaN,2)
#for(kk in 1:2){
# fit.single=iClusterPlus(dt1=Data1,dt2=Data2,type=c("gaussian","gaussian"),
# K=kk,maxiter=10)
# BIC[kk] <- fit.single$BIC
#}
#fit.single=iClusterPlus(dt1=Data1,dt2=Data2,type=c("gaussian","gaussian"),
# K=which.min(BIC),maxiter=10)
#cl_bayes <- fit.single$clusters
#BAY[[count]] <- c(BAY[[count]],clustComp(cl_bayes,truelab)$ARI) #rand.index(cl,truelab)
#BAY_sil[[count]] <- c(BAY_sil[[count]],sil[2])
## PINSPLUS
result <- try(SubtypingOmicsData(dataList = list(Data1,Data2), iterMin=20,verbose=FALSE),silent = TRUE)
if(class(result) == "try-error"){
cl_pins <- rep(1,dim(Data1)[1])
}else{
#print(result)
cl_pins <- result$cluster2
}
PINS[[count]] <- c(PINS[[count]],clustComp(cl_pins,truelab)$ARI)
#NEMO
#source("NEMO.R")
omics_list = list(as.data.frame(t(Data1)),as.data.frame(t(Data2)))
cl_nemo = nemo.clustering(omics_list,num.neighbors=7)
nag = nemo.affinity.graph(omics_list, k = 7)
sil_nemo <- calc.SIL(as.dist(1-nag),9)
X_NEMO[[count]] <- c(X_NEMO[[count]],clustComp(cl_nemo,truelab)$ARI)
NEMO_sil[[count]] <- c(NEMO_sil[[count]],sil_nemo[2])
####
}
}
#par(mfrow=c(2,4))
#boxplot(OURS, ylim=c(0,1), ylab="ARI", names=as.character(VVV), xlab="sd", main="HC-fused",col="grey")
#boxplot(SNF, ylim=c(0,1), ylab="ARI", names=as.character(VVV), xlab="sd", main="SNF",col="grey")
#boxplot(CON, ylim=c(0,1), ylab="ARI", names=as.character(VVV), xlab="sd", main="Concatenate",col="grey")
#boxplot(BAY, ylim=c(0,1), ylab="ARI", names=as.character(VVV), xlab="sd", main="iClusterPlus",col="grey")
#boxplot(OURS_sil, ylim=c(0,1), ylab="SIL (k=3)", names=as.character(VVV), xlab="sd", main="HC-fused",col="grey")
#boxplot(SNF_sil, ylim=c(0,1), ylab="SIL (k=3)", names=as.character(VVV), xlab="sd", main="SNF",col="grey")
#boxplot(CON_sil, ylim=c(0,1), ylab="SIL (k=3)", names=as.character(VVV), xlab="sd", main="Concatenate",col="grey")
#boxplot(BAY_sil, ylim=c(0,1), ylab="SIL (k=3)", names=as.character(VVV), xlab="sd", main="iClusterPlus",col="grey")
par(mfrow=c(2,3))
plot(sapply(OURS,mean), ylim=c(0,1), ylab="ARI", xaxt="n", xlab="variance", main="HC-fused",col="black", type="b", pch=19)
points(sapply(OURS,sd), col="dark grey", type="b", pch=19,lty=2)
axis(1,1:length(VVV),as.character(VVV))
plot(sapply(SNF,mean), ylim=c(0,1), ylab="ARI", xaxt="n", xlab="variance", main="SNF",col="blue", type="b", pch=9)
points(sapply(SNF,sd), col="dark grey", type="b", pch=19,lty=2)
axis(1,1:length(VVV),as.character(VVV))
#points(sapply(BAY,mean), ylim=c(0,1), ylab="ARI", xaxt="n", xlab="variance", main="",col="green", type="b", pch=10)
#axis(1,1:length(VVV),as.character(VVV))
plot(sapply(CON,mean), ylim=c(0,1), ylab="ARI", xaxt="n", xlab="variance", main="HC-concatenate",col="orange", type="b", pch=11)
points(sapply(CON,sd), col="dark grey", type="b", pch=19,lty=2)
axis(1,1:length(VVV),as.character(VVV))
plot(sapply(PINS,mean), ylim=c(0,1), ylab="ARI", xaxt="n", xlab="variance", main="PINSPLUS",col="red", type="b", pch=12)
points(sapply(PINS,sd), col="dark grey", type="b", pch=19,lty=2)
axis(1,1:length(VVV),as.character(VVV))
plot(sapply(X_NEMO,mean), ylim=c(0,1), ylab="ARI", xaxt="n", xlab="variance", main="NEMO",col="green", type="b", pch=5)
points(sapply(X_NEMO,sd), col="dark grey", type="b", pch=19,lty=2)
axis(1,1:length(VVV),as.character(VVV))
#plot(sapply(SPEC,mean), ylim=c(0,1), ylab="ARI", xaxt="n", xlab="variance", main="SPECTRUM",col="black", type="b", pch=19)
#axis(1,1:length(VVV),as.character(VVV))
#legend
#legend("bottomleft", c("HC-fused", "SNF","HC-concatenate","PINSPlus"),
#col = c("black", "blue", "orange","red"), text.col = "black",
#lty = c(1, 1, 1), pch = c(19, 9, 11, 12),
#merge = TRUE, bg = FALSE,bty = "n")
plot(sapply(OURS_sil,mean), ylim=c(0,1), ylab="SIL (k=3)", xaxt="n", xlab="variance", main="",col="black", type="b", pch=19)
#axis(1,1:length(VVV),as.character(VVV))
points(sapply(SNF_sil,mean), ylim=c(0,1), ylab="SIL (k=3)", xaxt="n", xlab="variance", main="",col="blue", type="b", pch=9)
#axis(1,1:length(VVV),as.character(VVV))
#points(sapply(BAY_sil,mean), ylim=c(0,1), ylab="SIL (k=3)", xaxt="n", xlab="variance", main="",col="black", type="b", pch=10)
#axis(1,1:length(VVV),as.character(VVV))
points(sapply(CON_sil,mean), ylim=c(0,1), ylab="SIL (k=3)", xaxt="n", xlab="variance", main="",col="orange", type="b", pch=11)
#axis(1,1:length(VVV),as.character(VVV))
points(sapply(NEMO_sil,mean), ylim=c(0,1), ylab="SIL (k=3)", xaxt="n", xlab="variance", main="",col="green", type="b", pch=5)
axis(1,1:length(VVV),as.character(VVV))
#legend
legend("topright", c("HC-fused", "SNF","HC-concatenate","NEMO"),
col = c("black", "blue", "orange","green"), text.col = "black",
lty = c(1, 1, 1,1), pch = c(19, 9, 11,5),
merge = TRUE, bg = FALSE,bty = "n")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.