|
2
|
1 import argparse
|
|
|
2 import os
|
|
|
3 import shutil
|
|
|
4 import subprocess
|
|
|
5 import sys
|
|
|
6 import tempfile
|
|
|
7
|
|
|
8 BUFF_SIZE = 1048576
|
|
|
9 DESIGN_FILE = 'design.tabular'
|
|
|
10
|
|
|
11 parser = argparse.ArgumentParser()
|
|
|
12 parser.add_argument('--all_events_table', dest='all_events_table', help='Output all events table file')
|
|
|
13 parser.add_argument('--alphascale', dest='alphascale', type=float, default=None, help='Alpha scaling factor')
|
|
|
14 parser.add_argument('--chrom_len_file', dest='chrom_len_file', help='File listing the lengths of all chromosomes')
|
|
|
15 parser.add_argument('--control', dest='control', default=None, help='Input control files and data formats')
|
|
|
16 parser.add_argument('--diffp', dest='diffp', type=float, default=None, help='Minimum p-value for reporting differential enrichment')
|
|
|
17 parser.add_argument('--edgerod', dest='edgerod', type=float, default=None, help='EdgeR over-dispersion parameter value')
|
|
|
18 parser.add_argument('--exclude', dest='exclude', default=None, help='File containing a set of regions to ignore during MultiGPS training')
|
|
|
19 parser.add_argument('--expt', dest='expt', default=None, help="Input expt files and data formats")
|
|
|
20 parser.add_argument('--eventsaretxt', dest='eventsaretxt', default=None, help='Append a .txt extension to the events file for browser rendering')
|
|
|
21 parser.add_argument('--fixedalpha', dest='fixedalpha', type=int, default=None, help='Impose this alpha')
|
|
|
22 parser.add_argument('--fixedmodelrange', dest='fixedmodelrange', default=None, help='Keep binding model range fixed to inital size')
|
|
|
23 parser.add_argument('--fixedpb', dest='fixedpb', type=int, default=None, help='Fixed per-base limit')
|
|
|
24 parser.add_argument('--fixedscaling', dest='fixedscaling', type=float, default=None, help='Multiply control counts by total tag count ratio and then by this factor')
|
|
|
25 parser.add_argument('--format', dest='format', default=None, help='Input expt file data format')
|
|
|
26 parser.add_argument('--gaussmodelsmoothing', dest='gaussmodelsmoothing', default=None, help='Use Gaussian model smoothing')
|
|
|
27 parser.add_argument('--gausssmoothparam', dest='gausssmoothparam', type=int, default=None, help='Smoothing factor')
|
|
4
|
28 parser.add_argument('--input_item', dest='input_items', action='append', nargs=8, default=None, help="Input files, attributes and options")
|
|
2
|
29 parser.add_argument('--jointinmodel', dest='jointinmodel', default=None, help='Allow joint events in model updates')
|
|
|
30 parser.add_argument('--mappability', dest='mappability', type=float, default=None, help='Fraction of the genome that is mappable for these experiments')
|
|
|
31 parser.add_argument('--maxtrainingrounds', dest='maxtrainingrounds', type=int, default=None, help='Maximum number of training rounds for updating binding event read distributions')
|
|
|
32 parser.add_argument('--medianscale', dest='medianscale', default=None, help='Use the median signal/control ratio as the scaling factor')
|
|
|
33 parser.add_argument('--meme1proc', dest='meme1proc', default=None, help='Do not run the parallel version of meme')
|
|
|
34 parser.add_argument('--mememaxw', dest='mememaxw', type=int, default=None, help='Maximum motif width for MEME')
|
|
|
35 parser.add_argument('--mememinw', dest='mememinw', type=int, default=None, help='Minimum motif width for MEME')
|
|
|
36 parser.add_argument('--memenmotifs', dest='memenmotifs', type=int, default=None, help='Number of motifs MEME should find for each condition')
|
|
|
37 parser.add_argument('--minfold', dest='minfold', type=float, default=None, help='Minimum event fold-change vs scaled control')
|
|
|
38 parser.add_argument('--minqvalue', dest='minqvalue', type=float, default=None, help='Minimum Q-value (corrected p-value) of reported binding events')
|
|
|
39 parser.add_argument('--minmodelupdateevents', dest='minmodelupdateevents', type=int, default=None, help='Minimum number of events to support an update of the read distribution')
|
|
|
40 parser.add_argument('--mlconfignotshared', dest='mlconfignotshared', default=None, help='Share component configs in the ML step')
|
|
|
41 parser.add_argument('--nocache', dest='nocache', default=None, help='Turn off caching of the entire set of experiments')
|
|
|
42 parser.add_argument('--nodifftests', dest='nodifftests', default=None, help='Run differential enrichment tests')
|
|
|
43 parser.add_argument('--nomodelsmoothing', dest='nomodelsmoothing', default=None, help='Perform binding model smoothing')
|
|
|
44 parser.add_argument('--nomodelupdate', dest='nomodelupdate', default=None, help='Perform binding model updates')
|
|
|
45 parser.add_argument('--nomotifprior', dest='nomotifprior', default=None, help='Perform motif-finding only')
|
|
|
46 parser.add_argument('--nomotifs', dest='nomotifs', default=None, help='Perform motif-finding and motif priors')
|
|
|
47 parser.add_argument('--nonunique', dest='nonunique', default=None, help='Use non-unique reads')
|
|
|
48 parser.add_argument('--noposprior', dest='noposprior', default=None, help='Perform inter-experiment positional prior')
|
|
|
49 parser.add_argument('--noscaling', dest='noscaling', default=None, help='Do not use signal vs control scaling')
|
|
|
50 parser.add_argument('--output_bed', dest='output_bed', help='Output bed results file')
|
|
|
51 parser.add_argument('--output_html', dest='output_html', help='Output html results file')
|
|
|
52 parser.add_argument('--output_html_files_path', dest='output_html_files_path', help='Output html extra files')
|
|
|
53 parser.add_argument('--plotscaling', dest='plotscaling', default=None, help='Plot diagnostic information for the chosen scaling method')
|
|
|
54 parser.add_argument('--poissongausspb', dest='poissongausspb', type=int, default=None, help='Poisson threshold for filtering per base')
|
|
|
55 parser.add_argument('--prlogconf', dest='prlogconf', type=int, default=None, help='Poisson log threshold for potential region scanning')
|
|
|
56 parser.add_argument('--probshared', dest='probshared', type=float, default=None, help='Probability that events are shared across conditions')
|
|
|
57 parser.add_argument('--readdistributionfile', dest='readdistributionfile', default=None, help='Optional binding event read distribution file for initializing models')
|
|
|
58 parser.add_argument('--regressionscale', dest='regressionscale', default=None, help='Use scaling by regression on binned tag counts')
|
|
|
59 parser.add_argument('--replicates_counts', dest='replicates_counts', help='Output replicates counts file')
|
|
|
60 parser.add_argument('--scalewin', dest='scalewin', type=int, default=None, help='Window size for estimating scaling ratios')
|
|
|
61 parser.add_argument('--seq', dest='seq', default=None, help='Reference genome path')
|
|
|
62 parser.add_argument('--sesscale', dest='sesscale', default=None, help='Estimate scaling factor by SES')
|
|
|
63 parser.add_argument('--splinesmoothparam', dest='splinesmoothparam', type=int, default=None, help='Spline smoothing parameter')
|
|
|
64 parser.add_argument('--threads', dest='threads', type=int, default=4, help='The number of threads to run')
|
|
|
65 args = parser.parse_args()
|
|
|
66
|
|
|
67
|
|
|
68 def generate_design_file(input_items):
|
|
4
|
69 control_file_names = []
|
|
2
|
70 design_file = open(DESIGN_FILE, 'w')
|
|
|
71 for item in input_items:
|
|
4
|
72 experiment_file_name, control_file_name, label, file_format, condition_name, replicate_name, experiment_type, fixed_read_count = item
|
|
|
73 if control_file_name not in ['None', None, '']:
|
|
|
74 control_file_names.append(control_file_name)
|
|
2
|
75 file_format = file_format.upper()
|
|
4
|
76 items = [experiment_file_name, label, file_format, condition_name]
|
|
2
|
77 if replicate_name not in ['None', None, '']:
|
|
|
78 items.append(replicate_name)
|
|
|
79 if experiment_type not in ['None', None, '']:
|
|
|
80 items.append(experiment_type)
|
|
|
81 if fixed_read_count not in ['None', None, '']:
|
|
|
82 items.append(fixed_read_count)
|
|
|
83 design_file.write('%s\n' % '\t'.join(items))
|
|
|
84 design_file.close()
|
|
4
|
85 return control_file_names
|
|
2
|
86
|
|
|
87
|
|
|
88 def get_file_with_extension(dir, ext):
|
|
|
89 file_list = [f for f in os.listdir(dir) if f.endswith(ext)]
|
|
|
90 if len(file_list) == 1:
|
|
|
91 return file_list[0]
|
|
|
92 stop_err('Error running MultiGPS: output file with extension "%s" not generated.' % ext)
|
|
|
93
|
|
|
94
|
|
|
95 def get_stderr_exception(tmp_err, tmp_stderr, tmp_out, tmp_stdout, include_stdout=False):
|
|
|
96 tmp_stderr.close()
|
|
|
97 # Get stderr, allowing for case where it's very large.
|
|
|
98 tmp_stderr = open(tmp_err, 'rb')
|
|
|
99 stderr_str = ''
|
|
|
100 buffsize = BUFF_SIZE
|
|
|
101 try:
|
|
|
102 while True:
|
|
|
103 stderr_str += tmp_stderr.read(buffsize)
|
|
|
104 if not stderr_str or len(stderr_str) % buffsize != 0:
|
|
|
105 break
|
|
|
106 except OverflowError:
|
|
|
107 pass
|
|
|
108 tmp_stderr.close()
|
|
|
109 if include_stdout:
|
|
|
110 tmp_stdout = open(tmp_out, 'rb')
|
|
|
111 stdout_str = ''
|
|
|
112 buffsize = BUFF_SIZE
|
|
|
113 try:
|
|
|
114 while True:
|
|
|
115 stdout_str += tmp_stdout.read(buffsize)
|
|
|
116 if not stdout_str or len(stdout_str) % buffsize != 0:
|
|
|
117 break
|
|
|
118 except OverflowError:
|
|
|
119 pass
|
|
|
120 tmp_stdout.close()
|
|
|
121 if include_stdout:
|
|
|
122 return 'STDOUT\n%s\n\nSTDERR\n%s\n' % (stdout_str, stderr_str)
|
|
|
123 return stderr_str
|
|
|
124
|
|
|
125
|
|
|
126 def stop_err(msg):
|
|
|
127 sys.stderr.write(msg)
|
|
|
128 sys.exit(1)
|
|
|
129
|
|
|
130
|
|
|
131 # Preparation.
|
|
|
132 tmp_dir = tempfile.mkdtemp(prefix="tmp-multigps-")
|
|
|
133 os.makedirs(args.output_html_files_path)
|
|
|
134 # Build the command line.
|
|
|
135 cmd = 'multigps'
|
|
|
136 # General options
|
|
|
137 cmd += ' --threads %s' % args.threads
|
|
|
138 if args.eventsaretxt is not None:
|
|
|
139 # Append .txt extensions to events hrefs
|
|
|
140 # in output dataset so files will render
|
|
|
141 # in the browser.
|
|
|
142 cmd += ' --eventsaretxt'
|
|
|
143 if args.meme1proc is not None:
|
|
|
144 # Do not run the parallel version of meme.
|
|
|
145 cmd += ' --meme1proc'
|
|
|
146 # Experiment.
|
|
|
147 if args.input_items is not None:
|
|
4
|
148 control_file_names = generate_design_file(args.input_items)
|
|
2
|
149 cmd += ' --design %s' % DESIGN_FILE
|
|
4
|
150 if len(control_file_names) > 0:
|
|
|
151 cmd += ' --ctrl %s' % ','.join(control_file_names)
|
|
2
|
152 else:
|
|
|
153 if args.expt is not None:
|
|
|
154 cmd += ' --expt %s' % args.expt
|
|
|
155 if args.format is not None:
|
|
|
156 cmd += ' --format %s' % args.format
|
|
|
157 if args.control is not None:
|
|
|
158 cmd += ' --ctrl %s' % args.control
|
|
|
159 cmd += ' --geninfo %s' % args.chrom_len_file
|
|
|
160 # Advanced options.
|
|
|
161 if args.seq is not None:
|
|
|
162 cmd += ' --seq %s' % args.seq
|
|
|
163 # Limits on how many reads
|
|
|
164 if args.fixedpb is not None:
|
|
|
165 cmd += ' --fixedpb %d' % args.fixedpb
|
|
|
166 if args.poissongausspb is not None:
|
|
|
167 cmd += ' --poissongausspb %d' % args.poissongausspb
|
|
|
168 if args.nonunique is not None:
|
|
|
169 cmd += ' --nonunique'
|
|
|
170 if args.mappability is not None:
|
|
|
171 cmd += ' --mappability %4f' % args.mappability
|
|
|
172 if args.nocache is not None:
|
|
|
173 cmd += ' --nocache'
|
|
|
174 # Scaling data.'
|
|
|
175 if args.noscaling is not None:
|
|
|
176 cmd += ' --noscaling %s' % args.noscaling
|
|
|
177 if args.medianscale is not None:
|
|
|
178 cmd += ' --medianscale %s' % args.medianscale
|
|
|
179 if args.regressionscale is not None:
|
|
|
180 cmd += ' --regressionscale %s' % args.regressionscale
|
|
|
181 if args.sesscale is not None:
|
|
|
182 cmd += ' --sesscale %s' % args.sesscale
|
|
|
183 if args.fixedscaling is not None:
|
|
|
184 cmd += ' --fixedscaling %4f' % args.fixedscaling
|
|
|
185 if args.scalewin is not None:
|
|
|
186 cmd += ' --scalewin %d' % args.scalewin
|
|
|
187 if args.plotscaling is not None:
|
|
|
188 cmd += ' --plotscaling %s' % args.plotscaling
|
|
|
189 # Running MultiGPS.
|
|
|
190 if args.readdistributionfile is not None:
|
|
|
191 cmd += ' --d %s' % args.readdistributionfile
|
|
|
192 if args.maxtrainingrounds is not None:
|
|
|
193 cmd += ' --r %s' % args.maxtrainingrounds
|
|
|
194 if args.nomodelupdate is not None:
|
|
|
195 cmd += ' --nomodelupdate'
|
|
|
196 if args.minmodelupdateevents is not None:
|
|
|
197 cmd += ' --minmodelupdateevents %d' % args.minmodelupdateevents
|
|
|
198 if args.nomodelsmoothing is not None:
|
|
|
199 cmd += ' --nomodelsmoothing'
|
|
|
200 if args.splinesmoothparam is not None:
|
|
|
201 cmd += ' --splinesmoothparam %d' % args.splinesmoothparam
|
|
|
202 if args.gaussmodelsmoothing is not None:
|
|
|
203 cmd += ' --gaussmodelsmoothing'
|
|
|
204 if args.gausssmoothparam is not None:
|
|
|
205 cmd += ' --gausssmoothparam %s' % args.gausssmoothparam
|
|
|
206 if args.jointinmodel is not None:
|
|
|
207 cmd += ' --jointinmodel'
|
|
|
208 if args.fixedmodelrange is not None:
|
|
|
209 cmd += ' --fixedmodelrange'
|
|
|
210 if args.prlogconf is not None:
|
|
|
211 cmd += ' --prlogconf %d' % args.prlogconf
|
|
|
212 if args.fixedalpha is not None:
|
|
|
213 cmd += ' --fixedalpha %d' % args.fixedalpha
|
|
|
214 if args.alphascale is not None:
|
|
|
215 cmd += ' --alphascale %4f' % args.alphascale
|
|
|
216 if args.mlconfignotshared is not None:
|
|
|
217 cmd += ' --mlconfignotshared'
|
|
|
218 if args.exclude not in [None, 'None']:
|
|
|
219 cmd += ' --exclude %s' % args.exclude_file
|
|
|
220 # MultiGPS priors
|
|
|
221 if args.noposprior is not None:
|
|
|
222 cmd += ' --noposprior'
|
|
|
223 if args.probshared is not None:
|
|
|
224 cmd += ' --probshared %4f' % args.probshared
|
|
|
225 if args.memenmotifs is not None:
|
|
|
226 cmd += ' --memenmotifs %d' % args.memenmotifs
|
|
|
227 if args.mememinw is not None:
|
|
|
228 cmd += ' --mememinw %d' % args.mememinw
|
|
|
229 if args.mememaxw is not None:
|
|
|
230 cmd += ' --mememaxw %d' % args.mememaxw
|
|
|
231 if args.nomotifs is not None:
|
|
|
232 cmd += ' --nomotifs'
|
|
|
233 if args.nomotifprior is not None:
|
|
|
234 cmd += ' --nomotifprior'
|
|
|
235 # Reporting binding events
|
|
|
236 if args.minqvalue is not None:
|
|
|
237 cmd += ' --q %4f' % args.minqvalue
|
|
|
238 if args.minfold is not None:
|
|
|
239 cmd += ' --minfold %4f' % args.minfold
|
|
|
240 if args.nodifftests is not None:
|
|
|
241 cmd += ' --nodifftests'
|
|
|
242 if args.edgerod is not None:
|
|
|
243 cmd += ' --edgerod %4f' % args.edgerod
|
|
|
244 if args.diffp is not None:
|
|
|
245 cmd += ' --diffp %4f' % args.diffp
|
|
|
246 # Output directory.
|
|
|
247 cmd += ' --out %s' % args.output_html_files_path
|
|
|
248 # Define command response buffers.
|
|
|
249 tmp_out = tempfile.NamedTemporaryFile(dir=tmp_dir).name
|
|
|
250 tmp_stdout = open(tmp_out, 'wb')
|
|
|
251 tmp_err = tempfile.NamedTemporaryFile(dir=tmp_dir).name
|
|
|
252 tmp_stderr = open(tmp_err, 'wb')
|
|
|
253 tmp_filename = tempfile.NamedTemporaryFile(dir=tmp_dir).name
|
|
|
254 # Execute the command.
|
|
|
255 proc = subprocess.Popen(args=cmd, stderr=tmp_stderr, stdout=tmp_stdout, shell=True)
|
|
|
256 rc = proc.wait()
|
|
|
257 if rc != 0:
|
|
|
258 error_message = get_stderr_exception(tmp_err, tmp_stderr, tmp_out, tmp_stdout)
|
|
|
259 stop_err(error_message)
|
|
|
260 # Move each output file to the approapriate output dataset path.
|
|
|
261 output_bed = get_file_with_extension(args.output_html_files_path, 'bed')
|
|
|
262 shutil.move(os.path.join(args.output_html_files_path, output_bed), args.output_bed)
|
|
|
263 output_html = get_file_with_extension(args.output_html_files_path, 'html')
|
|
|
264 shutil.move(os.path.join(args.output_html_files_path, output_html), args.output_html)
|
|
|
265 replicates_counts = get_file_with_extension(args.output_html_files_path, 'counts')
|
|
|
266 shutil.move(os.path.join(args.output_html_files_path, replicates_counts), args.replicates_counts)
|
|
|
267 all_events_table = get_file_with_extension(args.output_html_files_path, 'table.txt')
|
|
|
268 shutil.move(os.path.join(args.output_html_files_path, all_events_table), args.all_events_table)
|
|
|
269 # Clean up.
|
|
|
270 if os.path.exists(tmp_dir):
|
|
|
271 shutil.rmtree(tmp_dir)
|
