# Copyright (C) 2014 Alejandro Ochoa, Singh Research Group, Princeton University
# This file is part of dPUC.
# dPUC 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.
# dPUC 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 for more details.
# You should have received a copy of the GNU General Public License along with dPUC.  If not, see <http://www.gnu.org/licenses/>.

package DpucNet;
our $VERSION = 2.03;

# 2015-06-05 10:08:49 EDT - v2.01 
# - changed output so now list of nodes is also explicitly given as the first line (to catch version mismatches, which are otherwise catastrophic for performance by turning unobserved context pairs into negative context pairs)
# 2015-07-04 17:04:50 EDT - v2.01 still (haven't published)
# - parsing proteins removes version number (i.e. .1 in BLAH.1). For Pfam 28 memory usage exploded so this extra step is warranted.
# 2015-12-23 10:11:17 EST - v2.02
# - Pfam 29 introduced changes that affect my "Pfam-A.full" parser.  The file I actually want to parse from this version onward is Pfam-A.full.uniprot, and the lines I have to parse are different (some annotations I used to rely on [#=GS ... AC ...] are not present in this file, but the data is present in another form [the subsequence alignment]).
# - also added some sorting steps at output time, so the same input always generates the same output (previously, some things were ordered randomly, and it was a pain to compare outputs)
# 2019-08-23 15:23:11 EDT - v2.03
# - makeNet verbose now reports path of file being processed


use lib '.';
use FileGz;
use strict;

# extract a directed network of domain family pair counts as observed in Pfam-A.full

# only public function, does all the magic
sub makeNet {
    # a sub that puts all these parts together, to get output straight from input
    my ($fi, $fo, $comp, $verbose) = @_;
    
    print "Parsing Pfam-A.full ($fi)...\n" if $verbose;
    my ($prot2hits, $accs) = getDomsPfamAFull($fi);
    
    print "Extracting family pair counts from domains...\n" if $verbose;
    my $net = countsPairFromDoms($prot2hits);
    
    print "Saving context count network...\n" if $verbose;
    my $fho = FileGz::getOutFh($fo, $comp);
    # first line is list of nodes!
    # note PF prefix has to be added back before printing
    print $fho join("\t", sort map { 'PF'.$_ } @$accs)."\n";
    # add network edges
#    while (my ($pair, $c) = each %$net) {
    foreach my $pair (sort keys %$net) {
	my $c = $net->{$pair};
	print $fho "$pair\t$c\n";
    }
    close $fho;
}

### INTERNAL FUNCTIONS

sub getDomsPfamAFull {
    # parse architectures from Pfam-A.full
    # this is very minimalistic to increase speed and reduce memory usage as much as possible (so format is very weird, sorry)
    # main annoyance that this solves is that predictions are grouped by domain fmaily accession first, but we want things grouped by protein (this is where the architectures come from)
    my ($fi) = @_;
    
    # some data to initialize
    my %prot2hits; # build predictions per protein
    my @accs; # also keep list of accessions, it'll be useful!
    
    # temporary variables
    my $acc; # current domain we're parsing
    
    my $fhi = FileGz::getInFh($fi);
    while (<$fhi>) {
	# we don't care about the majority of the lines of this huge file...
	if (/^\#=GF AC   PF(\d{5})/) {
	    # note this won't keep the Pfam version, or the initial PF (will save us a bit of memory)
	    $acc = $1;
	    push @accs, $acc; # add to list of unique accessions
	}
	# OLD version, doesn't work with Pfam 29's Pfam-A.full.uniprot
#	elsif (/^\#=GS \w+\/(\d+)\-\d+ +AC (\w+)\.\d+/) { # start,prot # note "end" not needed for order
#	    $prot2hits{$2} .= $1.'-'.$acc.';'; # store as string, to reduce memory usage
	#	}
	# NEW version, for Pfam 29's Pfam-A.full.uniprot and backward compatible with Pfam-A.full of older Pfams
	# parsing is much more strict (prot IDs must be numbers or capital letters, rather than \w)
	elsif (/^([A-Z0-9_]+)[^\/]*\/(\d+)\-\d+/) { # prot,start # note "end" not needed for order
	    $prot2hits{$1} .= $2.'-'.$acc.';'; # store as string, to reduce memory usage
	}
    }
    close $fhi;
    
    return(\%prot2hits, \@accs);
}

sub countsPairFromDoms {
    # compute raw pair counts, no filters or normalizations
    # input must be the $prot2hits generated by getDomsPfamAFull, which is a low-memory domain structure (a hash that points to strings) in a custom format that none of my other code uses!
    my ($prot2hits) = @_;
    
    my %pair2c; # domain fam pairs to counts, for both pairs and repeats
    
    # for each protein: order hits, count fam pairs!
    while (my ($prot, $hits) = each %$prot2hits) {
	# Part 1: parse stringified domains, sort by start to have correct order
	my @hits; # the list of domains
	chop $hits; # the hits always end with one slash, remove it before splitting in the foreach loop below
	foreach my $hitTxt (split /;/, $hits) { # split by semicolon
	    my ($start,$acc) = split /-/, $hitTxt; # split by dash
	    push @hits, [$start, $acc]; # this structure will be easiest to sort
	}
	next unless @hits > 1; # skip protein if there was only one domain!
	# sort domains, ascending numerically by start coordinate (stored in [0])
	@hits = sort { $a->[0] <=> $b->[0] } @hits;
	# now we're only keeping accession (stored in [1]), so throw away start coordinate
	# also add all the missing PF prefixes, now things are cool
	@hits = map { 'PF'.$_->[1] } @hits;
	
	# Part 2: count domain family pairs!
	my $numHits = scalar @hits; # to set loops
	for (my $i = 0; $i < $numHits; $i++) { # navigate list of domains in order of protein sequence
	    my $acci = $hits[$i]; # get accession
	    for (my $j = $i+1; $j < $numHits; $j++) { # navigate second domain, not including first domain again
		$pair2c{$acci."\t".$hits[$j]}++; # increment pair count
	    }
	}
    }
    
    # return network of directed domain family pair counts, with edges as strings here
    return \%pair2c;
}

1;
