|
0
|
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!')
|
