alignTranscripts1.0

#!/usr/bin/env python

import os
import sys
import argparse
import string
import tempfile
import subprocess
from math import floor

LASTZ = 'lastz'
BEDTOOLS = 'bedtools'
FASTAFROMBED = 'fastaFromBed'
SHUFFLEBED = 'shuffleBed'
STOP_CODON = [
	'TAA',
	'TAG',
	'TGA'
	]
def checkDependencies():
	#check bedtools version
	cmd = [BEDTOOLS, '--version']
	try:
		out = subprocess.check_output(cmd)
	except:
		sys.exit("ERROR: bedtools not installed! You must have bedtools v2.17.0 or higher in your path! Exiting...")

	out = out.split()
	version = out[1][1:].split(".")
	if int(version[0]) > 2 or int(version[1]) < 17:
		sys.exit("ERROR: The version of bedtools you have installed is too old. You must have bedtools v2.17.0 or higher in your path! Exiting...")

	#check lastz
	cmd = [LASTZ, '--version']
	try:
		p = subprocess.Popen(cmd, stderr=subprocess.PIPE, stdout=subprocess.PIPE)
	except:
		sys.exit("ERROR: lastz not installed! You must have lastz in your path! Exiting...")

	out = p.stdout.read()
	out = out.split()
	if out[0].strip() != "lastz":
		sys.exit("ERROR: lastz not installed! You must have lastz in your path! Exiting...")

def intArray(array):
	for i in range(len(array)):
		if array[i] != "": array[i] = int(array[i])
	return array

def convertToGeneCoords(base, exonSizes, exonStarts, numExons):
	block = 0
	for i in range(numExons):
		if base <= block + exonSizes[i]:
			newBase = exonStarts[i] + base - block
			break
		block += exonSizes[i]
	return newBase

def inExon(base, start, exonStarts, exonSizes, numExons, strand):
	#print strand
	for i in range(numExons):
		if (base > start + exonStarts[i] and base <= (start + exonStarts[i] + exonSizes[i])):
			
			#get position for splicing, 1 for first base in exon; -1 for last base in exon
			if strand == "+":
				posAbs = base - (start + exonStarts[i])
			else:
				posAbs = (start + exonStarts[i] + exonSizes[i]) - base + 1

			if (posAbs == 1): pos = 1
			elif (posAbs == exonSizes[i]): pos = -1
			else: pos = 0

			

			if strand == "+": return i+1, pos
			else: return numExons-i, pos

	return -1, 0

def inIntron(base, exonStarts, exonSizes, numExons, strand, geneSize):
	#print strand
	for i in range(numExons-1):
		if (base > (exonStarts[i] + exonSizes[i]) and base < exonStarts[i+1]):
			return exonStarts[i+1] - (exonStarts[i] + exonSizes[i])
	return -1


def writePaddedBed(gene, pad, genomeDict):
	chr = gene[0]

	gene[1] = max(0, int(gene[1])-pad)
	gene[2] = min(int(gene[2])+pad, genomeDict[chr])
		
	tmpFd, tmpPath = tempfile.mkstemp()
	tmp = os.fdopen(tmpFd, 'w')
	tmp.write("%s\t%i\t%i\t%s\t%s\t%s\n" % (gene[0], gene[1], gene[2], gene[3], gene[4], gene[5]))
	return tmpPath

def writeBed(gene):
	tmpFd, tmpPath = tempfile.mkstemp()
	tmp = os.fdopen(tmpFd, 'w')
	tmp.write("%s" % gene[0])
	for i in range(11):
		tmp.write("\t%s" % gene[i+1])
	tmp.write("\n")
	return tmpPath

def getFa(bed, fasta):
	tmpFd, tmpPath = tempfile.mkstemp()
	#cmd = "%s -fi %s -bed %s -fo %s -name" % (FASTAFROMBED, fasta, bed, tmpPath)
	cmd = [FASTAFROMBED, "-s", "-fi", fasta, "-bed", bed, "-fo", tmpPath, "-name"]
	subprocess.check_call(cmd)
	return tmpPath

