Mercurial > repos > stemcellcommons > macs2
annotate macs2_wrapper.py @ 1:c05f607d116c draft default tip
Replace simplejson with json. Add option to call broad peaks.
author | stemcellcommons |
---|---|
date | Tue, 20 May 2014 12:10:38 -0400 |
parents | 642c0da30ca6 |
children |
rev | line source |
---|---|
1
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
1 # macs2 python wrapper |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
2 # based on http://toolshed.g2.bx.psu.edu/view/modencode-dcc/macs2 |
0 | 3 |
4 import sys, subprocess, tempfile, shutil, glob, os, os.path, gzip | |
5 from galaxy import eggs | |
1
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
6 import json |
0 | 7 |
8 CHUNK_SIZE = 1024 | |
9 | |
10 #========================================================================================== | |
11 #functions | |
12 #========================================================================================== | |
13 def gunzip_cat_glob_path( glob_path, target_filename, delete = False ): | |
14 out = open( target_filename, 'wb' ) | |
15 for filename in glob.glob( glob_path ): | |
16 fh = gzip.open( filename, 'rb' ) | |
17 while True: | |
18 data = fh.read( CHUNK_SIZE ) | |
19 if data: | |
20 out.write( data ) | |
21 else: | |
22 break | |
23 fh.close() | |
24 if delete: | |
25 os.unlink( filename ) | |
26 out.close() | |
27 | |
28 def xls_to_interval( xls_file, interval_file, header = None ): | |
29 out = open( interval_file, 'wb' ) | |
30 if header: | |
31 out.write( '#%s\n' % header ) | |
32 wrote_header = False | |
33 #From macs readme: Coordinates in XLS is 1-based which is different with BED format. | |
34 for line in open( xls_file ): | |
35 #keep all existing comment lines | |
36 if line.startswith( '#' ): | |
37 out.write( line ) | |
1
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
38 #added for macs2 since there is an extra newline |
0 | 39 elif line.startswith( '\n' ): |
40 out.write( line ) | |
41 elif not wrote_header: | |
42 out.write( '#%s' % line ) | |
1
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
43 print line |
0 | 44 wrote_header = True |
45 else: | |
46 fields = line.split( '\t' ) | |
47 if len( fields ) > 1: | |
48 fields[1] = str( int( fields[1] ) - 1 ) | |
49 out.write( '\t'.join( fields ) ) | |
50 out.close() | |
51 | |
52 #========================================================================================== | |
53 #main | |
54 #========================================================================================== | |
55 def main(): | |
56 #take in options file and output file names | |
1
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
57 options = json.load( open( sys.argv[1] ) ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
58 outputs = json.load( open( sys.argv[2] ) ) |
0 | 59 |
60 #================================================================================= | |
61 #parse options and execute macs2 | |
62 #================================================================================= | |
63 #default inputs that are in every major command | |
64 experiment_name = '_'.join( options['experiment_name'].split() ) #save experiment name here, it will be used by macs for some file names | |
65 cmdline = "macs2 %s -t %s" % ( options['command'], ",".join( options['input_chipseq'] ) ) | |
66 if options['input_control']: | |
67 cmdline = "%s -c %s" % ( cmdline, ",".join( options['input_control'] ) ) | |
68 | |
69 #================================================================================= | |
70 if (options['command'] == "callpeak"): | |
1
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
71 output_bed = outputs['output_bed_file'] |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
72 output_extra_html = outputs['output_extra_file'] |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
73 output_extra_path = outputs['output_extra_file_path'] |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
74 output_peaks = outputs['output_peaks_file'] |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
75 output_narrowpeaks = outputs['output_narrowpeaks_file'] |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
76 output_xls_to_interval_peaks_file = outputs['output_xls_to_interval_peaks_file'] |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
77 output_xls_to_interval_negative_peaks_file = outputs['output_xls_to_interval_negative_peaks_file'] |
0 | 78 |
1
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
79 if 'pvalue' in options: |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
80 cmdline = "%s --format='%s' --name='%s' --gsize='%s' --bw='%s' --pvalue='%s' --mfold %s %s %s %s" % ( cmdline, options['format'], experiment_name, options['gsize'], options['bw'], options['pvalue'], options['mfoldlo'], options['mfoldhi'], options['nolambda'], options['bdg'] ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
81 elif 'qvalue' in options: |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
82 cmdline = "%s --format='%s' --name='%s' --gsize='%s' --bw='%s' --qvalue='%s' --mfold %s %s %s %s" % ( cmdline, options['format'], experiment_name, options['gsize'], options['bw'], options['qvalue'], options['mfoldlo'], options['mfoldhi'], options['nolambda'], options['bdg'] ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
83 |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
84 if 'broad_cutoff' in options: |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
85 cmdline += " --broad --broad-cutoff=%s" % (options['broad_cutoff']) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
86 |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
87 if 'nomodel' in options: |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
88 cmdline = "%s --nomodel --shiftsize='%s'" % ( cmdline, options['nomodel'] ) |
0 | 89 #================================================================================= |
90 if (options['command'] == "bdgcmp"): | |
1
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
91 output_bdgcmp = outputs['output_bdgcmp_file'] |
0 | 92 |
1
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
93 cmdline = "%s -m %s -p %s -o bdgcmp_out.bdg" % ( cmdline, options['m'], options['pseudocount'] ) |
0 | 94 #================================================================================= |
95 | |
96 tmp_dir = tempfile.mkdtemp() #macs makes very messy output, need to contain it into a temp dir, then provide to user | |
97 stderr_name = tempfile.NamedTemporaryFile().name # redirect stderr here, macs provides lots of info via stderr, make it into a report | |
98 proc = subprocess.Popen( args=cmdline, shell=True, cwd=tmp_dir, stderr=open( stderr_name, 'wb' ) ) | |
99 proc.wait() | |
100 #We don't want to set tool run to error state if only warnings or info, e.g. mfold could be decreased to improve model, but let user view macs log | |
101 #Do not terminate if error code, allow dataset (e.g. log) creation and cleanup | |
102 if proc.returncode: | |
103 stderr_f = open( stderr_name ) | |
104 while True: | |
105 chunk = stderr_f.read( CHUNK_SIZE ) | |
106 if not chunk: | |
107 stderr_f.close() | |
108 break | |
109 sys.stderr.write( chunk ) | |
110 | |
111 #================================================================================= | |
112 #copy files created by macs2 to appripriate directory with the provided names | |
113 #================================================================================= | |
114 | |
115 #================================================================================= | |
116 #move files generated by callpeak command | |
117 if (options['command'] == "callpeak"): | |
1
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
118 #run R to create pdf from model script |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
119 if os.path.exists( os.path.join( tmp_dir, "%s_model.r" % experiment_name ) ): |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
120 cmdline = 'R --vanilla --slave < "%s_model.r" > "%s_model.r.log"' % ( experiment_name, experiment_name ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
121 proc = subprocess.Popen( args=cmdline, shell=True, cwd=tmp_dir ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
122 proc.wait() |
0 | 123 |
1
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
124 #move bed out to proper output file |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
125 created_bed_name = os.path.join( tmp_dir, "%s_peaks.bed" % experiment_name ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
126 if os.path.exists( created_bed_name ): |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
127 shutil.move( created_bed_name, output_bed ) |
0 | 128 |
1
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
129 #OICR peak_xls file |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
130 created_peak_xls_file = os.path.join( tmp_dir, "%s_peaks.xls" % experiment_name ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
131 if os.path.exists( created_peak_xls_file ): |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
132 # shutil.copy( created_peak_xls_file, os.path.join ( "/mnt/galaxyData/tmp/", "%s_peaks.xls" % ( os.path.basename(output_extra_path) ))) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
133 shutil.copyfile( created_peak_xls_file, output_peaks ) |
0 | 134 |
1
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
135 #peaks.encodepeaks (narrowpeaks) file |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
136 created_narrowpeak_file = os.path.join (tmp_dir, "%s_peaks.encodePeak" % experiment_name ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
137 if os.path.exists( created_narrowpeak_file ): |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
138 shutil.move (created_narrowpeak_file, output_narrowpeaks ) |
0 | 139 |
1
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
140 #parse xls files to interval files as needed |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
141 #if 'xls_to_interval' in options: |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
142 if (options['xls_to_interval'] == "True"): |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
143 create_peak_xls_file = os.path.join( tmp_dir, '%s_peaks.xls' % experiment_name ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
144 if os.path.exists( create_peak_xls_file ): |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
145 xls_to_interval( create_peak_xls_file, output_xls_to_interval_peaks_file, header = 'peaks file' ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
146 create_peak_xls_file = os.path.join( tmp_dir, '%s_negative_peaks.xls' % experiment_name ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
147 if os.path.exists( create_peak_xls_file ): |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
148 print "negative file exists" |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
149 xls_to_interval( create_peak_xls_file, output_xls_to_interval_negative_peaks_file, header = 'negative peaks file' ) |
0 | 150 |
1
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
151 #move all remaining files to extra files path of html file output to allow user download |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
152 out_html = open( output_extra_html, 'wb' ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
153 out_html.write( '<html><head><title>Additional output created by MACS (%s)</title></head><body><h3>Additional Files:</h3><p><ul>\n' % experiment_name ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
154 os.mkdir( output_extra_path ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
155 for filename in sorted( os.listdir( tmp_dir ) ): |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
156 shutil.move( os.path.join( tmp_dir, filename ), os.path.join( output_extra_path, filename ) ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
157 out_html.write( '<li><a href="%s">%s</a></li>\n' % ( filename, filename ) ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
158 #out_html.write( '<li><a href="%s">%s</a>peakxls %s SomethingDifferent tmp_dir %s path %s exp_name %s</li>\n' % ( created_peak_xls_file, filename, filename, tmp_dir, output_extra_path, experiment_name ) ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
159 out_html.write( '</ul></p>\n' ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
160 out_html.write( '<h3>Messages from MACS:</h3>\n<p><pre>%s</pre></p>\n' % open( stderr_name, 'rb' ).read() ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
161 out_html.write( '</body></html>\n' ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
162 out_html.close() |
0 | 163 |
164 #================================================================================= | |
165 #move files generated by bdgcmp command | |
166 if (options['command'] == "bdgcmp"): | |
1
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
167 created_bdgcmp_file = os.path.join (tmp_dir, "bdgcmp_out.bdg" ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
168 if os.path.exists( created_bdgcmp_file ): |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
169 shutil.move (created_bdgcmp_file, output_bdgcmp ) |
c05f607d116c
Replace simplejson with json. Add option to call broad peaks.
stemcellcommons
parents:
0
diff
changeset
|
170 |
0 | 171 #================================================================================= |
172 #cleanup | |
173 #================================================================================= | |
174 os.unlink( stderr_name ) | |
175 os.rmdir( tmp_dir ) | |
176 | |
177 if __name__ == "__main__": main() |