Mercurial > repos > biotechcoder > riboseqr_wrapper

comparison riboseqr/ribosome_profile.py @ 0:e01de823e919 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author biotechcoder
date Fri, 01 May 2015 05:41:51 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:e01de823e919
1 #!/usr/bin/env python
2 # TODO: check first if argument parameters are of the right type if given.
3 import os
4 import sys
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 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('library(riboSeqR)')
37 run_rscript('load("{}")'.format(rdata_load))
38
39 html = '<h2>Plot ribosome profile - results</h2><hr>'
40 if len(transcript_name):
41 cmd_args = (
42 '"{transcript_name}", main="{transcript_name}",'
43 'coordinates=ffCs@CDS, riboData=riboDat,'
44 'length={transcript_length}'.format(**options))
45 if transcript_cap:
46 cmd_args += ', cap={transcript_cap}'.format(**options)
47 cmd = 'plotTranscript({})'.format(cmd_args)
48
49 for fmat in ('pdf', 'png'):
50 run_rscript('{0}(file="{1}")'.format(
51 fmat, os.path.join(
52 output_path, 'Ribosome-profile-plot.{}'.format(fmat))))
53 run_rscript(cmd)
54 run_rscript('dev.off()')
55
56 html += ('<p>Selected ribosome footprint length: <strong>{0}</strong>'
57 '</p><p><img border="1" src="Ribosome-profile-plot.png" '
58 'alt="Ribosome profile plot" /><br>'
59 '<a href="Ribosome-profile-plot.pdf">PDF version</a>'
60 '</p>'.format(transcript_length))
61 else:
62 msg = 'No transcript name was provided. Did not generate plot.'
63 html += '<p>{}</p>'.format(msg)
64 logging.debug(msg)
65
66 logging.info('\n{:#^80}\n{}\n{:#^80}\n'.format(
67 ' R script for this session ', rscript, ' End R script '))
68
69 with open(os.path.join(output_path, 'ribosome-profile.R'), 'w') as r:
70 r.write(rscript)
71
72 html += ('<h4>R script for this session</h4>'
73 '<p><a href="ribosome-profile.R">ribosome-profile.R</a></p>')
74
75 with open(html_file, 'w') as f:
76 f.write(html)
77
78
79 if __name__ == '__main__':
80 parser = argparse.ArgumentParser(description='Plot Ribosome profile')
81
82 # required arguments
83 flags = parser.add_argument_group('required arguments')
84 flags.add_argument('--rdata_load', required=True,
85 help='Saved riboSeqR data from Step 2')
86 flags.add_argument('--transcript_name', required=True,
87 help='Name of the transcript to be plotted')
88 flags.add_argument(
89 '--transcript_length', required=True,
90 help='Size class of ribosome footprint data to be plotted',
91 default='27')
92 flags.add_argument(
93 '--transcript_cap', required=True,
94 help=('Cap on the largest value that will be plotted as an abundance '
95 'of the ribosome footprint data'))
96 parser.add_argument('--html_file', help='HTML file with reports')
97 parser.add_argument('--output_path', help='Directory to save output files')
98 parser.add_argument('--debug', help='Produce debug output',
99 action='store_true')
100
101 args = parser.parse_args()
102 if args.debug:
103 level = logging.DEBUG
104 else:
105 level = logging.INFO
106
107 logging.basicConfig(format='%(levelname)s - %(message)s', level=level,
108 stream=sys.stdout)
109 logging.debug('Supplied Arguments\n{}\n'.format(vars(args)))
110
111 if not os.path.exists(args.output_path):
112 os.mkdir(args.output_path)
113
114 plot_transcript(rdata_load=args.rdata_load,
115 transcript_name=args.transcript_name,
116 transcript_length=args.transcript_length,
117 transcript_cap=args.transcript_cap,
118 html_file=args.html_file, output_path=args.output_path)
119 logging.debug('Done!')