|
1
|
1 #!/usr/bin/env Rscript
|
|
|
2
|
|
14
|
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
|
|
1
|
14 suppressPackageStartupMessages(library("data.table"))
|
|
|
15 suppressPackageStartupMessages(library("optparse"))
|
|
|
16
|
|
|
17 option_list <- list(
|
|
12
|
18 make_option(c("--chrom_bed_input"), action="store", dest="chrom_bed_input", defaul=NULL, help="Chromosome windows positions file"),
|
|
14
|
19 make_option(c("--exclude_bed_input"), action="store", dest="exclude_bed_input", defaul=NULL, help="File(s) containing regions to exclude"),
|
|
3
|
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"),
|
|
5
|
22 make_option(c("--ideaspre_input_config"), action="store", dest="ideaspre_input_config", help="Preprocessing input config file"),
|
|
3
|
23 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"),
|
|
5
|
25 make_option(c("--output_hid"), action="store", dest="output_hid", help="Primary output dataset hid"),
|
|
3
|
26 make_option(c("--reads_per_bp"), action="store", dest="reads_per_bp", type="integer", help="Calculate the signal in each genomic window"),
|
|
|
27 make_option(c("--restrict_to_chroms"), action="store", dest="restrict_to_chroms", default=NULL, help="Restrict processing to specified chromosomes"),
|
|
|
28 make_option(c("--standardize_datasets"), action="store_true", dest="standardize_datasets", default=FALSE, help="Standardize all datasets"),
|
|
5
|
29 make_option(c("--chromosome_windows"), action="store", dest="chromosome_windows", default=NULL, help="Windows positions by chroms config file"),
|
|
3
|
30 make_option(c("--window_size"), action="store", dest="window_size", type="integer", default=NULL, help="Window size in base pairs")
|
|
1
|
31 )
|
|
|
32
|
|
|
33 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
|
|
|
34 args <- parse_args(parser, positional_arguments=TRUE)
|
|
|
35 opt <- args$options
|
|
|
36
|
|
5
|
37 create_primary_html = function(output, output_hid, output_files_path) {
|
|
|
38 files = list.files(path=output_files_path);
|
|
1
|
39 s <- paste('<html><head></head><body>', sep="\n");
|
|
5
|
40 s <- paste(s, '<h3>History item ', output_hid, ' files prepared for IDEAS</h3>\n', sep="");
|
|
1
|
41 s <- paste(s, '<ul>\n', sep="");
|
|
3
|
42 for (i in 1:length(files)) {
|
|
|
43 s <- paste(s, '<li><a href="', files[i], '">', files[i], '</a></li>\n', sep="");
|
|
1
|
44 }
|
|
|
45 s <- paste(s, '</ul>\n</body>\n</html>', sep="");
|
|
|
46 cat(s, file=output);
|
|
|
47 }
|
|
|
48
|
|
3
|
49 tmp_dir = "tmp";
|
|
|
50 output_tmp_dir = paste(opt$output_files_path, tmp_dir, sep="/");
|
|
|
51 dir.create(output_tmp_dir, showWarnings=FALSE);
|
|
1
|
52
|
|
5
|
53 # Read the ideaspre_input_config text file which has this format:
|
|
1
|
54 # "cell type name" "epigenetic factor name" "file path" "file name" "datatype"
|
|
5
|
55 ideaspre_input_config = as.matrix(read.table(opt$ideaspre_input_config));
|
|
14
|
56
|
|
|
57 # TODO: fix this
|
|
|
58 window_size = opt$window_size
|
|
|
59 if (is.null(opt$window_size)) {
|
|
|
60 window_size = 500;
|
|
|
61 }
|
|
1
|
62 # Process data to windows mean.
|
|
12
|
63 if (!is.null(opt$chrom_bed_input)) {
|
|
|
64 for (i in 1:dim(ideaspre_input_config)[1]) {
|
|
|
65 file_path = ideaspre_input_config[i, 3]
|
|
|
66 file_name = ideaspre_input_config[i, 4]
|
|
|
67 datatype = ideaspre_input_config[i, 5]
|
|
|
68 if (datatype == "bam") {
|
|
|
69 cmd = paste("samtools index", file_path);
|
|
|
70 system(cmd);
|
|
|
71 bigwig_file_name = paste(file_name, "bw", sep=".");
|
|
14
|
72
|
|
|
73 cmd = paste("bamCoverage --bam", file_path, "-o", bigwig_file_name, "--binSize", window_size);
|
|
12
|
74 system(cmd);
|
|
|
75 } else {
|
|
|
76 bigwig_file_name = file_path;
|
|
|
77 }
|
|
|
78 bed_file_name = paste(file_name, "bed", sep=".");
|
|
|
79 bed_file_path = paste("tmp", bed_file_name, sep="/");
|
|
|
80 cmd = paste("bigWigAverageOverBed", bigwig_file_name, opt$chrom_bed_input, "stdout | cut -f5 >", bed_file_path);
|
|
3
|
81 system(cmd);
|
|
12
|
82 cmd = paste("gzip -f", bed_file_path);
|
|
|
83 system(cmd);
|
|
1
|
84 }
|
|
|
85 }
|
|
|
86
|
|
9
|
87 # Create file1.txt.
|
|
5
|
88 cmd = paste("cut -d' '", opt$ideaspre_input_config, "-f1,2 > file1.txt", sep=" ");
|
|
3
|
89 system(cmd);
|
|
|
90 # Compress the bed files in the tmp directory.
|
|
|
91 tmp_gzipped_files = paste(tmp_dir, "*.bed.gz", sep="/");
|
|
9
|
92 # Create file2.txt.
|
|
3
|
93 cmd = paste("ls", tmp_gzipped_files, "> file2.txt", sep=" ");
|
|
|
94 system(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=" " );
|
|
3
|
98 system(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);
|
|
9
|
102 # Move the compressed bed files in the tmp
|
|
3
|
103 # directory to the output tmp directory.
|
|
|
104 tmp_files = list.files(path=tmp_dir);
|
|
|
105 for (i in 1:length(tmp_files)) {
|
|
|
106 from_path = paste(tmp_dir, tmp_files[i], sep="/");
|
|
|
107 to_path = paste(output_tmp_dir, tmp_files[i], sep="/");
|
|
|
108 file.rename(from_path, to_path);
|
|
|
109 }
|
|
12
|
110 if (!is.null(opt$chrom_bed_input) && !is.null(opt$chromosome_windows)) {
|
|
|
111 # Renane opt$chrom_bed_input to be chromosomes.bed
|
|
|
112 # and make a copy of it in the output directory.
|
|
|
113 to_path = paste(opt$output_files_path, "chromosomes.bed", sep="/");
|
|
15
|
114 file.copy(opt$chrom_bed_input, to_path);
|
|
12
|
115 # Move chromosome_windows.txt to the output directory.
|
|
5
|
116 to_path = paste(opt$output_files_path, opt$chromosome_windows, sep="/");
|
|
|
117 file.rename(opt$chromosome_windows, to_path);
|
|
3
|
118 }
|
|
1
|
119 # Create the primary HTML dataset.
|
|
9
|
120 create_primary_html(opt$output, opt$output_hid, opt$output_files_path);
|
