comparison cpt_dnaplotter.pl @ 1:d724f34e671d draft default tip

planemo upload commit 94b0cd1fff0826c6db3e7dc0c91c0c5a8be8bb0c
author cpt
date Mon, 05 Jun 2023 02:50:07 +0000
parents
children
comparison
equal deleted inserted replaced
0:e4de0a0e90c8 1:d724f34e671d
1 #!/usr/bin/perl
2 #
3 # Code written by Eric Rasche
4 # mailto:rasche.eric@yandex.ru
5 # tel: 404.692.2048
6 # http://eric.rasche.co.uk
7 # for
8 # Center for Phage Technology
9 #
10
11 use strict;
12 use warnings;
13
14 use CPT::GalaxyGetOpt;
15 use Data::Dumper;
16 use CPT::Bio;
17 use CPT::Circos::Conf;
18 my $bio = CPT::Bio->new();
19
20 my @argv_copy;
21 foreach(@ARGV){push(@argv_copy, "$_");}
22
23 my $ggo = CPT::GalaxyGetOpt->new();
24 my $options = $ggo->getOptions(
25 'options' => [
26 [ 'file|f', 'Input file', { validate => 'File/Input',
27 file_format => ['genbank', 'embl', 'txt'],
28 } ],
29 [],
30 ['Track Configuration'],
31 ['track_key', 'Key to select from genbank data', { validate => 'Genomic/Tag', required => 1, multiple => 1 } ],
32 ['track_feature_filter_invert', 'Should the qualifier search be inverted?', { validate => 'Option', options => { 'yes', 'Yes', 'no', 'No' } , multiple => 1 } ],
33 ['track_feature_filter_hastag', 'Select a tag which should be present in that qualifier (e.g., signal/tmhelix/pseudo)', { validate => 'String' , multiple => 1 } ],
34 ['track_feature_filter_textquery', 'Specify text which MUST be in that tag', { validate => 'String' , multiple => 1 } ],
35 ['track_feature_filter_strand', 'Which strand should the feature appear on?', { validate => 'Option', options => { 'f', 'Forward', 'r', 'Reverse', 'a', 'Any' } , multiple => 1 } ],
36 [],
37 ['enable_gc_skew_plot', 'Enable/Disable calculation of GC Skew Plot', { validate => 'Flag' } ],
38 ['gc_skew_plot_window_size', 'Window size for calculation of GC Skew', { validate => 'Int', min => 1000, default => 10000} ],
39 ['gc_skew_plot_step_size', 'Step size for calculation of GC Skew', { validate => 'Int', min => 200, default => 200 } ],
40 ],
41 'outputs' => [
42 [
43 'output_circos_confs',
44 'Circos Configuration Files',
45 {
46 validate => 'File/Output',
47 required => 1,
48 default => 'out',
49 data_format => 'archive',
50 default_format => 'zip',
51 }
52 ],
53 ],
54 'defaults' => [
55 'appid' => 'CircularDNAPlotter',
56 'appname' => 'Circos based DNAPlotter',
57 'appdesc' => 'plots genomes similar to Artemis\'s DNAPlotter',
58 'appvers' => '1.94.1',
59 ],
60 );
61
62 #perl cpt_dnaplotter.pl \
63 #-f ../t/test-files/moon.gbk \
64 #--track_key CDS --track_feature_filter_invert yes --track_feature_filter_hastag pseudo --track_feature_filter_strand f \
65 #--track_key CDS --track_feature_filter_invert yes --track_feature_filter_hastag pseudo --track_feature_filter_strand r \
66 #--track_key CDS --track_feature_filter_hastag pseudo --track_feature_filter_strand a \
67 #--track_key tRNA --track_feature_filter_strand a \
68 #--track_key CDS --track_feature_filter_hastag signal --track_feature_filter_strand a \
69 #--track_key CDS --track_feature_filter_hastag tmhelix --track_feature_filter_strand a
70
71
72 my @reorg_args = ();
73
74 my $cum_gc_ske_mean = 0;
75
76
77
78 my %latest = ();
79 for(my $i = 0; $i < scalar(@argv_copy); $i++){
80 my $c = $argv_copy[$i];
81 # We have entered a new one block
82 if($c eq '--track_key'){
83 # If we have loaded data
84 if(scalar(keys(%latest)) > 0){
85 my %copy;
86 foreach(keys(%latest)){
87 $copy{$_} = "" . $latest{$_};
88 }
89 push(@reorg_args, \%copy);
90 }
91
92 # Clean out latest to prep for new data
93 foreach(keys(%latest)){
94 delete $latest{$_};
95 }
96 }
97
98 if($c =~ /^--track_(.*)/){
99 $latest{$1} = $argv_copy[$i+1];
100 # Artificially bump so we don't bother looking at the answer to
101 # this question. We can "get away" with this because none of
102 # the options are flags. However, I've disabled it in the event
103 # that flags ARE introduced
104 #$i++;
105 }
106 }
107 push(@reorg_args, \%latest);
108 #$VAR1 = [
109 #{
110 #'feature_filter_invert' => 'yes',
111 #'feature_plot_color' => '005500',
112 #'feature_filter_strand' => 'f',
113 #'feature_filter_hastag' => 'pseudo',
114 #'key' => 'CDS'
115 #},
116 #{
117 #'feature_filter_strand' => 'a',
118 #'key' => 'RBS'
119 #}
120 #];
121 my @files = ();
122
123 my $number_of_tracks = 0;
124 sub register_track {
125 #my $r0 = ( 90 - (10 * $number_of_tracks - 1)/2) / 100;
126 #my $r1 = ( 90 - (10 * $number_of_tracks - 9)/2) / 100;
127 my $r0 = ( 90 - (10 * $number_of_tracks - 1)/1) / 100;
128 my $r1 = ( 90 - (10 * $number_of_tracks - 9)/1) / 100;
129 $number_of_tracks++;
130 return ($r0, $r1);
131 }
132
133 sub circosconf {
134 my $cc = CPT::Circos::Conf->new();
135 $cc->include('etc/colors_fonts_patterns.conf');
136 # Features to plot along the genome
137 $cc->include('ideogram.conf');
138 # markings indicating position along genome
139 $cc->include('ticks.conf');
140 # Genome data
141 $cc->set('karyotype', 'karyotype.conf');
142 # Default image params are fine
143 $cc->start_block('image');
144 $cc->include('etc/image.conf');
145 $cc->end_block();
146 # ???
147 $cc->set('chromosome_units', '1000');
148 $cc->set('chromosome_display_default', 'yes');
149 #$cc->include('highlights.conf');
150 $cc->include('plots.conf');
151
152 $cc->include('etc/housekeeping.conf');
153 my $result = $cc->finalize();
154 $cc = CPT::Circos::Conf->new();
155 return $result;
156 }
157 sub ideogramconf{
158 my $cc = CPT::Circos::Conf->new();
159 $cc->start_block('ideogram');
160 $cc->start_block('spacing');
161 $cc->set('default','0u');
162 $cc->set('break','0u');
163 $cc->end_block();
164
165 $cc->set('thickness', '20p');
166 $cc->set('stroke_thickness', '2');
167 $cc->set('stroke_color', 'black');
168 $cc->set('fill','no');
169 $cc->set('fill_color','black');
170 $cc->set('radius','0.85r');
171 $cc->set('show_label','yes');
172 $cc->set('label_font','default');
173 $cc->set('label_radius','dims(ideogram,radius) + 0.05');
174 $cc->set('label_size','36');
175 $cc->set('label_parallel','yes');
176 $cc->set('label_case','upper');
177
178 $cc->set('band_stroke_thickness','2');
179 $cc->set('show_bands','yes');
180 $cc->set('fill_bands','yes');
181 $cc->end_block();
182
183 return $cc->finalize();
184 }
185 sub generate_feature_table {
186 my ($filename, %filter) = @_;
187 print "Filtering on features\n";
188 print Dumper \%filter;
189 my $seqio_object = Bio::SeqIO->new(-file => $options->{file}, -format=>'genbank');
190 # Only handing first sequence.
191 my $seq_object = $seqio_object->next_seq;
192 my $parent = $seq_object->display_id();
193 # Feature data
194 my @features;
195 foreach my $feat($seq_object->get_SeqFeatures()){
196 if($feat->primary_tag() eq $filter{key}){
197 # If they've said "hastag" AND we do indeed have that tag AND we haven't inverted this filter.
198 if(
199 ($filter{feature_filter_hastag} && $feat->has_tag($filter{feature_filter_hastag}) && !$filter{feature_filter_invert})
200 ||
201 ($filter{feature_filter_hastag} && $filter{feature_filter_invert} && !$feat->has_tag($filter{feature_filter_hastag}))
202 ||
203 (! $filter{feature_filter_hastag})
204 ){
205 if(
206 !$filter{feature_filter_strand}
207 ||
208 ($feat->strand() == 1 && ( $filter{feature_filter_strand} eq 'f' || $filter{feature_filter_strand} eq 'a' ))
209 ||
210 ($feat->strand() == -1 && ( $filter{feature_filter_strand} eq 'r' || $filter{feature_filter_strand} eq 'a' ))
211 ||
212 ($feat->strand() == 0 && ( $filter{feature_filter_strand} eq 'a' ))
213 ){
214 push(@features, join(' ', $parent, $feat->start, $feat->end));
215 }
216 }
217 }
218 }
219 print "Found " . scalar @features . " features \n";
220 push(@files, [ 'data/' . $filename, join("\n", @features) ] );
221 }
222 sub plotsconf{
223 my $cc = CPT::Circos::Conf->new();
224
225 $cc->start_block('plots');
226
227
228 my $idx = 0;
229 foreach(@reorg_args){
230 my %filter = %{$_};
231 #{
232 #'feature_filter_invert' => 'yes',
233 #'feature_plot_color' => '005500',
234 #'feature_filter_strand' => 'f',
235 #'feature_filter_hastag' => 'pseudo',
236 #'key' => 'CDS'
237 #},
238 $idx++;
239 my $filename = sprintf('org.features.%s.txt', $idx);
240 generate_feature_table($filename, %filter);
241
242 my ($r0,$r1) = register_track();
243 $cc->start_block('plot');
244 $cc->set('type','tile');
245 $cc->set('file',$filename);
246 $cc->set('orientation', 'center');
247 $cc->set('thickness', '30');
248 $cc->set('r1', $r1 . 'r');# '0.78r');
249 $cc->set('r0', $r0 . 'r');# '0.72r');
250 $cc->set('layers','3');
251 $cc->set('layers_overflow','collapse');
252 $cc->set('color','paired-6-qual-' . $idx);
253 $cc->end_block();
254 }
255
256 if($options->{enable_gc_skew_plot}){
257 my ($r0,$r1) = register_track();
258 $cc->start_block('plot');
259 $cc->set('type','histogram');
260 $cc->set('file','gc.txt');
261 $cc->set('r1',$r1 . 'r');
262 $cc->set('r0',$r0 . 'r');
263 $cc->set('fill_color','purple');
264 $cc->set('orientation','in');
265 $cc->start_block('rules');
266 $cc->start_block('rule');
267 $cc->set('condition','var(value) < 0');
268 $cc->set('fill_color', 'green');
269 $cc->end_block();
270 $cc->end_block();
271 $cc->end_block();
272 }
273
274 #$cc->start_block('plot');
275 #$cc->set('type','histogram');
276 #$cc->set('file','gc_cumulative.txt');
277 #$cc->set('r1','0.6r');
278 #$cc->set('r0','0.55r');
279 #$cc->set('fill_color','purple');
280 #$cc->set('orientation','out');
281 #$cc->start_block('rules');
282 #$cc->start_block('rule');
283 #$cc->set('condition','var(value) < ' . $cum_gc_ske_mean);
284 #$cc->set('fill_color', 'green');
285 #$cc->end_block();
286 #$cc->end_block();
287 #$cc->end_block();
288
289 $cc->end_block();
290 return $cc->finalize();
291 }
292 sub ticksconf{
293 my $cc = CPT::Circos::Conf->new();
294
295 $cc->set('show_ticks','yes');
296 $cc->set('show_tick_labels','yes');
297 $cc->start_block('ticks');
298 $cc->set('radius','1r');
299 $cc->set('color','black');
300 $cc->set('thickness','2p');
301 $cc->set('multiplier','1e-3');
302 $cc->set('format','%d');
303
304 $cc->start_block('tick');
305 $cc->set('spacing','1000u');
306 $cc->set('size','10p');
307 $cc->end_block();
308
309 $cc->start_block('tick');
310 $cc->set('spacing','10000u');
311 $cc->set('size','15p');
312 $cc->set('show_label','yes');
313 $cc->set('label_size','20p');
314 $cc->set('label_offset','10p');
315 $cc->set('format','%d');
316 $cc->end_block();
317 $cc->end_block();
318 return $cc->finalize();
319 }
320 sub karyotype {
321 my @karyotype_data = ();
322 my $seqio_object = Bio::SeqIO->new(-file => $options->{file}, -format=>'genbank');
323 # Only handing first sequence.
324 my $seq_object = $seqio_object->next_seq;
325 # Main 'chromosome' data
326 push(@karyotype_data, join(' ', 'chr', '-',$seq_object->display_id(),$seq_object->display_id(),0, $seq_object->length(),'black'));
327
328 return join("\n", @karyotype_data);
329 }
330
331 sub gcgraph_cumulative {
332 my @gcdata = ();
333 my $seqio_object = Bio::SeqIO->new(-file => $options->{file}, -format=>'genbank');
334 # Only handing first sequence.
335 my $seq_object = $seqio_object->next_seq;
336
337 my $parent = $seq_object->display_id();
338
339 my $seq = $seq_object->seq();
340 my $sep = int($options->{gc_skew_plot_window_size}/2);
341 my $stepsep = int($options->{gc_skew_plot_step_size}/2);
342 my $cumulative_gc_skew = 0;
343 my @cumgc_vals;
344
345 my $count = 0;
346 foreach(my $i = $sep; $i < $seq_object->length() - $sep - $options->{gc_skew_plot_step_size}; $i += $options->{gc_skew_plot_step_size}){
347 $count++;
348 $cumulative_gc_skew += _calculate_gc_skew_for_seq(substr($seq,$i-$sep,$options->{gc_skew_plot_window_size})),
349 push(@cumgc_vals, $cumulative_gc_skew);
350 push(@gcdata, join(" ",
351 $parent,
352 $i - $stepsep,
353 $i + $stepsep,
354 $cumulative_gc_skew
355 ));
356 }
357
358 my $sum = 0;
359 foreach(@cumgc_vals){$sum += $_;}
360 $cum_gc_ske_mean = $sum / $count;
361
362 return join("\n", @gcdata);
363 # Main 'chromosome' data
364 }
365 sub gcgraph {
366 my @gcdata = ();
367 my $seqio_object = Bio::SeqIO->new(-file => $options->{file}, -format=>'genbank');
368 # Only handing first sequence.
369 my $seq_object = $seqio_object->next_seq;
370
371 my $parent = $seq_object->display_id();
372
373 my $seq = $seq_object->seq();
374 my $sep = int($options->{gc_skew_plot_window_size}/2);
375 my $stepsep = int($options->{gc_skew_plot_step_size}/2);
376 foreach(my $i = $sep; $i < $seq_object->length() - $sep - $options->{gc_skew_plot_step_size}; $i += $options->{gc_skew_plot_step_size}){
377 push(@gcdata, join(" ",
378 $parent,
379 $i - $stepsep,
380 $i + $stepsep,
381 _calculate_gc_skew_for_seq(substr($seq,$i-$sep,$options->{gc_skew_plot_window_size})),
382 ));
383 }
384 return join("\n", @gcdata);
385 # Main 'chromosome' data
386 }
387 sub _calculate_gc_skew_for_seq {
388 my ($seq) = @_;
389 $seq = uc($seq);
390 my %counts;
391 foreach(split //,$seq){
392 $counts{$_}++;
393 }
394 if($counts{G} + $counts{C} > 0){
395 return ($counts{G} - $counts{C}) / ($counts{G} + $counts{C});
396 }
397 return 0;
398 }
399
400 push(@files, [ 'etc/karyotype.conf', karyotype() ]);
401 push(@files, [ 'etc/circos.conf', circosconf() ]);
402 push(@files, [ 'etc/ideogram.conf', ideogramconf() ]);
403 if($options->{enable_gc_skew_plot}){
404 push(@files, [ 'data/gc.txt', gcgraph() ]);
405 push(@files, [ 'data/gc_cumulative.txt', gcgraph_cumulative() ]);
406 }
407
408 push(@files, [ 'etc/plots.conf', plotsconf() ]);
409 push(@files, [ 'etc/ticks.conf', ticksconf() ]);
410
411
412 use Archive::Any::Create;
413 my $archive = Archive::Any::Create->new();
414 $archive->container('conf');
415 foreach(@files){
416 my ($location, $content) = @{$_};
417 $archive->add_file($location, $content);
418 }
419
420 use CPT::OutputFiles;
421 my $crr_output = CPT::OutputFiles->new(
422 name => 'output_circos_confs',
423 GGO => $ggo,
424 );
425 $crr_output->CRR(data => $archive);
426
427 =head1 NAME
428
429 DNAPlotter
430
431 =head1 DESCRIPTION
432
433 Much like artemis's DNAPlotter, this tool plots genomes in a ciruclar fashion, and can plot gc-deviation tracks as well. The options are somewhat reduced compared to artemis, so if you need something that isn't available in this version please file a bug report.
434
435 =head1 USE
436
437 Each track has several parameters:
438
439 =over 4
440
441 =item track_key
442
443 This selects a set of features from a genbank file, e.g., CDSs or tRNAs
444
445 =item track_feature_filter_invert
446
447 This parameter will invert (negate) the results of whatever query parameters you specify after it.
448
449 =item track_feature_filter_hastag
450
451 Require that a feature has a specific tag....
452
453 =item track_feature_filter_textquery
454
455 ...with this specific text in it
456
457 =item track_feature_filter_strand
458
459 Which strand should the feature be on (not inverted)
460
461 =back
462
463 Additionally, users are able to enable/disable GC skew plots. it's suggested that these are generally left alone, as they can quickly increase runtime.
464
465 =cut