updated pseudobulk, pbfacet, quad split

R/plotting.R

@@ -32,28 +32,35 @@ savePlot <- function(p, pltprefix, w=7, h=7, dpi=600, ...){
   dev.off()
   print(pltprefix)
 }
+#' Heatmap sort matrix
+#'
+#' @param M input matrix
+#' @param cutoff Value cutoff for counting rows as significant to be counted
+#' @param ratio Ratio of of rows above cutoff to be counted
+#' @param method Method to be used for column/row ordering
+#' @param sort Whether to sort by rows, columns, or both
 #' @export
 htSortMatrix <- function(M, method="euclidean", ratio=0.5, cutoff=0.25, sort=c(1,2)) {
     if (is.logical(sort)) {
         if (sort) {
-            sort = c(1,2)
+            sort <- c(1,2)
         } else {
-            sort = c()
+            sort <- c()
         }
     }
-    rows = 1 %in% sort
-    columns = 2 %in% sort
+    rows <- 1 %in% sort
+    columns <- 2 %in% sort
     if (rows) {
-        M = order.tsp(M, rows=TRUE, method=method)
+        M <- order.tsp(M, rows=TRUE, method=method)
     }
     if (columns) {
-        M = order.tsp(M, rows=FALSE, method=method)
+        M <- order.tsp(M, rows=FALSE, method=method)
     }
     if (rows) {
-        M = diag.mat3(M, rows=TRUE, ratio=ratio, cutoff=cutoff)
+        M <- diag.mat3(M, rows=TRUE, ratio=ratio, cutoff=cutoff)
     }
     if (columns) {
-        M = diag.mat3(M, rows=FALSE, ratio=ratio, cutoff=cutoff)
+        M <- diag.mat3(M, rows=FALSE, ratio=ratio, cutoff=cutoff)
     }
     return(M)
 }
@@ -155,6 +162,32 @@ ht_triangle_split <- function(mat.ul, mat.lr, col.ul, col.lr, lwd=0, ...) {
     })
 }
 
+#' @export
+ht_quad_triangle_split <- function(mat.top, mat.right, mat.bot, mat.left, col.top, col.right, col.bot, col.left, lwd=0, ...) {
+    return(function(j, i, x, y, width, height, fill) {
+        ### start with UL, UR, LR, LL
+        corners = data.frame(x=as.numeric(width) * c(-1/2, 1/2, 1/2, -1/2, 0),
+                             y=as.numeric(height) * c(-1/2, -1/2, 1/2, 1/2, 0),
+                             row.names=c("LL", "UL", "UR", "LR", "C"))
+        ### top
+        grid::grid.polygon(x=as.numeric(x) + corners[c("UL", "UR", "C", "UL"),"x"],
+                           y=as.numeric(y) + corners[c("UL", "UR", "C", "UL"),"y"],
+                           gp=grid::gpar(fill=col.top(mat.top[i, j]), lwd=lwd, ...))
+        ### right
+        grid::grid.polygon(x=as.numeric(x) + corners[c("UR", "LR", "C", "UR"),"x"],
+                           y=as.numeric(y) + corners[c("UR", "LR", "C", "UR"),"y"],
+                           gp=grid::gpar(fill=col.right(mat.right[i, j]), lwd=lwd, ...))
+        ### bottom
+        grid::grid.polygon(x=as.numeric(x) + corners[c("LR", "LL", "C", "LR"),"x"],
+                           y=as.numeric(y) + corners[c("LR", "LL", "C", "LR"),"y"],
+                           gp=grid::gpar(fill=col.bot(mat.bot[i, j]), lwd=lwd, ...))
+        ### left
+        grid::grid.polygon(x=as.numeric(x) + corners[c("LL", "UL", "C", "LL"),"x"],
+                           y=as.numeric(y) + corners[c("LL", "UL", "C", "LL"),"y"],
+                           gp=grid::gpar(fill=col.left(mat.left[i, j]), lwd=lwd, ...))
+    })
+}
+
 #' Draw asterisks
 #' Pass cell_fun=ht_asterisks(...) to Heatmap()
 #' @param p.matrix P value matrix, unsorted
