Mercurial > repos > pjbriggs > macs21

comparison macs21_wrapper.py @ 2:00d73c812399 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Version 2.1.0-6: add sorting step in bigWig generation, and explicitly terminate tool on error from MACS2.
author pjbriggs
date Wed, 22 Mar 2017 11:36:07 -0400
parents 06cb587a5e87
children
comparison
equal deleted inserted replaced
1:02a01ea54722 2:00d73c812399
59 59
60 The protocol is: 60 The protocol is:
61 61
62 $ fetchChromSizes.sh mm9 > mm9.chrom.sizes 62 $ fetchChromSizes.sh mm9 > mm9.chrom.sizes
63 $ bedClip treat.bedgraph mm9.chrom.sizes treat.clipped 63 $ bedClip treat.bedgraph mm9.chrom.sizes treat.clipped
64 $ bedGraphToBigWig treat.clipped mm9.chrom.sizes treat.bw 64 $ bedSort treat.clipped treat.clipped.sorted
65 $ bedGraphToBigWig treat.clipped.sorted mm9.chrom.sizes treat.bw
65 66
66 Get the binaries from 67 Get the binaries from
67 http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/ 68 http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/
68 69
69 We skip the fetchChromSizes step if the 'chrom_sizes' 70 We skip the fetchChromSizes step if the 'chrom_sizes'
109 # Check that clipped file exists 110 # Check that clipped file exists
110 treat_clipped = os.path.join(working_dir,treat_clipped) 111 treat_clipped = os.path.join(working_dir,treat_clipped)
111 if not os.path.exists(treat_clipped): 112 if not os.path.exists(treat_clipped):
112 sys.stderr.write("Failed to create clipped bed file\n") 113 sys.stderr.write("Failed to create clipped bed file\n")
113 sys.exit(1) 114 sys.exit(1)
115 # Run bedSort
116 treat_clipped_sorted = "%s.sorted" % os.path.basename(treat_clipped)
117 cmd = "bedSort %s %s" % (treat_clipped,treat_clipped_sorted)
118 print "Running %s" % cmd
119 proc = subprocess.Popen(args=cmd,shell=True,cwd=working_dir)
120 proc.wait()
121 # Check that sorted file exists
122 treat_clipped_sorted = os.path.join(working_dir,treat_clipped_sorted)
123 if not os.path.exists(treat_clipped_sorted):
124 sys.stderr.write("Failed to create sorted clipped bed file\n")
125 sys.exit(1)
114 # Run bedGraphToBigWig 126 # Run bedGraphToBigWig
115 cmd = "bedGraphToBigWig %s %s %s" % (treat_clipped,chrom_sizes, 127 cmd = "bedGraphToBigWig %s %s %s" % (treat_clipped_sorted,
128 chrom_sizes,
116 bigwig_file) 129 bigwig_file)
117 print "Running %s" % cmd 130 print "Running %s" % cmd
118 proc = subprocess.Popen(args=cmd,shell=True,cwd=working_dir) 131 proc = subprocess.Popen(args=cmd,shell=True,cwd=working_dir)
119 proc.wait() 132 proc.wait()
120 # Clean up temporary chrom length file 133 # Clean up temporary chrom length file
204 # 217 #
205 # Run MACS2 218 # Run MACS2
206 proc = subprocess.Popen(args=cmdline,shell=True,cwd=working_dir, 219 proc = subprocess.Popen(args=cmdline,shell=True,cwd=working_dir,
207 stderr=open(stderr_filen,'wb')) 220 stderr=open(stderr_filen,'wb'))
208 proc.wait() 221 proc.wait()
222 exit_code = proc.returncode
223 if exit_code != 0:
224 sys.stderr.write(open(stderr_filen,'rb').read())
225 sys.exit(exit_code)
209 226
210 # Run R script to create PDF from model script 227 # Run R script to create PDF from model script
211 if os.path.exists(os.path.join(working_dir,"%s_model.r" % experiment_name)): 228 if os.path.exists(os.path.join(working_dir,"%s_model.r" % experiment_name)):
212 cmdline = 'R --vanilla --slave < "%s_model.r" > "%s_model.r.log"' % \ 229 cmdline = 'R --vanilla --slave < "%s_model.r" > "%s_model.r.log"' % \
213 (experiment_name, experiment_name) 230 (experiment_name, experiment_name)