|
1
|
1 #!/usr/bin/env Rscript
|
|
|
2
|
|
|
3 suppressPackageStartupMessages(library("data.table"))
|
|
|
4 suppressPackageStartupMessages(library("optparse"))
|
|
|
5
|
|
|
6 option_list <- list(
|
|
12
|
7 make_option(c("--chrom_bed_input"), action="store", dest="chrom_bed_input", defaul=NULL, help="Chromosome windows positions file"),
|
|
14
|
8 make_option(c("--exclude_bed_input"), action="store", dest="exclude_bed_input", defaul=NULL, help="File(s) containing regions to exclude"),
|
|
3
|
9 make_option(c("--chrom_len_file"), action="store", dest="chrom_len_file", default=NULL, help="Chromosome lengths file"),
|
|
5
|
10 make_option(c("--ideaspre_input_config"), action="store", dest="ideaspre_input_config", help="Preprocessing input config file"),
|
|
3
|
11 make_option(c("--output"), action="store", dest="output", help="Primary output dataset"),
|
|
|
12 make_option(c("--output_files_path"), action="store", dest="output_files_path", help="Primary output dataset extra files path"),
|
|
5
|
13 make_option(c("--chromosome_windows"), action="store", dest="chromosome_windows", default=NULL, help="Windows positions by chroms config file"),
|
|
3
|
14 make_option(c("--window_size"), action="store", dest="window_size", type="integer", default=NULL, help="Window size in base pairs")
|
|
1
|
15 )
|
|
|
16
|
|
|
17 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
|
|
|
18 args <- parse_args(parser, positional_arguments=TRUE)
|
|
|
19 opt <- args$options
|
|
|
20
|
|
32
|
21 run_cmd = function(cmd) {
|
|
|
22 rc = system(cmd);
|
|
|
23 if (rc != 0) {
|
|
|
24 quit(save="no", status=rc);
|
|
|
25 }
|
|
|
26 }
|
|
|
27
|
|
22
|
28 tmp_dir = "tmp";
|
|
27
|
29 cbi_file = "chrom_bed_input.bed";
|
|
22
|
30
|
|
27
|
31 if (is.null(opt$chrom_bed_input)) {
|
|
|
32 # Create a chromosome windows positions file
|
|
|
33 # using the received chromosome lengths file
|
|
|
34 # and the window size.
|
|
|
35 cmd = paste("bedtools makewindows -g", opt$chrom_len_file, "-w", opt$window_size, ">", cbi_file, sep=" ");
|
|
32
|
36 run_cmd(cmd);
|
|
27
|
37 } else {
|
|
|
38 if (!is.null(opt$exclude_bed_input)) {
|
|
|
39 # Copy the received chrom_bed_input
|
|
|
40 # since we will alter it.
|
|
|
41 file.copy(opt$chrom_bed_input, cbi_file);
|
|
|
42 } else {
|
|
|
43 cbi_file = opt$chrom_bed_input;
|
|
|
44 }
|
|
|
45 }
|
|
|
46 # Exclude regions if necessary.
|
|
|
47 if (!is.null(opt$exclude_bed_input)) {
|
|
|
48 exclude_bed_inputs = as.character(opt$exclude_bed_input);
|
|
|
49 exclude_bed_files = strsplit(exclude_bed_inputs, ",");
|
|
|
50 tmp_file = paste("tmp", cbi_file, sep="_");
|
|
|
51 for (exclude_bed_file in exclude_bed_files) {
|
|
|
52 cmd = paste("bedtools subtract -a", cbi_file, "-b", exclude_bed_file, ">", tmp_file, sep=" ");
|
|
32
|
53 run_cmd(cmd);
|
|
27
|
54 cmd = paste("mv", tmp_file, cbi_file, sep=" ");
|
|
32
|
55 run_cmd(cmd);
|
|
27
|
56 }
|
|
|
57 }
|
|
|
58 # Read the chromosome windows positions file
|
|
|
59 # to get the smallest window size in the file
|
|
|
60 # (i.e., the minimum of column 3 - column 2.
|
|
|
61 cbi = fread(cbi_file);
|
|
|
62 min_window_size = min(cbi[,3]-cbi[,2]);
|
|
5
|
63 # Read the ideaspre_input_config text file which has this format:
|
|
1
|
64 # "cell type name" "epigenetic factor name" "file path" "file name" "datatype"
|
|
5
|
65 ideaspre_input_config = as.matrix(read.table(opt$ideaspre_input_config));
|
|
27
|
66 # Process data to windows mean.
|
|
|
67 for (i in 1:dim(ideaspre_input_config)[1]) {
|
|
|
68 file_path = ideaspre_input_config[i, 3]
|
|
|
69 file_name = ideaspre_input_config[i, 4]
|
|
|
70 datatype = ideaspre_input_config[i, 5]
|
|
|
71 if (datatype == "bam") {
|
|
|
72 cmd = paste("samtools index", file_path);
|
|
32
|
73 run_cmd(cmd);
|
|
27
|
74 bigwig_file_name = paste(file_name, "bw", sep=".");
|
|
|
75 cmd = paste("bamCoverage --bam", file_path, "-o", bigwig_file_name, "--binSize", min_window_size);
|
|
32
|
76 run_cmd(cmd);
|
|
27
|
77 } else {
|
|
|
78 bigwig_file_name = file_path;
|
|
|
79 }
|
|
|
80 bed_file_name = paste(file_name, "bed", sep=".");
|
|
|
81 bed_file_path = paste("tmp", bed_file_name, sep="/");
|
|
|
82 cmd = paste("bigWigAverageOverBed", bigwig_file_name, opt$chrom_bed_input, "stdout | cut -f5 >", bed_file_path);
|
|
32
|
83 run_cmd(cmd);
|
|
27
|
84 cmd = paste("gzip -f", bed_file_path);
|
|
32
|
85 run_cmd(cmd);
|
|
14
|
86 }
|
|
9
|
87 # Create file1.txt.
|
|
5
|
88 cmd = paste("cut -d' '", opt$ideaspre_input_config, "-f1,2 > file1.txt", sep=" ");
|
|
32
|
89 run_cmd(cmd);
|
|
22
|
90 # Compress the bed files in the tmp directory.
|
|
|
91 tmp_gzipped_files = paste(tmp_dir, "*.bed.gz", sep="/");
|
|
|
92 # Create file2.txt.
|
|
|
93 cmd = paste("ls", tmp_gzipped_files, "> file2.txt", sep=" ");
|
|
32
|
94 run_cmd(cmd);
|
|
12
|
95 # Create IDEAS_input_config.txt with the format required by IDEAS.
|
|
|
96 ideas_input_config = "IDEAS_input_config.txt"
|
|
|
97 cmd = paste("paste -d' ' file1.txt file2.txt >", ideas_input_config, sep=" " );
|
|
32
|
98 run_cmd(cmd);
|
|
12
|
99 # Move IDEAS_input_config.txt to the output directory.
|
|
|
100 to_path = paste(opt$output_files_path, ideas_input_config, sep="/");
|
|
|
101 file.rename(ideas_input_config, to_path);
|
|
22
|
102 # Archive the tmp directory.
|
|
|
103 cmd = "tar -cvf tmp.tar tmp";
|
|
32
|
104 run_cmd(cmd);
|
|
25
|
105 # Compress the archive.
|
|
|
106 cmd = "gzip tmp.tar";
|
|
32
|
107 run_cmd(cmd);
|
|
22
|
108 # Move the tmp archive to the output directory.
|
|
25
|
109 to_path = paste(opt$output_files_path, "tmp.tar.gz", sep="/");
|
|
|
110 file.rename("tmp.tar.gz", to_path);
|
|
27
|
111 # Handle file names for display in the primary dataset if necessary.
|
|
29
|
112 to_path = paste(opt$output_files_path, "chromosomes.bed", sep="/");
|
|
30
|
113 if (is.null(opt$chrom_bed_input)) {
|
|
29
|
114 # Move cbi_file to the output directory,
|
|
|
115 # naming it chromosomes.bed.
|
|
|
116 file.rename(cbi_file, to_path);
|
|
|
117 } else {
|
|
|
118 # Copy opt$chrom_bed_input to the output
|
|
|
119 # directory, naming it chromosomes.bed.
|
|
15
|
120 file.copy(opt$chrom_bed_input, to_path);
|
|
29
|
121 }
|
|
|
122 if (!is.null(opt$chromosome_windows)) {
|
|
12
|
123 # Move chromosome_windows.txt to the output directory.
|
|
5
|
124 to_path = paste(opt$output_files_path, opt$chromosome_windows, sep="/");
|
|
|
125 file.rename(opt$chromosome_windows, to_path);
|
|
3
|
126 }
|
