circos 0.64-1 / usr / share / perl5 / Circos / Ideogram.pm

package Circos::Ideogram;

=pod

=head1 NAME

Circos::Ideogram - ideogram routines for Circos

=head1 SYNOPSIS

This module is not meant to be used directly.

=head1 DESCRIPTION

Circos is an application for the generation of publication-quality,
circularly composited renditions of genomic data and related
annotations.

Circos is particularly suited for visualizing alignments, conservation
and intra and inter-chromosomal relationships. However, Circos can be
used to plot any kind of 2D data in a circular layout - its use is not
limited to genomics. Circos' use of lines to relate position pairs
(ribbons add a thickness parameter to each end) is effective to
display relationships between objects or positions on one or more
scales.

All documentation is in the form of tutorials at L<http://www.circos.ca>.

=cut

# -------------------------------------------------------------------

use strict;
use warnings;

use base 'Exporter';
our @EXPORT = qw(
is_on_ideogram
);

use Carp qw( carp confess croak );
use Cwd;
use FindBin;
use File::Spec::Functions;
use Math::Round;
use Math::VecStat qw(max);
#use Memoize;
use Params::Validate qw(:all);
#use Regexp::Common qw(number);

use POSIX qw(floor ceil);

use lib "$FindBin::RealBin";
use lib "$FindBin::RealBin/../lib";
use lib "$FindBin::RealBin/lib";

use Circos::Configuration;
use Circos::Constants;
use Circos::Debug;
use Circos::Error;
use Circos::Utils;

################################################################
#
# Determine which chromosomes are going to be displayed. Several parameters
# are used together to determine the list and order of displayed chromosomes.
#
# - chromosomes
# - chromosomes_breaks
# - chromosomes_display_default
# - chromosomes_order_by_karyotype
#
# If chromosomes_display_default is set to 'yes', then any chromosomes that
# appear in the karyotype are appended to the 'chromosomes' parameter. 
# The order in which they are appended depends on the value of 'chromosomes_order_by_karyotype'.
# 
# If you want to display only those chromosomes that are mentioned in the 
# 'chromosomes' parameter, then set chromosomes_display_default=no.
#
# Both 'chromosomes' and 'chromosomes_breaks' define a list of chromosome regions
# to show, delimited by ;
#
# name{[tag]}{:runlist}
#
# e.g.   hs1
#        hs1[a]
#        hs1:0-100
#        hs1[a]:0-100
#        hs1[a]:0-100,150-200
#        hs1;hs2[a];hs3:0-100
#
# You can also use "=" as the field delimiter,
#
# e.g.   hs1=0-100
#        hs1[a]=0-100
#
# The functional role of 'chromosomes' and 'chromosomes_breaks' is the same. The latter
# gives an opportunity to separate definitions of regions which are not shown.
# 

