Mercurial > repos > vimalkumarvelayudhan > riboseqr_wrapper

annotate riboseqr/ribosome_profile.py @ 1:4419c2f77b60

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