← Back to branch summary

~ubuntu-branches/ubuntu/utopic/circos/utopic-proposed

« back to all changes in this revision

Viewing changes to lib/Circos/Ideogram.pm

  • Committer: Package Import Robot
  • Author(s): Olivier Sallou, Olivier Sallou, Charles Plessy, Andreas Tille
  • Date: 2012-06-14 12:56:33 UTC
  • mfrom: (1.1.1)
  • Revision ID: package-import@ubuntu.com-20120614125633-0wh7ovv69s5k1uiq
Tags: 0.61-1
[ Olivier Sallou ]
* New upstream release

[ Charles Plessy ]
* renamed debian/upstream-metadata.yaml to debian/upstream

[ Andreas Tille ]
* debian/upstream: enhanced citation information 

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/Circos/Ideogram.pm
 
1
package Circos::Ideogram;
 
2
 
 
3
=pod
 
4
 
 
5
=head1 NAME
 
6
 
 
7
Circos::Ideogram - ideogram routines for Circos
 
8
 
 
9
=head1 SYNOPSIS
 
10
 
 
11
This module is not meant to be used directly.
 
12
 
 
13
=head1 DESCRIPTION
 
14
 
 
15
Circos is an application for the generation of publication-quality,
 
16
circularly composited renditions of genomic data and related
 
17
annotations.
 
18
 
 
19
Circos is particularly suited for visualizing alignments, conservation
 
20
and intra and inter-chromosomal relationships. However, Circos can be
 
21
used to plot any kind of 2D data in a circular layout - its use is not
 
22
limited to genomics. Circos' use of lines to relate position pairs
 
23
(ribbons add a thickness parameter to each end) is effective to
 
24
display relationships between objects or positions on one or more
 
25
scales.
 
26
 
 
27
All documentation is in the form of tutorials at L<http://www.circos.ca>.
 
28
 
 
29
=cut
 
30
 
 
31
# -------------------------------------------------------------------
 
32
 
 
33
use strict;
 
34
use warnings;
 
35
 
 
36
use base 'Exporter';
 
37
our @EXPORT = qw(
 
38
 
 
39
);
 
40
 
 
41
use Carp qw( carp confess croak );
 
42
use Cwd;
 
43
use FindBin;
 
44
use File::Spec::Functions;
 
45
use Math::Round;
 
46
use Math::VecStat qw(max);
 
47
#use Memoize;
 
48
use Params::Validate qw(:all);
 
49
#use Regexp::Common qw(number);
 
50
 
 
51
use POSIX qw(floor ceil);
 
52
 
 
53
use lib "$FindBin::RealBin";
 
54
use lib "$FindBin::RealBin/../lib";
 
55
use lib "$FindBin::RealBin/lib";
 
56
 
 
57
use Circos::Configuration;
 
58
use Circos::Constants;
 
59
use Circos::Debug;
 
60
use Circos::Error;
 
61
use Circos::Utils;
 
62
 
 
63
################################################################
 
64
#
 
65
# Determine which chromosomes are going to be displayed. Several parameters
 
66
# are used together to determine the list and order of displayed chromosomes.
 
67
#
 
68
# - chromosomes
 
69
# - chromosomes_breaks
 
70
# - chromosomes_display_default
 
71
# - chromosomes_order_by_karyotype
 
72
#
 
73
# If chromosomes_display_default is set to 'yes', then any chromosomes that
 
74
# appear in the karyotype are appended to the 'chromosomes' parameter. 
 
75
# The order in which they are appended depends on the value of 'chromosomes_order_by_karyotype'.
 
76
# 
 
77
# If you want to display only those chromosomes that are mentioned in the 
 
78
# 'chromosomes' parameter, then set chromosomes_display_default=no.
 
79
#
 
80
# Both 'chromosomes' and 'chromosomes_breaks' define a list of chromosome regions
 
81
# to show, delimited by ;
 
82
#
 
