Mercurial > repos > stemcellcommons > macs2

comparison macs2_wrapper.py @ 1:c05f607d116c draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 )