-
Notifications
You must be signed in to change notification settings - Fork 5
/
PlotMBpca.R
104 lines (103 loc) · 3.03 KB
/
PlotMBpca.R
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
#' @title Plotting pseudo-time ordering or gene expression in Model-based
#' clustering in PCA
#' @description The PCA representation can either be used to show pseudo-time
#' ordering or the gene expression of a particular gene.
#' @param object \code{DISCBIO} class object.
#' @param type either `order` to plot pseudo-time ordering or `exp` to plot gene
#' expression
#' @param g Individual gene name or vector with a group of gene names
#' corresponding to a subset of valid row names of the \code{ndata} slot of
#' the \code{DISCBIO} object. Ignored if `type="order"`.
#' @param n String of characters representing the title of the plot. Default is
#' NULL and the first element of \code{g} is chosen. Ignored if
#' `type="order"`.
#' @importFrom RColorBrewer brewer.pal
#' @importFrom grDevices colorRampPalette
#' @importFrom graphics layout par image
#' @return A plot of the PCA.
#' @export
PlotMBpca <- function(object, type = "order", g = NULL, n = NULL) {
# ==========================================================================
# Validation
# ==========================================================================
data <- object@MBclusters
if (type == "exp") {
if (is.null(g)) {
stop('g must be provided if type="exp"')
}
if (length(intersect(g, rownames(object@ndata))) < length(unique(g))) {
stop(
"second argument does not correspond to set of rownames slot",
"ndata of SCseq object"
)
}
if (is.null(n)) {
n <- g[1]
}
logObj <- log(object@ndata)
l <- apply(logObj[g, ] - .1, 2, sum) + .1
x <- data$pcareduceres
} else if (type == "order") {
MBordertable <- cbind(data$pcareduceres, object@MBordering)
l <- MBordertable[, 3]
x <- MBordertable
} else {
stop("Invalid type. Valid alternatives as 'order' and 'exp'.")
}
# ==========================================================================
# Plotting
# ==========================================================================
mi <- min(l, na.rm = TRUE)
ma <- max(l, na.rm = TRUE)
ColorRamp <- colorRampPalette(rev(brewer.pal(n = 11, name = "RdYlBu")))(100)
ColorLevels <- seq(mi, ma, length = length(ColorRamp))
v <- round((l - mi) / (ma - mi) * 99 + 1, 0)
layout(
matrix(
data = c(1, 3, 2, 4),
nrow = 2,
ncol = 2
),
widths = c(5, 1, 5, 1),
heights = c(5, 1, 1, 1)
)
opar <- withr::local_par(mar = c(5, 5, 2.5, 2))
on.exit(withr::local_par(opar))
plot(
x[, 1],
x[, 2],
xlab = "PC1",
ylab = "PC2",
pch = 20,
cex = 0,
col = "grey",
las = 1,
main = n
)
for (k in seq_len(length(v))) {
points(
x[k, 1],
x[k, 2],
col = ColorRamp[v[k]],
pch = 20,
cex = 2
)
}
opar <- withr::local_par(mar = c(3, 2.5, 2.5, 2))
on.exit(withr::local_par(opar))
image(
1,
ColorLevels,
matrix(
data = ColorLevels,
ncol = length(ColorLevels),
nrow = 1
),
col = ColorRamp,
xlab = "",
ylab = "",
las = 2,
xaxt = "n"
)
layout(1)
}