R/hetset_vis.R

In ZytoHMGU/hetset: Identification of Heterogeneous Subsets in Data

# hetset visualization tools ==================================================

# plot heterogeneous networks -------------------------------------------------
plot_hetset <- function(H,plot_in = "no_file",file_name = "hetset_plot_temp"){
  if(is.na(H@metadata$slf[1])){
    H@metadata$slf <- H@NAMES
    H@metadata$color_type <- "raw"
  }
  plot_type <- as.character(min(3,sum(!is.na(H@metadata$slf))))
  if (sum(!is.na(H@metadata$slf))>10){
    plot_type <- "make_table"
  }
  off_flag <- F

  # plot on screen or file (select type)
  switch(plot_in,
         "png" = {
           off_flag <- T
           png(filename = paste0(file_name,".png"),width = 13,height = 13,
               units = "cm",res = 90)
         },
         "svg" = {
           off_flag <- T
           svg(filename = paste0(file_name,".svg"),width = 8,height = 8)
         },
         "pdf" = {
           off_flag <- T
           pdf(file = paste0(file_name,".pdf"),width = 8,height = 8)
         })

  # set colors for plotting
  color_sel <- colors()[c(130,35,17,50,500,578,610,506,640,636,654,96,24)]
  # format labelling and coloring information
  if (is.null(H@metadata$color_type)){ # case: "no type specified"
    if (!is.null(H$sample_color)){ #     subcase: "labels assigned manually"
      H@metadata$color_type <- "indiv"
      H$sample_color <- as.factor(H$sample_color)
    } else {
      if (is.null(H$prt)){         #     subcase: "raw data"
        H@metadata$color_type <- "raw"
        H$sample_color <- NA
      } else {                     #     subcase: "labels to subpopulations"
        H@metadata$color_type <- "subpop"
        H$sample_color <- H$prt
      }
    }
  } else {                            # case: "type specified"
    if (H@metadata$color_type %in% c("subpop","pop","prt","raw")){
      if (H@metadata$color_type == "raw"){ # subcase: "raw"
        H$sample_color <- rep(NA,ncol(H))
      } else {                        #      subcase: "subpop"
        H@metadata$color_type <- "subpop"
        H$sample_color <- H$prt
      }
    } else {
      if (H@metadata$color_type != "indiv"){
        # subcase: name of labeling factor in color_type
        label <- H@metadata$color_type
        eval(parse(text = c("H$sample_color <- H$",label)))
        H$sample_color <- as.factor(H$sample_color)
        rm(label)
        H@metadata$color_type <- "indiv"
      } else {
        H$sample_color <- as.factor(H$sample_color)
      }
    }
  }

  plot.info <- .prep_legend(H = H,color = color_sel)
  # type of plot
  switch(plot_type,
         "1" = {
           .plot_hetset_histogram(plot.info)
         },
         "2" = {
           .plot_hetset_scatterplot(plot.info)
         },
         "3" = {
           .plot_hetset_scattermatrix(plot.info)
         },
         "make_table" = {
           cat(" report for more than 10 selected features \n")
           cat("    ... not yet implemented ...  \n")
         })

  if (off_flag){
    dev.off()
  } else {
    par(mfrow = c(1,1),mar = c(5,4,4,2)+0.1,mgp = c(3,1,0))
  }
}

