refactoring; renamed reduce_bed to reduce_sequence

bin/reduce_bed.R → bin/reduce_sequence.R

@@ -9,32 +9,35 @@ option_list <- list(
   make_option(opt_str = c("-t", "--threads"), default = 1, help = "Number of threads to use. Use 0 for all available cores. Default = %default", metavar = "integer"),
   make_option(opt_str = c("-c", "--clean"), default = TRUE, help = "Delete all temporary files. Default = %default", metavar = "logical"),
   make_option(opt_str = c("-s", "--min_seq_length"), default = NULL, help = "Remove sequences below this length. Defaults to the maximum value of motif and kmer and can not be lower.", metavar = "integer", type = "integer"),
-  make_option(opt_str = c("-n", "--minoverlap_kmer"), default = NULL, help = "Minimum required overlap between kmer to merge kmer. Used to create reduced sequence ranges. Can not be greater than kmer length. Default = kmer - 1", metavar = "integer", type = "integer"),
+  make_option(opt_str = c("-n", "--minoverlap_kmer"), default = NULL, help = "Minimum required overlap between kmer. Used to create reduced sequence ranges out of merged kmer. Can not be greater than kmer length. Default = kmer - 1", metavar = "integer", type = "integer"),
   make_option(opt_str = c("-v", "--minoverlap_motif"), default = NULL, help = "Minimum required overlap between motif and kmer to consider kmer significant. Used for kmer cutoff calculation. Can not be greater than motif and kmer length. Default = ceiling(motif / 2)", metavar = "integer", type = "integer"),
-  make_option(opt_str = c("-f", "--motif_occurence"), default = 1, help = "Number of motifs per sequence any value above 0. Default = %default.", metavar = "double")
+  make_option(opt_str = c("-f", "--motif_occurence"), default = 1, help = "Define how many motifs are expected per sequence. This value is used during kmer cutoff calculation. Default = %default meaning that there should be approximately one motif per sequence.", metavar = "double")
 )
 
 opt_parser <- OptionParser(option_list = option_list, 
-                           description = "Reduce sequences to frequent regions.")
+                           description = "Reduces each sequence to its most frequent region.")
 
 opt <- parse_args(opt_parser)
 
-#' Reduce bed file to conserved regions
+#' Reduces each sequence to its most frequent region.
 #' 
-#' @param input bed file
-#' @param kmer Length of kmer.
-#' @param motif Estimated motif length.
-#' @param output Output file
-#' @param threads Number of threads. Default = 1. 0 for all cores.
+#' @param input Input bed-file. Last column must be sequences.
+#' @param kmer Kmer length. Default = 10
+#' @param motif Estimated motif length. Default = 10
+#' @param output Output file. Default = reduced.bed
+#' @param threads Number of threads to use. Default = 1. Use 0 for all cores.
 #' @param clean Delete all temporary files.
-#' @param minoverlap_kmer Minimum required overlap between kmer to merge kmer. Used to create reduced sequence ranges. Can not be greater than kmer length. Default = kmer - 1
+#' @param minoverlap_kmer Minimum required overlap between kmer. Used to create reduced sequence ranges out of merged kmer. Can not be greater than kmer length . Default = kmer - 1
 #' @param minoverlap_motif Minimum required overlap between motif and kmer to consider kmer significant. Used for kmer cutoff calculation. Can not be greater than motif and kmer length. Default = ceiling(motif / 2)
-#' @param min_seq_length Must be greater or equal to kmer and motif. Default = max(c(motif, kmer)).
+#' @param min_seq_length Remove sequences below this length. Defaults to the maximum value of motif and kmer and can not be lower.
 #' @param motif_occurence Define how many motifs are expected per sequence. This value is used during kmer cutoff calculation. Default = 1 meaning that there should be approximately one motif per sequence.
 #' 
-#' @return reduced bed
 #' TODO check whether jellyfish is installed
-reduce_bed <- function(input, kmer = 10, motif = 10, output = "reduced.bed", threads = NULL, clean = TRUE, minoverlap_kmer = kmer - 1, minoverlap_motif = ceiling(motif / 2), min_seq_length = max(c(motif, kmer)), motif_occurence = 1) {
+reduce_sequence <- function(input, kmer = 10, motif = 10, output = "reduced.bed", threads = NULL, clean = TRUE, minoverlap_kmer = kmer - 1, minoverlap_motif = ceiling(motif / 2), min_seq_length = max(c(motif, kmer)), motif_occurence = 1) {
+  if (missing(input)) {
+    stop("No input specified! Please forward a valid bed-file.")
+  }
+
   # get number of available cores
   if (threads == 0) {
     threads <- parallel::detectCores()
@@ -117,12 +120,12 @@ reduce_bed <- function(input, kmer = 10, motif = 10, output = "reduced.bed", thr
   merged[, sequence := stringr::str_sub(sequence, relative_start, relative_end)]  
 
   # bed files count from 0
-  merged[, `:=`(relative_start = relative_start - 1, relative_end = relative_end - 1)]
+  merged[, data.table::`:=`(relative_start = relative_start - 1, relative_end = relative_end - 1)]
   # change start end location
-  merged[, `:=`(start = start + relative_start, end = start + relative_end)]
+  merged[, data.table::`:=`(start = start + relative_start, end = start + relative_end)]
 
   # clean table
-  merged[, `:=`(relative_start = NULL, relative_end = NULL, width = NULL)]
+  merged[, data.table::`:=`(relative_start = NULL, relative_end = NULL, width = NULL)]
 
   if (clean) {
     file.remove(fasta_file, count_output_binary, mer_count_table)
@@ -171,6 +174,7 @@ significant_kmer <- function(bed, kmer, motif, minoverlap = ceiling(motif / 2),
 reduce_kmer <- function(kmer, significant) {
   data.table::setorderv(kmer, cols = names(kmer)[2], order = -1)
 
+  # TODO don't use 'V2'
   kmer[, cumsum := cumsum(V2)]
 
   return(kmer[cumsum <= significant])
@@ -255,5 +259,5 @@ if (!interactive()) {
   pbo <- pbapply::pboptions(type = "timer")
   # remove last parameter (help param)
   params <- opt[-length(opt)]
-  do.call(reduce_bed, args = params)  
+  do.call(reduce_sequence, args = params)  
 }