/
cor_adjust.R
67 lines (65 loc) · 1.98 KB
/
cor_adjust.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
#' Calculates gene-to-gene and cell-to-cell SAVER correlation
#'
#' Adjusts for SAVER estimation uncertainty by calculating and adjusting
#' gene-to-gene and cell-to-cell correlation matrices
#'
#' The SAVER estimates that are produced have varying levels of uncertainty
#' depending on the gene and the cell. These functions adjust the gene-to-gene
#' and cell-to-cell correlations of the SAVER estimates to reflect the
#' estimation uncertainty.
#'
#' @param x A \code{saver} object.
#'
#' @param cor.mat If a correlation matrix of the SAVER estimates was already
#' obtained, then it can be provided as an input to avoid recomputation.
#'
#' @return An adjusted correlation matrix.
#'
#' @examples
#' data("linnarsson_saver")
#'
#' gene.cor <- cor.genes(linnarsson_saver)
#'
#' @importFrom stats cor
#' @export
#' @rdname cor_adjust
cor.genes <- function(x, cor.mat = NULL) {
if (!is.null(x$alpha)) {
return(cor.genes.old(x, cor.mat))
}
if (is.null(cor.mat)) {
message("Calculating correlation matrix...")
cor.mat <- cor(t(x$estimate))
}
ngenes <- nrow(x$estimate)
adj.vec <- rep(0, ngenes)
for (i in 1:ngenes) {
adj.vec[i] <- sqrt(var(x$estimate[i, ], na.rm = TRUE)/
(var(x$estimate[i, ], na.rm = TRUE) +
mean(x$se[i, ]^2, na.rm = TRUE)))
}
adj.mat <- outer(adj.vec, adj.vec)
cor.adj <- adj.mat*cor.mat
return(cor.adj)
}
#' @export
#' @rdname cor_adjust
cor.cells <- function(x, cor.mat = NULL) {
if (!is.null(x$alpha)) {
return(cor.cells.old(x, cor.mat))
}
if (is.null(cor.mat)) {
message("Calculating correlation matrix...")
cor.mat <- cor(x$estimate)
}
ncells <- ncol(x$estimate)
adj.vec <- rep(0, ncells)
for (i in 1:ncells) {
adj.vec[i] <- sqrt(var(x$estimate[, i], na.rm = TRUE)/
(var(x$estimate[, i], na.rm = TRUE) +
mean(x$se[, i]^2, na.rm = TRUE)))
}
adj.mat <- outer(adj.vec, adj.vec)
cor.adj <- adj.mat*cor.mat
return(cor.adj)
}