# prepare legend, colors and plotting order -----------------------------------
.prep_legend <- function(H,color){
  na.group.term <- "unlabelled"
  # legend for 1 selected feature
  if (sum(!is.na(H@metadata$slf)) == 1){
    switch(H@metadata$color_type,
           "raw" = {
             legend_info <- data.frame("typ" = "H.A",
                                       "ID" = NA,
                                       "tex" = "all samples",
                                       "pch" = NA,
                                       "lwd" = NA,
                                       "col" = NA,
                                       "fil" = color[3],
                                       "lty" = NA,
                                       "bor" = "black")
             H$sample_color <- 1
           },
           "subpop" = {
             legend_text <- c("all samples",
                              "subpopulation A",
                              "subpopulation B",
                              na.group.term)
             legend_info <- data.frame("typ" = c("H.A","H.S","H.S","H.S"),
                                       "ID" = c(NA,"A","B",na.group.term),
                                       "tex" = legend_text,
                                       "pch" = NA,
                                       "lwd" = NA,
                                       "col" = c("white",color[1:3]),
                                       "fil" = c("white",color[1:3]),
                                       "lty" = NA,
                                       "bor" = "black")
             H$sample_color <- factor(x = H$sample_color,
                                      levels = c(levels(x = H$sample_color),
                                                 na.group.term))
             if (sum(is.na(H$sample_color))==0){
               legend_info <- legend_info[-4,]
             } else {
               H$sample_color[is.na(H$sample_color)] <- na.group.term
             }
             legend_info$ID <- droplevels(legend_info$ID)
             H$sample_color <- droplevels(H$sample_color)
           },
           "indiv" = {
             color <- color[-(1:2)]
             gr.distr <- table(H$sample_color)
             Gi <- names(-sort(-gr.distr))
             H$sample_color <- factor(x = H$sample_color,
                                      levels = c(levels(x = H$sample_color),
                                                 na.group.term))
             H$sample_color[is.na(H$sample_color)] <- na.group.term
             Gi <- c(Gi,na.group.term)
             legend_text <- c("all samples",Gi)

             legend_info <- data.frame("typ" = c("H.A",rep("H.S",length(Gi))),
                                       "ID" = c(NA,Gi),
                                       "tex" = legend_text,
                                       "pch" = NA,
                                       "lwd" = NA,
                                       "col" = NA,
                                       "fil" = c("white",color[2:length(Gi)],
                                                 color[1]),
                                       "lty" = NA,
                                       "bor" = "black")
             rm(gr.distr,color,Gi,legend_text)

             if (sum(H$sample_color==na.group.term)==0){
               legend_info <- legend_info[-nrow(legend_info),]
             }
             H$sample_color <- factor(H$sample_color,
                                      levels = legend_info$"ID"[2:nrow(legend_info)])
             legend_info$ID <- droplevels(legend_info$ID)
           })
  } else { # legend for scatterplots
    no.dens <- is.na(H@metadata$prm.A$mean[1])
    switch(H@metadata$color_type,
           "raw" = {
             legend_text <- c(" ","location subpop. A",
                              "location subpop. B","location all samples")
             legend_info <- data.frame("typ" = c("H.A","d","d","d"),
                                       "ID" = c(na.group.term,NA,NA,NA),
                                       "tex" = legend_text,
                                       "pch" = c(16,NA,NA,NA),
                                       "lwd" = c(NA,2,2,2),
                                       "col" = color[c(3,1,2,3)],
                                       "fil" = NA, #color[c(3,1,2,3)],
                                       "lty" = c(0,1,1,1),
                                       "bor" = NA)
             H$sample_color <- as.factor(na.group.term)
             if (no.dens){
               legend_info <- legend_info[1,]
             }
             legend_info$ID <- droplevels(legend_info$ID)
           },
           "subpop" = {
             legend_text <- c("subpopulation A",
                              "subpopulation B",
                              na.group.term,
                              "location subpop. A",
                              "location subpop. B",
                              "location all samples")
             legend_info <- data.frame("typ" = c("H.S","H.S","H.S","d","d","d"),
                                       "ID" = c("A","B",na.group.term,NA,NA,NA),
                                       "tex" = legend_text,
                                       "pch" = c(rep(16,3),NA,NA,NA),
                                       "lwd" = c(rep(NA,3),2,2,2),
                                       "col" = color[c(1,2,3,1,2,3)],
                                       "fil" = NA,
                                       "lty" = c(rep(NA,3),1,1,1),
                                       "bor" = NA)
             if (no.dens){
               legend_info <- legend_info[1:3,]
             }
             H$sample_color <- factor(x = H$sample_color,
                                      levels = c(levels(x = H$sample_color),
                                                 na.group.term))
             H$sample_color[is.na(H$sample_color)] <- na.group.term
             if (sum(H$sample_color==na.group.term)==0){
               legend_info <- legend_info[-3,]
             }
             H$sample_color <- factor(H$sample_color,
                                      levels = legend_info$"ID"[1:sum(legend_info$typ != "d")])
             legend_info$ID <- droplevels(legend_info$ID)
           },
           "indiv" = {
             color.loc <- color[1:3]
             color <- color[-(1:2)]
             gr.distr <- table(H$sample_color)
             Gi <- names(-sort(-gr.distr))
             H$sample_color <- factor(x = H$sample_color,
                                      levels = c(levels(x = H$sample_color),
                                                 na.group.term))
             H$sample_color[is.na(H$sample_color)] <- na.group.term
             Gi <- c(Gi,na.group.term)
             legend_text <- c(Gi,"location subpop. A","location subpop. B",
                              "location all samples")
             lGi <- length(Gi)

             legend_info <- data.frame("typ" = c(rep("H.S",lGi),"d","d","d"),
                                       "ID" = c(Gi,NA,NA,NA),
                                       "tex" = legend_text,
                                       "pch" = c(rep(16,lGi),NA,NA,NA),
                                       "lwd" = c(rep(NA,lGi),2,2,2),
                                       "col" = c(color[2:lGi],1,color.loc),
                                       "fil" = NA,
                                       "lty" = c(rep(0,lGi),1,1,1),
                                       "bor" = NA)
             rm(gr.distr,color,color.loc,Gi,legend_text)
             if (no.dens){
               legend_info <- legend_info[1:lGi,]
             }
             if (sum(H$sample_color==na.group.term)==0){
               legend_info <- legend_info[-lGi,]
             }
             H$sample_color <- factor(H$sample_color,
                                      levels = legend_info$"ID"[1:nrow(legend_info)])
             legend_info$ID <- droplevels(legend_info$ID)
           })
  }

  return(list("H" = H,"legend_info" = legend_info))
}

