parsemap

#!/bin/env perl

=pod

=head1 NAME

parsemap - create Circos segment and color files from brain region map

=head1 SYNOPSIS

  cat map.txt | ./parsemap [-links map.links.txt] [-confdir etc/] [-datadir data/]

  ./parsemap -map map.txt -links map.links.txt

  # help
  ./parsemap -h

  # man page
  ./parsemap -man

=head1 DESCRIPTION

Using a map of brain regions, codes and colors, this script generates Circos files to draw a connectome diagram. 

The format of the map input (passed via -map or STDIN) is

  # region parcelated_structure_code r g b

  Fro TrFPoG/S 255 153 153
  Fro FMarG/S 204 0 51
  ...

For details and examples, see

  http://www.circos.ca/documentation/tutorials/recipes/cortical_maps/

=head1 OPTIONS

=head2 -map FILE

Definitions of regions, parcelation structures and colors. 

  Fro TrFPoG/S 255 153 153
  Fro FMarG/S 204 0 51 
  Fro MFS 255 153 51 
  ...

Optionally, you can add any number of measures to each structure - these will be rendered as heatmaps.

  Fro TrFPoG/S 255 153 153 0.910094 0.265257 0.893188 0.220351 0.810623
  Fro FMarG/S 204 0 51 0.631798 0.571077 0.332158 0.104455 0.173531
  Fro MFS 255 153 51 0.502931 0.567394 0.854165 0.0401409 0.484983
  ...

=head2 -links FILE

Links between structures. The start and end of the link are defined by the structure name along with "r" or "l" to indicate the hemisphere. The values will be associated with the 'type' and 'score' parameters for each link, which can be referenced in the configuration file using var(type) and var(score).

  r InfFGOrp l PosCS 1 0.0229917613607071
  l BSt l SbOrS 1 0.213414893099078
  l TPl r Pu 0 0.26688626172767

=head2 -structuresize 250

Change the size of each parcelation unit in the image. 

You should not have to adjust this value, since all the regions have
the same size and the absolute value of this number is not relevant in
how the figure appears.

=head1 HISTORY

=over

=item * 21 June 2014



=item * 17 May 2013

Measures for parcelation regions no longer required.

Segment labels are created in C<DATADIR/structure.label.txt>.

=item * 14 Jun 2012

First version.

=back 

=head1 REFERENCES

Irimia, A., M. C. Chambers, et al. (2012). "Circular representation of human cortical networks for subject and population-level connectomic visualization." NeuroImage.

Irimia, A., M. C. Chambers, et al. (2012). "Patient-tailored connectomics visualization for the assessment of white matter atrophy in traumatic brain injury." Frontiers in Neurology 3.

Van Horn, J. D., A. Irimia, et al. (2012). "Mapping connectivity damage in the case of phineas gage." PLoS One 7(5): e37454.

=head1 AUTHOR

Martin Krzywinski

=head1 CONTACT

Martin Krzywinski
Genome Sciences Center
BC Cancer Research Center
100-570 W 7th Ave
Vancouver BC V5Z 4S6

mkweb.bcgsc.ca
martink@bcgsc.ca

=cut

use strict;
use warnings FATAL=>"all";

use Carp;
use Config::General;
use Cwd qw(getcwd abs_path);
use Data::Dumper;
use File::Basename;
use FindBin;
use Getopt::Long;
use Pod::Usage;
use lib "$FindBin::RealBin";
use lib "$FindBin::RealBin/../lib";
use lib "$FindBin::RealBin/lib";

our (%OPT,%CONF,$conf);
our @COMMAND_LINE = ("map=s",
										 "links=s",
										 "configfile=s",
										 "datadir=s",
										 "confdir=s",
										 "structuresize=i",
										 "help",
										 "man",
										 "debug+");
our $VERSION = 0.01;

# common and custom module imports below
#use Regexp::Common;
#use IO::File;
#use List::Util;
use List::MoreUtils qw(uniq);
#use Set::IntSpan;
#use Statistics::Descriptive;

# read and parse configuration file
_parse_config();
#_dump_config();

