[381c22]: / R / PlotMBpca.R

Download this file

105 lines (103 with data), 3.1 kB

  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)
}