def alignFa(paddedAFaPath, paddedBFaPath, unmask, gap_open, gap_extend):
	tmpFd, tmpPath = tempfile.mkstemp()
	if unmask:
		cmd = [LASTZ, '--strand=plus', paddedAFaPath+"[unmask]", paddedBFaPath+"[unmask]", '--format=maf', '--output=%s' % tmpPath, '--chain', '--gap=%d,%d' %(gap_open, gap_extend)]
	else:
		cmd = [LASTZ, '--strand=plus', paddedAFaPath, paddedBFaPath, '--format=maf', '--output=%s' % tmpPath, '--chain', '--gap=%d,%d' %(gap_open, gap_extend)]

	nullFd, nullPath = tempfile.mkstemp()
	null = open(nullPath, 'w')
	subprocess.check_call(cmd, stderr = null)
	null.close()
	cmd = ['rm', nullPath]
	subprocess.check_call(cmd)
	return tmpPath

code = {
	'TTT': 'Phe',
	'TTC': 'Phe',
	'TTA': 'Leu',
	'TTG': 'Leu',
	'TCT': 'Ser',
	'TCC': 'Ser',
	'TCA': 'Ser',
	'TCG': 'Ser',
	'TAT': 'Tyr',
	'TAC': 'Tyr',
	'TAA': 'Stop',
	'TAG': 'Stop',
	'TGT': 'Cys',
	'TGC': 'Cys',
	'TGA': 'Stop',
	'TGG': 'Trp',
	'CTT': 'Leu',
	'CTC': 'Leu',
	'CTA': 'Leu',
	'CTG': 'Leu',
	'CCT': 'Pro',
	'CCC': 'Pro',
	'CCA': 'Pro',
	'CCG': 'Pro',
	'CAT': 'His',
	'CAC': 'His',
	'CAA': 'Gln',
	'CAG': 'Gln',
	'CGT': 'Arg',
	'CGC': 'Arg',
	'CGA': 'Arg',
	'CGG': 'Arg',
	'ATT': 'Ile',
	'ATC': 'Ile',
	'ATA': 'Ile',
	'ATG': 'Met',
	'ACT': 'Thr',
	'ACC': 'Thr',
	'ACA': 'Thr',
	'ACG': 'Thr',
	'AAT': 'Asn',
	'AAC': 'Asn',
	'AAA': 'Lys',
	'AAG': 'Lys',
	'AGT': 'Ser',
	'AGC': 'Ser',
	'AGA': 'Arg',
	'AGG': 'Arg',
	'GTT': 'Val',
	'GTC': 'Val',
	'GTA': 'Val',
	'GTG': 'Val',
	'GCT': 'Ala',
	'GCC': 'Ala',
	'GCA': 'Ala',
	'GCG': 'Ala',
	'GAT': 'Asp',
	'GAC': 'Asp',
	'GAA': 'Glu',
	'GAG': 'Glu',
	'GGT': 'Gly',
	'GGC': 'Gly',
	'GGA': 'Gly',
	'GGG': 'Gly'
}

def kNkS(seqA, seqB):
	totalN = 0
	totalS = 0
	for i in range(3, len(seqB)-3):
		if i%3 == 0:
			codonA = seqA[i:i+3].upper()
			codonB = seqB[i:i+3].upper()
		elif i%3 == 1: 
			codonA = seqA[i-1:i+2].upper()
			codonB = seqB[i-1:i+2].upper()
		else: 
			codonA = seqA[i-2:i+1].upper()
			codonB = seqB[i-2:i+1].upper()

		if codonA.find("-") != -1 or codonB.find("-") != -1: continue
		
		for j in ['A', 'T', 'C', 'G']:
			if j == seqB[i].strip(): continue
			newcodon = list(codonB)
			newcodon[(i%3)] = j
			newcodon = "".join(newcodon)

			if code[codonB] == code[newcodon]: totalS += 1
			else: totalN += 1
	
	kN = 0
	kS = 0
	for i in range(3, len(seqB)-3):
		if seqB[i] != seqA[i]:
			if i%3 == 0:
				codonA = seqA[i:i+3].upper()
				codonB = seqB[i:i+3].upper()
			elif i%3 == 1: 
				codonA = seqA[i-1:i+2].upper()
				codonB = seqB[i-1:i+2].upper()
			else: 
				codonA = seqA[i-2:i+1].upper()
				codonB = seqB[i-2:i+1].upper()

			if codonA.find("-") != -1 or codonB.find("-") != -1: continue

			if code[codonA] == code[codonB]: kS += 1
			else: kN += 1
	
	kN = kN * 1.0 / totalN
	kS = kS * 1.0 / totalS
	
	return kN, kS

