Merge branches of RepliSeq package

.Rbuildignore

@@ -0,0 +1,3 @@
+^.*\.Rproj$
+^\.Rproj\.user$
+^LICENSE

DESCRIPTION

@@ -1,18 +1,42 @@
 Package: RepliSeq
 Type: Package
-Title: Repliseq
-Version: 0.0.0.9000
-Author: Sami EL HILALI
+Title: Process Repli-Seq data to analyze the DNA replication timing
+Version: 0.1.1
+Authors@R: c(person(given="Sami",
+ family ="El Hilali",
+ email = "elhilali.sami@gmail.com",
+ role = c("aut", "cre"),
+ comment = c(ORCID = "0000-0003-4417-8399")),
+ person(given="Chun-Long",
+ family="Chen",
+ role = "aut",
+ comment = c(ORCID = "0000-0002-4795-0295"),
+ person(given="Stefano",
+ family="Gnan",
+ role = "ctb",
+ comment = c(ORCID = "0000-0003-3314-4302"))
+  )
 Maintainer: Sami EL HILALI <elhilali.sami@gmail.com>
-Description: This Package contains a set of functions for the analysis of Repli-Seq data to study DNA replication timing program. These functions handles repliseq assays as data frames from count matrices (bedgraph format !!) and performs such things as rescaling, data normalization, replication timing calculation (S50), URI calculation.
-License: GPL-3
+Depends: R (>= 4.1.0)
+Imports: 
+    dplyr,
+    magrittr,
+    RColorBrewer,
+    ggplot2,
+    corrplot,
+    reshape2
+Description: This Package contains a set of functions for the analysis of Repli-Seq data to 
+ study DNA replication timing program. These functions handles repliseq assays as data frames 
+ from count matrices (Bedgraph format) and performs rescaling, data normalization, 
+ replication timing calculation (S50), URI calculation.
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 7.1.0
+License: GPL-3
+URL: https://github.com/SamiLhll/RepliSeq
+BugReports: https://github.com/SamiLhll/RepliSeq/issues
+Roxygen: list(markdown = TRUE)
+RoxygenNote: 7.2.1
 Suggests: 
     knitr,
     rmarkdown
 VignetteBuilder: knitr
-Imports: 
-    dplyr (>= 0.8.3),
-    magrittr (>= 1.5)

NAMESPACE

@@ -7,13 +7,25 @@ export(doubleXchr)
 export(normalizeRS)
 export(readRS)
 export(removeNoise)
+export(repliseqCorrPlot)
+export(repliseqDensityPlot)
 export(rescaleRS)
 export(smoothRS)
 export(weight_phases)
 export(writeBedgraph)
 export(writeBigwig)
+import(ggplot2)
 import(magrittr)
+importFrom(RColorBrewer,brewer.pal)
+importFrom(corrplot,corrplot)
+importFrom(corrplot,corrplot.mixed)
 importFrom(dplyr,funs)
 importFrom(dplyr,mutate_each)
 importFrom(dplyr,rename)
 importFrom(dplyr,select)
+importFrom(reshape2,melt)
+importFrom(stats,cor)
+importFrom(stats,cor.test)
+importFrom(stats,density)
+importFrom(utils,read.table)
+importFrom(utils,write.table)

NEWS.md

@@ -0,0 +1,17 @@
+# RepliSeq 0.1.1
+
+### New features :
+
+* Added the function repliseqCorrPlot() to draw correlograms of several Repli-Seq assays.   
+
+### Enhancements :
+
+* Added a `NEWS.md` file to track changes to the package.
+
+* Improved the documentation and vignette   
+
+* repliseqDensityPlot now allows to draw counts through the argument stat, and it now has arguments to display the labels.
+
+### Bug fixes :
+
+* Fixed a bug in rescaleRS()

R/calculateNoiseRatios.R

@@ -3,23 +3,32 @@
 #'
 #' @param rs_assay a Repli-seq assay (data.frame) loaded with readRS() and formatted as chr,start,stop,S1,...,Sx
 #' @param rs_control a Repli-seq control assay (data.frame) loaded with readRS() and formatted as chr,start,stop,S0
-#'
+#' @param noise_method either "max_density"(default) or "min_downslope"
 #' @return an array with the ratios in the order of rs_assay fractions
+#' 
+#' @importFrom stats density
+#' 
 #' @export
 #'
 
 
-calculateNoiseRatios <- function(rs_assay, rs_control) {
+calculateNoiseRatios <- function(rs_assay, rs_control, noise_method = "max_density") {
   # select control and fractions :
   temp_fractions <- rs_assay[,!(names(rs_assay) %in% c("chr","start","stop"))]
   temp_control <- rs_control[,!(names(rs_control) %in% c("chr","start","stop"))]
   # calculate ratios :
   to_return <- NULL
   for (fraction in names(temp_fractions)) {
     temp_fraction_vs_control <- round(rs_assay[,fraction] / temp_control, digits = 3)
-    temp_density <- density(temp_fraction_vs_control,na.rm = TRUE)
-    temp_max <- which.max(temp_density$y)
-    temp_ratio <- round(temp_density$x[temp_max],digits = 3)
+    temp_density <- stats::density(temp_fraction_vs_control,na.rm = TRUE)
+    if (noise_method == "min_downslope") {
+      temp_downslopes <- which(diff(temp_density$y) < 0)
+      temp_ratio <- temp_density$x[min(temp_downslopes)]
+    }
+    else if (noise_method == "max_density") {
+      temp_max <- which.max(temp_density$y)
+      temp_ratio <- round(temp_density$x[temp_max],digits = 3)
+    }
     to_return <- c(to_return,temp_ratio)
   }
   return(to_return)

R/readRS.R

@@ -1,31 +1,53 @@
-#' @title Read Repli-seq assay data (bedgraph format).
-#' @description Bedgraphs with a one line header and with similar organization (same length and sorted)
+#' @title Read Repli-Seq assay data (Bedgraph format).
+#' @description Read and integrate all Repli-Seq fractions in a single dataframe.
+#'  Fractions must be formatted as Bedgraph with a one line header and with similar resolutions.
 #'
-#' @param paths_data list of paths to repliseq fractions files
-#' @param fractions list of fractions names
+#' @param paths_data Array. List of paths to Repli-Seq fractions files
+#' @param fraction_names Array. List of fractions names (default = S1, ..., Sn)
 #'
 #' @return a dataframe composed of genomic coordinates plus all the specified fractions as for example : chr,start,stop,S1,S2,S3,S4,S5,S6
 #' @import magrittr
 #' @importFrom dplyr select
 #' @importFrom dplyr rename
+#' @importFrom utils read.table
+#' 
+#' @examples 
+#' # basic usage of read_RS :
+#' 
+#' S1_path <- system.file("extdata", "NT_chr22-s1.bdg", package="RepliSeq")
+#' S2_path <- system.file("extdata", "NT_chr22-s2.bdg", package="RepliSeq")
+#' 
+#' my_rs_assay <- readRS(paths_data = c(S1_path , S2_path),
+#'                       fraction_names = c("S1","S2"))
+#' 
 #' @export
 
 
-readRS <- function(paths_data,fractions) {
+readRS <- function(paths_data = NULL,
+                   fraction_names = NULL) {
+
+  V1 <- V2 <- V3 <- NULL
+
+  ### Create default faction_names
+  if (is.null(fraction_names)) {
+    amount_of_fractions <- length(paths_data)
+    fractions_names <- paste0("S",c(1:amount_of_fractions))
+  }
 
   #### Control that path data and fractions are the same length
-  if (length(paths_data) != length(fractions)) {
+  if (length(paths_data) != length(fraction_names)) {
     stop("arguments paths_data and fractions must have the same sizes")
   }
   ### read genomic coordinates from first fraction :
-  rs_temp <- read.table(paths_data[1] , header = FALSE , skip=1 , sep="") %>%
+  rs_temp <- utils::read.table(paths_data[1] , header = FALSE , skip=1 , sep="") %>%
     select(V1,V2,V3) %>%
     rename(chr = V1, start = V2, stop = V3)
+
   ### read other fractions and create the corresponding columns in rs_temp
   for (i in c(1:length(paths_data))) {
-    temp_fraction <- fractions[i]
+    temp_fraction <- fraction_names[i]
     temp_path <- paths_data[i]
-    temp_fraction_df <- read.table(temp_path,header=FALSE,skip=1,sep="")
+    temp_fraction_df <- utils::read.table(temp_path,header=FALSE,skip=1,sep="")
     rs_temp[temp_fraction] <- temp_fraction_df$V4
   }
 

R/repliseqCorrPlot.R

@@ -0,0 +1,63 @@
+#' @title Calculate and plot correlation plots.
+#' @description Calculate and draw a correlation plot from Repli-seq assay (data.frame) see http://www.sthda.com/french/wiki/visualiser-une-matrice-de-correlation-par-un-correlogramme
+#'
+#' @param rs_assay a Repli-seq assay (data.frame) loaded with readRS() and formatted as chr,start,stop,S1,...,Sx
+#' @param corr_method character in either "pearson","spearman" or "kendall" to use as a parameter for cor.test function
+#' @param ordering a character in either "original","AOE", "FPC", "hclust", "alphabet" to reorder the compared variables (to fill the corrplot() order parameter )
+#' @param plottype from corrplot type argument. either "upper", "lower" or "mixed". default = "mixed".
+#' 
+#' @return a mixed correlation plot
+#' @importFrom RColorBrewer brewer.pal
+#' @importFrom corrplot corrplot.mixed corrplot
+#' @importFrom stats cor.test cor
+#' 
+#' @export
+#'
+
+
+repliseqCorrPlot <- function(rs_assay,
+                             corr_method = "pearson",
+                             ordering = "original",
+                             plottype="mixed") {
+
+  temp_rs_assay <- rs_assay[,!(names(rs_assay) %in% c("chr","start","stop"))]
+  # ressource : http://www.sthda.com/french/wiki/visualiser-une-matrice-de-correlation-par-un-correlogramme
+  # Calculate correlation matrix
+  temp.cor <- stats::cor(temp_rs_assay,method = corr_method)
+  # Calculate p-values associated with each correlation
+  p.mat <- cor.mtest(temp.cor)
+  # plot the correlation matrix
+  if (plottype == "mixed") {
+    corrplot.mixed(temp.cor,is.corr = FALSE,
+                   order=ordering, tl.col= "black",addrect = 2,
+                   upper.col = brewer.pal(n = 10, name = "PuOr"),
+                   lower.col= brewer.pal(n = 10, name = "PuOr"),
+                   tl.cex=0.8,number.cex = 0.8,cl.cex = 0.8)
+  }
+  else {
+    corrplot(temp.cor,is.corr = FALSE, type = plottype,
+             order=ordering, tl.col="black",addrect = 2,
+             upper.col = brewer.pal(n = 10, name = "PuOr"),
+             lower.col= brewer.pal(n = 10, name = "PuOr"),
+             tl.cex=0.8,number.cex = 0.8,cl.cex = 0.8)
+  }
+
+}
+
+###########
+
+# Calculation of the p-values associated with a correlation matrix
+cor.mtest <- function(mat, ...) {
+  mat <- as.matrix(mat)
+  n <- ncol(mat)
+  p.mat<- matrix(NA, n, n)
+  diag(p.mat) <- 0
+  for (i in 1:(n - 1)) {
+    for (j in (i + 1):n) {
+      tmp <- stats::cor.test(mat[, i], mat[, j], ...)
+      p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
+    }
+  }
+  colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
+  p.mat
+}

R/repliseqDensityPlot.R

@@ -0,0 +1,50 @@
+#' @title Calculate and plot density plots.
+#' @description Calculate and draw a density plot from Repli-seq assay (data.frame)
+#'
+#' @param rs_assay a Repli-seq assay (data.frame) loaded with readRS() and formatted as chr,start,stop,S1,...,Sx
+#' @param grided boolean (default = F) if True then each fraction s displayed on a different panel, otherwise all the fractions are colored in a single panel
+#' @param stat "counts" (default) or "density" to use in the stat argument
+#' @param ncols (default = 3) amount of columns when grided == T (see ggplot2 facet_wrap ncol parameter) 
+#' @param rs_resolution resolution of the rs_assay to display in the labels in the plots. Default equals "50kb"
+#'
+#' @return a ggplot object displaying density for the provided repliseq assay
+#' @import ggplot2
+#' @importFrom reshape2 melt
+#' @export
+#'
+
+
+repliseqDensityPlot <- function (rs_assay,
+                                 stat = "counts",
+                                 grided = F,
+                                 ncols = 3,
+                                 rs_resolution = "50kb") {
+
+  value <- variable <- NULL
+  # extract and melt the fractions :
+  temp_fractions <- rs_assay[,!(names(rs_assay) %in% c("chr","start","stop"))]
+  temp_fractions_melted <- melt(temp_fractions)
+
+
+  label_x <- paste0("Counts of reads / ",rs_resolution, "windows")
+  label_y <- paste0(stat," of windows")
+
+  # grided plot :
+  if (grided) {
+    p <- ggplot(temp_fractions_melted,aes(x = value)) +
+      geom_line(stat=stat) +
+      facet_wrap(~ variable,ncol=ncols) +
+      theme_bw(20) +
+      labs(x = label_x, y = label_y)
+  }
+
+  # not grided plot :
+  else {
+    p <- ggplot(temp_fractions_melted,aes(x = value,color=variable)) +
+      geom_line(stat=stat) +
+      theme_bw(20) +
+      labs(x = label_x, y = label_y ,color="sample")
+  }
+
+  return(p)
+}

R/rescaleRS.R

@@ -13,15 +13,18 @@
 #' 
 
 
-rescaleRS <- function(rs_assay,scale_factor) {
+rescaleRS <- function(rs_assay,
+                      scale_factor) {
+
+  chr <- NULL
   # calculate_inital scale :
   initial_scale <- (rs_assay$stop - rs_assay$start)[1]
   # get names
   rs_names <- names(rs_assay)
   rs_fractions <- rs_names[rs_names != "chr" & rs_names != "start" & rs_names != "stop"]
   # initialize to_return variable
   to_return <- NULL
-  rs_chroms <- levels(unique(rs_assay$chr))
+  rs_chroms <- unique(rs_assay$chr)
   # iterate on chroms :
   for (i in rs_chroms) {
     # select data :

R/smoothRS.R

@@ -13,7 +13,11 @@
 #' 
 
 
-smoothRS <- function(rs_assay,smooth_factor) {
+smoothRS <- function(rs_assay,
+                     smooth_factor) {
+
+  chr <- NULL
+
   # calculate_inital scale :
   initial_scale <- (rs_assay$stop - rs_assay$start)[1]
   # get names

R/writeBedgraph.R

@@ -6,6 +6,8 @@
 #' @param sample_name a string for the file names
 #' @param system_separator default is "/" as dedicated for linux file system
 #'
+#' @importFrom utils write.table
+#'
 #' @return NULL
 #' @export
 #'
@@ -21,6 +23,6 @@ writeBedgraph <- function(rs_assay,path_file,sample_name,system_separator = "/")
     to_write <- temp_coords
     to_write[i] <- temp_fraction
     file_name <- paste0(path_file,system_separator,sample_name,"-",i,".bdg")
-    write.table(to_write,file = file_name,quote = FALSE, col.names = header,row.names = FALSE)
+    utils::write.table(to_write,file = file_name,quote = FALSE, col.names = header,row.names = FALSE)
   }
 }

R/writeBigWig.R

@@ -10,6 +10,8 @@
 #' @param wiggle_span a numeric value for wiggle header
 #' @param system_separator default is "/" as dedicated for linux file system
 #'
+#' @importFrom utils write.table
+#'
 #' @return NULL
 #' @export
 #' 
@@ -29,7 +31,7 @@ writeBigwig <- function(rs_assay,path_file,sample_name,chromsizes,wiggle_start,w
     }
     # write wiggle ; convert to BigWig and rm wiggle :
     file_name <- paste0(path_file,system_separator,sample_name,"-",i,".wig")
-    write.table(to_write,file = file_name,quote = FALSE, col.names = temp_header,row.names = FALSE)
+    utils::write.table(to_write,file = file_name,quote = FALSE, col.names = temp_header,row.names = FALSE)
     system(paste0("wigToBigWig -clip ",file_name," ",chromsizes," ",file_name,".bw"))
     system(paste0("rm ",file_name))
   }