annotate riboseqr/prepare.py @ 0:e01de823e919 draft default tip

Uploaded
author biotechcoder
date Fri, 01 May 2015 05:41:51 -0400
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
1 #!/usr/bin/env python
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
2 import os
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
3 import sys
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
4 import argparse
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
5 import logging
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
6 import rpy2.robjects as robjects
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
7
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
8 import utils
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
9
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
10 rscript = ''
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
11 R = robjects.r
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
12
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
13
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
14 def run_rscript(command=None):
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
15 """Run R command, log it, append to rscript"""
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
16 global rscript
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
17 if not command:
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
18 return
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
19 logging.debug(command)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
20 rscript += '{}\n'.format(command)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
21 R(command)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
22
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
23
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
24 def prep_riboseqr_input(sam_file, output_file):
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
25 """Generate input file for riboSeqR from SAM format file."""
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
26 with open(output_file, 'w') as f:
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
27 for line in open(sam_file):
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
28 if line.startswith('@'):
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
29 continue
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
30 line = line.split()
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
31 flag = line[1]
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
32
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
33 if flag != '0':
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
34 continue
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
35 # make start 0-indexed, sam alignments are 1-indexed
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
36 start = int(line[3]) - 1
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
37 (name, sequence) = (line[2], line[9])
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
38 f.write('"+"\t"{0}"\t{1}\t"{2}"\n'.format(name, start, sequence))
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
39
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
40
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
41 def batch_process(sam_files, seq_type, output_path):
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
42 """Batch process the conversion of SAM format files -> riboSeqR format
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
43 input files.
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
44
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
45 Files are saved with file names corresponding to their sequence type.
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
46
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
47 """
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
48 outputs = []
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
49 prefix = '{}'
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
50 if seq_type == 'riboseq':
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
51 prefix = 'RiboSeq file {}'
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
52 elif seq_type == 'rnaseq':
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
53 prefix = 'RNASeq file {}'
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
54
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
55 for count, fname in enumerate(sam_files):
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
56 count += 1
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
57 out_file = os.path.join(output_path, prefix.format(count))
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
58 logging.debug('Processing: {}'.format(fname))
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
59 logging.debug('Writing output to: {}'.format(out_file))
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
60 prep_riboseqr_input(fname, out_file)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
61 outputs.append(out_file)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
62 return outputs
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
63
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
64
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
65 def generate_ribodata(ribo_files='', rna_files='', replicate_names='',
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
66 seqnames='', rdata_save='Prepare.rda', sam_format=True,
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
67 html_file='Prepare-report.html', output_path=os.getcwd()):
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
68 """Prepares Ribo and RNA seq data in the format required for riboSeqR. Calls
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
69 the readRibodata function of riboSeqR and saves the result objects in an
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
70 R data file which can be used as input for the next step.
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
71
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
72 """
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
73 input_ribo_files = utils.process_args(ribo_files, ret_mode='list')
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
74 logging.debug('Found {} Ribo-Seq files'.format(len(input_ribo_files)))
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
75 logging.debug(input_ribo_files)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
76
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
77 input_rna_files = []
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
78 if rna_files:
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
79 input_rna_files = utils.process_args(rna_files, ret_mode='list')
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
80 logging.debug('Found {} RNA-Seq files'.format(len(input_rna_files)))
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
81 logging.debug(input_rna_files)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
82
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
83 replicates = utils.process_args(replicate_names, ret_mode='charvector')
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
84 logging.debug('Replicates: {}\n'.format(replicates))
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
85
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
86 if sam_format:
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
87 ribo_seq_files = batch_process(input_ribo_files, 'riboseq', output_path)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
88 else:
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
89 ribo_seq_files = input_ribo_files
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
90
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
91 html = '<h2>Prepare riboSeqR input - results</h2><hr>'
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
92 if len(ribo_seq_files):
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
93 html += '<h4>Generated riboSeqR format input files ' \
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
94 '<em>(RiboSeq)</em></h4><p>'
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
95 for fname in ribo_seq_files:
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
96 html += '<a href="{0}">{0}</a><br>'.format(
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
97 os.path.basename(fname))
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
98 html += '</p>'
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
99
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
100 rna_seq_files = []
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
101 if len(input_rna_files):
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
102 if sam_format:
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
103 rna_seq_files = batch_process(
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
104 input_rna_files, 'rnaseq', output_path)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
105 else:
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
106 rna_seq_files = input_rna_files
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
107
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
108 if len(rna_seq_files):
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
109 html += ('<h4>Generated riboSeqR format input files '
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
110 '<em>(RNASeq)</em></h4><p>')
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
111 for fname in rna_seq_files:
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
112 html += '<a href="{0}">{0}</a><br>'.format(
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
113 os.path.basename(fname))
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
114 html += '</p>'
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
115
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
116 input_seqnames = utils.process_args(seqnames, ret_mode='charvector')
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
117 options = {'ribo_seq_files': 'c({})'.format(str(ribo_seq_files)[1:-1]),
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
118 'rna_seq_files': 'c({})'.format(str(rna_seq_files)[1:-1]),
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
119 'input_replicates': replicates,
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
120 'input_seqnames': input_seqnames}
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
121
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
122 logging.debug('{}'.format(R('sessionInfo()')))
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
123
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
124 script = ''
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
125 cmd = 'library(riboSeqR)'
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
126 run_rscript(cmd)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
127 script += '{}\n'.format(cmd)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
128
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
129 if len(rna_seq_files):
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
130 cmd_args = ('riboFiles={ribo_seq_files}, '
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
131 'rnaFiles={rna_seq_files}'.format(**options))
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
132 else:
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
133 cmd_args = 'riboFiles={ribo_seq_files}'.format(**options)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
134
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
135 if input_seqnames:
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
136 cmd_args += ', seqnames={input_seqnames}'.format(**options)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
137 if replicates:
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
138 cmd_args += ', replicates={input_replicates}'.format(**options)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
139 else:
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
140 cmd_args += ', replicates=c("")'
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
141 cmd = 'riboDat <- readRibodata({0})'.format(cmd_args)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
142 run_rscript(cmd)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
143 script += '{}\n'.format(cmd)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
144
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
145 ribo_data = R['riboDat']
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
146 logging.debug('riboDat \n{}\n'.format(ribo_data))
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
147 cmd = 'save("riboDat", file="{}", compress=FALSE)'.format(rdata_save)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
148 run_rscript(cmd)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
149 script += '{}\n'.format(cmd)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
150
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
151 msg = '\n{:#^80}\n{}\n{:#^80}\n'.format(
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
152 ' R script for this session ', script, ' End R script ')
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
153 logging.debug(msg)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
154
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
155 with open(os.path.join(output_path, 'prepare.R'), 'w') as r:
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
156 r.write(script)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
157
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
158 html += ('<h4>R script for this session</h4>'
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
159 '<p><a href="prepare.R">prepare.R</a></p>'
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
160 '<p>Next step: <em>Triplet periodicity</em></p>')
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
161
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
162 with open(html_file, 'w') as f:
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
163 f.write(html)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
164
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
165 return ribo_data
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
166
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
167
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
168 if __name__ == '__main__':
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
169
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
170 description = (
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
171 'Prepare riboSeqR input file from SAM format RNA/Ribo-Seq alignment.')
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
172 parser = argparse.ArgumentParser(description=description)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
173
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
174 # required arguments
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
175 flags = parser.add_argument_group('required arguments')
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
176 flags.add_argument('--ribo_files', required=True,
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
177 help='List of Ribo-Seq files, comma-separated')
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
178 # optional arguments
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
179 parser.add_argument('--rna_files',
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
180 help='List of RNA-Seq files. Comma-separated')
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
181 parser.add_argument('--replicate_names',
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
182 help='Replicate names, comma-separated')
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
183 parser.add_argument('--seqnames',
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
184 help='Transcript (seqname) names to be read')
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
185 parser.add_argument('--rdata_save',
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
186 help='File to write RData to (default: %(default)s)',
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
187 default='Prepare.rda')
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
188 parser.add_argument(
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
189 '--sam_format',
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
190 help='Flag. Input is in SAM format', action='store_true')
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
191 parser.add_argument('--html_file', help='Output file for results (HTML)')
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
192 parser.add_argument('--output_path',
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
193 help='Files are saved in this directory')
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
194 parser.add_argument('--debug', help='Flag. Produce debug output',
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
195 action='store_true')
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
196 args = parser.parse_args()
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
197
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
198 if args.debug:
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
199 level = logging.DEBUG
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
200 else:
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
201 level = logging.INFO
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
202
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
203 logging.basicConfig(format='%(module)s: %(levelname)s - %(message)s',
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
204 level=level, stream=sys.stdout)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
205 logging.debug('Supplied Arguments: {}'.format(vars(args)))
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
206
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
207 if not os.path.exists(args.output_path):
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
208 os.mkdir(args.output_path)
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
209
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
210 generate_ribodata(
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
211 ribo_files=args.ribo_files, rna_files=args.rna_files,
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
212 replicate_names=args.replicate_names, seqnames=args.seqnames,
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
213 rdata_save=args.rdata_save,
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
214 sam_format=args.sam_format, html_file=args.html_file,
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
215 output_path=args.output_path
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
216 )
e01de823e919 Uploaded
biotechcoder
parents:
diff changeset
217 logging.debug('Done')