sub parse_chromosomes {
  my $karyotype = shift;

  my @chrs;
  
  # get the chromosomes in the karyotype
  my @chrs_in_k = keys %$karyotype;
  # sort them by their appearance in the file
  @chrs_in_k = sort { $karyotype->{$a}{chr}{display_order} <=> $karyotype->{$b}{chr}{display_order} } @chrs_in_k;
  printdebug_group("chrfilter","karyotypeorder",@chrs_in_k);
  # sort them by digit in chromosome name (e.g. chr1 before chr11 before chrx) (done in Karyotype::sort_karyotype)
  my @chrs_in_k_native_sort = sort { $karyotype->{$a}{chr}{sort_idx} <=> $karyotype->{$b}{chr}{sort_idx} } @chrs_in_k;
  printdebug_group("chrfilter","nativesort",@chrs_in_k_native_sort);

  if ( $CONF{chromosomes_display_default} ) {
      #
      # The default order for chromosomes is string-then-number if
      # chromosomes contain a number, and if not then asciibetic
      #
      # I used to have this based on the order in the KARYOTYPE (use
      # {CHR}{chr}{display_order} field) but decided to change it.
      #
      if ( $CONF{chromosomes_order_by_karyotype} ) {
	  # @chrs_in_k already ordered by appearance
	  printdebug_group("chrfilter","using karyotypeorder");
      } else {
	  # sort the chromosomes with digits in them
	  @chrs_in_k = @chrs_in_k_native_sort;
	  printdebug_group("chrfilter","using nativesort");
      }
      
      ################################################################
      # Reconstruct the $CONF{chromosomes} argument using 
      # chromosomes from karyotype and those in $CONF{chromosomes}
      my ($chrs_in_c,@chrs_ordered);
      
      # First, parse chromosomes and regular expressions in the chromosomes string
      if ( $CONF{chromosomes} ) {
	  $chrs_in_c = parse_chromosomes_string($CONF{chromosomes});
	  #printdumper($chrs_in_c);exit;
      }
    CHR_IN_K:
      for my $chr (@chrs_in_k) {
	  my $found;
	  if ($chrs_in_c) {
	      my $accept = 1;
	      my @chr_in_c_found;
	      for my $isrx (1,0) {
		  # for each chromosome in $CONF{chromosomes}
		CHR_IN_C:
		  for my $chr_in_c (@$chrs_in_c) {
		      next if $isrx && ! $chr_in_c->{rx};
		      next if ! $isrx && $chr_in_c->{rx};
		      my $reject = $chr_in_c->{reject};
		      my $match  = 0;
		      if ($isrx && $chr =~ $chr_in_c->{rx}) {
			  $match = 1;
		      }
		      if (! $isrx && $chr eq $chr_in_c->{chr}) {
			  $match = 1;
		      }
		      printdebug_group("chrfilter",$chr,"inchromosomes",$chr_in_c->{chr},"isrx",defined $chr_in_c->{rx},"match",$match,"reject",$reject);
		      if ($match) {
			  push @chr_in_c_found, $chr_in_c;
			  if ($reject) {
			      $accept = 0;
			  } else {
			      $accept = 1;
			  }
			  #last CHR_IN_C;
		      }
		  }
		  printdebug_group("chrfilter",$chr,"rx",$isrx,"accept",$accept);
	      }
	      if ($accept) {
		  if (@chr_in_c_found) {
		      #printdumper(@chr_in_c_found);
		      for my $c (@chr_in_c_found) {
			  next if $c->{reject};
			  if (! $c->{rx}) {
			      my $str = $c->{chr};
			      if ($c->{reject}) {
				  $str = "-$str";
			      }
			      if ($c->{tag}) {
				  $str .= sprintf("[%s]",$c->{tag});
			      }
			      if ($c->{runlist}) {
				  $str .= sprintf(":%s",$c->{runlist});
			      }
			      push @chrs_ordered, $str;
			  } else {
			      push @chrs_ordered, $chr;
			  }
		      }
		  } else {
		      push @chrs_ordered, $chr;
		  }
	      } else {
		  push @chrs_ordered, "-$chr";
	      }
	  } else {
	      push @chrs_ordered, $chr;
	  }
      }
      $CONF{chromosomes} = join( $SEMICOLON, @chrs_ordered );
  } 
  printdebug_group("chrfilter","effective 'chromosomes' parameter",$CONF{chromosomes});
  
  my %karyotype_chrs_seen;
  
  for my $isrx (1,0) {
      for my $pair ([$CONF{chromosomes},1],[$CONF{chromosomes_breaks},0]) {
	  my ($string,$accept_default) = @$pair;
	  my $chrstring_list = Circos::Configuration::make_parameter_list_array($string,qr/\s*;\s*/);
	  for my $chrstring (@$chrstring_list) {
	      my $chr_record = parse_chromosomes_record($chrstring);
	      my ($reject,$chr,$runlist,$tag,$chrrx) = @{$chr_record}{qw(reject chr runlist tag rx)};
	      $tag       = $EMPTY_STR if !defined $tag;
	      $chr       = $EMPTY_STR if !defined $chr;
	      $runlist   = $EMPTY_STR if !defined $runlist;
	      if ($chr eq $EMPTY_STR) {
		  fatal_error("ideogram","unparsable_def",$chrstring);
	      }
	      next if $isrx   && ! defined $chrrx;
	      next if ! $isrx && defined $chrrx;
	      # $accept identifies whether the regions indicate inclusions or exclusions
	      # $accept=1 this region is to be included
	      # $accept=0 this region is to be included (region prefixed by -)
	      my $accept = $accept_default;
	      $accept = 0 if $reject;
	      if ( $isrx && $tag) {
		  fatal_error("ideogram","regex_tag",$chrstring,$tag);
	      }
	      my $chrkey = make_key($chr,$tag);
	      #printinfo($isrx, $chrrx, $isrx ? "RX" : "NOTRX", !$isrx && defined $chrrx ? "next" : "accept");
	      if ( ! $isrx && ! defined $karyotype->{$chr}{chr} ) {
		  fatal_error("ideogram","use_undefined",$chrstring,$chr);
	      }
	      
	      my @chrs_to_store;
	      if ($isrx) {
		  for my $c (@chrs_in_k_native_sort) {
		      next if $accept && $karyotype_chrs_seen{ make_key($c,$tag) };
		      if ($c =~ /$chrrx/i) {
			  push @chrs_to_store, $c;
			  $karyotype_chrs_seen{ make_key($c,$tag) }++;
			  $karyotype_chrs_seen{ make_key($c,"") }++;
		      }
		  }
	      } else {
		  for my $c (@chrs_in_k_native_sort) {
		      next if $accept && $karyotype_chrs_seen{ make_key($c,$tag) };
		      if ($c eq $chr) {
			  push @chrs_to_store, $c;
			  $karyotype_chrs_seen{ make_key($c,$tag) }++;
			  $karyotype_chrs_seen{ make_key($c,"") }++;
		      }
		  }
	      }
	      #printdumper(\%karyotype_chrs_seen);
	      printdebug_group("chrfilter","chrrx",$chrstring,"rx?",$isrx,"accept",$accept,"tag",$tag || "-","chrs",@chrs_to_store);
	      next unless @chrs_to_store;
	      
	      sub make_key {
		  my ($chr,$tag) = @_;
		  $tag ||= "";
		  return sprintf("%s_%s",$chr,$tag);
	      }
	      
	      # all numbers in runlist are automatically multiplied by
	      # chromosomes_units value - this saves you from having to type
	      # a lot of zeroes in the runlists
	      
	      if ( $CONF{chromosomes_units} ) {
		  $runlist =~ s/([\.\d]+)/$1*$CONF{chromosomes_units}/eg;
	      }
	      
	      for my $c (@chrs_to_store) {
		  # are we trying to remove this chromosome?
		  printdebug_group("chrfilter","parsed chromosome range", $c, $runlist || $DASH );
		  my $set = $runlist ? Set::IntSpan->new($runlist) : $karyotype->{$c}{chr}{set};
		  
		  ################################################################
		  # uncertain - what was I trying to do here?
		  $set->remove($set->max) if $runlist;
		  if ( ! $accept ) {
		      $set->remove( $set->min ) if $set->min;
		      $set->remove( $set->max );
		  }
		  ################################################################
		  
		  if ($accept) {
		      push @chrs,
		      {
			  chr    => $c,
			  tag    => $tag || $c,
			  idx    => int(@chrs),
			  accept => $accept,
			  set    => $set
		      };
		      $karyotype->{$c}{chr}{display_region}{accept} ||= Set::IntSpan->new();
		      $karyotype->{$c}{chr}{display_region}{accept} = $karyotype->{$c}{chr}{display_region}{accept}->union($set);
		  } else {
		      if ($accept_default) {
			  @chrs = grep($_->{chr} ne $c && $_->{tag} ne $tag, @chrs);
		      }
		      $karyotype->{$c}{chr}{display_region}{reject} ||= Set::IntSpan->new();
		      $karyotype->{$c}{chr}{display_region}{reject} = $karyotype->{$c}{chr}{display_region}{reject}->union($set);
		  }
	      }
	  }
      }
  }
  
  if ( ! grep( $_->{accept}, @chrs ) ) {
      fatal_error("ideogram","no_ideograms_to_draw");
  }
  
  for my $c (@chrs) {
      printdebug_group("chrfilter","chrlist",sprintf("%2d %4s %4s %d %10s %10s",
						     $c->{idx},
						     $c->{chr},
						     $c->{tag}||$EMPTY_STR,
						     $c->{accept}||"-",
						     defined $c->{set}->min ? $c->{set}->min : "(-",
						     defined $c->{set}->max ? $c->{set}->max : "-)"));
  }
  return @chrs;
}