my $ih;
if (my $file = $CONF{map}) {
  die "No such file $file" unless -e $file;
  open(FILE,$file);
  $ih = \*FILE;
} else {
  $ih = \*STDIN;
}

my @brain  = parse_map($ih);

report_colors(@brain);
write_karyotype(@brain);
write_seg_order(@brain);
write_label(@brain);
write_heatmap(@brain);
write_links(@brain);

exit;

for my $i (0..100) {
	my ($n1,$n2) = sort {rand()<=>rand()} (0..@brain-1);
	my $s1    = rand() < 0.5 ? "l" : "r";
	my $s2    = rand() < 0.5 ? "l" : "r";
	my $type  = rand() < 0.5 ? 0 : 1;
	my $score = rand();
	printinfo(
						$s1,
						$brain[$n1]{structure},
						$s2,
						$brain[$n2]{structure},
						$type,
						$score
					 );
}

sub write_links {
	my @brain = @_;
	return unless $CONF{links};
	open(L,$CONF{links});
	my $file = make_filename($CONF{datadir},"links.txt");
	open(F,">$file");
	while (<L>) {
		print $_;
		chomp;
		next if /^\s*\#/;
		next if /^\s*$/;
		my ($side1,$s1,
				$side2,$s2,
				$type,$score) = split;
		$score ||= 0;
		my $s1obj = get_structure($s1,@brain);
		my $s2obj = get_structure($s2,@brain);
		my $link = sprintf("%s %d %d %s %d %d type=%d,score=%f",
											 make_region_name($s1obj->{region},$side1) ||
											 make_region_name($s1obj->{region},"l"),
											 $s1obj->{start},
											 $s1obj->{end},
											 make_region_name($s2obj->{region},$side2) ||
											 make_region_name($s2obj->{region},"l"),
											 $s2obj->{start},
											 $s2obj->{end},
											 $type,$score);
		printf F ("%s\n",$link);
		printdebug(1,"writing link",$link);
	}
	close(L);
	close(F);
	printdebug(1,"wrote file",$file);
}

sub write_heatmap {
	my @brain = @_;
	my @regions = get_regions(@brain);

	unlink <make_filename($CONF{datadir},"measure.*.txt")>;
	my $seenfile;

	for my $side (qw(l r)) {
		for my $region (@regions) {
			my $region_name = make_region_name($region,$side);
			next unless $region_name;
			my @structures = get_structures($region,@brain);
			for my $structure (@structures) {
				my $sobj = get_structure($structure,@brain);
				next unless int @{$sobj->{heatmap}};
				for my $i (0 .. @{$sobj->{heatmap}}-1) {
					my $file = make_filename($CONF{datadir},"measure.$i.txt");
					open(F,">>$file");
					my $value = $sobj->{heatmap}[$i];
					printf F ("%s %d %d %f\n",
										$region_name,
										@{$sobj}{qw(start end)},
										$value);
					close(F);
					printdebug(1,"wrote file",$file) if ! $seenfile->{$file}++;
				}
			}
		}
	}
	close(F);
}

sub report_colors {
	my @brain = @_;
	my $file = "color.brain.conf";
	$file = make_filename($CONF{confdir},$file);
	open(F,">$file");
	my %seen;
	for my $s (@brain) {
		my $color_name = $s->{color};
		next if $seen{ $color_name } ++;
		printdebug(1,"writing color",$color_name,@{$s->{rgb}});
		printf F ("%s = %d,%d,%d\n",
							$color_name,
							@{$s->{rgb}});
	}
	close(F);
	printdebug(1,"wrote file",$file);
}

sub write_seg_order {
	my @brain = @_;
	my @r = grep($_ ne "BSt", uniq (map {$_->{region}} @brain));
	my $file = "segment.order.conf";
	$file = make_filename($CONF{confdir},$file);
	open(F,">$file");
	printdebug(1,"writing segment order");
	if(grep($_->{region} eq "BSt", @brain)) {
		printf F ("chromosomes_order = %s,%s,%s\n",
							join(",",map { lc $_ . "-r" } @r ),
							"bst",
							join(",",map { lc $_ . "-l" } reverse @r ));
	} else {
		printf F ("chromosomes_order = %s,%s\n",
							join(",",map { lc $_ . "-r" } @r ),
							join(",",map { lc $_ . "-l" } reverse @r ));

	}
	close(F);
	printdebug(1,"wrote file",$file);
}

