#!/usr/bin/perl -w
###########################################################################
#
# Prepare liftover chains for a "fastLift.pl".  Namely, split apart
# overlapping chains and sort by chromosome and position.
#
# Copyright 2010 Philip Johnson.
#
#     This program is free software: you can redistribute it and/or modify
#     it under the terms of the GNU General Public License as published
#     by the Free Software Foundation, either version 3 of the License,
#     or (at your option) any later version.
#
#     This program is distributed in the hope that it will be useful,
#     but WITHOUT ANY WARRANTY; without even the implied warranty of
#     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#     GNU General Public License at <http://www.gnu.org/licenses/> for
#     more details.
#
###########################################################################

use strict;
use Getopt::Long;
use constant {
    TNAME => 2,
    TSTART => 5,
    TEND => 6,
    QSTART => 10,
    QEND => 11,
    ID => 12
};

$| = 1;

my ($chainFile);
if (!GetOptions('chain=s' => \$chainFile)  ||
    !defined($chainFile)) {
    print "usage: $0 -chain <chainFile | chainFile.gz>\n".
        "\n".

        "Use to generate input for fastLift.pl.  Very inefficent and\n".
        "slow, but only need to run once.  Output is sent to stdout,\n".
        "so you probably want to redirect to a file.\n".
        "\n";

    exit 1;
}

my @chainsData;
my $chainDat;

#load raw UCSC chain file
if ($chainFile =~ /.gz$/) {
    open CHAIN, "zcat $chainFile |" or die $!;
} else {
    open CHAIN, $chainFile or die $!;
}
while (<CHAIN>) {
    next if /^\#/;
    if (/^c/) {
        if (defined($chainDat)) {
            push @chainsData, $chainDat;
        }
        my @header = split /\s+/;
        my ($tName, $tStart, $tEnd, $id);
        $chainDat = [@header[TNAME, TSTART, TEND, ID], ''];
    }
    $chainDat->[4] .= $_ if !/^$/;
}
if (@{$chainDat}) {
    push @chainsData, $chainDat;
}
close CHAIN;

#sort
@chainsData = sort { $a->[0] cmp $b->[0]  ||  $a->[1] <=> $b->[1] } @chainsData;

#break into non-overlapping
my $changed;
my $c = 0;
do {
    $changed = 0;
    my $numChains = @chainsData;
    for (my $i = 0;  $i < $numChains-1;  ++$i) {
        if ($chainsData[$i+1]->[0] eq $chainsData[$i]->[0]  &&
            $chainsData[$i+1]->[1] < $chainsData[$i]->[2]) {
            #print "splitting at ", $chainsData[$i+1]->[0], ": ", $chainsData[$i+1]->[1], "\n";
            push @chainsData, SplitChain($chainsData[$i], $chainsData[$i+1]->[1]);
            ++$changed;
#            if ($changed < 10) {
#                print STDERR "changed id: ", $chainsData[$i]->[3], "\n";
#            }
        }
    }
    print STDERR "round we go ", ++$c," ($changed changed)\n";
    @chainsData = sort { $a->[0] cmp $b->[0]  ||  $a->[1] <=> $b->[1] } @chainsData;
} while ($changed);

#print
##WHY DID THIS USE TO BE @chainsData-1 ???
for (my $i = 0;  $i < @chainsData;  ++$i) {
    print $chainsData[$i]->[4], "\n";
}

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

sub SplitChain {
    my ($chain, $splitPos) = @_;

    my @blocks = split /\n/, $chain->[4];
    my @header = split /\s+/, $blocks[0];
    shift @blocks;
    my ($tStart, $tEnd, $qStart, $qEnd) = @header[TSTART,TEND, QSTART,QEND];
    my ($size, $dt, $dq, $prev_tEnd, $prev_qEnd);

    ($size, $dt, $dq) = split /\s+/, $blocks[0];
    $tEnd = $tStart + $size;
    $qEnd = $qStart + $size;

    my $i = 1;
    for (;  $i < @blocks  &&  $tEnd <= $splitPos;  ++$i) {
        ($prev_tEnd, $prev_qEnd) = ($tEnd, $qEnd);
        $tStart += $size + $dt;
        $qStart += $size + $dq;
        ($size, $dt, $dq) = split /\s+/, $blocks[$i];
        $tEnd = $tStart + $size;
        $qEnd = $qStart + $size;
    }

    if ($i == 1  || $tEnd <= $splitPos) { #should have split! {
        print $i, "\n";
        print $splitPos, "\t", $tEnd, "\n";
        print join(' ', @header), "\n";
        print join("\n", @blocks), "\n";       
        die "bogus";
    }
    --$i;

    my @newHeader = @header;
    $newHeader[TSTART] = $tStart;
    $newHeader[QSTART] = $qStart;
    if ($newHeader[ID] =~ /(\d+)\.(\d+)/) {
        $newHeader[ID] = $1.'.'.($2+1);
    } else {
        $newHeader[ID] .= '.1';
    }
    my @newChain = (@newHeader[TNAME,TSTART,TEND,ID],
                    join("\n", join(' ', @newHeader), @blocks[$i..$#blocks]));
                    
    $header[TEND] = $prev_tEnd;
    $header[QEND] = $prev_qEnd;
    ($size, undef, undef) = split /\s+/, $blocks[$i-1];
    $blocks[$i-1] = $size."\n";
    @$chain = (@header[TNAME,TSTART,TEND,ID],
               join("\n", join(' ', @header), @blocks[0..($i-1)]));
    return \@newChain;
}