# plot histogram --------------------------------------------------------------
.plot_hetset_histogram <- function(plot.info){
  H <- plot.info$H
  info <- plot.info$legend_info
  rm(plot.info)

  x <- t(assays(H[H@metadata$slf,])[[1]])
  x_min <- min(x,na.rm = T)
  x_max <- max(x,na.rm = T)
  x_breaks <- seq(x_min,x_max,(x_max-x_min)/20)
  hist_out <- hist(x,breaks = x_breaks,
                   main = paste0("Histogram of ",H@metadata$slf),
                   xlab = paste0("log(",H@metadata$slf,")"))

  legend_fill <- rgb(t(col2rgb("white")),maxColorValue = 255)
  if (nrow(info) > 1){
    for (i in 2:nrow(info)){
      Gh <- H$sample_color==info$ID[i]
      hist(x[Gh],breaks = x_breaks,add = T,
           col = rgb(t(col2rgb(info$col[i])),alpha = 125,maxColorValue = 255),
           ann=F,xaxt='n',yaxt='n')
      legend_fill <- c(legend_fill,rgb(t(col2rgb(info$col[i])),alpha = 125,
                                       maxColorValue = 255))
    }

  }
  legend("topright",legend = info$tex,bty = 'n',fill = legend_fill)
}

# make scatterplot ------------------------------------------------------------
.plot_hetset_scatterplot <- function(plot.info){
  H <- plot.info$H
  info <- plot.info$legend_info
  rm(plot.info)

  # set layout
  layout(matrix(c(2,0,1,3),2,2,byrow = T),widths = c(3,1),heights = c(1,3))
  par(mar=c(2,2,0,0)+0.1)
  # calculate adjustment-points for plot modules
  n_bars <- 20
  x_barwidth <- diff(range(assays(H[H@metadata$slf[1],])[[1]],na.rm = T))/n_bars
  x_lowadj <- min(assays(H[H@metadata$slf[1],])[[1]],na.rm = T)-x_barwidth/2
  x_uppadj <- max(assays(H[H@metadata$slf[1],])[[1]],na.rm = T)+x_barwidth/2
  y_barwidth <- diff(range(assays(H[H@metadata$slf[2],])[[1]],na.rm = T))/n_bars
  y_lowadj <- min(assays(H[H@metadata$slf[2],])[[1]],na.rm = T)-y_barwidth/2
  y_uppadj <- max(assays(H[H@metadata$slf[2],])[[1]],na.rm = T)+y_barwidth/2

  # scatterplot
  N.cols <- sum(info$typ!="d")
  plot(x = assays(H[H@metadata$slf[1],])[[1]],
       y = assays(H[H@metadata$slf[2],])[[1]],
       type = 'o',xlim = c(x_lowadj,x_uppadj),ylim = c(y_lowadj,y_uppadj),
       main = NULL,xlab = 'n',ylab = 'n',ann = F,
       pch = info$pch[1:N.cols],lty = 0,
       col = rgb(t(col2rgb(info$col[H$sample_color])),alpha = 125,
                 maxColorValue = 255))
  # add densities
  if (sum(info$typ=="d")>0){
    .add_ellipses(He = H,info$col[info$"typ"=="d"])
  }

  # histograms x- and y-data --------------------------------------------------
  .add_histogram(H,var_no = 1,adj = c(x_lowadj,x_uppadj),n_bars = n_bars,
                info = info)
  .add_histogram(H,var_no = 2,adj = c(y_lowadj,y_uppadj),n_bars = n_bars,
                info = info)
}