sub write_label {
	my @brain = @_;
	my $file = make_filename($CONF{datadir},"structure.label.txt");
	open(L,">$file");
	printdebug(1,"writing structure labels");
	for my $region (@brain) {
		for my $side (qw(l r)) {
			printf L ("%s-%s %d %d %s\n",
								lc $region->{region},
								$side,
								@{$region}{qw(start end structure)});
		}
	}
	close(L);
	printdebug(1,"wrote file",$file);
}


sub make_region_name {
	my ($region,$side) = @_;
	croak if ! defined $region;
	croak if ! defined $side;
	if ($region eq "BSt") {
		if ($side eq "l") {
			return sprintf("%s",lc $region);
		} else {
			return;
		}
	} else {
		return sprintf("%s-%s",lc $region,$side);
	}
}

sub write_karyotype {
	my @brain = @_;
	my $file = make_filename($CONF{datadir},"segments.txt");
	open(K,">$file");
	my @regions = get_regions(@brain);
	for my $side (qw(l r)) {
		for my $region (@regions) {
			my $region_name = make_region_name($region,$side);
			next unless $region_name;
			my @structures = get_structures($region,@brain);
			printdebug(1,"writing structure to karyotype file",$region_name,$region,$CONF{structure_size}*int(@structures)-1);
			printf K ("chr - %s %s 0 %d black\n",
								$region_name,
								$region,
								$CONF{structure_size} * int(@structures) - 1);
			for my $structure (@structures) {
				my $sobj = get_structure($structure,@brain);
				printdebug(1,"writing substructure to karyotype file",$region_name,@{$sobj}{qw(structure start end color)});
				printf K ("band %s %s %s %d %d %s\n",
									$region_name,
									$sobj->{structure},
									$sobj->{structure},
									$sobj->{start},
									$sobj->{end},
									$sobj->{color});
			}
		}
	}
	printdebug(1,"wrote file",$file);
}

sub parse_map {
	my $ih = shift;
	my @brain;
	my $region_count;
	while (<$ih>) {
		chomp;
		next if /^\s*\#/;
		next if /^\s*$/;
		my ($region,$structure,$r,$g,$b,@hmap) = split;
		if (! defined $b) {
			confess "Could not parse the line [$_]. It must have 5 fields: region, structure, r, g, b.";
		}
		my $color_name = lc $structure;
		$color_name =~ s/\W//g;
		printdebug(1,$region,$structure,$color_name,$r,$g,$b);
		$region_count->{$region} ||= 0;
		push @brain, {region=>$region,
									rgb=>[$r,$g,$b],
									structure=>$structure,
									start=> $CONF{structure_size} * $region_count->{$region},
									end=> $CONF{structure_size} * (1+$region_count->{$region}) - 1,
									heatmap=>\@hmap,
									color=>$color_name};
		$region_count->{$region}++;
	}
	return @brain;
}

sub get_regions {
	my @brain = @_;
	return uniq(map {$_->{region}} @brain);
}

sub get_structure {
	my ($str,@brain) = @_;
	my ($s) = grep($_->{structure} eq $str, @brain);
	return $s;
}

sub get_structures {
	my ($region,@brain) = @_;
	return map {$_->{structure}} grep($_->{region} eq $region, @brain);
}


sub make_filename {
	my ($dir,$name) = @_;
	if (defined $dir) {
		my $path = "$dir/$name";
		$path =~ s/\/+/\//g;
		return $path;
	} else {
		return $name;
	}
}

sub validateconfiguration {
	$CONF{confdir}        ||= "etc/";
	$CONF{datadir}        ||= "data/";
	$CONF{structure_size} ||= 100;
	if(! -d $CONF{confdir}) {
		die "You asked for configuration files to be written to [$CONF{confdir}] but this directory is not found. Use -confdir DIR to specify the directory.";
	}
	if(! -d $CONF{datadir}) {
		die "You asked for data files to be written to [$CONF{datadir}] but this directory is not found. Use -datadir DIR to specify the directory.";
	}

}

