## ----message=FALSE------------------------------------------------------------
library(PINSPlus)
data(AML2004)
data <- as.matrix(AML2004$Gene)

## -----------------------------------------------------------------------------
system.time(result <- PerturbationClustering(data = data, verbose = FALSE))

## ----eval=FALSE---------------------------------------------------------------
# result <- PerturbationClustering(data = data)

## -----------------------------------------------------------------------------
result$k

## -----------------------------------------------------------------------------
result$cluster

## -----------------------------------------------------------------------------
condition <- seq(unique(AML2004$Group[, 2]))
names(condition) = unique(AML2004$Group[, 2])
plot(prcomp(AML2004$Gene)$x, col = result$cluster, 
     pch = condition[AML2004$Group[, 2]], main = "AML2004")
legend("bottomright", legend = paste("Cluster ", sort(unique(result$cluster)), sep = ""),
        fill = sort(unique(result$cluster)))
legend("bottomleft", legend = names(condition), pch = condition)

## ----eval=FALSE---------------------------------------------------------------
# result <- PerturbationClustering(data = data, kMax = 5,
#                                  clusteringMethod = "kmeans")

## ----eval=FALSE---------------------------------------------------------------
# result <- PerturbationClustering(data = data, kMax = 5,
#                                  clusteringMethod = "pam")

## ----eval=FALSE---------------------------------------------------------------
# result <- PerturbationClustering(data = data, kMax = 5,
#                                  clusteringMethod = "hclust")

## ----eval=FALSE---------------------------------------------------------------
# result <- PerturbationClustering(
#     data = data,
#     clusteringMethod = "kmeans",
#     clusteringOptions = list(nstart = 100, iter.max = 500),
#     verbose = FALSE
# )

## ----eval=FALSE---------------------------------------------------------------
# result <- PerturbationClustering(data = data,
#     clusteringFunction = function(data, k){
#     # this function must return a vector of cluster
#     kmeans(x = data, centers = k, nstart = k*10, iter.max = 2000)$cluster
# })

## ----eval=FALSE---------------------------------------------------------------
# result <- PerturbationClustering(data = data,
#                                  perturbMethod = "noise",
#                                  perturbOptions = list(noise = 1.23))

## ----eval=FALSE---------------------------------------------------------------
# result <- PerturbationClustering(data = data,
#                                  perturbMethod = "noise",
#                                  perturbOptions = list(noisePercent = 10))

## ----eval=FALSE---------------------------------------------------------------
# result <- PerturbationClustering(data = data,
#                                  perturbMethod = "subsampling",
#                                  perturbOptions = list(percent = 80))

## ----eval=FALSE---------------------------------------------------------------
# result <- PerturbationClustering(data = data, perturbFunction = function(data){
#     rowNum <- nrow(data)
#     colNum <- ncol(data)
#     epsilon <-
#         matrix(
#             data = rnorm(rowNum * colNum, mean = 0, sd = 1.23456),
#             nrow = rowNum, ncol = colNum
#         )
# 
#     list(
#         data = data + epsilon,
#         ConnectivityMatrixHandler = function(connectivityMatrix, iter, k) {
#             connectivityMatrix
#         }
#     )
# })

## ----eval=FALSE---------------------------------------------------------------
# sampleNum <- 50000 # Number of samples
# geneNum <- 5000 # Number of genes
# subtypeNum <- 3 # Number of subtypes
# 
# # Generate expression matrix
# exprs <- matrix(rnorm(sampleNum*geneNum, 0, 1), nrow = sampleNum, ncol = geneNum)
# rownames(exprs) <- paste0("S", 1:sampleNum) # Assign unique names for samples
# 
# # Generate subtypes
# group <- sort(rep(1:subtypeNum, sampleNum/subtypeNum + 1)[1:sampleNum])
# names(group) <- rownames(exprs)
# 
# # Make subtypes separate
# for (i in 1:subtypeNum) {
#   exprs[group == i, 1:100 + 100*(i-1)] <- exprs[group == i, 1:100 + 100*(i-1)] + 2
# }
# 
# # Plot the data
# library(irlba)
# exprs.pca <- irlba::prcomp_irlba(exprs, n = 2)$x
# plot(exprs.pca, main = "PCA")

## ----eval=FALSE---------------------------------------------------------------
# set.seed(1)
# t1 <- Sys.time()
# result <- PerturbationClustering(data = exprs.pca, ncore = 1)
# t2 <- Sys.time()

## ----eval=FALSE---------------------------------------------------------------
# t2-t1

## ----eval=FALSE---------------------------------------------------------------
# result$k

## ----eval=FALSE---------------------------------------------------------------
# subtype <- result$cluster

## ----eval=FALSE---------------------------------------------------------------
# if (!require("mclust")) install.packages("mclust")
# library(mclust)
# ari <- mclust::adjustedRandIndex(subtype, group)
# 

## ----eval=FALSE---------------------------------------------------------------
# 
# colors <- as.numeric(as.character(factor(subtype)))
# 
# plot(exprs.pca, col = colors, main = "Cluster assigments for simulation data")
# 
# legend("topright", legend = paste("ARI:", ari))
# 
# legend("bottomright", fill = unique(colors),
#     legend = paste("Group ",
#                    levels(factor(subtype)), ": ",
#                    table(subtype)[levels(factor(subtype))], sep = "" )
#     )

## ----eval=FALSE---------------------------------------------------------------
# # Load the kidney cancer carcinoma data
# data(KIRC)
# # SubtypingOmicsData`'s input data must be a list of
# # numeric matrices that have the same number of rows:
# dataList <- list (as.matrix(KIRC$GE), as.matrix(KIRC$ME), as.matrix(KIRC$MI))
# names(dataList) <- c("GE", "ME", "MI")
# # Run `SubtypingOmicsData`:
# result <- SubtypingOmicsData(dataList = dataList)

## ----eval=FALSE---------------------------------------------------------------
# result <- SubtypingOmicsData(
#     dataList = dataList,
#     clusteringMethod = "kmeans",
#     clusteringOptions = list(nstart = 50)
# )

## ----eval=FALSE---------------------------------------------------------------
# library(survival)
# cluster1=result$cluster1;cluster2=result$cluster2
# a <- intersect(unique(cluster2), unique(cluster1))
# names(a) <- intersect(unique(cluster2), unique(cluster1))
# a[setdiff(unique(cluster2), unique(cluster1))] <-
#     seq(setdiff(unique(cluster2), unique(cluster1))) + max(cluster1)
# colors <- a[levels(factor(cluster2))]
# coxFit <- coxph(
#      Surv(time = Survival, event = Death) ~ as.factor(cluster2),
#      data = KIRC$survival,
#      ties = "exact"
# )
# mfit <- survfit(Surv(Survival, Death == 1) ~ as.factor(cluster2), data = KIRC$survival)
# plot(
#      mfit, col = colors, main = "Survival curves for KIRC, level 2",
#      xlab = "Days", ylab = "Survival",lwd = 2
# )
# legend("bottomright",
#     legend = paste(
#         "Cox p-value:", round(summary(coxFit)$sctest[3], digits = 5), sep = ""
#     )
# )
# legend(
#     "bottomleft",
#     fill = colors,
#     legend = paste("Group ", levels(factor(cluster2)), ": ",
#         table(cluster2)[levels(factor(cluster2))], sep =""
#     )
# )