# add histograms to hetset-scatterplot ----------------------------------------
.add_histogram <- function(H,var_no,adj,n_bars,info){
  if(var_no==1){
    par(mar=c(0,2,0.42,0)+0.1)
  } else {
    par(mar=c(2,0,0,0.42)+0.1)
  }
  m_breaks <- seq(adj[1],adj[2],diff(adj)/n_bars)
  hist_out <- hist(assays(H[H@metadata$slf[var_no],])[[1]],
                   breaks = m_breaks,plot=F)
  barplot(height = hist_out$counts,horiz = (var_no==2),space = 0,col="white",
          axes = F)
  title(main = H@metadata$slf[var_no],line = -1)

  for (i in 1:sum(info$typ=="H.S")){
    Gh <- H$sample_color==info$ID[i]
    barplot(height = hist(assays(H[H@metadata$slf[var_no],Gh])[[1]],
                          breaks = m_breaks,plot = F)$counts,
            horiz = (var_no==2),space = 0,axes = F,
            col=rgb(t(col2rgb(info$col[i])),alpha = 125,maxColorValue = 255),
            add=T,ann=F,xaxt='n',yaxt='n')
  }

  if(var_no==1){
    legend("topright",legend = info$tex,box.lwd = 0,
           bg = rgb(t(col2rgb("white")),alpha = 128,maxColorValue = 255),
           col = rgb(t(col2rgb(info$col)),alpha = 175,maxColorValue = 255),
           border = info$bor,
           pch = info$pch,lty = info$lty)
  }
}

# make matrix of scatterplots -------------------------------------------------
.plot_hetset_scattermatrix <- function(plot.info){
  H <- plot.info$H
  info <- plot.info$legend_info
  rm(plot.info)
  N.cols <- sum(info$typ!="d")

  # organize plot
  paare <- t(combn(x = H@metadata$slf,m = 2))
  n_paare <- dim(paare)[1]
  par(mfrow = c(rep(ceiling(sqrt(n_paare+1)),2)))
  par(mgp = c(-1.2,0,0))
  par(mar = rep(0.2,4))
  for(i in 1:n_paare){
    namen <- paare[i,]
    x_lowadj <- min(assays(H[namen[1],])[[1]],na.rm = T)
    x_uppadj <- max(assays(H[namen[1],])[[1]],na.rm = T)
    y_lowadj <- min(assays(H[namen[2],])[[1]],na.rm = T)
    y_uppadj <- max(assays(H[namen[2],])[[1]],na.rm = T)

    plot(x = assays(H[namen[1],])[[1]],
         y = assays(H[namen[2],])[[1]],
         type = 'o',xlim = c(x_lowadj,x_uppadj),ylim = c(y_lowadj,y_uppadj),
         main = NULL,xlab = namen[1],ylab = namen[2],ann = T,cex.lab = 1.2,
         xaxt = 'n',yaxt = 'n',pch = info$pch[1:N.cols],lty = 0,
         col = rgb(t(col2rgb(info$col[H$sample_color])),alpha = 125,
                   maxColorValue = 255))
    par(new=F)
    # add densities
    if (sum(info$typ=="d")>0){
      .add_ellipses(He = subset_hetset(H = H,keep.features = namen),
                   color = info$col[info$"typ"=="d"])
    }
  }
  plot(x = c(0,1),y = c(0,1),xaxt='n',yaxt='n',type = 'n',ann=F,axes=F)
  legend("topright",legend = info$tex,bty = 'n',
         col = rgb(t(col2rgb(info$col)),alpha = 175,maxColorValue = 255),
         border = info$bor,
         fill = as.character(info$fil),
         lty = info$lty,
         pch = info$pch)
}

# add ellipses and correlation values to scatterplots -------------------------
.add_ellipses <- function(He,color){
  subpops <- c("A","B")
  for (i in 1:length(subpops)){
    X <- assays(He[He@metadata$slf[1],He$prt==subpops[i]])[[1]]
    mX <- mean(X)
    Y <- assays(He[He@metadata$slf[2],He$prt==subpops[i]])[[1]]
    mY <- mean(Y)
    if ((sd(X,na.rm = T) > 0) & (sd(Y,na.rm = T) > 0)){
      cor_out <- cor.test(x = X,y = Y)
      if (cor_out$p.value<0.05){
        text(x = mX,y = mY,
             labels = round(cor_out$estimate,digits = 2),col = "black")
      }
      mixtools::ellipse(mu = c(mX,mY),
                  sigma = cov(x = cbind(t(X),t(Y)),use = "complete.obs"),
                  alpha = 0.5,newplot = F,draw = T,
                  col = as.character(color[i]),lwd=2)
    }
  }
  if ((sd(assays(He[He@metadata$slf[1],])[[1]]) > 0) &
      (sd(assays(He[He@metadata$slf[2],])[[1]]) > 0)){
    corF_out <- cor.test(x = assays(He[He@metadata$slf[1],])[[1]],
                         y = assays(He[He@metadata$slf[2],])[[1]])
    if (corF_out$p.value<0.05){
      text(x = He@metadata$prm.full$mean[1],y = He@metadata$prm.full$mean[2],
           labels = round(corF_out$estimate,digits = 2),col = "black")
    }
    mixtools::ellipse(mu = He@metadata$prm.full$mean,
                sigma = He@metadata$prm.full$cov,
                alpha = 0.5,newplot = F,draw = T,
                col = as.character(color[3]),lty=1,lwd=1)
  }
}