Switch to unified view

a b/R/DIscBIO-generic-PlotmclustMB.R
1
#' @title Plotting the Model-based clusters in PCA.
2
#' @description Plot the model-based clustering results
3
#' @param object \code{DISCBIO} class object.
4
#' @importFrom ggplot2 ggplot aes geom_point aes_string scale_colour_manual
5
#'   geom_text geom_segment guides guide_legend xlab ylab theme element_blank
6
#'   element_line unit element_text element_rect
7
#' @importFrom igraph get.edgelist degree get.shortest.paths
8
#' @return A plot of the PCA.
9
setGeneric("PlotmclustMB", function(object) {
10
  standardGeneric("PlotmclustMB")
11
})
12
13
#' @export
14
#' @rdname PlotmclustMB
15
setMethod(
16
  "PlotmclustMB",
17
  signature = "DISCBIO",
18
  definition = function(object) {
19
    if (length(object@MBclusters) == 0) {
20
      stop("run ExprmclustMB before PlotmclustMB")
21
    }
22
    total <- object@MBclusters
23
24
    Plotmclust <- function(mclustobj, x = 1, y = 2, MSTorder = NULL,
25
                           show_tree = TRUE, show_full_tree = FALSE,
26
                           show_cell_names = FALSE, cell_name_size = 3,
27
                           markerexpr = NULL,
28
                           showcluster = TRUE) {
29
      lib_info_with_pseudo <- data.frame(
30
        State = mclustobj$clusterid,
31
        sample_name = names(mclustobj$clusterid)
32
      )
33
      lib_info_with_pseudo$State <-
34
        factor(lib_info_with_pseudo$State)
35
      S_matrix <- mclustobj$pcareduceres
36
      pca_space_df <- data.frame(S_matrix[, c(x, y)])
37
      colnames(pca_space_df) <- c("pca_dim_1", "pca_dim_2")
38
      pca_space_df$sample_name <- row.names(pca_space_df)
39
      edge_df <- merge(
40
        pca_space_df, lib_info_with_pseudo,
41
        by.x = "sample_name", by.y = "sample_name"
42
      )
43
      edge_df$Marker <- markerexpr[edge_df$sample_name]
44
      if (!is.null(markerexpr)) {
45
        g <- ggplot(
46
          data = edge_df,
47
          aes(
48
            x = pca_space_df$pca_dim_1,
49
            y = pca_space_df$pca_dim_2,
50
            size = edge_df$Marker
51
          )
52
        )
53
        if (showcluster) {
54
          g <- g + geom_point(
55
            aes_string(color = lib_info_with_pseudo$State),
56
            na.rm = TRUE
57
          )
58
          g <- g + scale_colour_manual(
59
            values = c(
60
              "1" = "black",
61
              "2" = "blue",
62
              "3" = "green",
63
              "4" = "red",
64
              "5" = "yellow",
65
              "6" = "gray"
66
            )
67
          )
68
        } else {
69
          g <- g + geom_point(na.rm = TRUE, color = "green")
70
        }
71
      } else {
72
        g <-
73
          ggplot(
74
            data = edge_df,
75
            aes(
76
              x = pca_space_df$pca_dim_1,
77
              y = pca_space_df$pca_dim_2
78
            )
79
          )
80
        if (showcluster) {
81
          g <- g + geom_point(
82
            aes_string(color = lib_info_with_pseudo$State),
83
            na.rm = TRUE, size = 3
84
          )
85
          g <- g + scale_colour_manual(
86
            values = c(
87
              "1" = "black",
88
              "2" = "blue",
89
              "3" = "green",
90
              "4" = "red",
91
              "5" = "yellow",
92
              "6" = "gray"
93
            )
94
          )
95
        } else {
96
          g <- g + geom_point(na.rm = TRUE, size = 3)
97
        }
98
      }
99
      if (show_cell_names) {
100
        g <- g + geom_text(
101
          aes(label = names(mclustobj$clusterid)),
102
          size = cell_name_size
103
        )
104
      }
105
106
      if (show_tree) {
107
        clucenter <- mclustobj$clucenter[, c(x, y)]
108
        clulines <- NULL
109
        if (show_full_tree) {
110
          alledges <-
111
            as.data.frame(get.edgelist(mclustobj$MSTtree),
112
              stringsAsFactors = FALSE
113
            )
114
          alledges[, 1] <- as.numeric(alledges[, 1])
115
          alledges[, 2] <- as.numeric(alledges[, 2])
116
          for (i in seq_len(nrow(alledges))) {
117
            clulines <- rbind(
118
              clulines,
119
              c(
120
                clucenter[alledges[i, 1], ],
121
                clucenter[alledges[i, 2], ]
122
              )
123
            )
124
          }
125
        } else {
126
          if (is.null(MSTorder)) {
127
            clutable <- table(mclustobj$clusterid)
128
            alldeg <- degree(mclustobj$MSTtree)
129
            allcomb <- expand.grid(
130
              as.numeric(names(alldeg)[alldeg == 1]),
131
              as.numeric(names(alldeg)[alldeg == 1])
132
            )
133
            allcomb <-
134
              allcomb[allcomb[, 1] < allcomb[, 2], ]
135
            numres <- t(apply(allcomb, 1, function(i) {
136
              tmp <- as.vector(
137
                get.shortest.paths(
138
                  mclustobj$MSTtree,
139
                  i[1], i[2]
140
                )$vpath[[1]]
141
              )
142
              c(length(tmp), sum(clutable[tmp]))
143
            }))
144
            orderedRows <- order(
145
              numres[, 1], numres[, 2],
146
              decreasing = TRUE
147
            )[1]
148
            optcomb <- allcomb[orderedRows, ]
149
            MSTorder <- get.shortest.paths(
150
              mclustobj$MSTtree, optcomb[1],
151
              optcomb[2]
152
            )$vpath[[1]]
153
          }
154
          for (i in 1:(length(MSTorder) - 1)) {
155
            clulines <- rbind(
156
              clulines,
157
              c(
158
                clucenter[MSTorder[i], ],
159
                clucenter[MSTorder[i + 1], ]
160
              )
161
            )
162
          }
163
        }
164
        clulines <- data.frame(
165
          x = clulines[, 1],
166
          xend = clulines[, 3],
167
          y = clulines[, 2],
168
          yend = clulines[, 4]
169
        )
170
        g <- g + geom_segment(
171
          aes_string(
172
            x = "x",
173
            xend = "xend",
174
            y = "y",
175
            yend = "yend",
176
            size = NULL
177
          ),
178
          data = clulines,
179
          size = 1,
180
          color = "orange"
181
        )
182
        clucenter <- data.frame(
183
          x = clucenter[, 1],
184
          y = clucenter[, 2],
185
          id = seq_len(nrow(clucenter))
186
        )
187
        g <- g + geom_text(
188
          aes_string(
189
            label = "id",
190
            x = "x",
191
            y = "y",
192
            size = NULL
193
          ),
194
          data = clucenter,
195
          size = 10,
196
          color = "orange"
197
        )
198
      }
199
      g <-
200
        g +
201
        guides(
202
          colour = guide_legend(override.aes = list(size = 5))
203
        ) +
204
        xlab(paste0("PCA_dimension_", x)) +
205
        ylab(paste0("PCA_dimension_", y)) +
206
        theme(
207
          panel.border = element_blank(),
208
          axis.line = element_line()
209
        ) +
210
        theme(
211
          panel.grid.minor.x = element_blank(),
212
          panel.grid.minor.y = element_blank()
213
        ) +
214
        theme(
215
          panel.grid.major.x = element_blank(),
216
          panel.grid.major.y = element_blank()
217
        ) +
218
        theme(
219
          legend.position = "top",
220
          legend.key.size = unit(0.3, "in"),
221
          legend.text = element_text(size = 20),
222
          legend.title = element_text(size = 20),
223
          legend.box = "vertical"
224
        ) +
225
        theme(legend.key = element_blank()) +
226
        theme(panel.background = element_rect(fill = "white")) +
227
        theme(
228
          axis.text.x = element_text(size = 17, color = "black"),
229
          axis.text.y = element_text(size = 17, color = "black"),
230
          axis.title.x = element_text(size = 20, vjust = -1),
231
          axis.title.y = element_text(size = 20, vjust = 1),
232
          plot.margin = unit(c(1, 1, 1, 1), "cm")
233
        )
234
      g
235
    }
236
    Plotmclust(total)
237
  }
238
)