from Bio import motifs
from Bio import SeqIO 
from Bio.Seq import Seq
from Bio.Alphabet import IUPAC, generic_dna, generic_protein
import pandas as pd
import os, sys


"""
Explanation: Locating TFBS
I will be using Biopython's Bio.motifs package.  I cannot get the python 3 kernel working with the biopython package.  This needs to be fixed.
End of explanation
"""


for seq_record in SeqIO.parse("../data/fasta/ludwig_eve-striped-2.fasta", "fasta"):
    print(seq_record.id)
    print(repr(seq_record.seq))
    print(len(seq_record))

records = list(SeqIO.parse("../data/fasta/ludwig_eve-striped-2.fasta", "fasta"))

# Use only the first record
seqTest = records[0]
print(seqTest.seq)
print(type(seqTest))

str(seqTest.id) 

# Why can't I put this first? Because I am writing over?
# Why do I have to put str in front of Seq
seqTest = Seq(str(seqTest.seq), IUPAC.IUPACUnambiguousDNA())
print(seqTest)
print(type(seqTest))

"""
Explanation: Sequences to Query
Reading in ES2 sequences from 11 species to test.
End of explanation
"""


bcd = motifs.read(open("../data/PWM/transpose_fm/bcd_FlyReg.fm"),"pfm")

print(bcd.counts)

# If you want to normalize, but has to be based off of actual sequences.
# the problem is that you have to normalize to get pwm and you need a pwm object 
# to get the Position Specific scoring Matrix. 
# Therefore I am just normalizing to zero.
# I may need to adjust psuedocounts for what the sequence really is or to the Drosophila genome??

pwm = bcd.counts.normalize(pseudocounts=0.0)
print(pwm)
print(type(pwm))
print(vars(pwm))
dir(pwm)
sys.getsizeof(pwm)


"""
Explanation: Motifs
End of explanation
"""


pssm = pwm.log_odds()

print(id(pssm.alphabet)) # id() gives a unique hash.  Good to see when data changes.
print(id(pssm))

print(pssm)
print(type(pssm))
print(pssm.alphabet)
#pssm #If you ask for the raw 

"""
Explanation: Postion Specific scoring matrix
pssm is a specific type of object class. More info here.
End of explanation
"""


# Forward
pssm.calculate(seqTest) 

# Reverse
rpssm = pssm.reverse_complement()
# rpssm.calculate(seqTest) #You can see the calculations

print("Background: %f" % bcd.pssm.mean(bcd.background))

distribution = pssm.distribution(background=bcd.background, precision=10**4)
#print(type(distribution)) #check
#print(distribution) #check
dir(distribution) #there is a threshold for patser. 

# that is what we should use. 

threshold = distribution.threshold_patser()
print("Patser Threshold %5.3f" % threshold)

for position, score in pssm.search(seqTest, threshold=6):
    print("Position %d: score = %5.3f" % (position, score))

"""
Explanation: Searching for TFBS
Selecting Threshold
Make sure that the alphabet of both pssm and and seqTest are the same. It won't work if they are different.
End of explanation
"""


#initiate 
allPos = []
print(seqTest)
#how do I output a vector?
for position, score in pssm.search(seqTest, threshold=3.262): 
    allPos.append(position)

print(type(allPos))

#length of motif
m = 7
print(len(seqTest))

print(seqTest[0:10]) # first 10 nucleotides
print(seqTest[885:898]) # last 10 nucleotides
print(seqTest[-13:-1]) # with negative numbers, which go from last position. 
# But you cannot use zero. How would you identify the last nucleotide?

print(seqTest[-1]) #first, can be called here though?
print(seqTest[-2])
print(seqTest[-3])

for pos in allPos:
    print(seqTest[pos:pos+m]) # this will find and print the motifs

"""
Explanation: Test Motifs present
We need to output a vector of all the positions from the search function, so we can can output all the motifs found for testing if the search threshold is valuable.
End of explanation
"""