# HOUSEKEEPING ###############################################################

sub _dump_config {
	printdumper(\%OPT,\%CONF);
}

sub _parse_config {
  my $dump_debug_level = 3;
  GetOptions(\%OPT,@COMMAND_LINE);
  pod2usage() if $OPT{help};
  pod2usage(-verbose=>2) if $OPT{man};
  loadconfiguration($OPT{configfile});
  populateconfiguration();	# copy command line options to config hash
  validateconfiguration(); 
  if (defined $CONF{debug} && $CONF{debug} == $dump_debug_level) {
    $Data::Dumper::Indent    = 2;
    $Data::Dumper::Quotekeys = 0;
    $Data::Dumper::Terse     = 0;
    $Data::Dumper::Sortkeys  = 1;
    $Data::Dumper::Varname = "OPT";
    printdumper(\%OPT);
    $Data::Dumper::Varname = "CONF";
    printdumper(\%CONF);
    exit;
  }
}

sub populateconfiguration {
  for my $var (keys %OPT) {
    $CONF{$var} = $OPT{$var};
  }
  repopulateconfiguration(\%CONF);
}

sub repopulateconfiguration {
  my $root     = shift;
  return unless ref($root) eq "HASH";
  for my $key (keys %$root) {
		my $value = $root->{$key};
		if (ref($value) eq "HASH") {
			repopulateconfiguration($value);
		} elsif (ref($value) eq "ARRAY") {
			for my $item (@$value) {
	      repopulateconfiguration($item);
			}
		} elsif (defined $value) {
			while ($value =~ /__([^_].+?)__/g) {
	      my $source = "__" . $1 . "__";
	      my $target = eval $1;
	      $value =~ s/\Q$source\E/$target/g;
			}
			$root->{$key} = $value;
		}
  }
}

################################################################
#
#

sub loadconfiguration {
  my $file = shift;
  if (defined $file) {
    if (-e $file && -r _) {
      # provided configuration file exists and can be read
      $file = abs_path($file);
    } else {
      confess "The configuration file [$file] passed with -configfile does not exist or cannot be read.";
    }
  } else {
    # otherwise, try to automatically find a configuration file
    my ($scriptname,$path,$suffix) = fileparse($0);
    my $cwd     = getcwd();
    my $bindir  = $FindBin::RealBin;
    my $userdir = $ENV{HOME};
    my @candidate_files = (
													 "$cwd/$scriptname.conf",
													 "$cwd/etc/$scriptname.conf",
													 "$cwd/../etc/$scriptname.conf",
													 "$bindir/$scriptname.conf",
													 "$bindir/etc/$scriptname.conf",
													 "$bindir/../etc/$scriptname.conf",
													 "$userdir/.$scriptname.conf",
													);
    my @additional_files = (
	
													 );
    for my $candidate_file (@additional_files,@candidate_files) {
			#printinfo("configsearch",$candidate_file);
			if (-e $candidate_file && -r _) {
				$file = $candidate_file;
				#printinfo("configfound",$candidate_file);
				last;
			}
    }
  }
  if (defined $file) {
    $OPT{configfile} = $file;
    $conf = new Config::General(
																-ConfigFile=>$file,
																-IncludeRelative=>1,
																-IncludeAgain=>1,
																-ExtendedAccess=>1,
																-AllowMultiOptions=>"yes",
																-LowerCaseNames=>1,
																-AutoTrue=>1
															 );
    %CONF = $conf->getall;
  }
}

sub printdebug {
	my ($level,@msg) = @_;
	my $prefix = "debug";
	if (defined $CONF{debug} && $CONF{debug} >= $level) {
		printinfo(sprintf("%s[%d]",$prefix,$level),@msg);
	}
}

sub printinfo {
	print join(" ",@_),"\n";
}

sub printdumper {
	print Dumper(@_);
}