r-rsc-statistics-data-mining/kmeans-clustering.R at master · sing-group/r-rsc-statistics-data-mining · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
imagesDirectory <- "images/clustering/kmeans/"

dir.create(imagesDirectory,
           recursive = TRUE,
           showWarnings = FALSE)

source("load-cancer.R")

#	---------------------------------------------------------------------------
#
#	1 K-means clustering
#
#	---------------------------------------------------------------------------

# Since k-means starts its centroids randomnly, we set a random seed in order
# to reproduce the same results

set.seed(2019)

#	---------------------------------------------------------------------------
#
#	1.1 K-means using the replicates matrix (binnedPeaksMatrix)
#
#	---------------------------------------------------------------------------

clustering.kmeans <- kmeans(binnedPeaksMatrix, centers = 3)
clustering.kmeans$cluster

#	---------------------------------------------------------------------------
#
#	1.2 K-means using the samples matrix (consensusBinnedPeaksMatrix)
#
#	---------------------------------------------------------------------------

clustering.kmeans.consensus <-
  kmeans(consensusBinnedPeaksMatrix, centers = 3)
clustering.kmeans.consensus$cluster

#	---------------------------------------------------------------------------
#
#	1.3 K-means using the replicates matrix (binnedPeaksMatrix) and using presence/absence
#
#	---------------------------------------------------------------------------

presenceBinnedPeaksMatrix <- asPresenceMatrix(binnedPeaksMatrix)
clustering.kmeans.presence <- kmeans(binnedPeaksMatrix, centers = 3)
clustering.kmeans.presence$cluster

#	---------------------------------------------------------------------------
#
#	1.4 K-means using the samples matrix (consensusBinnedPeaksMatrix) and using presence/abscence
#
#	---------------------------------------------------------------------------

presenceConsensusBinnedPeaksMatrix <-
  asPresenceMatrix(consensusBinnedPeaksMatrix)
clustering.kmeans.presence.consensus <-
  kmeans(consensusBinnedPeaksMatrix, centers = 3)
clustering.kmeans.presence.consensus$cluster


#	---------------------------------------------------------------------------
#
#	1.5 Try different number of clusters and plot the total withinss (calculated
#       with sum(withinss) or using tot.withinss
#
#	---------------------------------------------------------------------------

clustering.kmeans.k_withinss <- vector()
for (c in 1:10) {
  set.seed(2019)
  clust <- kmeans(binnedPeaksMatrix, centers = c)
  clustering.kmeans.k_withinss[c] <- sum(clust$withinss)
}
png(
  paste0(imagesDirectory, "clustering-kmeans-k-withinss.png"),
  width = 640,
  height = 480
)
plot(
  clustering.kmeans.k_withinss,
  type = "b",
  xlab = "Number of clusters",
  ylab = "total within_SS",
  main = "K-means quality at different k (replicates matrix)"
)
dev.off()

clustering.kmeans.k_withinss <- vector()
for (c in 1:10) {
  set.seed(2019)
  clust <- kmeans(consensusBinnedPeaksMatrix, centers = c)
  clustering.kmeans.k_withinss[c] <- sum(clust$withinss)
}
png(
  paste0(imagesDirectory, "clustering-kmeans-samples-k-withinss.png"),
  width = 640,
  height = 480
)
plot(
  clustering.kmeans.k_withinss,
  type = "b",
  xlab = "Number of clusters",
  ylab = "total within_SS",
  main = "K-means quality at different k (samples matrix)"
)
dev.off()

#	---------------------------------------------------------------------------
#
#	1.6 Get a consensus clustering with ConsensusClusterPlus
#
#	---------------------------------------------------------------------------

library("ConsensusClusterPlus")

pdf(paste0(imagesDirectory, "consensus-cluster-plus.pdf"))
ccp <-
  ConsensusClusterPlus(
    t(binnedPeaksMatrix),
    maxK = 3,
    distance = "euclidean",
    clusterAlg = "km",
    reps = 50,
    seed = 2019
  )
dev.off()
# show clustering results for k=3
ccp[[3]][["consensusClass"]]