Mercurial > repos > pjbriggs > macs21

comparison macs21_wrapper.py @ 1:d0986d2be693 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Substantial reimplementation of internals, also renamed id and version.
author pjbriggs
date Thu, 29 Jan 2015 11:11:21 -0500
parents fdad0c8c0957
children 15889783e759
comparison
equal deleted inserted replaced
0:fdad0c8c0957 1:d0986d2be693
1 #purpose: macs2 python wrapper 1 #!/bin/env python
2 #author: Ziru Zhou 2 #
3 #date: November, 2012 3 # Galaxy wrapper to run MACS 2.1
4 #
5 # Completely rewritten from the original macs2 wrapped by Ziru Zhou
6 # taken from http://toolshed.g2.bx.psu.edu/view/modencode-dcc/macs2
4 7
5 import sys, subprocess, tempfile, shutil, glob, os, os.path, gzip 8 import sys
6 from galaxy import eggs 9 import os
7 import pkg_resources 10 import subprocess
8 pkg_resources.require( "simplejson" ) 11 import tempfile
9 import simplejson 12 import shutil
10 13
11 CHUNK_SIZE = 1024 14 def move_file(working_dir,name,destination):
15 """Move a file 'name' from 'working_dir' to 'destination'
12 16
13 #========================================================================================== 17 """
14 #functions 18 if destination is None:
15 #========================================================================================== 19 # Nothing to do
16 def gunzip_cat_glob_path( glob_path, target_filename, delete = False ): 20 return
17 out = open( target_filename, 'wb' ) 21 source = os.path.join(working_dir,name)
18 for filename in glob.glob( glob_path ): 22 if os.path.exists(source):
19 fh = gzip.open( filename, 'rb' ) 23 shutil.move(source,destination)
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 24
31 def xls_to_interval( xls_file, interval_file, header = None ): 25 def convert_xls_to_interval(xls_file,interval_file,header=None):
32 out = open( interval_file, 'wb' ) 26 """Convert MACS XLS file to interval
27
28 From the MACS readme: "Coordinates in XLS is 1-based which is different with
29 BED format."
30
31 However this function no longer performs any coordinate conversions, it
32 simply ensures that any blank or non-data lines are commented out
33
34 """
35 fp = open(interval_file,'wb')
33 if header: 36 if header:
34 out.write( '#%s\n' % header ) 37 fp.write('#%s\n' % header)
35 wrote_header = False 38 for line in open(xls_file):
36 #From macs readme: Coordinates in XLS is 1-based which is different with BED format. 39 # Keep all existing comment lines
37 #PJB: updated to keep original 'start' coordinate i.e. do not convert to BED 40 if line.startswith('#'):
38 for line in open( xls_file ): 41 fp.write(line)
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: 42 else:
50 fields = line.split( '\t' ) 43 # Split line into fields and test to see if
51 if len( fields ) > 1: 44 # the 'start' field is actually an integer
45 fields = line.split('\t')
46 if len(fields) > 1:
52 try: 47 try:
53 # Try to convert 'start' to int and shift 48 int(fields[1])
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: 49 except ValueError:
58 # Integer conversion failed so comment out 50 # Integer conversion failed so comment out
59 # "bad" line instead 51 # "bad" line instead
60 fields[0] = "#%s" % fields[0] 52 fields[0] = "#%s" % fields[0]
61 out.write( '\t'.join( fields ) ) 53 fp.write( '\t'.join( fields ) )
62 out.close() 54 fp.close()
63 55
64 #========================================================================================== 56 if __name__ == "__main__":
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 57
72 #================================================================================= 58 # Echo the command line
73 #parse options and execute macs2 59 print ' '.join(sys.argv)
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 60
81 #================================================================================= 61 # Initialise output files - values are set by reading from
82 if (options['command'] == "callpeak"): 62 # the command line supplied by the Galaxy wrapper
83 output_summits_bed = outputs['output_summits_bed_file'] 63 output_extra_html = None
84 output_extra_html = outputs['output_extra_file'] 64 output_extra_path = None
85 output_extra_path = outputs['output_extra_file_path'] 65 output_broadpeaks = None
86 output_peaks = outputs['output_peaks_file'] 66 output_gappedpeaks = None
87 output_narrowpeaks = outputs['output_narrowpeaks_file'] 67 output_narrowpeaks = None
88 output_broadpeaks = outputs['output_broadpeaks_file'] 68 output_treat_pileup = None
89 output_gappedpeaks = outputs['output_gappedpeaks_file'] 69 output_lambda_bedgraph = None
90 output_xls_to_interval_peaks_file = outputs['output_xls_to_interval_peaks_file'] 70 output_xls_to_interval_peaks_file = None
91 output_treat_pileup = outputs['output_treat_pileup_file'] 71 output_peaks = None
92 output_lambda_bedgraph = outputs['output_lambda_bedgraph_file'] 72 output_bdgcmp = None
93 73
94 cmdline = "%s --format='%s' --name='%s' --gsize='%s' --bw='%s'" % ( cmdline, options['format'], experiment_name, options['gsize'], options['bw'] ) 74 # Build the MACS 2.1 command line
95 if (options['broad'] == 'broad'): 75 # Initial arguments are always the same: command & input ChIP-seq file name
96 cmdline = "%s --broad --broad-cutoff='%s'" % ( cmdline, options['broad_cutoff'] ) 76 cmdline = ["macs2 %s -t %s" % (sys.argv[1],sys.argv[2])]
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 77
109 cmdline = "%s -m %s -p %s -o bdgcmp_out.bdg" % ( cmdline, options['m'], options['pseudocount'] ) 78 # Process remaining args
110 #================================================================================= 79 for arg in sys.argv[3:]:
80 if arg.startswith('--format='):
81 # Convert format to uppercase
82 format_ = arg.split('=')[1].upper()
83 cmdline.append("--format=%s" % format_)
84 elif arg.startswith('--name='):
85 # Replace whitespace in name with underscores
86 experiment_name = '_'.join(arg.split('=')[1].split())
87 cmdline.append("--name=%s" % experiment_name)
88 elif arg.startswith('--output-'):
89 # Handle destinations for output files
90 arg0,filen = arg.split('=')
91 if arg0 == '--output-summits':
92 output_summits = filen
93 elif arg0 == '--output-extra-files':
94 output_extra_html = filen
95 elif arg0 == '--output-extra-files-path':
96 output_extra_path = filen
97 elif arg0 == '--output-broadpeaks':
98 output_broadpeaks = filen
99 elif arg0 == '--output-gappedpeaks':
100 output_gappedpeaks = filen
101 elif arg0 == '--output-narrowpeaks':
102 output_narrowpeaks = filen
103 elif arg0 == '--output-pileup':
104 output_treat_pileup = filen
105 elif arg0 == '--output-lambda-bedgraph':
106 output_lambda_bedgraph = filen
107 elif arg0 == '--output-xls-to-interval':
108 output_xls_to_interval_peaks_file = filen
109 elif arg0 == '--output-peaks':
110 output_peaks = filen
111 else:
112 # Pass remaining args directly to MACS
113 # command line
114 cmdline.append(arg)
115
116 cmdline = ' '.join(cmdline)
117 print "Generated command line:\n%s" % cmdline
111 118
112 tmp_dir = tempfile.mkdtemp() #macs makes very messy output, need to contain it into a temp dir, then provide to user 119 # Execute MACS2
113 stderr_name = tempfile.NamedTemporaryFile().name # redirect stderr here, macs provides lots of info via stderr, make it into a report 120 #
114 proc = subprocess.Popen( args=cmdline, shell=True, cwd=tmp_dir, stderr=open( stderr_name, 'wb' ) ) 121 # Make a working directory
122 working_dir = tempfile.mkdtemp()
123 #
124 # Collect stderr in a file for reporting later
125 stderr_filen = tempfile.NamedTemporaryFile().name
126 #
127 # Run MACS2
128 proc = subprocess.Popen(args=cmdline,shell=True,cwd=working_dir,
129 stderr=open(stderr_filen,'wb'))
115 proc.wait() 130 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 131
127 #================================================================================= 132 # Run R script to create PDF from model script
128 #copy files created by macs2 to appripriate directory with the provided names 133 if os.path.exists(os.path.join(working_dir,"%s_model.r" % experiment_name)):
129 #================================================================================= 134 cmdline = 'R --vanilla --slave < "%s_model.r" > "%s_model.r.log"' % \
135 (experiment_name, experiment_name)
136 proc = subprocess.Popen(args=cmdline,shell=True,cwd=working_dir)
137 proc.wait()
130 138
131 #================================================================================= 139 # Convert XLS to interval, if requested
132 #move files generated by callpeak command 140 if output_xls_to_interval_peaks_file is not None:
133 if (options['command'] == "callpeak"): 141 peaks_xls_file = os.path.join(working_dir,'%s_peaks.xls' % experiment_name )
134 #run R to create pdf from model script 142 if os.path.exists(peaks_xls_file):
135 if os.path.exists( os.path.join( tmp_dir, "%s_model.r" % experiment_name ) ): 143 convert_xls_to_interval(peaks_xls_file,output_xls_to_interval_peaks_file,
136 cmdline = 'R --vanilla --slave < "%s_model.r" > "%s_model.r.log"' % ( experiment_name, experiment_name ) 144 header='peaks file')
137 proc = subprocess.Popen( args=cmdline, shell=True, cwd=tmp_dir )
138 proc.wait()
139 145
140 #move bed out to proper output file 146 # Move MACS2 output files from working dir to their final destinations
141 created_bed_name = os.path.join( tmp_dir, "%s_summits.bed" % experiment_name ) 147 move_file(working_dir,"%s_summits.bed" % experiment_name,output_summits)
142 if os.path.exists( created_bed_name ): 148 move_file(working_dir,"%s_peaks.xls" % experiment_name,output_peaks)
143 shutil.move( created_bed_name, output_summits_bed ) 149 move_file(working_dir,"%s_peaks.narrowPeak" % experiment_name,output_narrowpeaks)
150 move_file(working_dir,"%s_peaks.broadPeak" % experiment_name,output_broadpeaks)
151 move_file(working_dir,"%s_peaks.gappedPeak" % experiment_name,output_gappedpeaks)
152 move_file(working_dir,"%s_treat_pileup.bdg" % experiment_name,output_treat_pileup)
153 move_file(working_dir,"%s_control_lambda.bdg" % experiment_name,output_lambda_bedgraph)
154 move_file(working_dir,"bdgcmp_out.bdg",output_bdgcmp)
144 155
145 #OICR peak_xls file 156 # Move remaining file to the 'extra files' path and link from the HTML
146 created_peak_xls_file = os.path.join( tmp_dir, "%s_peaks.xls" % experiment_name ) 157 # file to allow user to access them from within Galaxy
147 if os.path.exists( created_peak_xls_file ): 158 html_file = open(output_extra_html,'wb')
148 # shutil.copy( created_peak_xls_file, os.path.join ( "/mnt/galaxyData/tmp/", "%s_peaks.xls" % ( os.path.basename(output_extra_path) ))) 159 html_file.write('<html><head><title>Additional output created by MACS (%s)</title></head><body><h3>Additional Files:</h3><p><ul>\n' % experiment_name)
149 shutil.copyfile( created_peak_xls_file, output_peaks ) 160 # Make the 'extra files' directory
161 os.mkdir(output_extra_path)
162 # Move the files
163 for filen in sorted(os.listdir(working_dir)):
164 shutil.move(os.path.join(working_dir,filen),
165 os.path.join(output_extra_path,filen))
166 html_file.write( '<li><a href="%s">%s</a></li>\n' % (filen,filen))
167 # All files moved, close out HTML
168 html_file.write( '</ul></p>\n' )
169 # Append any stderr output
170 html_file.write('<h3>Messages from MACS:</h3>\n<p><pre>%s</pre></p>\n' %
171 open(stderr_filen,'rb').read())
172 html_file.write('</body></html>\n')
173 html_file.close()
150 174
151 #peaks.encodepeaks (narrowpeaks) file 175 # Clean up the working directory and files
152 created_narrowpeak_file = os.path.join (tmp_dir, "%s_peaks.narrowPeak" % experiment_name ) 176 os.unlink(stderr_filen)
153 if os.path.exists( created_narrowpeak_file ): 177 os.rmdir(working_dir)
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()