Mercurial > repos > greg > ideas_preprocessor
changeset 27:19881f817d25 draft
Uploaded
author | greg |
---|---|
date | Mon, 05 Feb 2018 13:42:11 -0500 |
parents | 06e63b5cb025 |
children | 55d264cd49e8 |
files | ideas_preprocessor.R |
diffstat | 1 files changed, 54 insertions(+), 47 deletions(-) [+] |
line wrap: on
line diff
--- a/ideas_preprocessor.R Mon Feb 05 13:42:03 2018 -0500 +++ b/ideas_preprocessor.R Mon Feb 05 13:42:11 2018 -0500 @@ -1,30 +1,15 @@ #!/usr/bin/env Rscript -# TODO: implement support for the following: -# 1. Scenario where user did not select chrom_bed_input -# 2. --exclude_bed_input -# 3. --bychr -# 4. --chrom_len_file -# 5. --reads_per_bp -# 6. --restrict_to_chroms -# 7. --standardize_datasets -# 8. Scenario where --window_size is NULL and need to handle bamCoverage - see TODO near line # 57. - - suppressPackageStartupMessages(library("data.table")) suppressPackageStartupMessages(library("optparse")) option_list <- list( make_option(c("--chrom_bed_input"), action="store", dest="chrom_bed_input", defaul=NULL, help="Chromosome windows positions file"), make_option(c("--exclude_bed_input"), action="store", dest="exclude_bed_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("--ideaspre_input_config"), action="store", dest="ideaspre_input_config", help="Preprocessing input config 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("--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("--chromosome_windows"), action="store", dest="chromosome_windows", default=NULL, help="Windows positions by chroms config file"), make_option(c("--window_size"), action="store", dest="window_size", type="integer", default=NULL, help="Window size in base pairs") ) @@ -34,41 +19,64 @@ opt <- args$options tmp_dir = "tmp"; +cbi_file = "chrom_bed_input.bed"; +if (is.null(opt$chrom_bed_input)) { + # Create a chromosome windows positions file + # using the received chromosome lengths file + # and the window size. + cmd = paste("bedtools makewindows -g", opt$chrom_len_file, "-w", opt$window_size, ">", cbi_file, sep=" "); + system(cmd); +} else { + if (!is.null(opt$exclude_bed_input)) { + # Copy the received chrom_bed_input + # since we will alter it. + file.copy(opt$chrom_bed_input, cbi_file); + } else { + cbi_file = opt$chrom_bed_input; + } +} +# Exclude regions if necessary. +if (!is.null(opt$exclude_bed_input)) { + exclude_bed_inputs = as.character(opt$exclude_bed_input); + exclude_bed_files = strsplit(exclude_bed_inputs, ","); + tmp_file = paste("tmp", cbi_file, sep="_"); + for (exclude_bed_file in exclude_bed_files) { + cmd = paste("bedtools subtract -a", cbi_file, "-b", exclude_bed_file, ">", tmp_file, sep=" "); + system(cmd); + cmd = paste("mv", tmp_file, cbi_file, sep=" "); + system(cmd); + } +} +# Read the chromosome windows positions file +# to get the smallest window size in the file +# (i.e., the minimum of column 3 - column 2. +cbi = fread(cbi_file); +min_window_size = min(cbi[,3]-cbi[,2]); # Read the ideaspre_input_config text file which has this format: # "cell type name" "epigenetic factor name" "file path" "file name" "datatype" ideaspre_input_config = as.matrix(read.table(opt$ideaspre_input_config)); - -# TODO: fix this -window_size = opt$window_size -if (is.null(opt$window_size)) { - window_size = 500; +# Process data to windows mean. +for (i in 1:dim(ideaspre_input_config)[1]) { + file_path = ideaspre_input_config[i, 3] + file_name = ideaspre_input_config[i, 4] + datatype = ideaspre_input_config[i, 5] + if (datatype == "bam") { + cmd = paste("samtools index", file_path); + system(cmd); + bigwig_file_name = paste(file_name, "bw", sep="."); + cmd = paste("bamCoverage --bam", file_path, "-o", bigwig_file_name, "--binSize", min_window_size); + system(cmd); + } else { + bigwig_file_name = file_path; + } + bed_file_name = paste(file_name, "bed", sep="."); + bed_file_path = paste("tmp", bed_file_name, sep="/"); + cmd = paste("bigWigAverageOverBed", bigwig_file_name, opt$chrom_bed_input, "stdout | cut -f5 >", bed_file_path); + system(cmd); + cmd = paste("gzip -f", bed_file_path); + system(cmd); } -# Process data to windows mean. -if (!is.null(opt$chrom_bed_input)) { - for (i in 1:dim(ideaspre_input_config)[1]) { - file_path = ideaspre_input_config[i, 3] - file_name = ideaspre_input_config[i, 4] - datatype = ideaspre_input_config[i, 5] - if (datatype == "bam") { - cmd = paste("samtools index", file_path); - system(cmd); - bigwig_file_name = paste(file_name, "bw", sep="."); - - cmd = paste("bamCoverage --bam", file_path, "-o", bigwig_file_name, "--binSize", window_size); - system(cmd); - } else { - bigwig_file_name = file_path; - } - bed_file_name = paste(file_name, "bed", sep="."); - bed_file_path = paste("tmp", bed_file_name, sep="/"); - cmd = paste("bigWigAverageOverBed", bigwig_file_name, opt$chrom_bed_input, "stdout | cut -f5 >", bed_file_path); - system(cmd); - cmd = paste("gzip -f", bed_file_path); - system(cmd); - } -} - # Create file1.txt. cmd = paste("cut -d' '", opt$ideaspre_input_config, "-f1,2 > file1.txt", sep=" "); system(cmd); @@ -93,7 +101,7 @@ # Move the tmp archive to the output directory. to_path = paste(opt$output_files_path, "tmp.tar.gz", sep="/"); file.rename("tmp.tar.gz", to_path); - +# Handle file names for display in the primary dataset if necessary. if (!is.null(opt$chrom_bed_input) && !is.null(opt$chromosome_windows)) { # Renane opt$chrom_bed_input to be chromosomes.bed # and make a copy of it in the output directory. @@ -103,4 +111,3 @@ to_path = paste(opt$output_files_path, opt$chromosome_windows, sep="/"); file.rename(opt$chromosome_windows, to_path); } -