version 0.0.2

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: JunJunZai
 Title: A Package to Visualize Homer Output Results
-Version: 0.0.1
+Version: 0.0.2
 Author: Jun Zhang
 Maintainer: Jun Zhang <[email protected]>
 BugReports: https://github.com/junjunlab/JunJunZai/issues
@@ -10,16 +10,21 @@ Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.2.3
 Imports: 
+    ChIPpeakAnno,
     ComplexHeatmap,
     dplyr,
     ggplot2,
     ggseqlogo,
     grid,
     gtable,
+    IRanges,
     magrittr,
     methods,
     monaLisa,
+    rtracklayer,
     seqLogo,
     stats,
-    TFBSTools
+    stringr,
+    TFBSTools,
+    tidyr
 biocViews: TFBSTools,monaLisa

NAMESPACE

@@ -2,14 +2,18 @@
 
 export("%>%")
 export(facet_sub)
+export(footprintPlot)
+export(getSeqFrombed)
 export(gglogo)
 export(homerResult)
 export(loadHomerRes)
 export(parseHomer)
 export(parseHomerMotif)
 export(plotMotifHeatmap)
 export(preparePFMmat)
+export(seq2pwmMatrix)
 export(seqLogoGrob2)
+export(zoom_motif_grob)
 exportClasses(homerResult)
 exportMethods(gglogo)
 exportMethods(plotMotifHeatmap)
@@ -20,12 +24,18 @@ import(ggseqlogo)
 import(grid)
 import(methods)
 import(utils)
+importFrom(ChIPpeakAnno,getAllPeakSequence)
+importFrom(IRanges,resize)
 importFrom(TFBSTools,PFMatrixList)
 importFrom(TFBSTools,PWMatrix)
 importFrom(dplyr,filter)
 importFrom(dplyr,if_else)
 importFrom(dplyr,mutate)
 importFrom(dplyr,select)
 importFrom(magrittr,"%>%")
+importFrom(rtracklayer,import.bed)
 importFrom(stats,as.dist)
 importFrom(stats,hclust)
+importFrom(stringr,str_detect)
+importFrom(stringr,str_replace)
+importFrom(tidyr,pivot_longer)

R/footprint_motif_vis.R