sub parse_chromosomes_string {
    my $str = shift;
    my $data;
    my $delim = ";";
    for my $record ( @{Circos::Configuration::make_parameter_list_array($str,qr/\s*$delim\s*/)} ) {
	push @$data, parse_chromosomes_record($record);
    }
    return $data;
}

#  hs1
# -hs1
#  hs1:1-10
#  hs1[a]:1-10
# -hs1:1-10
#  /h/
# -/h/
#
sub parse_chromosomes_record {
  my $str = shift;
  my $default_delim  = "[:=]";
  my ($chr,$runlist) = split(Circos::Configuration::fetch_configuration("list_field_delim") || $default_delim ,$str);
  my ($tag,$reject,$rx);
  ( $reject, $chr, $tag ) = $chr =~ /([!-])?(.+?)(?:\[([^\[\]]+)\])?$/;
  $reject = $reject ? 1 : 0;
  $rx = parse_as_rx($chr) || undef;
  my $isrx = $rx ? 1 : 0;
  printdebug_group("chrfilter","parsed chr record",$str,"chr",$chr,"tag",$tag,"reject",$reject,"runlist",$runlist,"rx",$rx,"rx?",$isrx);
  return { chr=>$chr,
	   rx=>$rx,
	   isrx=>$isrx,
	   reject=> $reject ? 1 : 0,
	   accept=> $reject ? 0 : 1,
	   tag=>$tag || undef,
	   runlist=>$runlist };
}

