Mercurial > repos > greg > ideas_preprocessor

comparison ideas_preprocessor.R @ 32:ea739695619d draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Fri, 16 Feb 2018 11:17:00 -0500
parents 0cb6ff8ba6df
children
comparison
equal deleted inserted replaced
31:c9f28cfea770 32:ea739695619d
16 16
17 parser <- OptionParser(usage="%prog [options] file", option_list=option_list) 17 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
18 args <- parse_args(parser, positional_arguments=TRUE) 18 args <- parse_args(parser, positional_arguments=TRUE)
19 opt <- args$options 19 opt <- args$options
20 20
21 run_cmd = function(cmd) {
22 rc = system(cmd);
23 if (rc != 0) {
24 quit(save="no", status=rc);
25 }
26 }
27
21 tmp_dir = "tmp"; 28 tmp_dir = "tmp";
22 cbi_file = "chrom_bed_input.bed"; 29 cbi_file = "chrom_bed_input.bed";
23 30
24 if (is.null(opt$chrom_bed_input)) { 31 if (is.null(opt$chrom_bed_input)) {
25 # Create a chromosome windows positions file 32 # Create a chromosome windows positions file
26 # using the received chromosome lengths file 33 # using the received chromosome lengths file
27 # and the window size. 34 # and the window size.
28 cmd = paste("bedtools makewindows -g", opt$chrom_len_file, "-w", opt$window_size, ">", cbi_file, sep=" "); 35 cmd = paste("bedtools makewindows -g", opt$chrom_len_file, "-w", opt$window_size, ">", cbi_file, sep=" ");
29 system(cmd); 36 run_cmd(cmd);
30 } else { 37 } else {
31 if (!is.null(opt$exclude_bed_input)) { 38 if (!is.null(opt$exclude_bed_input)) {
32 # Copy the received chrom_bed_input 39 # Copy the received chrom_bed_input
33 # since we will alter it. 40 # since we will alter it.
34 file.copy(opt$chrom_bed_input, cbi_file); 41 file.copy(opt$chrom_bed_input, cbi_file);
41 exclude_bed_inputs = as.character(opt$exclude_bed_input); 48 exclude_bed_inputs = as.character(opt$exclude_bed_input);
42 exclude_bed_files = strsplit(exclude_bed_inputs, ","); 49 exclude_bed_files = strsplit(exclude_bed_inputs, ",");
43 tmp_file = paste("tmp", cbi_file, sep="_"); 50 tmp_file = paste("tmp", cbi_file, sep="_");
44 for (exclude_bed_file in exclude_bed_files) { 51 for (exclude_bed_file in exclude_bed_files) {
45 cmd = paste("bedtools subtract -a", cbi_file, "-b", exclude_bed_file, ">", tmp_file, sep=" "); 52 cmd = paste("bedtools subtract -a", cbi_file, "-b", exclude_bed_file, ">", tmp_file, sep=" ");
46 system(cmd); 53 run_cmd(cmd);
47 cmd = paste("mv", tmp_file, cbi_file, sep=" "); 54 cmd = paste("mv", tmp_file, cbi_file, sep=" ");
48 system(cmd); 55 run_cmd(cmd);
49 } 56 }
50 } 57 }
51 # Read the chromosome windows positions file 58 # Read the chromosome windows positions file
52 # to get the smallest window size in the file 59 # to get the smallest window size in the file
53 # (i.e., the minimum of column 3 - column 2. 60 # (i.e., the minimum of column 3 - column 2.
61 file_path = ideaspre_input_config[i, 3] 68 file_path = ideaspre_input_config[i, 3]
62 file_name = ideaspre_input_config[i, 4] 69 file_name = ideaspre_input_config[i, 4]
63 datatype = ideaspre_input_config[i, 5] 70 datatype = ideaspre_input_config[i, 5]
64 if (datatype == "bam") { 71 if (datatype == "bam") {
65 cmd = paste("samtools index", file_path); 72 cmd = paste("samtools index", file_path);
66 system(cmd); 73 run_cmd(cmd);
67 bigwig_file_name = paste(file_name, "bw", sep="."); 74 bigwig_file_name = paste(file_name, "bw", sep=".");
68 cmd = paste("bamCoverage --bam", file_path, "-o", bigwig_file_name, "--binSize", min_window_size); 75 cmd = paste("bamCoverage --bam", file_path, "-o", bigwig_file_name, "--binSize", min_window_size);
69 system(cmd); 76 run_cmd(cmd);
70 } else { 77 } else {
71 bigwig_file_name = file_path; 78 bigwig_file_name = file_path;
72 } 79 }
73 bed_file_name = paste(file_name, "bed", sep="."); 80 bed_file_name = paste(file_name, "bed", sep=".");
74 bed_file_path = paste("tmp", bed_file_name, sep="/"); 81 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); 82 cmd = paste("bigWigAverageOverBed", bigwig_file_name, opt$chrom_bed_input, "stdout | cut -f5 >", bed_file_path);
76 system(cmd); 83 run_cmd(cmd);
77 cmd = paste("gzip -f", bed_file_path); 84 cmd = paste("gzip -f", bed_file_path);
78 system(cmd); 85 run_cmd(cmd);
79 } 86 }
80 # Create file1.txt. 87 # Create file1.txt.
81 cmd = paste("cut -d' '", opt$ideaspre_input_config, "-f1,2 > file1.txt", sep=" "); 88 cmd = paste("cut -d' '", opt$ideaspre_input_config, "-f1,2 > file1.txt", sep=" ");
82 system(cmd); 89 run_cmd(cmd);
83 # Compress the bed files in the tmp directory. 90 # Compress the bed files in the tmp directory.
84 tmp_gzipped_files = paste(tmp_dir, "*.bed.gz", sep="/"); 91 tmp_gzipped_files = paste(tmp_dir, "*.bed.gz", sep="/");
85 # Create file2.txt. 92 # Create file2.txt.
86 cmd = paste("ls", tmp_gzipped_files, "> file2.txt", sep=" "); 93 cmd = paste("ls", tmp_gzipped_files, "> file2.txt", sep=" ");
87 system(cmd); 94 run_cmd(cmd);
88 # Create IDEAS_input_config.txt with the format required by IDEAS. 95 # Create IDEAS_input_config.txt with the format required by IDEAS.
89 ideas_input_config = "IDEAS_input_config.txt" 96 ideas_input_config = "IDEAS_input_config.txt"
90 cmd = paste("paste -d' ' file1.txt file2.txt >", ideas_input_config, sep=" " ); 97 cmd = paste("paste -d' ' file1.txt file2.txt >", ideas_input_config, sep=" " );
91 system(cmd); 98 run_cmd(cmd);
92 # Move IDEAS_input_config.txt to the output directory. 99 # Move IDEAS_input_config.txt to the output directory.
93 to_path = paste(opt$output_files_path, ideas_input_config, sep="/"); 100 to_path = paste(opt$output_files_path, ideas_input_config, sep="/");
94 file.rename(ideas_input_config, to_path); 101 file.rename(ideas_input_config, to_path);
95 # Archive the tmp directory. 102 # Archive the tmp directory.
96 cmd = "tar -cvf tmp.tar tmp"; 103 cmd = "tar -cvf tmp.tar tmp";
97 system(cmd); 104 run_cmd(cmd);
98 # Compress the archive. 105 # Compress the archive.
99 cmd = "gzip tmp.tar"; 106 cmd = "gzip tmp.tar";
100 system(cmd); 107 run_cmd(cmd);
101 # Move the tmp archive to the output directory. 108 # Move the tmp archive to the output directory.
102 to_path = paste(opt$output_files_path, "tmp.tar.gz", sep="/"); 109 to_path = paste(opt$output_files_path, "tmp.tar.gz", sep="/");
103 file.rename("tmp.tar.gz", to_path); 110 file.rename("tmp.tar.gz", to_path);
104 # Handle file names for display in the primary dataset if necessary. 111 # Handle file names for display in the primary dataset if necessary.
105 to_path = paste(opt$output_files_path, "chromosomes.bed", sep="/"); 112 to_path = paste(opt$output_files_path, "chromosomes.bed", sep="/");