Mercurial > repos > greg > ideas_preprocessor

comparison ideas_preprocessor.R @ 27:19881f817d25 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Mon, 05 Feb 2018 13:42:11 -0500
parents f7563bb242fc
children 3b3001355f44
comparison
equal deleted inserted replaced
26:06e63b5cb025 27:19881f817d25
1 #!/usr/bin/env Rscript 1 #!/usr/bin/env Rscript
2
3 # TODO: implement support for the following:
4 # 1. Scenario where user did not select chrom_bed_input
5 # 2. --exclude_bed_input
6 # 3. --bychr
7 # 4. --chrom_len_file
8 # 5. --reads_per_bp
9 # 6. --restrict_to_chroms
10 # 7. --standardize_datasets
11 # 8. Scenario where --window_size is NULL and need to handle bamCoverage - see TODO near line # 57.
12
13 2
14 suppressPackageStartupMessages(library("data.table")) 3 suppressPackageStartupMessages(library("data.table"))
15 suppressPackageStartupMessages(library("optparse")) 4 suppressPackageStartupMessages(library("optparse"))
16 5
17 option_list <- list( 6 option_list <- list(
18 make_option(c("--chrom_bed_input"), action="store", dest="chrom_bed_input", defaul=NULL, help="Chromosome windows positions file"), 7 make_option(c("--chrom_bed_input"), action="store", dest="chrom_bed_input", defaul=NULL, help="Chromosome windows positions file"),
19 make_option(c("--exclude_bed_input"), action="store", dest="exclude_bed_input", defaul=NULL, help="File(s) containing regions to exclude"), 8 make_option(c("--exclude_bed_input"), action="store", dest="exclude_bed_input", defaul=NULL, help="File(s) containing regions to exclude"),
20 make_option(c("--bychr"), action="store_true", dest="bychr", defaul=FALSE, help="Separate files by chromosome"),
21 make_option(c("--chrom_len_file"), action="store", dest="chrom_len_file", default=NULL, help="Chromosome lengths file"), 9 make_option(c("--chrom_len_file"), action="store", dest="chrom_len_file", default=NULL, help="Chromosome lengths file"),
22 make_option(c("--ideaspre_input_config"), action="store", dest="ideaspre_input_config", help="Preprocessing input config file"), 10 make_option(c("--ideaspre_input_config"), action="store", dest="ideaspre_input_config", help="Preprocessing input config file"),
23 make_option(c("--output"), action="store", dest="output", help="Primary output dataset"), 11 make_option(c("--output"), action="store", dest="output", help="Primary output dataset"),
24 make_option(c("--output_files_path"), action="store", dest="output_files_path", help="Primary output dataset extra files path"), 12 make_option(c("--output_files_path"), action="store", dest="output_files_path", help="Primary output dataset extra files path"),
25 make_option(c("--reads_per_bp"), action="store", dest="reads_per_bp", type="integer", help="Calculate the signal in each genomic window"),
26 make_option(c("--restrict_to_chroms"), action="store", dest="restrict_to_chroms", default=NULL, help="Restrict processing to specified chromosomes"),
27 make_option(c("--standardize_datasets"), action="store_true", dest="standardize_datasets", default=FALSE, help="Standardize all datasets"),
28 make_option(c("--chromosome_windows"), action="store", dest="chromosome_windows", default=NULL, help="Windows positions by chroms config file"), 13 make_option(c("--chromosome_windows"), action="store", dest="chromosome_windows", default=NULL, help="Windows positions by chroms config file"),
29 make_option(c("--window_size"), action="store", dest="window_size", type="integer", default=NULL, help="Window size in base pairs") 14 make_option(c("--window_size"), action="store", dest="window_size", type="integer", default=NULL, help="Window size in base pairs")
30 ) 15 )
31 16
32 parser <- OptionParser(usage="%prog [options] file", option_list=option_list) 17 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
33 args <- parse_args(parser, positional_arguments=TRUE) 18 args <- parse_args(parser, positional_arguments=TRUE)
34 opt <- args$options 19 opt <- args$options
35 20
36 tmp_dir = "tmp"; 21 tmp_dir = "tmp";
22 cbi_file = "chrom_bed_input.bed";
37 23
24 if (is.null(opt$chrom_bed_input)) {
25 # Create a chromosome windows positions file
26 # using the received chromosome lengths file
27 # and the window size.
28 cmd = paste("bedtools makewindows -g", opt$chrom_len_file, "-w", opt$window_size, ">", cbi_file, sep=" ");
29 system(cmd);
30 } else {
31 if (!is.null(opt$exclude_bed_input)) {
32 # Copy the received chrom_bed_input
33 # since we will alter it.
34 file.copy(opt$chrom_bed_input, cbi_file);
35 } else {
36 cbi_file = opt$chrom_bed_input;
37 }
38 }
39 # Exclude regions if necessary.
40 if (!is.null(opt$exclude_bed_input)) {
41 exclude_bed_inputs = as.character(opt$exclude_bed_input);
42 exclude_bed_files = strsplit(exclude_bed_inputs, ",");
43 tmp_file = paste("tmp", cbi_file, sep="_");
44 for (exclude_bed_file in exclude_bed_files) {
45 cmd = paste("bedtools subtract -a", cbi_file, "-b", exclude_bed_file, ">", tmp_file, sep=" ");
46 system(cmd);
47 cmd = paste("mv", tmp_file, cbi_file, sep=" ");
48 system(cmd);
49 }
50 }
51 # Read the chromosome windows positions file
52 # to get the smallest window size in the file
53 # (i.e., the minimum of column 3 - column 2.
54 cbi = fread(cbi_file);
55 min_window_size = min(cbi[,3]-cbi[,2]);
38 # Read the ideaspre_input_config text file which has this format: 56 # Read the ideaspre_input_config text file which has this format:
39 # "cell type name" "epigenetic factor name" "file path" "file name" "datatype" 57 # "cell type name" "epigenetic factor name" "file path" "file name" "datatype"
40 ideaspre_input_config = as.matrix(read.table(opt$ideaspre_input_config)); 58 ideaspre_input_config = as.matrix(read.table(opt$ideaspre_input_config));
41 59 # Process data to windows mean.
42 # TODO: fix this 60 for (i in 1:dim(ideaspre_input_config)[1]) {
43 window_size = opt$window_size 61 file_path = ideaspre_input_config[i, 3]
44 if (is.null(opt$window_size)) { 62 file_name = ideaspre_input_config[i, 4]
45 window_size = 500; 63 datatype = ideaspre_input_config[i, 5]
64 if (datatype == "bam") {
65 cmd = paste("samtools index", file_path);
66 system(cmd);
67 bigwig_file_name = paste(file_name, "bw", sep=".");
68 cmd = paste("bamCoverage --bam", file_path, "-o", bigwig_file_name, "--binSize", min_window_size);
69 system(cmd);
70 } else {
71 bigwig_file_name = file_path;
72 }
73 bed_file_name = paste(file_name, "bed", sep=".");
74 bed_file_path = paste("tmp", bed_file_name, sep="/");
75 cmd = paste("bigWigAverageOverBed", bigwig_file_name, opt$chrom_bed_input, "stdout | cut -f5 >", bed_file_path);
76 system(cmd);
77 cmd = paste("gzip -f", bed_file_path);
78 system(cmd);
46 } 79 }
47 # Process data to windows mean.
48 if (!is.null(opt$chrom_bed_input)) {
49 for (i in 1:dim(ideaspre_input_config)[1]) {
50 file_path = ideaspre_input_config[i, 3]
51 file_name = ideaspre_input_config[i, 4]
52 datatype = ideaspre_input_config[i, 5]
53 if (datatype == "bam") {
54 cmd = paste("samtools index", file_path);
55 system(cmd);
56 bigwig_file_name = paste(file_name, "bw", sep=".");
57
58 cmd = paste("bamCoverage --bam", file_path, "-o", bigwig_file_name, "--binSize", window_size);
59 system(cmd);
60 } else {
61 bigwig_file_name = file_path;
62 }
63 bed_file_name = paste(file_name, "bed", sep=".");
64 bed_file_path = paste("tmp", bed_file_name, sep="/");
65 cmd = paste("bigWigAverageOverBed", bigwig_file_name, opt$chrom_bed_input, "stdout | cut -f5 >", bed_file_path);
66 system(cmd);
67 cmd = paste("gzip -f", bed_file_path);
68 system(cmd);
69 }
70 }
71
72 # Create file1.txt. 80 # Create file1.txt.
73 cmd = paste("cut -d' '", opt$ideaspre_input_config, "-f1,2 > file1.txt", sep=" "); 81 cmd = paste("cut -d' '", opt$ideaspre_input_config, "-f1,2 > file1.txt", sep=" ");
74 system(cmd); 82 system(cmd);
75 # Compress the bed files in the tmp directory. 83 # Compress the bed files in the tmp directory.
76 tmp_gzipped_files = paste(tmp_dir, "*.bed.gz", sep="/"); 84 tmp_gzipped_files = paste(tmp_dir, "*.bed.gz", sep="/");
91 cmd = "gzip tmp.tar"; 99 cmd = "gzip tmp.tar";
92 system(cmd); 100 system(cmd);
93 # Move the tmp archive to the output directory. 101 # Move the tmp archive to the output directory.
94 to_path = paste(opt$output_files_path, "tmp.tar.gz", sep="/"); 102 to_path = paste(opt$output_files_path, "tmp.tar.gz", sep="/");
95 file.rename("tmp.tar.gz", to_path); 103 file.rename("tmp.tar.gz", to_path);
96 104 # Handle file names for display in the primary dataset if necessary.
97 if (!is.null(opt$chrom_bed_input) && !is.null(opt$chromosome_windows)) { 105 if (!is.null(opt$chrom_bed_input) && !is.null(opt$chromosome_windows)) {
98 # Renane opt$chrom_bed_input to be chromosomes.bed 106 # Renane opt$chrom_bed_input to be chromosomes.bed
99 # and make a copy of it in the output directory. 107 # and make a copy of it in the output directory.
100 to_path = paste(opt$output_files_path, "chromosomes.bed", sep="/"); 108 to_path = paste(opt$output_files_path, "chromosomes.bed", sep="/");
101 file.copy(opt$chrom_bed_input, to_path); 109 file.copy(opt$chrom_bed_input, to_path);
102 # Move chromosome_windows.txt to the output directory. 110 # Move chromosome_windows.txt to the output directory.
103 to_path = paste(opt$output_files_path, opt$chromosome_windows, sep="/"); 111 to_path = paste(opt$output_files_path, opt$chromosome_windows, sep="/");
104 file.rename(opt$chromosome_windows, to_path); 112 file.rename(opt$chromosome_windows, to_path);
105 } 113 }
106