@@ -0,0 +1,141 @@
+globalVariables(c("Occupancy", "TF", "group", "name", "pos"))
+
+makePFM <- getFromNamespace("makePFM","ggseqlogo")
+
+# ==============================================================================
+# footprint plot
+# ==============================================================================
+
+#' plot footprint line plot from rgt-hint output
+#'
+#' @param linePath the rgt-hint output directory, default NULL.
+#' @param TFname the transcription factor name to be plotted, default NULL.
+#'
+#' @return a ggplot
+#' @export
+#' @importFrom stringr str_detect str_replace
+#' @importFrom tidyr pivot_longer
+#' @import ggplot2
+footprintPlot <- function(linePath = NULL,TFname = NULL){
+  fp_file <- list.files(path = linePath,pattern = "txt",full.names = T)
+  fp_file_name <- list.files(path = linePath,pattern = "txt")
+  fp_file_name <- sapply(strsplit(fp_file_name,split = ".txt"), "[",1)
+
+  # loop load footprinting
+  # x = 719
+  lapply(seq_along(fp_file),function(x){
+    # check names
+
+    var <- stringr::str_detect(fp_file_name[x],pattern = "var")
+
+    if(var){
+      tmp_name <- stringr::str_replace(fp_file_name[x],pattern = "_",replacement = "(")
+      tmp_name <- paste0(tmp_name,")")
+    }else{
+      tmp_name <- fp_file_name[x]
+    }
+
+    # load data
+    ft_tmp <- read.delim(fp_file[x]) %>%
+      mutate(pos = -100:99,TF = tmp_name)
+    ft <- ft_tmp %>%
+      tidyr::pivot_longer(cols = colnames(ft_tmp)[1:(ncol(ft_tmp) - 2)],
+                          names_to = "group",
+                          values_to = "Occupancy")
+
+  }) %>% do.call("rbind",.) -> ft_df
+
+  # footprinting line plot
+  footprint <-
+    ggplot(ft_df %>%
+             dplyr::filter(TF == TFname),
+           aes(x = pos,y = Occupancy,color = group)) +
+    geom_line() +
+    theme_bw() +
+    theme(panel.grid = element_blank(),
+          strip.background = element_rect(fill = "grey90"),
+          legend.background = element_blank(),
+          strip.text = element_text(face = "bold.italic",size = rel(1)),
+          axis.text = element_text(colour = "black")) +
+    facet_wrap(~TF)
+
+  footprint
+}
+
+
+# ==============================================================================
+# get TF sequence
+# ==============================================================================
+
+
+#' extract sequence from bed
+#'
+#' @param mbpsPath the mbps file from rgt-motifanalysis output directory, default NULL.
+#' @param TFname the transcription factor name to be plotted, default NULL.
+#' @param genome the BSgenome object for genome sequence, default NULL.
+#'
+#' @return sequence
+#' @export
+#' @importFrom rtracklayer import.bed
+#' @importFrom IRanges resize
+#' @importFrom ChIPpeakAnno getAllPeakSequence
+getSeqFrombed <- function(mbpsPath = NULL,TFname = NULL,genome = NULL){
+  mbps <- list.files(path = mbpsPath,pattern = "mpbs.bed",full.names = T)
+  lapply(mbps,function(x){
+    rtracklayer::import.bed(x) %>% subset(.,name == TFname)
+  }) %>% do.call("c",.) -> TFtarget
+
+  TFtarget_ed <- IRanges::resize(TFtarget,width = 200,fix = "center")
+
+  # subtract sequence for TF
+  mtf_seq <- ChIPpeakAnno::getAllPeakSequence(myPeakList = TFtarget_ed,
+                                              upstream = 0,downstream = 0,
+                                              genome = genome)
+
+  return(mtf_seq$sequence)
+}
+
+
+
+# ==============================================================================
+# sequence to pwmMatrix
+# ==============================================================================
+
+
+#' create a pwmMatrix
+#'
+#' @param seqs sequence to make pwmMatrix, default NULL.
+#' @param shift a number of shift to zoom the motif, default c(-10,10).
+#'
+#' @return A list pwmMatrix
+#' @export
+#' @importFrom TFBSTools PWMatrix
+seq2pwmMatrix <- function(seqs = NULL,shift = c(-10,10)){
+  pfm <- makePFM(seqs = seqs,seq_type = "dna")
+
+  pwm_mat <- TFBSTools::PWMatrix(ID = "Unknown", name = "Unknown", matrixClass = "Unknown",
+                                 strand = "+",
+                                 bg = c(A=0.25, C=0.25, G=0.25, T=0.25),
+                                 tags = list(),
+                                 profileMatrix = matrix(pfm[c(1,4,3,2),],
+                                                        byrow = FALSE, nrow = 4,
+                                                        dimnames=list(c("A", "C", "G", "T"))),
+                                 pseudocounts = numeric())
+
+  # get zoom range
+  mid <- ncol(pfm)/2
+  zoom_rg <- c((mid + shift[1]):(mid + shift[2]))
+
+  pwm_mat_sub <- TFBSTools::PWMatrix(ID = "Unknown", name = "Unknown", matrixClass = "Unknown",
+                                     strand = "+",
+                                     bg = c(A=0.25, C=0.25, G=0.25, T=0.25),
+                                     tags = list(),
+                                     profileMatrix = matrix(pfm[c(1,4,3,2),zoom_rg],
+                                                            byrow = FALSE, nrow = 4,
+                                                            dimnames=list(c("A", "C", "G", "T"))),
+                                     pseudocounts = numeric())
+
+  return(list(pwm_mat = pwm_mat,
+              pwm_mat_sub = pwm_mat_sub,
+              shift = shift))
+}

R/zoom_motif_grob.R

