diff ideas_preprocessor.R @ 1:91c5dbb14a13 draft

Uploaded
author greg
date Mon, 22 Jan 2018 14:35:28 -0500
parents
children e97851e8951b
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ideas_preprocessor.R	Mon Jan 22 14:35:28 2018 -0500
@@ -0,0 +1,62 @@
+#!/usr/bin/env Rscript
+
+suppressPackageStartupMessages(library("data.table"))
+suppressPackageStartupMessages(library("optparse"))
+
+option_list <- list(
+            make_option(c("--bed_input"), action="store", dest="bed_input", defaul=NULL, help="Chromosome windows positions file"),
+            make_option(c("--exclude_input"), action="store", dest="exclude_input", defaul=NULL, help="File(s) containing regions to exclude"),
+            make_option(c("--bychr"), action="store_true", dest="bychr", defaul=FALSE, help="Separate files by chromosome"),
+            make_option(c("--chrom_len_file"), action="store", dest="chrom_len_file", default=NULL, help="Chromosome lengths file"),
+            make_option(c("--output"), action="store", dest="output", help="Primary output dataset"),
+            make_option(c("--output_files_path"), action="store", dest="output_files_path", help="Primary output dataset extra files path"),
+            make_option(c("--prep_input_config"), action="store", dest="prep_input_config", help="Preprocessing input config file"),
+            make_option(c("--reads_per_bp"), action="store", dest="reads_per_bp", type="integer", help="Calculate the signal in each genomic window"),
+            make_option(c("--restrict_to_chroms"), action="store", dest="restrict_to_chroms", default=NULL, help="Restrict processing to specified chromosomes"),
+            make_option(c("--standardize_datasets"), action="store_true", dest="standardize_datasets", default=FALSE, help="Standardize all datasets"),
+            make_option(c("--window_size"), action="store", dest="window_size", type="integer", default=NULL, help="Window size in base pairs")
+)
+
+parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
+args <- parse_args(parser, positional_arguments=TRUE)
+opt <- args$options
+
+create_primary_html = function(output, output_files_path) {
+    tmp_files = list.files(path=output_files_path);
+    s <- paste('<html><head></head><body>', sep="\n");
+    s <- paste(s, '<h3>Contents of directory required by IDEAS</h3>\n', sep="");
+    s <- paste(s, '<ul>\n', sep="");
+    s <- paste(s, '<li><a href="prep_output_config.txt">', prep_output_config.txt, '</a></li>\n', sep="");
+    for (i in 1:length(tmp_files)) {
+        s <- paste(s, '<li><a href="', 'tmp/', tmp_files[i], '">', tmp_files[i], '</a></li>\n', sep="");
+    }
+    s <- paste(s, '</ul>\n</body>\n</html>', sep="");
+    cat(s, file=output);
+}
+
+# Create the directories that will contain all of the output files.
+dir.create(opt$output_files_path, showWarnings=FALSE);
+dir.create(paste(opt$output_files_path, "tmp", sep="/"), showWarnings=FALSE);
+
+# Read the prep_input_config text file which has this format:
+# "cell type name" "epigenetic factor name" "file path" "file name" "datatype"
+prep_input_config_matrix = as.matrix(read.table(opt$prep_input_config));
+# Process data to windows mean.
+for (i in 1:dim(prep_input_config_matrix)[1]) {
+    file_path = prep_input_config_matrix[i, 3]
+    file_name = prep_input_config_matrix[i, 4]
+    datatype = prep_input_config_matrix[i, 5]
+    if (datatype == "bam") {
+        system(paste("samtools index", file_path));
+        bw = paste(file_name, "bw", sep=".");
+        system(paste("bamCoverage --bam", file_path, "-o", bw, "--binSize", opt$window_size));
+    } else {
+        bw = file_path;
+    }
+    bd = paste(output_files_path, "tmp", file_name, ".bed", sep="/");
+    system(paste("bigWigAverageOverBed", bw, opt$bed_input, "stdout | cut -f5 >", bd));
+    system(paste("gzip -f", bd));
+}
+
+# Create the primary HTML dataset.
+create_primary_html(opt$output, opt$output_files_path);