# -------------------------------------------------------------------
sub report_chromosomes {
  my $karyotype = shift;
  for my $chr (
	       sort {
		 $karyotype->{$a}{chr}{display_order} <=> $karyotype->{$b}{chr}{display_order}
	       } keys %$karyotype
	      ) {
    next unless $karyotype->{$chr}{chr}{display};
    
    printinfo(
	      $chr,
	      $karyotype->{$chr}{chr}{display_order},
	      $karyotype->{$chr}{chr}{scale},
	      $karyotype->{$chr}{chr}{display_region}
	      ? $karyotype->{$chr}{chr}{display_region}->run_list
	      : $DASH,
	      $karyotype->{$chr}{chr}{length_cumul}
	     );
  }
}

# Quickly check whether a point appears on a drawn portion of
# a chromosome. An ideogram must exist and have a cover that
# intersects with the data point coordinate. Not using Set::IntSpan
# here because it is too slow.
sub is_on_ideogram {
	my $datum       = shift;
	my $on_ideogram = 1;
	for my $point ( @{$datum->{data}}) {
		my $chr       = $point->{chr};
		my $ideograms = Circos::get_ideograms_by_name($chr);
		# If the chromosome is not being displayed, it has no ideogram.
		return if ! defined $ideograms;
		my ($start,$end) = ($point->{start},$point->{end});
		for my $ideogram (@$ideograms) {
			for my $cover (@{$ideogram->{covers}}) {
				my ($cstart,$cend) = ($cover->{set}->min,$cover->{set}->max);
				if($cstart > $end || $cend < $start) {
					# not on ideogram
				} else {
					return 1;
				}
			}
		}
	}
	# none of the points coordinates fell inside any ideogram cover
	return;
}

1;