--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/bed2bigwig.xml	Tue Apr 30 08:27:27 2024 +0000
@@ -0,0 +1,32 @@
+<tool id="calculate_density" name="Calculate Density" version="1.0.0">
+    <description>Calculate coverage density from BED/GFF files</description>
+    <requirements>
+        <requirement type="package" version="1.7.3">r-optparse</requirement>
+        <requirement type="package" version="1.54.0">bioconductor-rtracklayer</requirement>
+    </requirements>
+    <required_files>
+        <include type="literal" path="calculate_density.R"/>
+    </required_files>
+    <command detect_errors="exit_code"><![CDATA[
+        Rscript '$__tool_directory__/calculate_density.R'
+        --bed '$input_bed'
+        --window '$window_size'
+        --output '$output_bigwig'
+        --format '$input_bed.ext'
+        #if $merge:
+          --merge
+        #end if
+    ]]></command>
+    <inputs>
+        <param type="data" format="bed,gff" name="input_bed" label="Input BED/GFF file" help="BED or GFF file containing genomic regions" />
+        <param type="integer" name="window_size" label="Window size" value="1000000" min="1" help="Size of the window for density calculation" />
+        <param type="boolean" name="merge" label="Merge overlapping regions" truevalue="true" falsevalue="false" checked="true" help="Merge overlapping regions before calculating density" />
+    </inputs>
+    <outputs>
+        <data format="bigwig" name="output_bigwig" label="${tool.name} on ${on_string}: BigWig output"/>
+    </outputs>
+    <help><![CDATA[
+        This tool calculates the coverage density from a provided BED or GFF file over specified window sizes and outputs the results in BigWig format.
+    ]]></help>
+
+</tool>
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/calculate_density.R	Tue Apr 30 08:27:27 2024 +0000
@@ -0,0 +1,47 @@
+#!/usr/bin/env Rscript
+library(optparse)
+
+get_density <- function(x, chr_size=NULL, tw=1000000){
+  cvg <- coverage(x)
+  bins <- tileGenome(chr_size, tilewidth = tw)
+  d <- binnedAverage(unlist(bins), cvg, "score")
+  d
+}
+max_chr_length <- function(g){
+  x <- split(g, ~seqnames)
+  L <- sapply(x, function(x)max(end(x), na.rm=TRUE))
+  L
+}
+
+# get input arguments
+# bed file
+# window size
+# output bigwig file
+option_list <- list(
+  make_option(c("-b", "--bed"), type="character", default=NULL, help="BED or GFF file"),
+  make_option(c("-w", "--window"), type="integer", default=1000000, help="Window size"),
+  make_option(c("-o", "--output"), type="character", default=NULL, help="Output BigWig file"),
+  make_option(c("-f", "--format"), type="character", default="gff3", help="Input format (gff3 or bed)"),
+  make_option(c("-m", "--merge"), type="logical", action="store_true", default=FALSE, help="Merge overlapping regions")
+
+
+)
+
+opt_parser <- OptionParser(option_list=option_list)
+opt <- parse_args(opt_parser)
+
+# check mandatory arguments
+if (is.null(opt$bed) || is.null(opt$output)){
+  stop("Please provide bed file and output file")
+}
+
+suppressPackageStartupMessages(library(rtracklayer))
+print(opt)
+g <- import(opt$bed, format=opt$format)
+if (opt$merge){
+  g <- reduce(g)
+}
+chr_size <- max_chr_length(g)
+d <- get_density(g, chr_size, opt$window)
+
+export(d, opt$output, format="bigwig")