83
# name{[tag]}{:runlist}
 
84
#
 
85
# e.g.   hs1
 
86
#        hs1[a]
 
87
#        hs1:0-100
 
88
#        hs1[a]:0-100
 
89
#        hs1[a]:0-100,150-200
 
90
#        hs1;hs2[a];hs3:0-100
 
91
#
 
92
# You can also use "=" as the field delimiter,
 
93
#
 
94
# e.g.   hs1=0-100
 
95
#        hs1[a]=0-100
 
96
#
 
97
# The functional role of 'chromosomes' and 'chromosomes_breaks' is the same. The latter
 
98
# gives an opportunity to separate definitions of regions which are not shown.
 
99
# 
 
100
 
 
101
sub parse_chromosomes {
 
102
  my $karyotype = shift;
 
103
 
 
104
  my @chrs;
 
105
  
 
106
  # get the chromosomes in the karyotype
 
107
  my @chrs_in_k = keys %$karyotype;
 
108
  # sort them by their appearance in the file
 
109
  @chrs_in_k = sort { $karyotype->{$a}{chr}{display_order} <=> $karyotype->{$b}{chr}{display_order} } @chrs_in_k;
 
110
  printdebug_group("chrfilter","karyotypeorder",@chrs_in_k);
 
111
  # sort them by digit in chromosome name (e.g. chr1 before chr11 before chrx) (done in Karyotype::sort_karyotype)
 
112
  my @chrs_in_k_native_sort = sort { $karyotype->{$a}{chr}{sort_idx} <=> $karyotype->{$b}{chr}{sort_idx} } @chrs_in_k;
 
113
  printdebug_group("chrfilter","nativesort",@chrs_in_k_native_sort);
 
114
 
 
115
  if ( $CONF{chromosomes_display_default} ) {
 
116
      #
 
117
      # The default order for chromosomes is string-then-number if
 
118
      # chromosomes contain a number, and if not then asciibetic
 
119
      #
 
120
      # I used to have this based on the order in the KARYOTYPE (use
 
121
      # {CHR}{chr}{display_order} field) but decided to change it.
 
122
      #
 
123
      if ( $CONF{chromosomes_order_by_karyotype} ) {
 
124
          # @chrs_in_k already ordered by appearance
 
125
          printdebug_group("chrfilter","using karyotypeorder");
 
126
      } else {
 
127
          # sort the chromosomes with digits in them
 
128
          @chrs_in_k = @chrs_in_k_native_sort;
 
129
          printdebug_group("chrfilter","using nativesort");
 
130
      }
 
131
      
 
132
      ################################################################
 
133
      # Reconstruct the $CONF{chromosomes} argument using 
 
134
      # chromosomes from karyotype and those in $CONF{chromosomes}
 
135
      my ($chrs_in_c,@chrs_ordered);
 
136
      
 
137
      # First, parse chromosomes and regular expressions in the chromosomes string
 
138
      if ( $CONF{chromosomes} ) {
 
139
          $chrs_in_c = parse_chromosomes_string($CONF{chromosomes});
 
140
          #printdumper($chrs_in_c);exit;
 
141
      }
 
142
    CHR_IN_K:
 
143
      for my $chr (@chrs_in_k) {
 
144
          my $found;
 
145
          if ($chrs_in_c) {
 
146
              my $accept = 1;
 
147
              my @chr_in_c_found;
 
148
              for my $isrx (1,0) {
 
149
                  # for each chromosome in $CONF{chromosomes}
 
150
                CHR_IN_C:
 
151
                  for my $chr_in_c (@$chrs_in_c) {
 
152
                      next if $isrx && ! $chr_in_c->{rx};
 
153
                      next if ! $isrx && $chr_in_c->{rx};
 
154
                      my $reject = $chr_in_c->{reject};
 
155
                      my $match  = 0;
 
156
                      if ($isrx && $chr =~ $chr_in_c->{rx}) {
 
157
                          $match = 1;
 
158
                      }
 
159
                      if (! $isrx && $chr eq $chr_in_c->{chr}) {
 
160
                          $match = 1;
 
161
                      }
 
162
                      printdebug_group("chrfilter",$chr,"inchromosomes",$chr_in_c->{chr},"isrx",defined $chr_in_c->{rx},"match",$match,"reject",$reject);
 
163
                      if ($match) {
 
164
                          push @chr_in_c_found, $chr_in_c;
 
165
                          if ($reject) {
 
166
                              $accept = 0;
 
167
                          } else {
 
168
                              $accept = 1;
 
169
                          }
 
170
                          #last CHR_IN_C;
 
171
                      }
 
172
                  }
 
173
                  printdebug_group("chrfilter",$chr,"rx",$isrx,"accept",$accept);
 
174
              }
 
175
              if ($accept) {
 
176
                  if (@chr_in_c_found) {
 
177
                      #printdumper(@chr_in_c_found);
 
178
                      for my $c (@chr_in_c_found) {
 
179
                          next if $c->{reject};
 
180
                          if (! $c->{rx}) {
 
181
                              my $str = $c->{chr};
 
182
                              if ($c->{reject}) {
 
183
                                  $str = "-$str";
 
184
                              }
 
185
                              if ($c->{tag}) {
 
186
                                  $str .= sprintf("[%s]",$c->{tag});
 
187
                              }
 
188
                              if ($c->{runlist}) {
 
189
                                  $str .= sprintf(":%s",$c->{runlist});
 
190
                              }
 
191
                              push @chrs_ordered, $str;
 
192
                          } else {
 
193
                              push @chrs_ordered, $chr;
 
194
                          }
 
195
                      }
 
196
                  } else {
 
197
                      push @chrs_ordered, $chr;
 
198
                  }
 
199
              } else {
 
200
                  push @chrs_ordered, "-$chr";
 
201
              }
 
202
          } else {
 
203
              push @chrs_ordered, $chr;
 
204
          }
 
205
      }
 
206
      $CONF{chromosomes} = join( $SEMICOLON, @chrs_ordered );
 
207
  } 
 
208
  printdebug_group("chrfilter","effective 'chromosomes' parameter",$CONF{chromosomes});
 
209
  
 
210
  my %karyotype_chrs_seen;
 
211
  
 
212
  for my $isrx (1,0) {
 
213
      for my $pair ([$CONF{chromosomes},1],[$CONF{chromosomes_breaks},0]) {
 
214
          my ($string,$accept_default) = @$pair;
 
215
          my $chrstring_list = Circos::Configuration::make_parameter_list_array($string,qr/\s*;\s*/);
 
216
          for my $chrstring (@$chrstring_list) {
 
217
              my $chr_record = parse_chromosomes_record($chrstring);
 
218
              my ($reject,$chr,$runlist,$tag,$chrrx) = @{$chr_record}{qw(reject chr runlist tag rx)};
 
219
              $tag       = $EMPTY_STR if !defined $tag;
 
220
              $chr       = $EMPTY_STR if !defined $chr;
 
221
              $runlist   = $EMPTY_STR if !defined $runlist;
 
222
              if ($chr eq $EMPTY_STR) {
 
223
                  fatal_error("ideogram","unparsable_def",$chrstring);
 
224
              }
 
225
              next if $isrx   && ! defined $chrrx;
 
226
              next if ! $isrx && defined $chrrx;
 
227
              # $accept identifies whether the regions indicate inclusions or exclusions
 
228
              # $accept=1 this region is to be included
 
229
              # $accept=0 this region is to be included (region prefixed by -)
 
230
              my $accept = $accept_default;
 
231
              $accept = 0 if $reject;
 
232
              if ( $isrx && $tag) {
 
233
                  fatal_error("ideogram","regex_tag",$chrstring,$tag);
 
234
              }
 
235
              my $chrkey = make_key($chr,$tag);
 
236
              #printinfo($isrx, $chrrx, $isrx ? "RX" : "NOTRX", !$isrx && defined $chrrx ? "next" : "accept");
 
237
              if ( ! $isrx && ! defined $karyotype->{$chr}{chr} ) {
 
238
                  fatal_error("ideogram","use_undefined",$chrstring,$chr);
 
239
              }
 
240
              
 
241
              my @chrs_to_store;
 
242
              if ($isrx) {
 
243
                  for my $c (@chrs_in_k_native_sort) {
 
244
                      next if $accept && $karyotype_chrs_seen{ make_key($c,$tag) };
 
245
                      if ($c =~ /$chrrx/i) {
 
246
                          push @chrs_to_store, $c;
 
247
                          $karyotype_chrs_seen{ make_key($c,$tag) }++;
 
248
                          $karyotype_chrs_seen{ make_key($c,"") }++;
 
249
                      }
 
250
                  }
 
251
              } else {
 
252
                  for my $c (@chrs_in_k_native_sort) {
 
253
                      next if $accept && $karyotype_chrs_seen{ make_key($c,$tag) };
 
254
                      if ($c eq $chr) {
 
255
                          push @chrs_to_store, $c;
 
256
                          $karyotype_chrs_seen{ make_key($c,$tag) }++;
 
257
                          $karyotype_chrs_seen{ make_key($c,"") }++;
 
258
                      }
 
259
                  }
 
260
              }
 
261
              #printdumper(\%karyotype_chrs_seen);
 
262
              printdebug_group("chrfilter","chrrx",$chrstring,"rx?",$isrx,"accept",$accept,"tag",$tag || "-","chrs",@chrs_to_store);
 
263
              next unless @chrs_to_store;
 
264
              
 
265
              sub make_key {
 
266
                  my ($chr,$tag) = @_;
 
267
                  $tag ||= "";
 
268
                  return sprintf("%s_%s",$chr,$tag);
 
269
              }
 
270
              
 
271
              # all numbers in runlist are automatically multiplied by
 
272
              # chromosomes_units value - this saves you from having to type
 
273
              # a lot of zeroes in the runlists
 
274
              
 
275
              if ( $CONF{chromosomes_units} ) {
 
276
                  $runlist =~ s/([\.\d]+)/$1*$CONF{chromosomes_units}/eg;
 
277
              }
 
278
              
 
279
              for my $c (@chrs_to_store) {
 
280
                  # are we trying to remove this chromosome?
 
281
                  printdebug_group("chrfilter","parsed chromosome range", $c, $runlist || $DASH );
 
282
                  my $set = $runlist ? Set::IntSpan->new($runlist) : $karyotype->{$c}{chr}{set};
 
283
                  
 
284
                  ################################################################
 
285
                  # uncertain - what was I trying to do here?
 
286
                  $set->remove($set->max) if $runlist;
 
287
                  if ( ! $accept ) {
 
288
                      $set->remove( $set->min ) if $set->min;
 
289
                      $set->remove( $set->max );
 
290
                  }
 
291
                  ################################################################
 
292
                  
 
293
                  if ($accept) {
 
294
                      push @chrs,
 
295
                      {
 
296
                          chr    => $c,
 
297
                          tag    => $tag || $c,
 
298
                          idx    => int(@chrs),
 
299
                          accept => $accept,
 
300
                          set    => $set
 
301
                      };
 
302
                      $karyotype->{$c}{chr}{display_region}{accept} ||= Set::IntSpan->new();
 
303
                      $karyotype->{$c}{chr}{display_region}{accept} = $karyotype->{$c}{chr}{display_region}{accept}->union($set);
 
304
                  } else {
 
305
                      if ($accept_default) {
 
306
                          @chrs = grep($_->{chr} ne $c && $_->{tag} ne $tag, @chrs);
 
307
                      }
 
308
                      $karyotype->{$c}{chr}{display_region}{reject} ||= Set::IntSpan->new();
 
309
                      $karyotype->{$c}{chr}{display_region}{reject} = $karyotype->{$c}{chr}{display_region}{reject}->union($set);
 
310
                  }
 
311
              }
 
312
          }
 
313
      }
 
314
  }
 
315
  
 
316
  if ( ! grep( $_->{accept}, @chrs ) ) {
 
317
      fatal_error("ideogram","no_ideograms_to_draw");
 
318
  }
 
319
  
 
320
  for my $c (@chrs) {
 
321
      printdebug_group("chrfilter","chrlist",sprintf("%2d %4s %4s %d %10s %10s",
 
322
                                                     $c->{idx},
 
323
                                                     $c->{chr},
 
324
                                                     $c->{tag}||$EMPTY_STR,
 
325
                                                     $c->{accept}||"-",
 
326
                                                     defined $c->{set}->min ? $c->{set}->min : "(-",
 
327
                                                     defined $c->{set}->max ? $c->{set}->max : "-)"));
 
328
  }
 
329
  return @chrs;
 
330
}
 
