Mercurial > repos > petr-novak > various_galaxy_tools

comparison calculate_density.R @ 3:1069776f7ae2 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/kavonrtep/galaxy_packages commit 3b9f93ed06cc32dbfa271789739e7a1e8fac528c
author petr-novak
date Tue, 30 Apr 2024 08:27:27 +0000
parents
children
comparison
equal deleted inserted replaced
2:90c46e40d222 3:1069776f7ae2
1 #!/usr/bin/env Rscript
2 library(optparse)
3
4 get_density <- function(x, chr_size=NULL, tw=1000000){
5 cvg <- coverage(x)
6 bins <- tileGenome(chr_size, tilewidth = tw)
7 d <- binnedAverage(unlist(bins), cvg, "score")
8 d
9 }
10 max_chr_length <- function(g){
11 x <- split(g, ~seqnames)
12 L <- sapply(x, function(x)max(end(x), na.rm=TRUE))
13 L
14 }
15
16 # get input arguments
17 # bed file
18 # window size
19 # output bigwig file
20 option_list <- list(
21 make_option(c("-b", "--bed"), type="character", default=NULL, help="BED or GFF file"),
22 make_option(c("-w", "--window"), type="integer", default=1000000, help="Window size"),
23 make_option(c("-o", "--output"), type="character", default=NULL, help="Output BigWig file"),
24 make_option(c("-f", "--format"), type="character", default="gff3", help="Input format (gff3 or bed)"),
25 make_option(c("-m", "--merge"), type="logical", action="store_true", default=FALSE, help="Merge overlapping regions")
26
27
28 )
29
30 opt_parser <- OptionParser(option_list=option_list)
31 opt <- parse_args(opt_parser)
32
33 # check mandatory arguments
34 if (is.null(opt$bed) || is.null(opt$output)){
35 stop("Please provide bed file and output file")
36 }
37
38 suppressPackageStartupMessages(library(rtracklayer))
39 print(opt)
40 g <- import(opt$bed, format=opt$format)
41 if (opt$merge){
42 g <- reduce(g)
43 }
44 chr_size <- max_chr_length(g)
45 d <- get_density(g, chr_size, opt$window)
46
47 export(d, opt$output, format="bigwig")