@@ -0,0 +1,109 @@
+# ==============================================================================
+# construct grob
+# ==============================================================================
+#' construct zoom_motif_grob
+#'
+#' @param pwmMatrix a pwmMatrix from seq2pwmMatrix function, default NULL.
+#' @param rel_height the relative height for each part, default c(0.4,0.2,0.4).
+#' @param link_gp_col the color of link, default "white".
+#' @param link_gp_fill the fill of link, default "grey90".
+#' @param border_col the border color, default "white".
+#' @param name name
+#' @param gp gp
+#' @param vp vp
+#'
+#' @return a grob
+#' @export
+#' @import grid
+zoom_motif_grob <- function(pwmMatrix = NULL,
+                            rel_height = c(0.4,0.2,0.4),
+                            link_gp_col = "white",
+                            link_gp_fill = "grey90",
+                            border_col = "white",
+                            name = NULL,gp = NULL,vp = NULL){
+  gTree(pwmMatrix = pwmMatrix,
+        rel_height = rel_height,
+        link_gp_col = link_gp_col,
+        link_gp_fill = link_gp_fill,
+        border_col = border_col,
+        name = name,gp = gp,vp = vp,
+        cl = "zoom_motif_grob")
+}
+
+
+
+makeContent.zoom_motif_grob <- function(x){
+  # get relative panel y position
+  rel_y = cumsum(x$rel_height)
+
+  # get zoom range
+  mid <- ncol(x$pwmMatrix$pwm_mat)/2
+  n_col <- ncol(x$pwmMatrix$pwm_mat)
+  zoom_rg <- c((mid + x$pwmMatrix$shift[1]):(mid + x$pwmMatrix$shift[2]))
+
+  # ============================================================================
+  # first layer
+  # ============================================================================
+  # grid.newpage()
+  # pushViewport(viewport(y = unit(rel_y[1],"npc"),height = unit(rel_height[1],"npc"),
+  #                       just = "top",
+  #                       xscale = c(1,ncol(pfm)))
+  # )
+  vp1 <- viewport(y = unit(rel_y[3],"npc"),height = unit(x$rel_height[1],"npc"),
+                  just = "top",
+                  xscale = c(1,n_col))
+
+  zoom_rect1 <-
+    rectGrob(x = 100,
+             width = abs(range(zoom_rg)[2] - range(zoom_rg)[1]),
+             gp = gpar(fill = x$link_gp_fill,col = x$link_gp_col),
+             default.units = "native",
+             vp = vp1)
+  # grid::grid.draw(JunJunZai::seqLogoGrob2(x = pwmMatrix$pwm_mat))
+  logo1 <- JunJunZai::seqLogoGrob2(x = x$pwmMatrix$pwm_mat,vp = vp1)
+  panel_rect1 <- rectGrob(gp = gpar(fill = "transparent",col = x$border_col),vp = vp1)
+  # popViewport()
+
+  # ============================================================================
+  # second layer
+  # ============================================================================
+  # pushViewport(viewport(y = unit(rel_y[2],"npc"),height = unit(rel_height[2],"npc"),
+  #                       just = "top",
+  #                       xscale = c(1,ncol(pfm)))
+  # )
+  vp2 <- viewport(y = unit(rel_y[2],"npc"),height = unit(x$rel_height[2],"npc"),
+                  just = "top",
+                  xscale = c(1,n_col))
+
+  zoom_rect2 <-
+    polygonGrob(x = c(1,range(zoom_rg),n_col),
+                y = c(0,1,1,0),
+                gp = gpar(fill = x$link_gp_fill,col = x$link_gp_col),
+                default.units = "native",
+                vp = vp2)
+  panel_rect2 <- rectGrob(gp = gpar(fill = "transparent",col = x$border_col),vp = vp2)
+  # popViewport()
+
+  # ============================================================================
+  # third layer
+  # ============================================================================
+  # pushViewport(viewport(y = unit(rel_y[3],"npc"),height = unit(rel_height[3],"npc"),
+  #                       just = "top",
+  #                       xscale = c(1,ncol(pfm)))
+  # )
+  vp3 <- viewport(y = unit(rel_y[1],"npc"),height = unit(x$rel_height[3],"npc"),
+                  just = "top",
+                  xscale = c(1,n_col))
+
+  # grid::grid.draw(JunJunZai::seqLogoGrob2(x = pwmMatrix$pwm_mat_sub))
+  logo2 <- JunJunZai::seqLogoGrob2(x = x$pwmMatrix$pwm_mat_sub,vp = vp3)
+  panel_rect3 <- rectGrob(gp = gpar(fill = "transparent",col = x$border_col),vp = vp3)
+  # popViewport()
+
+
+  # return
+  setChildren(x,gList(zoom_rect1,logo1,panel_rect1,
+                      zoom_rect2,panel_rect2,
+                      logo2,panel_rect3))
+
+}