Mercurial > repos > vimalkumarvelayudhan > riboseqr_wrapper

annotate riboseqr/ribosome_profile.py @ 6:5a242f289347 default tip

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