R/visualization.R

Defines functions plot_indel_binding makepfull makemotif makepwm makevector rep.col revert.columns

Documented in plot_indel_binding

revert.columns <- function(mat) {
  mat[, rev(seq(ncol(mat)))]
}

rep.col <- function(x, n) {
  matrix(rep(x, each = n), ncol = n, byrow = TRUE)
}

makevector <- function(x, strand_sign) {
  if (strand_sign == 1) {
    if (x == 1) {
      a <- c(1, 0, 0, 0)
    } else if (x == 2) {
      a <- c(0, 1, 0, 0)
    } else if (x == 3) {
      a <- c(0, 0, 1, 0)
    } else{
      a <- c(0, 0, 0, 1)
    }
  } else{
    if (x == 1) {
      a <- c(0, 0, 0, 1)
    } else if (x == 2) {
      a <- c(0, 0, 1, 0)
    } else if (x == 3) {
      a <- c(0, 1, 0, 0)
    } else{
      a <- c(1, 0, 0, 0)
    }
  }
  a
}

makepwm <- function(x, strand_sign) {
  a <- mapply(makevector, x, strand_sign)
  a <- unlist(a)
  b <- matrix(a, ncol = length(x))
  b
}

makemotif <- function(x, strand_sign) {
  p <- t(x[[1]])
  if (strand_sign == 1) {
    rownames(p) <- c("A", "C", "G", "T")
  } else{
    p <- revert.columns(p)
    rownames(p) <- c("A", "C", "G", "T")
  }
  p
}

makepfull <- function(p, left, right) {
  if (left > 0 & right > 0) {
    p1 <-
      cbind(rep.col(c(0.25, 0.25, 0.25, 0.25), left), p, rep.col(c(0.25, 0.25, 0.25, 0.25), right))
    p1
  } else if (left > 0 & right == 0) {
    p1 <- cbind(rep.col(c(0.25, 0.25, 0.25, 0.25), left), p)
    p1
  } else{
    p1 <- cbind(p, rep.col(c(0.25, 0.25, 0.25, 0.25), right))
    p1
  }
}


