#!/usr/bin/env python3

import argparse
import sys
import json

import grimoire.genome as genome
from grimoire.feature import FeatureTable

## Command line stuff ##

extended_help = """

%(prog)s is used for comparing gene sequence features between two
annotation sources. In a gene prediction setting this would be
predictions vs. reference, but it could also be two genome annotation
pipelines or two versions produced at different times. All of the
features must reference the same exact same chromosome sequences.

%(prog)s operates under some very simplistic assumptions:

	1. Each nt of the sequence can be uniquely labeled
	2. Strand doesn't matter

In cases where a nt has more than one labeling due to multiple isoforms,
priority is given in the order of CDS > exon > intron > intergenic.

In 'cds' mode, only CDS features are considered.
"""

parser = argparse.ArgumentParser(
	description='Compares two sets of genome annotations.',
	formatter_class=argparse.RawDescriptionHelpFormatter,
	epilog=extended_help)
parser.add_argument('--fasta', required=True, type=str,
	metavar='<path>', help='path to fasta file')
parser.add_argument('--file1', required=True, type=str,
	metavar='<path>', help='path to file1 in gff')
parser.add_argument('--file2', required=True, type=str,
	metavar='<path>', help='path to file2 in gff')
parser.add_argument('--cds', action='store_true',
	help='sets the comparision to CDS only')
arg = parser.parse_args()

class MorlisError(Exception):
	pass

def label_sequence(chrom):
	genes = chrom.ftable.build_genes()
	aseq = []
	for i in range(len(chrom.seq)):
		aseq.append('n')
		
	for gene in genes:
		for tx in gene.transcripts():
			for cds in tx.cdss:
				for i in range(cds.beg, cds.end +1):
					aseq[i-1] = 'c'
			if arg.cds: continue
			for intron in tx.introns:
				for i in range(intron.beg, intron.end +1):
					if aseq[i-1] == 'n':
						aseq[i-1] = 'i'
			for exon in tx.exons:
				for i in range(exon.beg, exon.end +1):
					if aseq[i-1] == 'n' or aseq[i-1] == 'i':
						aseq[i-1] = 'e'
	return aseq

if __name__ == '__main__':

	g1 = genome.Reader(fasta=arg.fasta, gff=arg.file1)
	g2 = genome.Reader(fasta=arg.fasta, gff=arg.file2)
	
	same, diff = 0, 0
	for c1, c2 in zip(g1, g2):
		s1 = label_sequence(c1)
		s2 = label_sequence(c2)
		s, d = 0, 0
		for i in range(len(s1)):
			if s1[i] == s2[i]: s += 1
			else:              d += 1
	#	sys.stderr.write('{} {}\n'.format(s, d))
		same += s
		diff += d
	
	print('{}\t{}\t{:.3f}'.format(same, diff, same / (same + diff)))
	
	