def checkForOrfs(stringA, stringB, baseA, baseB, pad, geneAStart, geneAEnd, exonStartsA, exonSizesA, numExonsA, strandA, geneBStart, geneBEnd, exonStartsB, exonSizesB, numExonsB, strandB):
	orfs = []

	starts = []
	startBaseA = []
	startBaseB = []
	
	if strandA == "+": baseAGen = geneAStart + baseA - pad + 1
	else: baseAGen = geneAEnd - (baseA - pad) 

	if strandB == "+": baseBGen = geneBStart - pad + baseB + 1
	else: baseBGen = geneBEnd - (baseB - pad)
	
	#get all possible start codons
	for i in range(len(stringA)):
		exonA, pos = inExon(baseAGen, geneAStart, exonStartsA, exonSizesA, numExonsA, strandA) 
		exonB, pos = inExon(baseBGen, geneBStart, exonStartsB, exonSizesB, numExonsB, strandB) 

		
		if exonA != -1 and exonB != -1:
			if stringA[i:i+3].upper() == "ATG" and stringB[i:i+3].upper() == "ATG":
				#print baseAGen
				starts.append(i)
				startBaseA.append(baseAGen)
				startBaseB.append(baseBGen)

		if stringA[i] !="-": 
			if strandA == "+": baseAGen += 1
			else: baseAGen -= 1

		if stringB[i] !="-": 
			if strandB == "+": baseBGen += 1
			else: baseBGen -= 1

	for i in range(len(starts)):
		valid = True
		baseAGen = startBaseA[i]
		baseBGen = startBaseB[i]
		orfA = ""
		orfB = ""

		#extract orf sequences
		for j in range(starts[i], len(stringA)):

			exonA, pos = inExon(baseAGen, geneAStart, exonStartsA, exonSizesA, numExonsA, strandA) 
			exonB, pos = inExon(baseBGen, geneBStart, exonStartsB, exonSizesB, numExonsB, strandB)
		
			if exonA != -1 and exonB != -1: 
				orfA += stringA[j]
				orfB += stringB[j]
			
			if stringA[j] !="-": 
				if strandA == "+": baseAGen += 1
				else: baseAGen -= 1

			if stringB[j] !="-":
				if strandB == "+" : baseBGen += 1
				else: baseBGen -= 1

		
		#truncate at first stop codon

		truncOrfA = ""
		truncOrfB = ""
		
		
		for j in range(min(len(orfA), len(orfB))):

			if len(truncOrfA.replace("-", "")) % 3 == 0 and len(truncOrfB.replace("-", "")) %3 == 0 and orfA[j:j+3] in STOP_CODON and orfB[j:j+3] in STOP_CODON:
					truncOrfA += orfA[j:j+3]
					truncOrfB += orfB[j:j+3]
					break
			else:
				truncOrfA += orfA[j]
				truncOrfB += orfB[j]

		if truncOrfA[-3:] not in STOP_CODON and truncOrfB[-3:] not in STOP_CODON: valid = False

	
		#check indels
		if valid:
			indel = False
			indelSize = 0
			for j in range(len(truncOrfA)):
				if truncOrfA[j] == "-":
					if not indel: indel = True
					indelSize += 1
				else:
					if indel: indel = False
					if indelSize %3 != 0: 
						valid = False
						break
			indel = False
			indelSize = 0
			for j in range(len(truncOrfB)):
				if truncOrfB[j] == "-":
					if not indel: indel = True
					indelSize += 1
				else:
					if indel: indel = False
					if indelSize %3 != 0: 
						valid = False
						break
		
		if valid and len(truncOrfA.replace("-", "")) > 29 and len(truncOrfB.replace("-", ""))> 29:
			kN, kS = kNkS(truncOrfA, truncOrfB)
			
			orfs.append([truncOrfA, truncOrfB, baseAGen, baseAGen+len(truncOrfA.replace("-", "")), len(truncOrfA.replace("-", "")), kN, kS])
	
	sorted(orfs, key=lambda x: x[4])
	return orfs
