Mercurial > repos > vimalkumarvelayudhan > riboseqr_wrapper
comparison riboseqr/ribosome_profile.py @ 0:c34c364ce75d
First commit
| author | Vimalkumar Velayudhan <vimal@biotechcoder.com> |
|---|---|
| date | Tue, 21 Jul 2015 14:48:39 +0100 |
| parents | |
| children | 423ad61697c4 |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:c34c364ce75d |
|---|---|
| 1 #!/usr/bin/env python | |
| 2 import os | |
| 3 import sys | |
| 4 import glob | |
| 5 import argparse | |
| 6 import logging | |
| 7 import rpy2.robjects as robjects | |
| 8 import utils | |
| 9 | |
| 10 rscript = '' | |
| 11 R = robjects.r | |
| 12 | |
| 13 | |
| 14 def run_rscript(command=None): | |
| 15 """Run R command, log it, append to rscript""" | |
| 16 global rscript | |
| 17 if not command: | |
| 18 return | |
| 19 logging.debug(command) | |
| 20 rscript += '{}\n'.format(command) | |
| 21 msg = R(command) | |
| 22 | |
| 23 | |
| 24 def plot_transcript(rdata_load='Metagene.rda', transcript_name='', | |
| 25 transcript_length='27', transcript_cap='', | |
| 26 html_file='Plot-ribosome-profile.html', | |
| 27 output_path=os.getcwd()): | |
| 28 """Plot ribosome profile for a given transcript. """ | |
| 29 options = {} | |
| 30 for key, value, rtype, rmode in ( | |
| 31 ('transcript_name', transcript_name, 'str', None), | |
| 32 ('transcript_length', transcript_length, 'int', 'charvector'), | |
| 33 ('transcript_cap', transcript_cap, 'int', None)): | |
| 34 options[key] = utils.process_args(value, ret_type=rtype, ret_mode=rmode) | |
| 35 | |
| 36 run_rscript('suppressMessages(library(riboSeqR))') | |
| 37 run_rscript('load("{}")'.format(rdata_load)) | |
| 38 | |
| 39 html = """<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 3.2//EN"> | |
| 40 <html> | |
| 41 <head> | |
| 42 <title>Ribosome Profile Plot - Report</title> | |
| 43 </head> | |
| 44 <body> | |
| 45 """ | |
| 46 html += '<h2>Plot ribosome profile - results</h2>\n<hr>\n' | |
| 47 if len(transcript_name): | |
| 48 cmd_args = ( | |
| 49 '"{transcript_name}", main="{transcript_name}",' | |
| 50 'coordinates=ffCs@CDS, riboData=riboDat,' | |
| 51 'length={transcript_length}'.format(**options)) | |
| 52 if transcript_cap: | |
| 53 cmd_args += ', cap={transcript_cap}'.format(**options) | |
| 54 plot_file = os.path.join(output_path, 'Ribosome-profile-plot') | |
| 55 for fmat in ('pdf', 'png'): | |
| 56 if fmat == 'png': | |
| 57 cmd = 'png(file="{}_%1d.png", type="cairo")'.format(plot_file) | |
| 58 else: | |
| 59 cmd = 'pdf(file="{}.pdf")'.format(plot_file) | |
| 60 run_rscript(cmd) | |
| 61 cmd = 'plotTranscript({})'.format(cmd_args) | |
| 62 run_rscript(cmd) | |
| 63 run_rscript('dev.off()') | |
| 64 | |
| 65 html += ('<p>Selected ribosome footprint length: ' | |
| 66 '<strong>{0}</strong>\n'.format(transcript_length)) | |
| 67 | |
| 68 for image in sorted(glob.glob('{}_*.png'.format(plot_file))): | |
| 69 html += '<p><img border="1" src="{0}" alt="{0}"></p>\n'.format( | |
| 70 os.path.basename(image)) | |
| 71 html += '<p><a href="Ribosome-profile-plot.pdf">PDF version</a></p>\n' | |
| 72 else: | |
| 73 msg = 'No transcript name was provided. Did not generate plot.' | |
| 74 html += '<p>{}</p>'.format(msg) | |
| 75 logging.debug(msg) | |
| 76 | |
| 77 logging.debug('\n{:#^80}\n{}\n{:#^80}\n'.format( | |
| 78 ' R script for this session ', rscript, ' End R script ')) | |
| 79 | |
| 80 with open(os.path.join(output_path, 'ribosome-profile.R'), 'w') as r: | |
| 81 r.write(rscript) | |
| 82 | |
| 83 html += ('<h4>R script for this session</h4>\n' | |
| 84 '<p><a href="ribosome-profile.R">ribosome-profile.R</a></p>\n' | |
| 85 '</body>\n</html>\n') | |
| 86 | |
| 87 with open(html_file, 'w') as f: | |
| 88 f.write(html) | |
| 89 | |
| 90 | |
| 91 if __name__ == '__main__': | |
| 92 parser = argparse.ArgumentParser(description='Plot Ribosome profile') | |
| 93 | |
| 94 # required arguments | |
| 95 flags = parser.add_argument_group('required arguments') | |
| 96 flags.add_argument('--rdata_load', required=True, | |
| 97 help='Saved riboSeqR data from Step 2') | |
| 98 flags.add_argument('--transcript_name', required=True, | |
| 99 help='Name of the transcript to be plotted') | |
| 100 flags.add_argument( | |
| 101 '--transcript_length', required=True, | |
| 102 help='Size class of ribosome footprint data to be plotted', | |
| 103 default='27') | |
| 104 flags.add_argument( | |
| 105 '--transcript_cap', required=True, | |
| 106 help=('Cap on the largest value that will be plotted as an abundance ' | |
| 107 'of the ribosome footprint data')) | |
| 108 parser.add_argument('--html_file', help='HTML file with reports') | |
| 109 parser.add_argument('--output_path', help='Directory to save output files') | |
| 110 parser.add_argument('--debug', help='Produce debug output', | |
| 111 action='store_true') | |
| 112 | |
| 113 args = parser.parse_args() | |
| 114 if args.debug: | |
| 115 logging.basicConfig(format='%(levelname)s - %(message)s', | |
| 116 level=logging.DEBUG, stream=sys.stdout) | |
| 117 logging.debug('Supplied Arguments\n{}\n'.format(vars(args))) | |
| 118 | |
| 119 if not os.path.exists(args.output_path): | |
| 120 os.mkdir(args.output_path) | |
| 121 | |
| 122 plot_transcript(rdata_load=args.rdata_load, | |
| 123 transcript_name=args.transcript_name, | |
| 124 transcript_length=args.transcript_length, | |
| 125 transcript_cap=args.transcript_cap, | |
| 126 html_file=args.html_file, output_path=args.output_path) | |
| 127 logging.debug('Done!') |