--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ribowaltz.R	Thu Sep 22 20:30:20 2022 +0000
@@ -0,0 +1,131 @@
+# setup R error handling to go to stderr
+options(show.error.messages = FALSE, error = function() {
+  cat(geterrmessage(), file = stderr())
+  q("no", 1, FALSE)
+})
+
+# we need that to not crash galaxy with an UTF8 error on German LC settings.
+loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
+library("getopt")
+options(stringAsFactors = FALSE, useFancyQuotes = FALSE)
+args <- commandArgs(trailingOnly = TRUE)
+
+# get options, using the spec as defined by the enclosed list.
+# we read the options from the default: commandArgs(TRUE).
+spec <- matrix(c(
+  "quiet", "q", 0, "logical",
+  "help", "h", 0, "logical",
+  "bamdir", "b", 1, "character",
+  "gtffile", "g", 1, "character",
+  "codon_coverage_info", "Y", 1, "character",
+  "cds_coverage_info", "Z", 1, "character",
+  "psite_info_rdata", "O", 0, "character",
+  "refseq_sep", "s", 0, "character",
+  "params_duplicate_filterting", "d", 0, "character",
+  "params_peridiocity_filterting", "l", 0, "character",
+  "params_custom_filterting", "c", 0, "character",
+  "params_psite_additional", "p", 0, "character",
+  "params_coverage_additional", "C", 0, "character"
+), byrow = TRUE, ncol = 4)
+opt <- getopt(spec)
+
+# if help was asked for print a friendly message
+# and exit with a non-zero error code
+if (!is.null(opt$help)) {
+  cat(getopt(spec, usage = TRUE))
+  q(status = 1)
+}
+
+verbose <- is.null(opt$quiet)
+
+library("riboWaltz")
+
+# create annotation data table
+annotation_dt <- create_annotation(opt$gtffile)
+
+sep <- opt$refseq_sep
+if (opt$refseq_sep == "") {
+  sep <- NULL
+}
+# convert alignments in BAM files into list of data tables
+reads_list <- bamtolist(bamfolder = opt$bamdir, annotation = annotation_dt, refseq_sep = sep)
+
+library("jsonlite")
+# remove duplicate reads
+if (!is.null(opt$params_duplicate_filterting)) {
+  json_duplicate_filterting <- fromJSON(opt$params_duplicate_filterting)
+  reads_list <- duplicates_filter(
+    data = reads_list,
+    extremity = json_duplicate_filterting$extremity,
+    keep = json_duplicate_filterting$keep
+  )
+}
+
+# selection of read lengths - periodicity filtering
+if (!is.null(opt$params_peridiocity_filterting)) {
+    json_peridiocity_filterting <- fromJSON(opt$params_peridiocity_filterting)
+    reads_list <- length_filter(
+      data = reads_list,
+      length_filter_mode = "periodicity",
+      periodicity_threshold = json_peridiocity_filterting$periodicity_threshold
+    )
+}
+# selection of read lengths - length range filtering
+if (!is.null(opt$params_custom_filterting)) {
+    json_custom_filterting <- fromJSON(opt$params_custom_filterting)
+    reads_list <- length_filter(
+      data = reads_list,
+      length_filter_mode = "custom",
+      length_range = json_custom_filterting$length_range
+    )
+}
+
+# compute P-site offset
+json_psite_additional <- fromJSON(opt$params_psite_additional)
+psite_offset <- psite(
+  reads_list,
+  start = json_psite_additional$use_start,
+  flanking = json_psite_additional$flanking,
+  extremity = json_psite_additional$psite_extrimity,
+  plot = TRUE,
+  cl = json_psite_additional$cl,
+  plot_format = "pdf",
+  plot_dir = "plots"
+)
+psite_offset
+
+reads_psite_list <- psite_info(reads_list, psite_offset)
+reads_psite_list
+# write a separate P-site offset info table for each sample
+for (sample in names(reads_psite_list)) {
+  write.table(
+    reads_psite_list[[sample]],
+    file = paste(sample, "psite_info.tsv",  sep = "_"),
+    sep = "\t", row.names = FALSE, quote = FALSE
+  )
+  print(paste(sample, "psite_info.tsv",  sep = "_"))
+}
+
+# write R object to a file
+if (!is.null(opt$psite_info_rdata)) {
+  save(reads_psite_list, annotation_dt, file = opt$psite_info_rdata)
+}
+
+json_coverage_additional <- fromJSON(opt$params_coverage_additional)
+# codon coverage
+codon_coverage_out <- codon_coverage(
+  reads_psite_list,
+  annotation_dt,
+  psite = json_coverage_additional$psites_per_region,
+  min_overlap = json_coverage_additional$min_overlap
+)
+write.table(codon_coverage_out, file = opt$codon_coverage_info, sep = "\t", row.names = FALSE, quote = FALSE)
+
+# CDS coverage
+cds_coverage_out <- cds_coverage(
+  reads_psite_list,
+  annotation_dt,
+  start_nts = json_coverage_additional$start_nts,
+  stop_nts = json_coverage_additional$stop_nts
+)
+write.table(cds_coverage_out, file = opt$cds_coverage_info, sep = "\t", row.names = FALSE, quote = FALSE)