#' @name plot_indel_binding
#' @title Plot the indel binding for both reference and mutated alleles.
#' @param indel_info A list of length 1 containing the indel information. See
#' \link{indel_info} for details and an example.
#' @param motif_scores A list of the same structure as the '$list' field in the
#' output of \code{\link{indel_motif_scores}}.
#' @param motif A list of length 1 containing the motif PWM. The
#' name of this element is the motif name.
#' @examples
#' \dontrun{data(example)
#' motif_scores <- indel_motif_scores(motif_lib, indel_info)$list
#' plot_indel_binding(
#'    indel_info[1],
#'    motif_scores,
#'    motif_lib[1]
#' )
#' }
#' @import grid
#' @importFrom motifStack plotMotifLogo
#' @author Qinyi Zhou \email{qinyi.zhou@utdallas.edu},
#' Sunyoung Shin \email{sunyoung.shin@@utdallas.edu}
#' @export
plot_indel_binding <- function(indel_info, motif_scores, motif) {
  if (length(indel_info) != 1) {
    stop("indel_info must have length 1.")
  }
  validate_motif_scores(motif_scores, names(motif), names(indel_info))
  j1 <- which(names(motif_scores$ref) == names(indel_info))
  j2 <- which(motif_scores$motif == names(motif))
  if (length(j1) != 1) {
    stop("Cannot find the scores corresponding to 'indel_info' in 'motif_scores'.")
  }
  if (length(j2) != 1) {
    stop("Cannot find the scores corresponding to 'motif' in 'motif_scores'.")
  }
  insertion <- indel_info[[1]]$insertion
  m <- indel_info[[1]]$insertion_len
  long_best_match <- motif_scores$match_pos_long[j1, j2]
  long_strand_sign <- sign(long_best_match)
  seque <- indel_info[[1]]$inserted_sequence
  n <- length(seque)
  n1 <- (n - m) / 2
  lm <- nrow(motif[[1]])
  s1 <- n1 - lm + 2
  e1 <- n1 + m + lm - 1
  longer <- seque[s1:e1]
  nlong <- length(longer)
  short_best_match <- motif_scores$match_pos_short[j1, j2]
  short_strand_sign <- sign(short_best_match)
  s2 <- n1
  e2 <- n1 + m + 1
  shorter <- seque[c(s1:s2, e2:e1)]
  nshort <- length(shorter)
  plot.new()
  if (insertion) {
    if (long_strand_sign == 1) {
      left <- long_best_match - s1
      right <- nlong - left - lm
      a <- makepwm(longer, long_strand_sign)
      rownames(a) <- c("A", "C", "G", "T")
      p <- makemotif(motif, long_strand_sign)
      p_full <- makepfull(p, left, right)
      pushViewport(viewport(y = 0.125, height = 0.25))
      motifStack::plotMotifLogo(
        p_full,
        yaxis = FALSE,
        xaxis = FALSE,
        xlab = "",
        ylab = "",
        newpage = FALSE,
        margins = c(2, 2.5, 1.5, 2.5)
      )
      grid.lines(
        x = c(
          convertUnit(unit(3, "lines"),
                      "npc", valueOnly = TRUE),
          1 - convertUnit(unit(2,
                               "lines"), "npc", valueOnly = TRUE)
        ),
        y = unit(1, "lines"),
        gp = gpar(
          col = "blue",
          lwd = 1.5,
          xpd = NA
        ),
        arrow = arrow(
          length = unit(0.1,
                        "inches"),
          angle = 15,
          ends = "last"
        )
      )
      grid.text(
        "3'",
        x = unit(1, "npc") -
          unit(1, "lines"),
        gp = gpar(col = "blue",
                  cex = 1),
        y = unit(0.5, "lines")
      )
      grid.text(
        "5'",
        x = unit(2, "lines"),
        gp = gpar(col = "blue", cex = 1),
        y = unit(0.5,
                 "lines")
      )
      popViewport()
      pushViewport(viewport(y = 0.325, height = 0.25))
      motifStack::plotMotifLogo(
        a,
        xaxis = FALSE,
        yaxis = FALSE,
        xlab = "",
        ylab = "(+)",
        newpage = FALSE,
        margins = c(2, 2, 1.5, 2)
      )
      grid.text(
        "3'",
        x = unit(1, "npc") -
          unit(1, "lines"),
        gp = gpar(col = "blue",
                  cex = 1),
        y = unit(0.5, "lines")
      )
      grid.text(
        "5'",
        x = unit(2, "lines"),
        gp = gpar(col = "blue", cex = 1),
        y = unit(0.5,
                 "lines")
      )
      grid.rect(
        just = "center",
        width = (m) / (nlong +
                         2),
        height = unit(0.6, "npc"),
        gp = gpar(
          col = "blue",
          lty = 3,
          lwd = 2,
          fill = NA
        )
      )
      if (short_strand_sign == 1) {
        a <- makepwm(shorter, short_strand_sign)
        rownames(a) <- c("A", "C", "G", "T")
        l1 <- floor(m / 2)
        r1 <- m - l1
        left <- short_best_match - s1 + l1
        right <- nlong - left - lm
        if (l1 == 0) {
          a_full <- cbind(a, rep.col(c(0.25, 0.25, 0.25,
                                       0.25), r1))
        }
        else {
          a_full <- cbind(rep.col(c(0.25, 0.25, 0.25,
                                    0.25), l1), a, rep.col(c(0.25, 0.25, 0.25,
                                                             0.25), r1))
        }
        p <- makemotif(motif, short_strand_sign)
        p_full <- makepfull(p, left, right)
        popViewport()
        pushViewport(viewport(y = 0.575, height = 0.25))
        motifStack::plotMotifLogo(
          a_full,
          xaxis = FALSE,
          yaxis = FALSE,
          xlab = "",
          ylab = "(+)",
          newpage = FALSE,
          margins = c(2, 2.5, 1.5, 2.5)
        )
        grid.text(
          "3'",
          x = unit(1, "npc") -
            unit(1, "lines"),
          gp = gpar(col = "blue",
                    cex = 1),
          y = unit(0.5, "lines")
        )
        grid.text(
          "5'",
          x = unit(2, "lines"),
          gp = gpar(col = "blue", cex = 1),
          y = unit(0.5,
                   "lines")
        )
        popViewport()
        pushViewport(viewport(y = 0.825, height = 0.25))
        motifStack::plotMotifLogo(
          p_full,
          paste(
            names(motif),
            "binding change by InDel",
            names(indel_info)
          ),
          font = "Helvetica-Bold",
          ncex = 1.2,
          yaxis = FALSE,
          xaxis = FALSE,
          xlab = "",
          ylab = "",
          newpage = FALSE,
          margins = c(2,
                      2.5, 1.5, 2.5)
        )
        grid.lines(
          x = c(
            convertUnit(unit(3, "lines"),
                        "npc", valueOnly = TRUE),
            1 - convertUnit(unit(2,
                                 "lines"), "npc", valueOnly = TRUE)
          ),
          y = unit(1, "lines"),
          gp = gpar(
            col = "blue",
            lwd = 1.5,
            xpd = NA
          ),
          arrow = arrow(
            length = unit(0.1,
                          "inches"),
            angle = 15,
            ends = "last"
          )
        )
        grid.text(
          "3'",
          x = unit(1, "npc") -
            unit(1, "lines"),
          gp = gpar(col = "blue",
                    cex = 1),
          y = unit(0.5, "lines")
        )
        grid.text(
          "5'",
          x = unit(2, "lines"),
          gp = gpar(col = "blue", cex = 1),
          y = unit(0.5,
                   "lines")
        )
        popViewport()
        pushViewport(viewport(y = 0.375, height = 0.25))
        grid.lines(
          x = unit(c((l1 + lm) / (nlong + 2),
                     (lm) / (nlong + 2)
          ), "npc"),
          y = unit(c(1.2,
                     0.75), "npc"),
          gp = gpar(fill = "black"),
          arrow = arrow(
            length = unit(0.05, "inches"),
            ends = "last",
            type = "closed"
          )
        )
        grid.lines(
          x = unit(c((l1 + lm) / (nlong + 2),
                     (lm + m) / (nlong + 2)
          ), "npc"),
          y = unit(c(1.2,
                     0.75), "npc"),
          gp = gpar(fill = "black"),
          arrow = arrow(
            length = unit(0.05, "inches"),
            ends = "last",
            type = "closed"
          )
        )
        grid.lines(
          x = unit(c((l1 + lm) / (nlong + 2),
                     (lm + l1) / (nlong + 2)
          ), "npc"),
          y = unit(c(2,
                     1.2), "npc"),
          gp = gpar(fill = "black")
        )
      }
      else {
        a <- makepwm(shorter, short_strand_sign)
        rownames(a) <- c("A", "C", "G",
                         "T")
        l1 <- floor(m / 2)
        r1 <- m - l1
        right <-
          (nlong - m) / 2 - (short_best_match + n1 + 1) + ceiling(m / 2)
        left <- nlong - lm - right
        if (l1 == 0) {
          a_full <- cbind(a, rep.col(c(0.25, 0.25, 0.25,
                                       0.25), r1))
        }
        else {
          a_full <- cbind(rep.col(c(0.25, 0.25, 0.25,
                                    0.25), l1), a, rep.col(c(0.25, 0.25, 0.25,
                                                             0.25), r1))
        }
        p <- makemotif(motif, short_strand_sign)
        p_full <- makepfull(p, left, right)
        title <-
          ifelse(
            insertion,
            "Best match to the reference sequence m=",
            "Best match to the mutation sequence m="
          )
        popViewport()
        pushViewport(viewport(y = 0.575, height = 0.25))
        motifStack::plotMotifLogo(
          a_full,
          xaxis = FALSE,
          yaxis = FALSE,
          xlab = "",
          ylab = "(-)",
          newpage = FALSE,
          margins = c(2, 2.5, 1.5, 2.5)
        )
        grid.text(
          "5'",
          x = unit(1, "npc") -
            unit(1, "lines"),
          gp = gpar(col = "blue",
                    cex = 1),
          y = unit(2.5, "lines")
        )
        grid.text(
          "3'",
          x = unit(2, "lines"),
          gp = gpar(col = "blue", cex = 1),
          y = unit(2.5,
                   "lines")
        )
        popViewport()
        pushViewport(viewport(y = 0.825, height = 0.25))
        motifStack::plotMotifLogo(
          p_full,
          paste(
            names(motif),
            "binding change by InDel",
            names(indel_info)
          ),
          font = "Helvetica-Bold",
          ncex = 1.2,
          yaxis = FALSE,
          xaxis = FALSE,
          xlab = "",
          ylab = "",
          newpage = FALSE,
          margins = c(2,
                      2.5, 1.5, 2.5)
        )
        grid.text(
          "5'",
          x = unit(1, "npc") -
            unit(1, "lines"),
          gp = gpar(col = "blue",
                    cex = 1),
          y = unit(2.5, "lines")
        )
        grid.text(
          "3'",
          x = unit(2, "lines"),
          gp = gpar(col = "blue", cex = 1),
          y = unit(2.5,
                   "lines")
        )
        grid.lines(
          x = c(
            convertUnit(unit(3, "lines"),
                        "npc", valueOnly = TRUE),
            1 - convertUnit(unit(2,
                                 "lines"), "npc", valueOnly = TRUE)
          ),
          y = unit(1, "lines"),
          gp = gpar(
            col = "blue",
            lwd = 1.5,
            xpd = NA
          ),
          arrow = arrow(
            length = unit(0.1,
                          "inches"),
            angle = 15,
            ends = "first"
          )
        )
        popViewport()
        pushViewport(viewport(y = 0.325, height = 0.25))
        grid.lines(
          x = unit(c((l1 + lm) / (nlong + 2),
                     (lm) / (nlong + 2)
          ), "npc"),
          y = unit(c(1.2,
                     0.75), "npc"),
          gp = gpar(fill = "black"),
          arrow = arrow(
            length = unit(0.05, "inches"),
            ends = "last",
            type = "closed"
          )
        )
        grid.lines(
          x = unit(c((l1 + lm) / (nlong + 2),
                     (lm + m) / (nlong + 2)
          ), "npc"),
          y = unit(c(1.2,
                     0.75), "npc"),
          gp = gpar(fill = "black"),
          arrow = arrow(
            length = unit(0.05, "inches"),
            ends = "last",
            type = "closed"
          )
        )
        grid.lines(
          x = unit(c((l1 + lm) / (nlong + 2),
                     (lm + l1) / (nlong + 2)
          ), "npc"),
          y = unit(c(2,
                     1.2), "npc"),
          gp = gpar(fill = "black")
        )
      }
    }
    else if (long_strand_sign == -1) {
      right <- (nlong - m) / 2 - (long_best_match + n1 +
                                    1)
      left <- nlong - lm - right
      a <- makepwm(longer, long_strand_sign)
      rownames(a) <- c("A", "C", "G",
                       "T")
      p <- makemotif(motif, long_strand_sign)
      p_full <- makepfull(p, left, right)
      title <-
        ifelse(
          insertion,
          "Best match to the mutation sequence m=",
          "Best match to the reference sequence m="
        )
      pushViewport(viewport(y = 0.125, height = 0.25))
      motifStack::plotMotifLogo(
        p_full,
        yaxis = FALSE,
        xaxis = FALSE,
        xlab = "",
        ylab = "",
        newpage = FALSE,
        margins = c(2, 2.5, 1.5, 2.5)
      )
      grid.text(
        "5'",
        x = unit(1, "npc") -
          unit(1, "lines"),
        gp = gpar(col = "blue",
                  cex = 1),
        y = unit(2.5, "lines")
      )
      grid.text(
        "3'",
        x = unit(2, "lines"),
        gp = gpar(col = "blue", cex = 1),
        y = unit(2.5,
                 "lines")
      )
      grid.lines(
        x = c(
          convertUnit(unit(3, "lines"),
                      "npc", valueOnly = TRUE),
          1 - convertUnit(unit(2,
                               "lines"), "npc", valueOnly = TRUE)
        ),
        y = unit(1, "lines"),
        gp = gpar(
          col = "blue",
          lwd = 1.5,
          xpd = NA
        ),
        arrow = arrow(
          length = unit(0.1,
                        "inches"),
          angle = 15,
          ends = "first"
        )
      )
      popViewport()
      pushViewport(viewport(y = 0.325, height = 0.25))
      motifStack::plotMotifLogo(
        a,
        xaxis = FALSE,
        yaxis = FALSE,
        xlab = "",
        ylab = "(-)",
        newpage = FALSE,
        margins = c(2, 2.5, 1.5, 2.5)
      )
      grid.text(
        "5'",
        x = unit(1, "npc") -
          unit(1, "lines"),
        gp = gpar(col = "blue",
                  cex = 1),
        y = unit(2.5, "lines")
      )
      grid.text(
        "3'",
        x = unit(2, "lines"),
        gp = gpar(col = "blue", cex = 1),
        y = unit(2.5,
                 "lines")
      )
      grid.rect(
        just = "center",
        width = (m) / (nlong +
                         2),
        height = unit(0.6, "npc"),
        gp = gpar(
          col = "blue",
          lty = 3,
          lwd = 2,
          fill = NA
        )
      )
      if (short_strand_sign == 1) {
        a <- makepwm(shorter, short_strand_sign)
        rownames(a) <- c("A", "C", "G",
                         "T")
        l1 <- floor(m / 2)
        r1 <- m - l1
        left <- short_best_match - s1 + l1
        right <- nlong - left - lm
        if (l1 == 0) {
          a_full <- cbind(a, rep.col(c(0.25, 0.25, 0.25,
                                       0.25), r1))
        }
        else {
          a_full <- cbind(rep.col(c(0.25, 0.25, 0.25,
                                    0.25), l1), a, rep.col(c(0.25, 0.25, 0.25,
                                                             0.25), r1))
        }
        p <- makemotif(motif, short_strand_sign)
        p_full <- makepfull(p, left, right)
        popViewport()
        pushViewport(viewport(y = 0.575, height = 0.25))
        motifStack::plotMotifLogo(
          a_full,
          xaxis = FALSE,
          yaxis = FALSE,
          xlab = "",
          ylab = "(+)",
          newpage = FALSE,
          margins = c(2, 2.5, 1.5, 2.5)
        )
        grid.text(
          "3'",
          x = unit(1, "npc") -
            unit(1, "lines"),
          gp = gpar(col = "blue",
                    cex = 1),
          y = unit(0.5, "lines")
        )
        grid.text(
          "5'",
          x = unit(2, "lines"),
          gp = gpar(col = "blue", cex = 1),
          y = unit(0.5,
                   "lines")
        )
        popViewport()
        pushViewport(viewport(y = 0.825, height = 0.25))
        motifStack::plotMotifLogo(
          p_full,
          paste(
            names(motif),
            "binding change by InDel",
            names(indel_info)
          ),
          font = "Helvetica-Bold",
          ncex = 1.2,
          yaxis = FALSE,
          xaxis = FALSE,
          xlab = "",
          ylab = "",
          newpage = FALSE,
          margins = c(2,
                      2.5, 1.5, 2.5)
        )
        grid.lines(
          x = c(
            convertUnit(unit(3, "lines"),
                        "npc", valueOnly = TRUE),
            1 - convertUnit(unit(2,
                                 "lines"), "npc", valueOnly = TRUE)
          ),
          y = unit(1, "lines"),
          gp = gpar(
            col = "blue",
            lwd = 1.5,
            xpd = NA
          ),
          arrow = arrow(
            length = unit(0.1,
                          "inches"),
            angle = 15,
            ends = "last"
          )
        )
        grid.text(
          "3'",
          x = unit(1, "npc") -
            unit(1, "lines"),
          gp = gpar(col = "blue",
                    cex = 1),
          y = unit(0.5, "lines")
        )
        grid.text(
          "5'",
          x = unit(2, "lines"),
          gp = gpar(col = "blue", cex = 1),
          y = unit(0.5,
                   "lines")
        )
        popViewport()
        pushViewport(viewport(y = 0.325, height = 0.25))
        grid.lines(
          x = unit(c((l1 + lm) / (nlong + 2),
                     (lm) / (nlong + 2)
          ), "npc"),
          y = unit(c(1.2,
                     0.75), "npc"),
          gp = gpar(fill = "black"),
          arrow = arrow(
            length = unit(0.05, "inches"),
            ends = "last",
            type = "closed"
          )
        )
        grid.lines(
          x = unit(c((l1 + lm) / (nlong + 2),
                     (lm + m) / (nlong + 2)
          ), "npc"),
          y = unit(c(1.2,
                     0.75), "npc"),
          gp = gpar(fill = "black"),
          arrow = arrow(
            length = unit(0.05, "inches"),
            ends = "last",
            type = "closed"
          )
        )
        grid.lines(
          x = unit(c((l1 + lm) / (nlong + 2),
                     (lm + l1) / (nlong + 2)
          ), "npc"),
          y = unit(c(2,
                     1.2), "npc"),
          gp = gpar(fill = "black")
        )
      }
      else {
        a <- makepwm(shorter, short_strand_sign)
        rownames(a) <- c("A", "C", "G",
                         "T")
        l1 <- floor(m / 2)
        r1 <- m - l1
        right <-
          (nlong - m) / 2 - (short_best_match + n1 + 1) + ceiling(m / 2)
        left <- nlong - lm - right
        if (l1 == 0) {
          a_full <- cbind(a, rep.col(c(0.25, 0.25, 0.25,
                                       0.25), r1))
        }
        else {
          a_full <- cbind(rep.col(c(0.25, 0.25, 0.25,
                                    0.25), l1), a, rep.col(c(0.25, 0.25, 0.25,
                                                             0.25), r1))
        }
        p <- makemotif(motif, short_strand_sign)
        p_full <- makepfull(p, left, right)
        title <-
          ifelse(
            insertion,
            "Best match to the reference sequence m=",
            "Best match to the mutation sequence m="
          )
        popViewport()
        pushViewport(viewport(y = 0.575, height = 0.25))
        motifStack::plotMotifLogo(
          a_full,
          xaxis = FALSE,
          yaxis = FALSE,
          xlab = "",
          ylab = "(-)",
          newpage = FALSE,
          margins = c(2, 2.5, 1.5, 2.5)
        )
        grid.text(
          "5'",
          x = unit(1, "npc") -
            unit(1, "lines"),
          gp = gpar(col = "blue",
                    cex = 1),
          y = unit(2.5, "lines")
        )
        grid.text(
          "3'",
          x = unit(2, "lines"),
          gp = gpar(col = "blue", cex = 1),
          y = unit(2.5,
                   "lines")
        )
        popViewport()
        pushViewport(viewport(y = 0.825, height = 0.25))
        motifStack::plotMotifLogo(
          p_full,
          paste(
            names(motif),
            "binding change by InDel",
            names(indel_info)
          ),
          font = "Helvetica-Bold",
          ncex = 1.2,
          yaxis = FALSE,
          xaxis = FALSE,
          xlab = "",
          ylab = "",
          newpage = FALSE,
          margins = c(2,
                      2.5, 1.5, 2.5)
        )
        grid.text(
          "5'",
          x = unit(1, "npc") -
            unit(1, "lines"),
          gp = gpar(col = "blue",
                    cex = 1),
          y = unit(2.5, "lines")
        )
        grid.text(
          "3'",
          x = unit(2, "lines"),
          gp = gpar(col = "blue", cex = 1),
          y = unit(2.5,
                   "lines")
        )
        grid.lines(
          x = c(
            convertUnit(unit(3, "lines"),
                        "npc", valueOnly = TRUE),
            1 - convertUnit(unit(2,
                                 "lines"), "npc", valueOnly = TRUE)
          ),
          y = unit(1, "lines"),
          gp = gpar(
            col = "blue",
            lwd = 1.5,
            xpd = NA
          ),
          arrow = arrow(
            length = unit(0.1,
                          "inches"),
            angle = 15,
            ends = "first"
          )
        )
        popViewport()
        pushViewport(viewport(y = 0.325, height = 0.25))
        grid.lines(
          x = unit(c((l1 + lm) / (nlong + 2),
                     (lm) / (nlong + 2)
          ), "npc"),
          y = unit(c(1.2,
                     0.75), "npc"),
          gp = gpar(fill = "black"),
          arrow = arrow(
            length = unit(0.05, "inches"),
            ends = "last",
            type = "closed"
          )
        )
        grid.lines(
          x = unit(c((l1 + lm) / (nlong + 2),
                     (lm + m) / (nlong + 2)
          ), "npc"),
          y = unit(c(1.2,
                     0.75), "npc"),
          gp = gpar(fill = "black"),
          arrow = arrow(
            length = unit(0.05, "inches"),
            ends = "last",
            type = "closed"
          )
        )
        grid.lines(
          x = unit(c((l1 + lm) / (nlong + 2),
                     (lm + l1) / (nlong + 2)
          ), "npc"),
          y = unit(c(2,
                     1.2), "npc"),
          gp = gpar(fill = "black")
        )
      }
    }
  }
  else if (insertion == 0) {
    if (long_strand_sign == 1) {
      left <- long_best_match - s1
      right <- nlong - left - lm
      a <- makepwm(longer, long_strand_sign)
      rownames(a) <- c("A", "C", "G",
                       "T")
      p <- makemotif(motif, long_strand_sign)
      p_full <- makepfull(p, left, right)
      pushViewport(viewport(y = 0.825, height = 0.25))
      motifStack::plotMotifLogo(
        p_full,
        paste(
          names(motif),
          "binding change by InDel",
          names(indel_info)
        ),
        font = "Helvetica-Bold",
        ncex = 1.2,
        yaxis = FALSE,
        xaxis = FALSE,
        xlab = "",
        ylab = "",
        newpage = FALSE,
        margins = c(2, 2.5, 1.5, 2.5)
      )
      grid.lines(
        x = c(
          convertUnit(unit(3, "lines"),
                      "npc", valueOnly = TRUE),
          1 - convertUnit(unit(2,
                               "lines"), "npc", valueOnly = TRUE)
        ),
        y = unit(1, "lines"),
        gp = gpar(
          col = "blue",
          lwd = 1.5,
          xpd = NA
        ),
        arrow = arrow(
          length = unit(0.1,
                        "inches"),
          angle = 15,
          ends = "last"
        )
      )
      grid.text(
        "3'",
        x = unit(1, "npc") -
          unit(1, "lines"),
        gp = gpar(col = "blue",
                  cex = 1),
        y = unit(0.5, "lines")
      )
      grid.text(
        "5'",
        x = unit(2, "lines"),
        gp = gpar(col = "blue", cex = 1),
        y = unit(0.5,
                 "lines")
      )
      popViewport()
      pushViewport(viewport(y = 0.575, height = 0.25))
      motifStack::plotMotifLogo(
        a,
        xaxis = FALSE,
        yaxis = FALSE,
        xlab = "",
        ylab = "(+)",
        newpage = FALSE,
        margins = c(2, 2, 1.5, 2)
      )
      grid.text(
        "3'",
        x = unit(1, "npc") -
          unit(1, "lines"),
        gp = gpar(col = "blue",
                  cex = 1),
        y = unit(0.5, "lines")
      )
      grid.text(
        "5'",
        x = unit(2, "lines"),
        gp = gpar(col = "blue", cex = 1),
        y = unit(0.5,
                 "lines")
      )
      grid.rect(
        just = "center",
        width = (m) / (nlong +
                         2),
        height = unit(0.6, "npc"),
        gp = gpar(
          col = "blue",
          lty = 3,
          lwd = 2,
          fill = NA
        )
      )
      if (short_strand_sign == 1) {
        left <- short_best_match - s1
        right <- nlong - left - lm
        a <- makepwm(shorter, short_strand_sign)
        rownames(a) <- c("A", "C", "G",
                         "T")
        a_full <- cbind(a[, 1:(lm - 1)], rep.col(c(0.25,
                                                   0.25, 0.25, 0.25), m), a[, lm:(2 * lm - 2)])
        p <- makemotif(motif, short_strand_sign)
        p_full <- makepfull(p, left, right)
        popViewport()
        pushViewport(viewport(y = 0.325, height = 0.25))
        motifStack::plotMotifLogo(
          a_full,
          xaxis = FALSE,
          yaxis = FALSE,
          xlab = "",
          ylab = "(+)",
          newpage = FALSE,
          margins = c(2, 2.5, 1.5, 2.5)
        )
        grid.text(
          "3'",
          x = unit(1, "npc") -
            unit(1, "lines"),
          gp = gpar(col = "blue",
                    cex = 1),
          y = unit(0.5, "lines")
        )
        grid.text(
          "5'",
          x = unit(2, "lines"),
          gp = gpar(col = "blue", cex = 1),
          y = unit(0.5,
                   "lines")
        )
        popViewport()
        pushViewport(viewport(y = 0.125, height = 0.25))
        motifStack::plotMotifLogo(
          p_full,
          yaxis = FALSE,
          xaxis = FALSE,
          xlab = "",
          ylab = "",
          newpage = FALSE,
          margins = c(2, 2.5, 1.5, 2.5)
        )
        grid.lines(
          x = c(
            convertUnit(unit(3, "lines"),
                        "npc", valueOnly = TRUE),
            1 - convertUnit(unit(2,
                                 "lines"), "npc", valueOnly = TRUE)
          ),
          y = unit(1, "lines"),
          gp = gpar(
            col = "blue",
            lwd = 1.5,
            xpd = NA
          ),
          arrow = arrow(
            length = unit(0.1,
                          "inches"),
            angle = 15,
            ends = "last"
          )
        )
        grid.text(
          "3'",
          x = unit(1, "npc") -
            unit(1, "lines"),
          gp = gpar(col = "blue",
                    cex = 1),
          y = unit(0.5, "lines")
        )
        grid.text(
          "5'",
          x = unit(2, "lines"),
          gp = gpar(col = "blue", cex = 1),
          y = unit(0.5,
                   "lines")
        )
        popViewport()
        pushViewport(viewport(y = 0.125, height = 0.25))
        grid.lines(
          x = unit(c(lm / (nlong + 2), lm / (nlong +
                                               2)), "npc"),
          y = unit(c(1.2, 0.75), "npc"),
          gp = gpar(fill = "black"),
          arrow = arrow(
            length = unit(0.05,
                          "inches"),
            ends = "last",
            type = "closed"
          )
        )
        grid.lines(
          x = unit(c((lm + m) / (nlong + 2), lm / (nlong +
                                                     2)), "npc"),
          y = unit(c(1.2, 0.75), "npc"),
          gp = gpar(fill = "black"),
          arrow = arrow(
            length = unit(0.05,
                          "inches"),
            ends = "last",
            type = "closed"
          )
        )
      }
      else {
        right <- (nlong - m) / 2 - (short_best_match +
                                      n1 + 1) + m
        left <- nlong - lm - right
        a <- makepwm(shorter, short_strand_sign)
        rownames(a) <- c("A", "C", "G",
                         "T")
        a_full <- cbind(a[, 1:(lm - 1)], rep.col(c(0.25,
                                                   0.25, 0.25, 0.25), m), a[, lm:(2 * lm - 2)])
        p <- makemotif(motif, short_strand_sign)
        p_full <- makepfull(p, left, right)
        title <-
          ifelse(
            insertion,
            "Best match to the reference sequence m=",
            "Best match to the mutation sequence m="
          )
        popViewport()
        pushViewport(viewport(y = 0.325, height = 0.25))
        motifStack::plotMotifLogo(
          a_full,
          xaxis = FALSE,
          yaxis = FALSE,
          xlab = "",
          ylab = "(-)",
          newpage = FALSE,
          margins = c(2, 2.5, 1.5, 2.5)
        )
        grid.text(
          "5'",
          x = unit(1, "npc") -
            unit(1, "lines"),
          gp = gpar(col = "blue",
                    cex = 1),
          y = unit(2.5, "lines")
        )
        grid.text(
          "3'",
          x = unit(2, "lines"),
          gp = gpar(col = "blue", cex = 1),
          y = unit(2.5,
                   "lines")
        )
        popViewport()
        pushViewport(viewport(y = 0.125, height = 0.25))
        motifStack::plotMotifLogo(
          p_full,
          yaxis = FALSE,
          xaxis = FALSE,
          xlab = "",
          ylab = "",
          newpage = FALSE,
          margins = c(2, 2.5, 1.5, 2.5)
        )
        grid.text(
          "5'",
          x = unit(1, "npc") -
            unit(1, "lines"),
          gp = gpar(col = "blue",
                    cex = 1),
          y = unit(2.5, "lines")
        )
        grid.text(
          "3'",
          x = unit(2, "lines"),
          gp = gpar(col = "blue", cex = 1),
          y = unit(2.5,
                   "lines")
        )
        grid.lines(
          x = c(
            convertUnit(unit(3, "lines"),
                        "npc", valueOnly = TRUE),
            1 - convertUnit(unit(2,
                                 "lines"), "npc", valueOnly = TRUE)
          ),
          y = unit(1, "lines"),
          gp = gpar(
            col = "blue",
            lwd = 1.5,
            xpd = NA
          ),
          arrow = arrow(
            length = unit(0.1,
                          "inches"),
            angle = 15,
            ends = "first"
          )
        )
        popViewport()
        pushViewport(viewport(y = 0.125, height = 0.25))
        grid.lines(
          x = unit(c(lm / (nlong + 2), lm / (nlong +
                                               2)), "npc"),
          y = unit(c(1.2, 0.75), "npc"),
          gp = gpar(fill = "black"),
          arrow = arrow(
            length = unit(0.05,
                          "inches"),
            ends = "last",
            type = "closed"
          )
        )
        grid.lines(
          x = unit(c((lm + m) / (nlong + 2), lm / (nlong +
                                                     2)), "npc"),
          y = unit(c(1.2, 0.75), "npc"),
          gp = gpar(fill = "black"),
          arrow = arrow(
            length = unit(0.05,
                          "inches"),
            ends = "last",
            type = "closed"
          )
        )
      }
    }
    else if (long_strand_sign == -1) {
      right <- (nlong - m) / 2 - (long_best_match + n1 +
                                    1)
      left <- nlong - lm - right
      a <- makepwm(longer, long_strand_sign)
      rownames(a) <- c("A", "C", "G",
                       "T")
      p <- makemotif(motif, long_strand_sign)
      p_full <- makepfull(p, left, right)
      title <-
        ifelse(
          insertion,
          "Best match to the mutation sequence m=",
          "Best match to the reference sequence m="
        )
      pushViewport(viewport(y = 0.825, height = 0.25))
      motifStack::plotMotifLogo(
        p_full,
        paste(
          names(motif),
          "binding change by InDel",
          names(indel_info)
        ),
        font = "Helvetica-Bold",
        ncex = 1.2,
        yaxis = FALSE,
        xaxis = FALSE,
        xlab = "",
        ylab = "",
        newpage = FALSE,
        margins = c(2, 2.5, 1.5, 2.5)
      )
      grid.text(
        "5'",
        x = unit(1, "npc") -
          unit(1, "lines"),
        gp = gpar(col = "blue",
                  cex = 1),
        y = unit(2.5, "lines")
      )
      grid.text(
        "3'",
        x = unit(2, "lines"),
        gp = gpar(col = "blue", cex = 1),
        y = unit(2.5,
                 "lines")
      )
      grid.lines(
        x = c(
          convertUnit(unit(3, "lines"),
                      "npc", valueOnly = TRUE),
          1 - convertUnit(unit(2,
                               "lines"), "npc", valueOnly = TRUE)
        ),
        y = unit(1, "lines"),
        gp = gpar(
          col = "blue",
          lwd = 1.5,
          xpd = NA
        ),
        arrow = arrow(
          length = unit(0.1,
                        "inches"),
          angle = 15,
          ends = "first"
        )
      )
      popViewport()
      pushViewport(viewport(y = 0.575, height = 0.25))
      motifStack::plotMotifLogo(
        a,
        font = "mono,Courier",
        xaxis = FALSE,
        yaxis = FALSE,
        xlab = "",
        ylab = "(-)",
        newpage = FALSE,
        margins = c(2,
                    2.5, 1.5, 2.5)
      )
      grid.text(
        "5'",
        x = unit(1, "npc") -
          unit(1, "lines"),
        gp = gpar(col = "blue",
                  cex = 1),
        y = unit(2.5, "lines")
      )
      grid.text(
        "3'",
        x = unit(2, "lines"),
        gp = gpar(col = "blue", cex = 1),
        y = unit(2.5,
                 "lines")
      )
      grid.rect(
        just = "center",
        width = (m) / (nlong +
                         2),
        height = unit(0.6, "npc"),
        gp = gpar(
          col = "blue",
          lty = 3,
          lwd = 2,
          fill = NA
        )
      )
      if (short_strand_sign == 1) {
        left <- short_best_match - s1
        right <- nlong - left - lm
        a <- makepwm(shorter, short_strand_sign)
        rownames(a) <- c("A", "C", "G",
                         "T")
        a_full <- cbind(a[, 1:(lm - 1)], rep.col(c(0.25,
                                                   0.25, 0.25, 0.25), m), a[, lm:(2 * lm - 2)])
        p <- makemotif(motif, short_strand_sign)
        p_full <- cbind(rep.col(c(0.25, 0.25, 0.25, 0.25),
                                left), p, rep.col(c(0.25, 0.25, 0.25, 0.25),
                                                  right))
        title <-
          ifelse(
            insertion,
            "Best match to the reference sequence m=",
            "Best match to the mutation sequence m="
          )
        popViewport()
        pushViewport(viewport(y = 0.325, height = 0.25))
        motifStack::plotMotifLogo(
          a_full,
          font = "mono,Courier",
          xaxis = FALSE,
          yaxis = FALSE,
          xlab = "",
          ylab = "(+)",
          newpage = FALSE,
          margins = c(2,
                      2.5, 1.5, 2.5)
        )
        grid.text(
          "3'",
          x = unit(1, "npc") -
            unit(1, "lines"),
          gp = gpar(col = "blue",
                    cex = 1),
          y = unit(0.5, "lines")
        )
        grid.text(
          "5'",
          x = unit(2, "lines"),
          gp = gpar(col = "blue", cex = 1),
          y = unit(0.5,
                   "lines")
        )
        popViewport()
        pushViewport(viewport(y = 0.125, height = 0.25))
        motifStack::plotMotifLogo(
          p_full,
          yaxis = FALSE,
          xaxis = FALSE,
          xlab = "",
          ylab = "",
          newpage = FALSE,
          margins = c(2, 2.5, 1.5, 2.5)
        )
        grid.lines(
          x = c(
            convertUnit(unit(3, "lines"),
                        "npc", valueOnly = TRUE),
            1 - convertUnit(unit(2,
                                 "lines"), "npc", valueOnly = TRUE)
          ),
          y = unit(1, "lines"),
          gp = gpar(
            col = "blue",
            lwd = 1.5,
            xpd = NA
          ),
          arrow = arrow(
            length = unit(0.1,
                          "inches"),
            angle = 15,
            ends = "last"
          )
        )
        grid.text(
          "3'",
          x = unit(1, "npc") -
            unit(1, "lines"),
          gp = gpar(col = "blue",
                    cex = 1),
          y = unit(0.5, "lines")
        )
        grid.text(
          "5'",
          x = unit(2, "lines"),
          gp = gpar(col = "blue", cex = 1),
          y = unit(0.5,
                   "lines")
        )
        popViewport()
        pushViewport(viewport(y = 0.125, height = 0.25))
        grid.lines(
          x = unit(c(lm / (nlong + 2), lm / (nlong +
                                               2)), "npc"),
          y = unit(c(1.2, 0.75), "npc"),
          gp = gpar(fill = "black"),
          arrow = arrow(
            length = unit(0.05,
                          "inches"),
            ends = "last",
            type = "closed"
          )
        )
        grid.lines(
          x = unit(c((lm + m) / (nlong + 2), lm / (nlong +
                                                     2)), "npc"),
          y = unit(c(1.2, 0.75), "npc"),
          gp = gpar(fill = "black"),
          arrow = arrow(
            length = unit(0.05,
                          "inches"),
            ends = "last",
            type = "closed"
          )
        )
      }
      else {
        right <- (nlong - m) / 2 - (short_best_match +
                                      n1 + 1) + m
        left <- nlong - lm - right
        a <- makepwm(shorter, short_strand_sign)
        rownames(a) <- c("A", "C", "G",
                         "T")
        a_full <- cbind(a[, 1:(lm - 1)], rep.col(c(0.25,
                                                   0.25, 0.25, 0.25), m), a[, lm:(2 * lm - 2)])
        p <- makemotif(motif, short_strand_sign)
        p_full <- makepfull(p, left, right)
        title <-
          ifelse(
            insertion,
            "Best match to the reference sequence m=",
            "Best match to the mutation sequence m="
          )
        popViewport()
        pushViewport(viewport(y = 0.375, height = 0.25))
        motifStack::plotMotifLogo(
          a_full,
          xaxis = FALSE,
          yaxis = FALSE,
          xlab = "",
          ylab = "(-)",
          newpage = FALSE,
          margins = c(2, 2.5, 1.5, 2.5)
        )
        grid.text(
          "5'",
          x = unit(1, "npc") -
            unit(1, "lines"),
          gp = gpar(col = "blue",
                    cex = 1),
          y = unit(2.5, "lines")
        )
        grid.text(
          "3'",
          x = unit(2, "lines"),
          gp = gpar(col = "blue", cex = 1),
          y = unit(2.5,
                   "lines")
        )
        popViewport()
        pushViewport(viewport(y = 0.125, height = 0.25))
        motifStack::plotMotifLogo(
          p_full,
          yaxis = FALSE,
          xaxis = FALSE,
          xlab = "",
          ylab = "",
          newpage = FALSE,
          margins = c(2, 2.5, 1.5, 2.5)
        )
        grid.text(
          "5'",
          x = unit(1, "npc") -
            unit(1, "lines"),
          gp = gpar(col = "blue",
                    cex = 1),
          y = unit(2.5, "lines")
        )
        grid.text(
          "3'",
          x = unit(2, "lines"),
          gp = gpar(col = "blue", cex = 1),
          y = unit(2.5,
                   "lines")
        )
        grid.lines(
          x = c(
            convertUnit(unit(3, "lines"),
                        "npc", valueOnly = TRUE),
            1 - convertUnit(unit(2,
                                 "lines"), "npc", valueOnly = TRUE)
          ),
          y = unit(1, "lines"),
          gp = gpar(
            col = "blue",
            lwd = 1.5,
            xpd = NA
          ),
          arrow = arrow(
            length = unit(0.1,
                          "inches"),
            angle = 15,
            ends = "first"
          )
        )
        popViewport()
        pushViewport(viewport(y = 0.125, height = 0.25))
        grid.lines(
          x = unit(c(lm / (nlong + 2), lm / (nlong +
                                               2)), "npc"),
          y = unit(c(1.2, 0.75), "npc"),
          gp = gpar(fill = "black"),
          arrow = arrow(
            length = unit(0.05,
                          "inches"),
            ends = "last",
            type = "closed"
          )
        )
        grid.lines(
          x = unit(c((lm + m) / (nlong + 2), lm / (nlong +
                                                     2)), "npc"),
          y = unit(c(1.2, 0.75), "npc"),
          gp = gpar(fill = "black"),
          arrow = arrow(
            length = unit(0.05,
                          "inches"),
            ends = "last",
            type = "closed"
          )
        )
      }
    }
  }
}
chandlerzuo/atIndel documentation built on Jan. 20, 2024, 4:10 a.m.