Mercurial > repos > stemcellcommons > macs2
comparison 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 |
comparison
equal
deleted
inserted
replaced
0:642c0da30ca6 | 1:c05f607d116c |
---|---|
1 #purpose: macs2 python wrapper | 1 # macs2 python wrapper |
2 #author: Ziru Zhou | 2 # based on http://toolshed.g2.bx.psu.edu/view/modencode-dcc/macs2 |
3 #date: November, 2012 | |
4 | 3 |
5 import sys, subprocess, tempfile, shutil, glob, os, os.path, gzip | 4 import sys, subprocess, tempfile, shutil, glob, os, os.path, gzip |
6 from galaxy import eggs | 5 from galaxy import eggs |
7 import pkg_resources | 6 import json |
8 pkg_resources.require( "simplejson" ) | |
9 import simplejson | |
10 | 7 |
11 CHUNK_SIZE = 1024 | 8 CHUNK_SIZE = 1024 |
12 | 9 |
13 #========================================================================================== | 10 #========================================================================================== |
14 #functions | 11 #functions |
36 #From macs readme: Coordinates in XLS is 1-based which is different with BED format. | 33 #From macs readme: Coordinates in XLS is 1-based which is different with BED format. |
37 for line in open( xls_file ): | 34 for line in open( xls_file ): |
38 #keep all existing comment lines | 35 #keep all existing comment lines |
39 if line.startswith( '#' ): | 36 if line.startswith( '#' ): |
40 out.write( line ) | 37 out.write( line ) |
41 #added for macs2 since there is an extra newline | 38 #added for macs2 since there is an extra newline |
42 elif line.startswith( '\n' ): | 39 elif line.startswith( '\n' ): |
43 out.write( line ) | 40 out.write( line ) |
44 elif not wrote_header: | 41 elif not wrote_header: |
45 out.write( '#%s' % line ) | 42 out.write( '#%s' % line ) |
46 print line | 43 print line |
47 wrote_header = True | 44 wrote_header = True |
48 else: | 45 else: |
49 fields = line.split( '\t' ) | 46 fields = line.split( '\t' ) |
50 if len( fields ) > 1: | 47 if len( fields ) > 1: |
51 fields[1] = str( int( fields[1] ) - 1 ) | 48 fields[1] = str( int( fields[1] ) - 1 ) |
55 #========================================================================================== | 52 #========================================================================================== |
56 #main | 53 #main |
57 #========================================================================================== | 54 #========================================================================================== |
58 def main(): | 55 def main(): |
59 #take in options file and output file names | 56 #take in options file and output file names |
60 options = simplejson.load( open( sys.argv[1] ) ) | 57 options = json.load( open( sys.argv[1] ) ) |
61 outputs = simplejson.load( open( sys.argv[2] ) ) | 58 outputs = json.load( open( sys.argv[2] ) ) |
62 | 59 |
63 #================================================================================= | 60 #================================================================================= |
64 #parse options and execute macs2 | 61 #parse options and execute macs2 |
65 #================================================================================= | 62 #================================================================================= |
66 #default inputs that are in every major command | 63 #default inputs that are in every major command |
69 if options['input_control']: | 66 if options['input_control']: |
70 cmdline = "%s -c %s" % ( cmdline, ",".join( options['input_control'] ) ) | 67 cmdline = "%s -c %s" % ( cmdline, ",".join( options['input_control'] ) ) |
71 | 68 |
72 #================================================================================= | 69 #================================================================================= |
73 if (options['command'] == "callpeak"): | 70 if (options['command'] == "callpeak"): |
74 output_bed = outputs['output_bed_file'] | 71 output_bed = outputs['output_bed_file'] |
75 output_extra_html = outputs['output_extra_file'] | 72 output_extra_html = outputs['output_extra_file'] |
76 output_extra_path = outputs['output_extra_file_path'] | 73 output_extra_path = outputs['output_extra_file_path'] |
77 output_peaks = outputs['output_peaks_file'] | 74 output_peaks = outputs['output_peaks_file'] |
78 output_narrowpeaks = outputs['output_narrowpeaks_file'] | 75 output_narrowpeaks = outputs['output_narrowpeaks_file'] |
79 output_xls_to_interval_peaks_file = outputs['output_xls_to_interval_peaks_file'] | 76 output_xls_to_interval_peaks_file = outputs['output_xls_to_interval_peaks_file'] |
80 output_xls_to_interval_negative_peaks_file = outputs['output_xls_to_interval_negative_peaks_file'] | 77 output_xls_to_interval_negative_peaks_file = outputs['output_xls_to_interval_negative_peaks_file'] |
81 | 78 |
82 if 'pvalue' in options: | 79 if 'pvalue' in options: |
83 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'] ) | 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'] ) |
84 elif 'qvalue' in options: | 81 elif 'qvalue' in options: |
85 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'] ) | 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'] ) |
86 | 83 |
87 if 'nomodel' in options: | 84 if 'broad_cutoff' in options: |
88 cmdline = "%s --nomodel --shiftsize='%s'" % ( cmdline, options['nomodel'] ) | 85 cmdline += " --broad --broad-cutoff=%s" % (options['broad_cutoff']) |
86 | |
87 if 'nomodel' in options: | |
88 cmdline = "%s --nomodel --shiftsize='%s'" % ( cmdline, options['nomodel'] ) | |
89 #================================================================================= | 89 #================================================================================= |
90 if (options['command'] == "bdgcmp"): | 90 if (options['command'] == "bdgcmp"): |
91 output_bdgcmp = outputs['output_bdgcmp_file'] | 91 output_bdgcmp = outputs['output_bdgcmp_file'] |
92 | 92 |
93 cmdline = "%s -m %s -p %s -o bdgcmp_out.bdg" % ( cmdline, options['m'], options['pseudocount'] ) | 93 cmdline = "%s -m %s -p %s -o bdgcmp_out.bdg" % ( cmdline, options['m'], options['pseudocount'] ) |
94 #================================================================================= | 94 #================================================================================= |
95 | 95 |
96 tmp_dir = tempfile.mkdtemp() #macs makes very messy output, need to contain it into a temp dir, then provide to user | 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 | 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' ) ) | 98 proc = subprocess.Popen( args=cmdline, shell=True, cwd=tmp_dir, stderr=open( stderr_name, 'wb' ) ) |
113 #================================================================================= | 113 #================================================================================= |
114 | 114 |
115 #================================================================================= | 115 #================================================================================= |
116 #move files generated by callpeak command | 116 #move files generated by callpeak command |
117 if (options['command'] == "callpeak"): | 117 if (options['command'] == "callpeak"): |
118 #run R to create pdf from model script | 118 #run R to create pdf from model script |
119 if os.path.exists( os.path.join( tmp_dir, "%s_model.r" % experiment_name ) ): | 119 if os.path.exists( os.path.join( tmp_dir, "%s_model.r" % experiment_name ) ): |
120 cmdline = 'R --vanilla --slave < "%s_model.r" > "%s_model.r.log"' % ( experiment_name, experiment_name ) | 120 cmdline = 'R --vanilla --slave < "%s_model.r" > "%s_model.r.log"' % ( experiment_name, experiment_name ) |
121 proc = subprocess.Popen( args=cmdline, shell=True, cwd=tmp_dir ) | 121 proc = subprocess.Popen( args=cmdline, shell=True, cwd=tmp_dir ) |
122 proc.wait() | 122 proc.wait() |
123 | 123 |
124 #move bed out to proper output file | 124 #move bed out to proper output file |
125 created_bed_name = os.path.join( tmp_dir, "%s_peaks.bed" % experiment_name ) | 125 created_bed_name = os.path.join( tmp_dir, "%s_peaks.bed" % experiment_name ) |
126 if os.path.exists( created_bed_name ): | 126 if os.path.exists( created_bed_name ): |
127 shutil.move( created_bed_name, output_bed ) | 127 shutil.move( created_bed_name, output_bed ) |
128 | 128 |
129 #OICR peak_xls file | 129 #OICR peak_xls file |
130 created_peak_xls_file = os.path.join( tmp_dir, "%s_peaks.xls" % experiment_name ) | 130 created_peak_xls_file = os.path.join( tmp_dir, "%s_peaks.xls" % experiment_name ) |
131 if os.path.exists( created_peak_xls_file ): | 131 if os.path.exists( created_peak_xls_file ): |
132 # shutil.copy( created_peak_xls_file, os.path.join ( "/mnt/galaxyData/tmp/", "%s_peaks.xls" % ( os.path.basename(output_extra_path) ))) | 132 # shutil.copy( created_peak_xls_file, os.path.join ( "/mnt/galaxyData/tmp/", "%s_peaks.xls" % ( os.path.basename(output_extra_path) ))) |
133 shutil.copyfile( created_peak_xls_file, output_peaks ) | 133 shutil.copyfile( created_peak_xls_file, output_peaks ) |
134 | 134 |
135 #peaks.encodepeaks (narrowpeaks) file | 135 #peaks.encodepeaks (narrowpeaks) file |
136 created_narrowpeak_file = os.path.join (tmp_dir, "%s_peaks.encodePeak" % experiment_name ) | 136 created_narrowpeak_file = os.path.join (tmp_dir, "%s_peaks.encodePeak" % experiment_name ) |
137 if os.path.exists( created_narrowpeak_file ): | 137 if os.path.exists( created_narrowpeak_file ): |
138 shutil.move (created_narrowpeak_file, output_narrowpeaks ) | 138 shutil.move (created_narrowpeak_file, output_narrowpeaks ) |
139 | 139 |
140 #parse xls files to interval files as needed | 140 #parse xls files to interval files as needed |
141 #if 'xls_to_interval' in options: | 141 #if 'xls_to_interval' in options: |
142 if (options['xls_to_interval'] == "True"): | 142 if (options['xls_to_interval'] == "True"): |
143 create_peak_xls_file = os.path.join( tmp_dir, '%s_peaks.xls' % experiment_name ) | 143 create_peak_xls_file = os.path.join( tmp_dir, '%s_peaks.xls' % experiment_name ) |
144 if os.path.exists( create_peak_xls_file ): | 144 if os.path.exists( create_peak_xls_file ): |
145 xls_to_interval( create_peak_xls_file, output_xls_to_interval_peaks_file, header = 'peaks file' ) | 145 xls_to_interval( create_peak_xls_file, output_xls_to_interval_peaks_file, header = 'peaks file' ) |
146 create_peak_xls_file = os.path.join( tmp_dir, '%s_negative_peaks.xls' % experiment_name ) | 146 create_peak_xls_file = os.path.join( tmp_dir, '%s_negative_peaks.xls' % experiment_name ) |
147 if os.path.exists( create_peak_xls_file ): | 147 if os.path.exists( create_peak_xls_file ): |
148 print "negative file exists" | 148 print "negative file exists" |
149 xls_to_interval( create_peak_xls_file, output_xls_to_interval_negative_peaks_file, header = 'negative peaks file' ) | 149 xls_to_interval( create_peak_xls_file, output_xls_to_interval_negative_peaks_file, header = 'negative peaks file' ) |
150 | 150 |
151 #move all remaining files to extra files path of html file output to allow user download | 151 #move all remaining files to extra files path of html file output to allow user download |
152 out_html = open( output_extra_html, 'wb' ) | 152 out_html = open( output_extra_html, 'wb' ) |
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 ) | 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 ) |
154 os.mkdir( output_extra_path ) | 154 os.mkdir( output_extra_path ) |
155 for filename in sorted( os.listdir( tmp_dir ) ): | 155 for filename in sorted( os.listdir( tmp_dir ) ): |
156 shutil.move( os.path.join( tmp_dir, filename ), os.path.join( output_extra_path, filename ) ) | 156 shutil.move( os.path.join( tmp_dir, filename ), os.path.join( output_extra_path, filename ) ) |
157 out_html.write( '<li><a href="%s">%s</a></li>\n' % ( filename, filename ) ) | 157 out_html.write( '<li><a href="%s">%s</a></li>\n' % ( filename, filename ) ) |
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 ) ) | 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 ) ) |
159 out_html.write( '</ul></p>\n' ) | 159 out_html.write( '</ul></p>\n' ) |
160 out_html.write( '<h3>Messages from MACS:</h3>\n<p><pre>%s</pre></p>\n' % open( stderr_name, 'rb' ).read() ) | 160 out_html.write( '<h3>Messages from MACS:</h3>\n<p><pre>%s</pre></p>\n' % open( stderr_name, 'rb' ).read() ) |
161 out_html.write( '</body></html>\n' ) | 161 out_html.write( '</body></html>\n' ) |
162 out_html.close() | 162 out_html.close() |
163 | 163 |
164 #================================================================================= | 164 #================================================================================= |
165 #move files generated by bdgcmp command | 165 #move files generated by bdgcmp command |
166 if (options['command'] == "bdgcmp"): | 166 if (options['command'] == "bdgcmp"): |
167 created_bdgcmp_file = os.path.join (tmp_dir, "bdgcmp_out.bdg" ) | 167 created_bdgcmp_file = os.path.join (tmp_dir, "bdgcmp_out.bdg" ) |
168 if os.path.exists( created_bdgcmp_file ): | 168 if os.path.exists( created_bdgcmp_file ): |
169 shutil.move (created_bdgcmp_file, output_bdgcmp ) | 169 shutil.move (created_bdgcmp_file, output_bdgcmp ) |
170 | 170 |
171 #================================================================================= | 171 #================================================================================= |
172 #cleanup | 172 #cleanup |
173 #================================================================================= | 173 #================================================================================= |
174 os.unlink( stderr_name ) | 174 os.unlink( stderr_name ) |
175 os.rmdir( tmp_dir ) | 175 os.rmdir( tmp_dir ) |