def main():
	parser = argparse.ArgumentParser(description='''
		Wrapper script for lastz to align two transcripts and report back exonic and sequence identity.\n
		If aligning non-coding genes, set --noncoding flag for more sensitive alignment.
	''')
	parser.add_argument('bedA', type=file,  help='bed file A')
	parser.add_argument('genomeFastaA', type=str, help='fasta file A')
	parser.add_argument('bedB', type=file, help='bed file B')
	parser.add_argument('genomeFastaB', type=str, help='fasta file B')
	parser.add_argument('out_prefix', type=str, help='out prefix')
	parser.add_argument('--geneA', type=str, help='if bedA has more than one entry, specify which gene to align')
	parser.add_argument('--geneB', type=str, help='if bedB has more than one entry, specify which gene to align')
	parser.add_argument('--pad', type=int, default=15000)
	parser.add_argument('--gap_open', type=int, default=200)
	parser.add_argument('--gap_extend', type=int, default=40)
	parser.add_argument('--unmask', action='store_true', help='unmask repeats when aligning')
	parser.add_argument('--orf', action='store_true', help='flag for checking for orfs in alignment')
	parser.add_argument('--bedtools_path', type=str)
	parser.add_argument('--lastz_path', type=str)
	parser.add_argument('--shuffle_bg', action='store_true', help='Aligns genes to random background and removes nonsignificant alignments')
	parser.add_argument('--pvalue', type=float, default=0.05)
	args = parser.parse_args()

	if args.bedtools_path is not None:
		global BEDTOOLS
		global FASTAFROMBED
		global SHUFFLEBED
		BEDTOOLS  = args.bedtools_path+"/bedtools"
		FASTAFROMBED = args.bedtools_path+"/fastaFromBed"
		SHUFFLEBED = args.bedtools_path+"shuffleBed"

	if args.lastz_path is not None:
		global LASTZ
		LASTZ = args.lastz_path+"/lastz"
	
	checkDependencies()

	genomeA = {}
	if os.path.exists(args.genomeFastaA+".fai"):
		genome = open(args.genomeFastaA+".fai", 'r')
		for line in genome.readlines():
			line = line.split()
			genomeA[line[0].strip()] = int(line[1])
	else:
		print "Fasta index for %s does not exist.  Please run samtools faidx %s. Exiting..." % (args.genomeFastaA, args.genomeFastaA)
		sys.exit(1)

	genomeB = {}
	if os.path.exists(args.genomeFastaB+".fai"):
		genome = open(args.genomeFastaB+".fai", 'r')
		for line in genome.readlines():
			line = line.split()
			genomeB[line[0].strip()] = int(line[1])
	else:
		print "Fasta index for %s does not exist.  Please run samtools faidx %s. Exiting..." % (args.genomeFastaB, args.genomeFastaB)
		sys.exit(1)

	#read in bed file A into array geneA
	geneA = []
	if args.geneA is None:
		geneA = args.bedA.readline().split()
	else:
		for line in args.bedA.readlines():
			if line.strip() == "": continue
			line = line.split()
			if (line[3].strip() == args.geneA.strip()): 
				geneA = line
				break
		if (len(geneA)==0):
			print "geneA %s not found! exiting..." % args.geneA
			sys.exit(1)

	#read in bed file B
	geneB = []
	if args.geneB is None:
		geneB = args.bedB.readline().split()
	else:
		for line in args.bedB.readlines():
			if line.strip() == "": continue
			line = line.split()
			if (line[3].strip() == args.geneB.strip()): 
				geneB = line
				break
		if (len(geneB)==0):
			print "geneB %s not found! exiting..." % args.geneB
			sys.exit(1)

	paddedABedPath = writePaddedBed(geneA, args.pad, genomeA)
	paddedBBedPath = writePaddedBed(geneB, args.pad, genomeB)

	paddedAFaPath = getFa(paddedABedPath, args.genomeFastaA)
	paddedBFaPath = getFa(paddedBBedPath, args.genomeFastaB)

	#ALIGN GENES OF INTEREST

	alignments = []
	mafs = []
	mafHeader = []
	finalOrfs = []

	mafPath = alignFa(paddedAFaPath, paddedBFaPath, args.unmask, args.gap_open, args.gap_extend)
	mafFile = open(mafPath, 'r')
	maf = mafFile.readlines()
	mafFile.close()



	if len(maf) > 14:
		chrA = geneA[0].strip()
		geneAStart = int(geneA[1])+args.pad
		geneAEnd = int(geneA[2])-args.pad
		geneAName = geneA[3]
		strandA = geneA[5].strip()
		exonStartsA = intArray(geneA[11].split(','))
		exonSizesA = intArray(geneA[10].split(','))
		numExonsA = int(geneA[9])
		geneASize = geneAEnd - geneAStart
		exonASize = 0
		for i in range(numExonsA):
			exonASize += exonSizesA[i]
		intronASize = 0
		for i in range(numExonsA):
			intronASize += exonStartsA[i]

		chrB = geneB[0].strip()
		geneBStart = int(geneB[1])+args.pad
		geneBEnd = int(geneB[2])-args.pad
		geneBName = geneB[3]
		strandB = geneB[5].strip()
		if strandB == "*": strandB = "+"
		exonStartsB = intArray(geneB[11].split(','))
		exonSizesB = intArray(geneB[10].split(','))
		numExonsB = int(geneB[9])
		geneBSize = geneBEnd - geneBStart
		exonBSize = 0
		for i in range(numExonsB):
			exonBSize += exonSizesB[i]
		
		#CRAWL THROUGH ALIGNMENTS
		
		id = 0
		exonId = 0
		intronId = 0
		maxIntronId = 0
		maxIntronNum = 0
		stringB = ""
		stringA = ""
		baseB = -1
		baseA = -1
		truncA = ""
		truncB = ""
		truncAStart = -1
		truncBStart = -1
		exonsAlignedA = []
		exonsAlignedB = []
		spliceStart = [0]*(numExonsA-1)
		spliceEnd = [0]*(numExonsA-1)
		counter = 0
		indelRate = 0
		indelRateIntron = 0
		inIndel = False
		inIndelIntron = False
		seen = 0
		seenIntron = 0
		for line in maf:
			line = line.strip()
			if line.startswith("#"): 
				mafHeader.append(line)

			if line.startswith("a"):
				counter = counter + 1
				a = line
				score = int(float(line.split("=")[1]))
			if line.startswith("s"):
				if stringA == "":
					sA = line.split()
					stringA = sA[6].strip()
					baseA = int(sA[2])
				else:
					sB = line.split()
					stringB = sB[6].strip()
					baseB = int(sB[2])

			#at line break, reset
			if line == "" and stringA !="":
				if seen == 0: indelRateNorm = 'NA'
				else: indelRateNorm = indelRate / (seen * 1.0)
				
				if seenIntron == 0: indelRateIntronNorm = 'NA'
				else: indelRateIntronNorm = indelRateIntron / (seenIntron * 1.0)
			
				#print counter	
				#print spliceStart
				#print spliceEnd
				
				numSpliceConserved = 0.0
				for i in range(numExonsA-1):
					if spliceStart[i] == spliceEnd[i] - 1: numSpliceConserved += 1
					elif spliceStart[i]> 0 or spliceEnd[i] > 0: numSpliceConserved += 0.5
				#print "numSpliceConserved: ", numSpliceConserved
				alignments.append([counter, geneAName, geneBName, exonId / (exonASize*1.0), id / (geneASize*1.0), exonId / (exonBSize*1.0), id / (geneBSize*1.0), indelRateNorm, indelRateIntronNorm, exonsAlignedA, exonsAlignedB, score, numSpliceConserved])
					

				#write maf
				if (strandA == "-"): truncAStart = genomeA[chrA] - truncAStart
				if (strandB == "-"): truncBStart = genomeB[chrB] - truncBStart
						
						
				mafStr = a+"\n"
				mafStr += "s %s %d %d %s %d %s\n" % (chrA, truncAStart, len(truncA.replace("-", "")), strandA, genomeA[chrA], truncA)
				mafStr += "s %s %d %d %s %d %s\n" % (chrB, truncBStart, len(truncB.replace("-", "")), strandB, genomeB[chrB], truncB)
				mafStr += ("\n")	
				mafs.append(mafStr)
					

				#check for small orfs
				if args.orf and exonId != 0:
					orfs = checkForOrfs(stringA, stringB, baseA, baseB, args.pad, geneAStart, geneAEnd, exonStartsA, exonSizesA, numExonsA, strandA, geneBStart, geneBEnd, exonStartsB, exonSizesB, numExonsB, strandB)
					for i in range(len(orfs)):
						eclipsed = False
						for j in range(i):
							if orfs[i][2] >= orfs[j][2] and orfs[i][2] <= orfs[j][3] and orfs[i][2] % 3 == orfs[j][2] % 3:
								eclipsed = True
								break
						if not eclipsed:
							curOrf = orfs[i]
							curOrf.append(score)
							curOrf.append(geneAName)
							curOrf.append(geneBName)
							finalOrfs.append(curOrf)

				#reset stats
				exonId = 0
				intronId = 0
				maxIntronId = 0
				maxIntronNum = 0
				id = 0
				score = 0
				stringB = ""
				stringA = ""
				truncA = ""
				truncB = ""
				truncAStart = -1
				truncBStart = -1
				exonsAlignedA = []
				exonsAlignedB = []
				spliceStart = [0]*(numExonsA-1)
				spliceEnd = [0]*(numExonsA-1)
				indelRate = 0
				indelRateIntron = 0
				inIndel = False	
				inIndelIntron = False
				seen = 0
				seenIntron = 0

			#print "*****"
			#if stringA and stringB populated, crawl alignments
			if stringB != "" and stringA !="":
				if strandA == "+": baseAGen = geneAStart + baseA - args.pad + 1

				else: baseAGen = geneAEnd - (baseA - args.pad) 

				if strandB == "+": baseBGen = geneBStart - args.pad + baseB + 1
				else: baseBGen = geneBEnd - (baseB - args.pad)
				
				for i in range(len(stringA)):
					
					if baseAGen >= geneAStart - 1000 and baseAGen <= geneAEnd+1000 and baseBGen >=geneBStart-1000 and baseBGen <= geneBEnd+1000:
						if truncA == "":
							truncAStart = baseAGen
							truncBStart = baseBGen
						truncA += stringA[i]
						truncB += stringB[i]


					exonA, exonAPos = inExon(baseAGen, geneAStart, exonStartsA, exonSizesA, numExonsA, strandA) 
					exonB, exonBPos = inExon(baseBGen, geneBStart, exonStartsB, exonSizesB, numExonsB, strandB)
					
					#if (exonAPos == -1 or exonAPos == 1): print exonA, exonAPos, stringA[i]
					#splice start
					if (exonAPos == -1 and exonBPos == -1): 
						#print (exonA-0), "splice start", exonAPos, exonBPos, stringA[i-20:i+20], stringB[i-20:i+20]
						if (exonA < numExonsA): spliceStart[exonA-1] = exonB
					
					if (exonAPos == 1 and exonBPos == 1): 
						#print (exonA-1), "splice end", exonAPos, exonBPos, stringA[i-20:i+20], stringB[i-20:i+20]
						if (exonA > 1): spliceEnd[exonA-2] = exonB

					#if its inside gene
					if baseAGen >= geneAStart and baseAGen < geneAEnd and baseBGen >= geneBStart and baseBGen < geneBEnd:
						if stringB[i].upper()==stringA[i].upper():
							id += 1

						#if it's inside exons
						if exonA != -1 and exonB != -1:
							seen += 1
							
							#reset inIndelIntron
							if inIndelIntron: inIndelIntron = False

							#if in exons and '-' at alignment
							if stringA[i] == "-" or stringB[i] == "-":
								if not inIndel:
									inIndel = True
									indelRate += 1
							else:
								inIndel = False
							#if stringA == stringB, increase exonID
							if stringB[i].upper()==stringA[i].upper():
								exonId += 1
								if exonA not in exonsAlignedA: exonsAlignedA.append(exonA)
								if exonB not in exonsAlignedB: exonsAlignedB.append(exonB)


						#else if inside intron-intron alignment
						elif exonA == -1 and exonB == -1:
							seenIntron += 1

							#reset inIndel(exons)
							if inIndel: inIndel = False
							
							#if in introns and '-' at alignment
							if stringA[i] == "-" or stringB[i] == "-":
								if not inIndelIntron:
									inIndelIntron = True
									indelRateIntron += 1
							else:
								inIndelIntron = False
						
						#else if at exon/intron alignment, just reset inIndel
						else:
							if inIndel: inIndel = False
							if inIndelIntron: inIndelIntron = False

					if stringA[i] !="-": 
						if strandA == "+": baseAGen += 1
						else: baseAGen -= 1

					if stringB[i] !="-":
						if strandB == "+" : baseBGen += 1
						else: baseBGen -= 1
	
	#if len(alignments)==0 : print "no alignment!"
	#else: print alignments

	#MAKE BACKGROUND
	if args.shuffle_bg and len(alignments)>0:
		shuffleScores = []
		bedB = writeBed(geneB)
		for i in range(25):
			#shuffle geneA
			cmd = [SHUFFLEBED, '-i', bedB, '-g', args.genomeFastaB+".fai"]
			shuffle = subprocess.check_output(cmd)
			shuffle = shuffle.split()
			paddedShuffleBed = writePaddedBed(shuffle, args.pad, genomeB)
			paddedShuffleFa = getFa(paddedShuffleBed, args.genomeFastaB)
			
			#align geneA to shuffled
			shuffleMafPath = alignFa(paddedAFaPath, paddedShuffleFa, args.unmask, args.gap_open, args.gap_extend)
			mafFile = open(shuffleMafPath, 'r')
			maf = mafFile.readlines()
			if len(maf) <=14: 
				shuffleScores.append(0)
			else:
				for line in maf:
					if line.startswith("a"): 
						shuffleScores.append(int(line.split("=")[1]))
			cmd = ['rm', shuffleMafPath, paddedShuffleBed, paddedShuffleFa]
			subprocess.check_call(cmd)

		bedA = writeBed(geneA)
		for i in range(25):
			#shuffle geneB
			cmd = [SHUFFLEBED, '-i', bedA, '-g', args.genomeFastaA+".fai"]
			shuffle = subprocess.check_output(cmd)
			shuffle = shuffle.split()
			paddedShuffleBed = writePaddedBed(shuffle, args.pad, genomeA)
			paddedShuffleFa = getFa(paddedShuffleBed, args.genomeFastaA)
			
			#align geneB to shuffled
			shuffleMafPath = alignFa(paddedBFaPath, paddedShuffleFa, args.unmask, args.gap_open, args.gap_extend)
			mafFile = open(shuffleMafPath, 'r')
			maf = mafFile.readlines()
			if len(maf) <=1: shuffleScores.append(0)
			else:
				for line in maf:
					if line.startswith("a"): shuffleScores.append(int(line.split("=")[1]))
			cmd = ['rm', shuffleMafPath, paddedShuffleBed, paddedShuffleFa]
			subprocess.check_call(cmd)

		cmd = ['rm', bedB, bedA]
		subprocess.check_call(cmd)
		
		shuffleScores.sort()
		maxShuffleScore = shuffleScores[int(floor(len(shuffleScores)* (1.0-args.pvalue)))]
		maxScore = maxShuffleScore
		for a in alignments:
			if a[11] > maxScore: 
				maxScore = a[11]
				break
	
	cmd = ['rm', paddedABedPath, paddedBBedPath, paddedAFaPath, paddedBFaPath, mafPath]
	subprocess.check_call(cmd)	
	
	if len(alignments) > 0 and (not args.shuffle_bg or maxScore > maxShuffleScore):
		alignment = open(args.out_prefix+".alignment_identity.txt", 'w')
		mafout = open(args.out_prefix+".maf", 'w')
		if len(finalOrfs) > 0: 
			orfout = open(args.out_prefix+".orfs.txt", 'w')
			orfout.write("#alignmentID\tgeneA\tgeneB\tlengthOrfA\tlengthOrfB\tkN\tkS\tkN/kS\torfA\torfB\n")

			for orf in finalOrfs:
				if orf[6] != 0: orfout.write("%d\t%s\t%s\t%d\t%d\t%.2f\t%.2f\t%.3f\t%s\t%s\n" % (orf[7], orf[8], orf[9], len(orf[0].replace("-", "")), len(orf[1].replace("-", "")), orf[5], orf[6], orf[5] / orf[6], orf[0].upper(), orf[1].upper()))
				else: orfout.write("%d\t%s\t%s\t%d\t%d\t%.2f\t%.2f\tinf\t%s\t%s\n" % (orf[7], orf[8], orf[9], len(orf[0].replace("-", "")), len(orf[1].replace("-", "")), orf[5], orf[6],orf[0].upper(), orf[1].upper()))
		
		for line in mafHeader:
			mafout.write(line+"\n")
		

		alignment.write("#alignmentID\talignmentScore\tgeneA\tgeneB\texonID_A\tgeneID_A\texonID_B\tgeneID_B\tindelRate(exons)\tindelRate(introns)\texonsAlignedA\texonsAlignedB\tnumSpliceSitesAligned\tnumTotalSpliceSites\n")
		counter = -1
		

		for i in range(len(alignments)):
			alignments[i].append(mafs[i])
		alignments.sort(key = lambda x: (x[3], x[4], x[10]), reverse=True)
		
		bestExonID = alignments[0][3]
		bestLociID = alignments[0][4]
		written = False
		for a in alignments:
			counter += 1
			if (not args.shuffle_bg or a[11] > maxShuffleScore):
				if (a[4] == 0.0 and not written):
					written = True
					alignment.write("%d\t%d\t%s\t%s\t%0.3f\t%0.3f\t%0.3f\t%0.3f\t" % (a[0],a[11], a[1], a[2], a[3], a[4], a[5], a[6]))
					
					#write indelRate
					if str(a[7]) != "NA": alignment.write("%0.3f\t" % a[7])
					else: alignment.write("NA\t")

					if str(a[8]) != "NA": alignment.write("%0.3f\t" % a[8])
					else: alignment.write("NA\t")

					a[9].sort()
					a[10].sort()
					if len(a[9]) > 0:
						for exon in a[9]:
							alignment.write("%d," % exon)
					else:
						alignment.write("NA")
					alignment.write("\t")
					if len(a[10]) > 0:
						for exon in a[10]:
							alignment.write("%d," % exon)
					else:
						alignment.write("NA")

					alignment.write("\t%.1f\t%d" % (a[12], min(numExonsA, numExonsB)-1))
					alignment.write("\n")
					if (float(a[4] > 0.0)):
						mafout.write(a[13])
				elif (a[4] != 0.0):
					if ((bestExonID > 0.0 and a[3] > 0.0) or (bestExonID == 0)):
						written = True
						alignment.write("%d\t%d\t%s\t%s\t%0.3f\t%0.3f\t%0.3f\t%0.3f\t" % (a[0], a[11], a[1], a[2], a[3], a[4], a[5], a[6]))
					
						#write indelRate
						if str(a[7]) != "NA": alignment.write("%0.3f\t" % a[7])
						else: alignment.write("NA\t")

						if str(a[8]) != "NA": alignment.write("%0.3f\t" % a[8])
						else: alignment.write("NA\t")
	
						a[9].sort()
						a[10].sort()
						if len(a[9]) > 0:
							for exon in a[9]:
								alignment.write("%d," % exon)
						else:
							alignment.write("NA")
						alignment.write("\t")
						if len(a[10]) > 0:
							for exon in a[10]:
								alignment.write("%d," % exon)
						else:
							alignment.write("NA")

						alignment.write("\t%.1f\t%d" % (a[12], min(numExonsA, numExonsB)-1))
						alignment.write("\n")
						if (float(a[4] > 0.0)):
							mafout.write(a[13])

		alignment.close()

if __name__ == "__main__":
	main()