331
 
 
332
sub parse_chromosomes_string {
 
333
    my $str = shift;
 
334
    my $data;
 
335
    my $delim = ";";
 
336
    for my $record ( @{Circos::Configuration::make_parameter_list_array($str,qr/\s*$delim\s*/)} ) {
 
337
        push @$data, parse_chromosomes_record($record);
 
338
    }
 
339
    return $data;
 
340
}
 
341
 
 
342
#  hs1
 
343
# -hs1
 
344
#  hs1:1-10
 
345
#  hs1[a]:1-10
 
346
# -hs1:1-10
 
347
#  /h/
 
348
# -/h/
 
349
#
 
350
sub parse_chromosomes_record {
 
351
  my $str = shift;
 
352
  my $default_delim  = "[:=]";
 
353
  my ($chr,$runlist) = split(Circos::Configuration::fetch_configuration("list_field_delim") || $default_delim ,$str);
 
354
  my ($tag,$reject,$rx);
 
355
  ( $reject, $chr, $tag ) = $chr =~ /([!-])?(.+?)(?:\[([^\[\]]+)\])?$/;
 
356
  $reject = $reject ? 1 : 0;
 
357
  $rx = parse_as_rx($chr) || undef;
 
358
  my $isrx = $rx ? 1 : 0;
 
359
  printdebug_group("chrfilter","parsed chr record",$str,"chr",$chr,"tag",$tag,"reject",$reject,"runlist",$runlist,"rx",$rx,"rx?",$isrx);
 
360
  return { chr=>$chr,
 
361
           rx=>$rx,
 
362
           isrx=>$isrx,
 
363
           reject=> $reject ? 1 : 0,
 
364
           accept=> $reject ? 0 : 1,
 
365
           tag=>$tag || undef,
 
366
           runlist=>$runlist };
 
367
}
 
368
 
 
369
# -------------------------------------------------------------------
 
370
sub report_chromosomes {
 
371
  my $karyotype = shift;
 
372
  for my $chr (
 
373
               sort {
 
374
                 $karyotype->{$a}{chr}{display_order} <=> $karyotype->{$b}{chr}{display_order}
 
375
               } keys %$karyotype
 
376
              ) {
 
377
    next unless $karyotype->{$chr}{chr}{display};
 
378
    
 
379
    printinfo(
 
380
              $chr,
 
381
              $karyotype->{$chr}{chr}{display_order},
 
382
              $karyotype->{$chr}{chr}{scale},
 
383
              $karyotype->{$chr}{chr}{display_region}
 
384
              ? $karyotype->{$chr}{chr}{display_region}->run_list
 
385
              : $DASH,
 
386
              $karyotype->{$chr}{chr}{length_cumul}
 
387
             );
 
388
  }
 
389
}
 
390
 
 
391
1;