#
# Vorlesung "Vorlesung mit \"Ubungen: Bioinformatik
#
# 4. \"Ubung -- Needleman-Wunsch.py
#
#

# Einbinden der Matrixklasse aus dem Numeric-Modul
#import Numeric
import numarray as Numeric
#import Numpy as Numeric
import sys

# Funktion zum Lesen der Sequenz
def getSequence(fastaName):
   # Variable, die die Sequenz speichert
   sequence = ''

   # oeffne die Datei
   myFasta = open(fastaName, 'r')

   # lese den Header, ohne ihn auszuwerten
   line = myFasta.readline();

   # lese nun die Sequenz-Dateien mit einer
   # while-Schleife ein
   while line != "":
      line = myFasta.readline()
      line = line.rstrip("\n")
      # Verketten der Strings
      sequence += line
   # Ende der while-Schleife

   # Schliesse die Fasta-Datei
   myFasta.close()

   # gebe die Sequenz zurueck
   return sequence
# Ende der Funktion zum Lesen der Sequenz

# Eine einfache Scorefunktion/Aehnlichkeitsfunktion: Z\"ahlen der Mismatches.
# Achtung: Reihenfolge der if-elif-else Abfragen wichtig, wenn man nicht "unnoetige" Abfragen
#          aufschreiben will, denn wenn ich erst "a == '-' sowie "b == '-' teste, habe ich das
#          Gap in einer der beiden Sequenzen schon abgefangen
def s(a,b): 
   if a == '-':
      return -1
   elif b == '-':
      return -1
   elif a != b:
#      if b != '-': # ist somit unnoetig (Erklaerung oben)
#         return 0
      return 0
   elif a == b:
      return 1

# Verbesserung: definiere Usage Methode (Funktion)
def usage(name):
   sys.stderr.write("\n")
   sys.stderr.write("Usage:" +  name + " FastaEingabeDatei1 FastaEingabeDatei2")
   sys.stderr.write("\n")
   sys.stderr.write("\n")

# Die zu alignierenden Sequenzen.
#A = 'ACGTATGTAC'
#B = 'GATTTTATACATAC'

# sys.argv[0] traegt den Namen des aufrufenden Skriptes
programName = sys.argv[0]

# len-Funktion gibt die Laenge eines Strings zurueck
lengthParam = len(sys.argv)

# werden 3 Argumente beim Aufruf angeben, d.h Programm-Name Sequenz1 Sequenz2?
if lengthParam != 3:
   usage(programName)
   sys.exit(1)

# Die Dateinamen werden eingelesen
fastaName1 = sys.argv[1] 
fastaName2 = sys.argv[2]

# oeffne und lese Sequenzen aus den Dateien
A = getSequence(fastaName1)
B = getSequence(fastaName2)

# Die L\"angen der Sequenzen
m = len(A)
n = len(B)

# Initialisiere die DP-Matrix und die Traceback-Matrix.
# Beide haben Dimension (m+1)x(n+1)
S = Numeric.zeros([m+1, n+1])
T = Numeric.zeros([m+1, n+1])

# Initialisiere die erste Zeile und erste Spalte
# der beiden Matrizen.
# ''vertikal" bekommt ''-1'', horizontal bekommt ''1''
for i in range(m+1):
   S[i, 0] = i * s("-", "A")
for j in range(n+1):
   S[0, j] = j * s("-", "A")
for i in range(m+1):
   T[i,0] = -1 
for j in range(n+1):
   T[0,j] =  1

# Dynamische Programmierung: f\"ulle DP-Matrix.
for i in range(1,m+1):
   for j in range(1,n+1):
      # Berechne Distanzen
      s_horiz = S[i-1,j] + s(A[i-1],"-")
      s_diag  = S[i-1,j-1] + s(A[i-1],B[j-1])
      s_vert  = S[i,j-1] + s("-", B[j-1])

      # Speichere Maximum in DP-Matrix S
      S[i,j] = max([s_horiz, s_diag, s_vert])

      # Merke Richtung in Traceback-Matrix T
      if s_diag >= s_horiz and s_diag >= s_vert:
         T[i,j] = 0
      elif s_horiz >= s_vert:
         T[i,j] = -1
      else:
         T[i,j] = +1

# Die alignierten Strings (mit Gaps)
# nennen wir X und Y (A', B' sind keine
# zu\"assigen Variablennamen).
X = ""
Y = ""

# Traceback: beginne in der unteren rechten
# Zelle der Matrix und laufe solange bis oben
# links angekommen.
i = m
j = n
while i > 0 or j > 0:
   if T[i,j] == 0:   # diagonal: keine Gaps
      X = A[i-1] + X
      Y = B[j-1] + Y
      (i,j) = (i-1, j-1)
   elif T[i,j] == -1: # nach links (horizontal): Gap in B
      X = A[i-1] + X
      Y = "-"    + Y
      i = i - 1
   else:              # nach oben (vertical): Gap in A
      X = "-"    + X
      Y = B[j-1] + Y
      j = j -1

# Gib das Alignment aus.
print X
print Y