@@ -211,6 +244,37 @@ volcano <- function(df, label="gene", title="Volcano plot of",
     return(g)
 }
 
+#' @export
+pbviolin_facet <- function(sce, gene, groupby, facet.by, facet.keep=NULL, prefix="Gene", sample.col="Sample", target_sum=10000, max_pb=TRUE) {
+    require(ggpubr)
+    if (!is.null(facet.keep)) {
+        sce <- sce[,sce[[facet.by]] %in% facet.keep]
+    }
+    X = SummarizedExperiment::assays(sce)$counts
+    cd = as.data.frame(SummarizedExperiment::colData(sce))
+    if ("total_counts" %in% colnames(cd)) {
+        D = Matrix::Diagonal(x=target_sum / cd$total_counts)
+    } else {
+        D = Matrix::Diagonal(x=target_sum / Matrix::colSums(X))
+    }
+    X = log1p(as.matrix(X[gene,,drop=FALSE] %*% D))
+    uid <- paste0(cd[[groupby]], "#", cd[[sample.col]], "#", cd[[facet.by]])
+    P <- X %*% make_average(uid)
+    pd <- cd[!duplicated(uid), c(groupby, facet.by, sample.col)]
+    rownames(pd) <- uid[!duplicated(uid)]
+    L <- list()
+    for (g in rownames(X)) {
+        cd$gene <- X[g,]
+        pd$gene <- P[g, rownames(pd)]
+        G <- ggviolin(cd, x=groupby, y="gene", fill=groupby, ylim=c(min(P[g,]), max(P[g,]))) + geom_point(data=pd, position=position_jitter(width=0.2))
+        L[[g]] <- facet(G, facet.by=facet.by, nrow=1) + ggtitle(paste0(prefix, " ", g)) + ylab(g)
+        ## if (max_pb) {
+        ##     L[[g]] <- ggpar(L[[g]], ylim=c(min(P), max(P)))
+        ## }
+    }
+    return(L)
+}
+
 #' @export
 pbviolin <- function(sce, gene, groupby, sample.col="Sample", target_sum=10000, plot=TRUE) {
     require(ggplot2)

R/sc.R

@@ -7,7 +7,7 @@
 #' @param se SummarizedExperiment
 #' @param on Column in colData() to split upon
 #' @export
-se_make_pseudobulk <- function(se, on, missing.levels=FALSE, unlevel=TRUE) {
+se_make_pseudobulk <- function(se, on, missing.levels=FALSE, unlevel=TRUE, ncell_col="ncell") {
     cd = SummarizedExperiment::colData(se)
     rd = SummarizedExperiment::rowData(se)
     A = SummarizedExperiment::assays(se)
@@ -17,6 +17,7 @@ se_make_pseudobulk <- function(se, on, missing.levels=FALSE, unlevel=TRUE) {
     })
     cd[[on]] = factor(as.character(cd[[on]]), levels=colnames(to_P))
     pcd = data.frame(row.names=levels(cd[[on]]))
+    pcd[[ncell_col]] <- Matrix::colSums(to_P)
     for (cn in colnames(cd)) {
 ### for each "on", ensure all equal
         u = lapply(split(cd[[cn]], cd[[on]]), unique)
@@ -271,7 +272,7 @@ se_as_Seurat <- function(sce) {
     C = SummarizedExperiment::assays(sce)$counts
     cd = SummarizedExperiment::colData(sce)
     s.obj = Seurat::CreateSeuratObject(counts=C, meta.data=as.data.frame(cd))
-    s.obj[["X_umap"]] <- Seurat::CreateDimReducObject(
+    s.obj[["umap"]] <- Seurat::CreateDimReducObject(
          embeddings = SingleCellExperiment::reducedDims(sce)$X_umap,
          key = "UMAP_", assay = "RNA")
     s.obj[["X_pca"]] <- Seurat::CreateDimReducObject(
@@ -280,3 +281,9 @@ se_as_Seurat <- function(sce) {
     s.obj = Seurat::NormalizeData(s.obj)
     return(s.obj)
 }
+
+ad_as_Seurat <- function(h5ad) {
+  s.obj = anndataR::read_h5ad(h5ad, to="Seurat")
+  obsm = rhdf5::h5read(h5ad, "obsm")
+  s.obj[["X_umap"]] <- Seurat::CreateDimReducObject(embeddings=t(obsm$X_umap), key="UMAP_", assay="RNA")
+}

R/util_matrix.R

@@ -97,11 +97,11 @@ make_average <- function(x, u=NULL, unlevel=FALSE) {
 #' @export
 order.tsp <- function(mat, rows=TRUE, method="euclidean") {
     if (rows) {
-        tsp = seriation::seriate(dist(mat, method=method))
+        tsp <- seriation::seriate(dist(mat, method=method))
     } else {
-        tsp = seriation::seriate(dist(t(mat), method=method))
+        tsp <- seriation::seriate(dist(t(mat), method=method))
     }
-    ord = seriation::get_order(tsp, 1)
+    ord <- seriation::get_order(tsp, 1)
     if (rows) {
         return(mat[ord,])
     } else {
@@ -111,43 +111,43 @@ order.tsp <- function(mat, rows=TRUE, method="euclidean") {
 
 #' Diagonally cluster the matrix
 #'
-#' This function clusters a matrix
+#' This function clusters a matrix, by either rows or by columns
 #'
 #' @param mat Input matrix
 #' @param ratio Ratio
 #' @param cutoff Cutoff for value being counted
 #' @return A list of the sorted matrix and associated columns
 #' @export
 diag.mat3 <- function(mat, ratio=0.5, cutoff=0.25, rows=TRUE) {
-    prev.cutoff = attr(mat, "cutoff")
+    prev.cutoff <- attr(mat, "cutoff")
     if (rows) {
-        ord = order(rowSums(mat > cutoff)  > ratio * ncol(mat),
+        ord <- order(rowSums(mat > cutoff)  > ratio * ncol(mat),
                     apply(mat,1,which.max), decreasing=TRUE)
-        mat = mat[ord,]
-        cto = apply(mat, 1, which.max)
-        idx = rowSums(mat > cutoff)  > ratio * ncol(mat)
-        cto[idx] = 0
+        mat <- mat[ord,]
+        cto <- apply(mat, 1, which.max)
+        idx <- rowSums(mat > cutoff)  > ratio * ncol(mat)
+        cto[idx] <- 0
     } else {
-        ord = order(colSums(mat > cutoff)  > ratio * nrow(mat),
+        ord <- order(colSums(mat > cutoff)  > ratio * nrow(mat),
                     apply(mat,2,which.max), decreasing=TRUE)
-        mat = mat[,ord]
-        cto = apply(mat, 2, which.max)
-        idx = colSums(mat > cutoff)  > ratio * nrow(mat)
-        cto[idx] = 0
+        mat <- mat[,ord]
+        cto <- apply(mat, 2, which.max)
+        idx <- colSums(mat > cutoff)  > ratio * nrow(mat)
+        cto[idx] <- 0
 
     }
-    attr(mat, "cutoff") = prev.cutoff
+    attr(mat, "cutoff") <- prev.cutoff
     if (is.null(attr(mat, "cutoff"))) {
         if (rows) {
-            attr(mat, "cutoff") = list(row=cto, col=NULL)
+            attr(mat, "cutoff") <- list(row=cto, col=NULL)
         } else {
-            attr(mat, "cutoff") = list(row=NULL, col=cto)
+            attr(mat, "cutoff") <- list(row=NULL, col=cto)
         }
     } else {
         if (rows) {
-            attr(mat, "cutoff")$row = cto
+            attr(mat, "cutoff")$row <- cto
         } else {
-            attr(mat, "cutoff")$col = cto
+            attr(mat, "cutoff")$col <- cto
         }
     }
     return(mat)