Mercurial > repos > stemcellcommons > macs2
comparison macs2_wrapper.py @ 21:eed6c8dcef5d draft default tip
Add option to call broad peaks.
| author | stemcellcommons |
|---|---|
| date | Wed, 16 Apr 2014 15:30:35 -0400 |
| parents | db2805f9df66 |
| children |
comparison
equal
deleted
inserted
replaced
| 20:db2805f9df66 | 21:eed6c8dcef5d |
|---|---|
| 33 #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. |
| 34 for line in open( xls_file ): | 34 for line in open( xls_file ): |
| 35 #keep all existing comment lines | 35 #keep all existing comment lines |
| 36 if line.startswith( '#' ): | 36 if line.startswith( '#' ): |
| 37 out.write( line ) | 37 out.write( line ) |
| 38 #added for macs2 since there is an extra newline | 38 #added for macs2 since there is an extra newline |
| 39 elif line.startswith( '\n' ): | 39 elif line.startswith( '\n' ): |
| 40 out.write( line ) | 40 out.write( line ) |
| 41 elif not wrote_header: | 41 elif not wrote_header: |
| 42 out.write( '#%s' % line ) | 42 out.write( '#%s' % line ) |
| 43 print line | 43 print line |
| 44 wrote_header = True | 44 wrote_header = True |
| 45 else: | 45 else: |
| 46 fields = line.split( '\t' ) | 46 fields = line.split( '\t' ) |
| 47 if len( fields ) > 1: | 47 if len( fields ) > 1: |
| 48 fields[1] = str( int( fields[1] ) - 1 ) | 48 fields[1] = str( int( fields[1] ) - 1 ) |
| 66 if options['input_control']: | 66 if options['input_control']: |
| 67 cmdline = "%s -c %s" % ( cmdline, ",".join( options['input_control'] ) ) | 67 cmdline = "%s -c %s" % ( cmdline, ",".join( options['input_control'] ) ) |
| 68 | 68 |
| 69 #================================================================================= | 69 #================================================================================= |
| 70 if (options['command'] == "callpeak"): | 70 if (options['command'] == "callpeak"): |
| 71 output_bed = outputs['output_bed_file'] | 71 output_bed = outputs['output_bed_file'] |
| 72 output_extra_html = outputs['output_extra_file'] | 72 output_extra_html = outputs['output_extra_file'] |
| 73 output_extra_path = outputs['output_extra_file_path'] | 73 output_extra_path = outputs['output_extra_file_path'] |
| 74 output_peaks = outputs['output_peaks_file'] | 74 output_peaks = outputs['output_peaks_file'] |
| 75 output_narrowpeaks = outputs['output_narrowpeaks_file'] | 75 output_narrowpeaks = outputs['output_narrowpeaks_file'] |
| 76 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'] |
| 77 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'] |
| 78 | 78 |
| 79 if 'pvalue' in options: | 79 if 'pvalue' in options: |
| 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'] ) | 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'] ) |
| 81 elif 'qvalue' in options: | 81 elif 'qvalue' in options: |
| 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'] ) | 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'] ) |
| 83 | 83 |
| 84 if 'nomodel' in options: | 84 if 'broad_cutoff' in options: |
| 85 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'] ) | |
| 86 #================================================================================= | 89 #================================================================================= |
| 87 if (options['command'] == "bdgcmp"): | 90 if (options['command'] == "bdgcmp"): |
| 88 output_bdgcmp = outputs['output_bdgcmp_file'] | 91 output_bdgcmp = outputs['output_bdgcmp_file'] |
| 89 | 92 |
| 90 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'] ) |
| 91 #================================================================================= | 94 #================================================================================= |
| 92 | 95 |
| 93 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 |
| 94 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 |
| 95 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' ) ) |
| 110 #================================================================================= | 113 #================================================================================= |
| 111 | 114 |
| 112 #================================================================================= | 115 #================================================================================= |
| 113 #move files generated by callpeak command | 116 #move files generated by callpeak command |
| 114 if (options['command'] == "callpeak"): | 117 if (options['command'] == "callpeak"): |
| 115 #run R to create pdf from model script | 118 #run R to create pdf from model script |
| 116 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 ) ): |
| 117 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 ) |
| 118 proc = subprocess.Popen( args=cmdline, shell=True, cwd=tmp_dir ) | 121 proc = subprocess.Popen( args=cmdline, shell=True, cwd=tmp_dir ) |
| 119 proc.wait() | 122 proc.wait() |
| 120 | 123 |
| 121 #move bed out to proper output file | 124 #move bed out to proper output file |
| 122 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 ) |
| 123 if os.path.exists( created_bed_name ): | 126 if os.path.exists( created_bed_name ): |
| 124 shutil.move( created_bed_name, output_bed ) | 127 shutil.move( created_bed_name, output_bed ) |
| 125 | 128 |
| 126 #OICR peak_xls file | 129 #OICR peak_xls file |
| 127 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 ) |
| 128 if os.path.exists( created_peak_xls_file ): | 131 if os.path.exists( created_peak_xls_file ): |
| 129 # 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) ))) |
| 130 shutil.copyfile( created_peak_xls_file, output_peaks ) | 133 shutil.copyfile( created_peak_xls_file, output_peaks ) |
| 131 | 134 |
| 132 #peaks.encodepeaks (narrowpeaks) file | 135 #peaks.encodepeaks (narrowpeaks) file |
| 133 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 ) |
| 134 if os.path.exists( created_narrowpeak_file ): | 137 if os.path.exists( created_narrowpeak_file ): |
| 135 shutil.move (created_narrowpeak_file, output_narrowpeaks ) | 138 shutil.move (created_narrowpeak_file, output_narrowpeaks ) |
| 136 | 139 |
| 137 #parse xls files to interval files as needed | 140 #parse xls files to interval files as needed |
| 138 #if 'xls_to_interval' in options: | 141 #if 'xls_to_interval' in options: |
| 139 if (options['xls_to_interval'] == "True"): | 142 if (options['xls_to_interval'] == "True"): |
| 140 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 ) |
| 141 if os.path.exists( create_peak_xls_file ): | 144 if os.path.exists( create_peak_xls_file ): |
| 142 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' ) |
| 143 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 ) |
| 144 if os.path.exists( create_peak_xls_file ): | 147 if os.path.exists( create_peak_xls_file ): |
| 145 print "negative file exists" | 148 print "negative file exists" |
| 146 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' ) |
| 147 | 150 |
| 148 #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 |
| 149 out_html = open( output_extra_html, 'wb' ) | 152 out_html = open( output_extra_html, 'wb' ) |
| 150 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 ) |
| 151 os.mkdir( output_extra_path ) | 154 os.mkdir( output_extra_path ) |
| 152 for filename in sorted( os.listdir( tmp_dir ) ): | 155 for filename in sorted( os.listdir( tmp_dir ) ): |
| 153 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 ) ) |
| 154 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 ) ) |
| 155 #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 ) ) |
| 156 out_html.write( '</ul></p>\n' ) | 159 out_html.write( '</ul></p>\n' ) |
| 157 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() ) |
| 158 out_html.write( '</body></html>\n' ) | 161 out_html.write( '</body></html>\n' ) |
| 159 out_html.close() | 162 out_html.close() |
| 160 | 163 |
| 161 #================================================================================= | 164 #================================================================================= |
| 162 #move files generated by bdgcmp command | 165 #move files generated by bdgcmp command |
| 163 if (options['command'] == "bdgcmp"): | 166 if (options['command'] == "bdgcmp"): |
| 164 created_bdgcmp_file = os.path.join (tmp_dir, "bdgcmp_out.bdg" ) | 167 created_bdgcmp_file = os.path.join (tmp_dir, "bdgcmp_out.bdg" ) |
| 165 if os.path.exists( created_bdgcmp_file ): | 168 if os.path.exists( created_bdgcmp_file ): |
| 166 shutil.move (created_bdgcmp_file, output_bdgcmp ) | 169 shutil.move (created_bdgcmp_file, output_bdgcmp ) |
| 167 | 170 |
| 168 #================================================================================= | 171 #================================================================================= |
| 169 #cleanup | 172 #cleanup |
| 170 #================================================================================= | 173 #================================================================================= |
| 171 os.unlink( stderr_name ) | 174 os.unlink( stderr_name ) |
| 172 os.rmdir( tmp_dir ) | 175 os.rmdir( tmp_dir ) |