Mercurial > repos > pjbriggs > macs21

comparison macs21_wrapper.py @ 0:fdad0c8c0957 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded initial version to test toolshed.
author pjbriggs
date Wed, 21 Jan 2015 11:07:37 -0500
parents
children d0986d2be693
comparison
equal deleted inserted replaced
-1:000000000000 0:fdad0c8c0957
1 #purpose: macs2 python wrapper
2 #author: Ziru Zhou
3 #date: November, 2012
4
5 import sys, subprocess, tempfile, shutil, glob, os, os.path, gzip
6 from galaxy import eggs
7 import pkg_resources
8 pkg_resources.require( "simplejson" )
9 import simplejson
10
11 CHUNK_SIZE = 1024
12
13 #==========================================================================================
14 #functions
15 #==========================================================================================
16 def gunzip_cat_glob_path( glob_path, target_filename, delete = False ):
17 out = open( target_filename, 'wb' )
18 for filename in glob.glob( glob_path ):
19 fh = gzip.open( filename, 'rb' )
20 while True:
21 data = fh.read( CHUNK_SIZE )
22 if data:
23 out.write( data )
24 else:
25 break
26 fh.close()
27 if delete:
28 os.unlink( filename )
29 out.close()
30
31 def xls_to_interval( xls_file, interval_file, header = None ):
32 out = open( interval_file, 'wb' )
33 if header:
34 out.write( '#%s\n' % header )
35 wrote_header = False
36 #From macs readme: Coordinates in XLS is 1-based which is different with BED format.
37 #PJB: updated to keep original 'start' coordinate i.e. do not convert to BED
38 for line in open( xls_file ):
39 #keep all existing comment lines
40 if line.startswith( '#' ):
41 out.write( line )
42 #added for macs2 since there is an extra newline
43 ##elif line.startswith( '\n' ):
44 ## out.write( line )
45 elif not wrote_header:
46 out.write( '#%s' % line )
47 print line
48 wrote_header = True
49 else:
50 fields = line.split( '\t' )
51 if len( fields ) > 1:
52 try:
53 # Try to convert 'start' to int and shift
54 ##fields[1] = str( int( fields[1] ) - 1 )
55 # PJB don't correct, just test to see if it's an integer
56 int( fields[1] )
57 except ValueError:
58 # Integer conversion failed so comment out
59 # "bad" line instead
60 fields[0] = "#%s" % fields[0]
61 out.write( '\t'.join( fields ) )
62 out.close()
63
64 #==========================================================================================
65 #main
66 #==========================================================================================
67 def main():
68 #take in options file and output file names
69 options = simplejson.load( open( sys.argv[1] ) )
70 outputs = simplejson.load( open( sys.argv[2] ) )
71
72 #=================================================================================
73 #parse options and execute macs2
74 #=================================================================================
75 #default inputs that are in every major command
76 experiment_name = '_'.join( options['experiment_name'].split() ) #save experiment name here, it will be used by macs for some file names
77 cmdline = "macs2 %s -t %s" % ( options['command'], ",".join( options['input_chipseq'] ) )
78 if options['input_control']:
79 cmdline = "%s -c %s" % ( cmdline, ",".join( options['input_control'] ) )
80
81 #=================================================================================
82 if (options['command'] == "callpeak"):
83 output_summits_bed = outputs['output_summits_bed_file']
84 output_extra_html = outputs['output_extra_file']
85 output_extra_path = outputs['output_extra_file_path']
86 output_peaks = outputs['output_peaks_file']
87 output_narrowpeaks = outputs['output_narrowpeaks_file']
88 output_broadpeaks = outputs['output_broadpeaks_file']
89 output_gappedpeaks = outputs['output_gappedpeaks_file']
90 output_xls_to_interval_peaks_file = outputs['output_xls_to_interval_peaks_file']
91 output_treat_pileup = outputs['output_treat_pileup_file']
92 output_lambda_bedgraph = outputs['output_lambda_bedgraph_file']
93
94 cmdline = "%s --format='%s' --name='%s' --gsize='%s' --bw='%s'" % ( cmdline, options['format'], experiment_name, options['gsize'], options['bw'] )
95 if (options['broad'] == 'broad'):
96 cmdline = "%s --broad --broad-cutoff='%s'" % ( cmdline, options['broad_cutoff'] )
97 if 'pvalue' in options:
98 cmdline = "%s --pvalue='%s'" % ( cmdline, options['pvalue'] )
99 elif 'qvalue' in options:
100 cmdline = "%s --qvalue='%s'" % ( cmdline, options['qvalue'] )
101 cmdline = "%s --mfold %s %s %s %s %s %s --keep-dup %s" % (cmdline, options['mfoldlo'], options['mfoldhi'], options['nolambda'], options['bdg'], options['spmr'], options['call_summits'], options['keep_dup'] )
102
103 if 'nomodel' in options:
104 cmdline = "%s --nomodel --extsize='%s'" % ( cmdline, options['nomodel'] )
105 #=================================================================================
106 if (options['command'] == "bdgcmp"):
107 output_bdgcmp = outputs['output_bdgcmp_file']
108
109 cmdline = "%s -m %s -p %s -o bdgcmp_out.bdg" % ( cmdline, options['m'], options['pseudocount'] )
110 #=================================================================================
111
112 tmp_dir = tempfile.mkdtemp() #macs makes very messy output, need to contain it into a temp dir, then provide to user
113 stderr_name = tempfile.NamedTemporaryFile().name # redirect stderr here, macs provides lots of info via stderr, make it into a report
114 proc = subprocess.Popen( args=cmdline, shell=True, cwd=tmp_dir, stderr=open( stderr_name, 'wb' ) )
115 proc.wait()
116 #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
117 #Do not terminate if error code, allow dataset (e.g. log) creation and cleanup
118 if proc.returncode:
119 stderr_f = open( stderr_name )
120 while True:
121 chunk = stderr_f.read( CHUNK_SIZE )
122 if not chunk:
123 stderr_f.close()
124 break
125 sys.stderr.write( chunk )
126
127 #=================================================================================
128 #copy files created by macs2 to appripriate directory with the provided names
129 #=================================================================================
130
131 #=================================================================================
132 #move files generated by callpeak command
133 if (options['command'] == "callpeak"):
134 #run R to create pdf from model script
135 if os.path.exists( os.path.join( tmp_dir, "%s_model.r" % experiment_name ) ):
136 cmdline = 'R --vanilla --slave < "%s_model.r" > "%s_model.r.log"' % ( experiment_name, experiment_name )
137 proc = subprocess.Popen( args=cmdline, shell=True, cwd=tmp_dir )
138 proc.wait()
139
140 #move bed out to proper output file
141 created_bed_name = os.path.join( tmp_dir, "%s_summits.bed" % experiment_name )
142 if os.path.exists( created_bed_name ):
143 shutil.move( created_bed_name, output_summits_bed )
144
145 #OICR peak_xls file
146 created_peak_xls_file = os.path.join( tmp_dir, "%s_peaks.xls" % experiment_name )
147 if os.path.exists( created_peak_xls_file ):
148 # shutil.copy( created_peak_xls_file, os.path.join ( "/mnt/galaxyData/tmp/", "%s_peaks.xls" % ( os.path.basename(output_extra_path) )))
149 shutil.copyfile( created_peak_xls_file, output_peaks )
150
151 #peaks.encodepeaks (narrowpeaks) file
152 created_narrowpeak_file = os.path.join (tmp_dir, "%s_peaks.narrowPeak" % experiment_name )
153 if os.path.exists( created_narrowpeak_file ):
154 shutil.move (created_narrowpeak_file, output_narrowpeaks )
155
156 #peaks.broadpeaks file
157 created_broadpeak_file = os.path.join (tmp_dir, "%s_peaks.broadPeak" % experiment_name )
158 if os.path.exists( created_broadpeak_file ):
159 shutil.move (created_broadpeak_file, output_broadpeaks )
160
161 #peaks.gappedpeaks file
162 created_gappedpeak_file = os.path.join (tmp_dir, "%s_peaks.gappedPeak" % experiment_name )
163 if os.path.exists( created_gappedpeak_file ):
164 shutil.move (created_gappedpeak_file, output_gappedpeaks )
165
166 #parse xls files to interval files as needed
167 #if 'xls_to_interval' in options:
168 if (options['xls_to_interval'] == "True"):
169 create_peak_xls_file = os.path.join( tmp_dir, '%s_peaks.xls' % experiment_name )
170 if os.path.exists( create_peak_xls_file ):
171 xls_to_interval( create_peak_xls_file, output_xls_to_interval_peaks_file, header = 'peaks file' )
172
173 #bedgraph bedgraph files
174 if (options['bdg']):
175 created_treat_pileup_file = os.path.join (tmp_dir, "%s_treat_pileup.bdg" % experiment_name )
176 if os.path.exists( created_treat_pileup_file ):
177 shutil.move (created_treat_pileup_file, output_treat_pileup )
178 created_lambda_bedgraph_file = os.path.join (tmp_dir, "%s_control_lambda.bdg" % experiment_name )
179 if os.path.exists( created_lambda_bedgraph_file ):
180 shutil.move (created_lambda_bedgraph_file, output_lambda_bedgraph )
181
182 #move all remaining files to extra files path of html file output to allow user download
183 out_html = open( output_extra_html, 'wb' )
184 out_html.write( '<html><head><title>Additional output created by MACS (%s)</title></head><body><h3>Additional Files:</h3><p><ul>\n' % experiment_name )
185 os.mkdir( output_extra_path )
186 for filename in sorted( os.listdir( tmp_dir ) ):
187 shutil.move( os.path.join( tmp_dir, filename ), os.path.join( output_extra_path, filename ) )
188 out_html.write( '<li><a href="%s">%s</a></li>\n' % ( filename, filename ) )
189 #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 ) )
190 out_html.write( '</ul></p>\n' )
191 out_html.write( '<h3>Messages from MACS:</h3>\n<p><pre>%s</pre></p>\n' % open( stderr_name, 'rb' ).read() )
192 out_html.write( '</body></html>\n' )
193 out_html.close()
194
195 #=================================================================================
196 #move files generated by bdgcmp command
197 if (options['command'] == "bdgcmp"):
198 created_bdgcmp_file = os.path.join (tmp_dir, "bdgcmp_out.bdg" )
199 if os.path.exists( created_bdgcmp_file ):
200 shutil.move (created_bdgcmp_file, output_bdgcmp )
201
202 #=================================================================================
203 #cleanup
204 #=================================================================================
205 os.unlink( stderr_name )
206 os.rmdir( tmp_dir )
207
208 if __name__ == "__main__": main()