from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import single_letter_alphabet, DNAAlphabet, IUPAC
from Bio import pairwise2
from Bio.SubsMat.MatrixInfo import blosum62
from Bio.pairwise2 import format_alignment



#create a sequence object
my_seq = Seq('CATGTAGACTAGCATGATACACGTAATCGTGGCTATTACTGGGATGGAGGTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACTGGCACCGCCGCGCCACCCGCGACCACGGCTGGTGGAAACAACATGAC')

my_seq2 = Seq('CATGTAGACTAGCATGATACACGTAATCGTGGCTAGTGTGTGTGTAATATGTAGTAGATAGATAAAAAAAAAACTGGCACCGCCGCGCCACCCGCGACCACGGCTGGTGGAAACAACATGAC')


my_dna = Seq("ATCTATCGTACTCGCACTTTCCCTGGTTCTGGTCGCTCCCA", IUPAC.ambiguous_dna)




alignments = pairwise2.align.globalxx(my_seq, my_seq2)




print(format_alignment(*alignments[0]))
for a in pairwise2.align.localds("CATGTAGACTAGCATGATACACGTAATCGTGGCTATTACTGGGATGGAGGTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACTGGCACCGCCGCGCCACCCGCGACCACGGCTGGTGGAAACAACATGAC", "CATGTAGACTAGCATGATACACGTAATCGTGGCTAGTGTGTGTGTAATATGTAGTAGATAGATAAAAAAAAAACTGGCACCGCCGCGCCACCCGCGACCACGGCTGGTGGAAACAACATGAC", blosum62, -10, -1):
    print(format_alignment(*a))





#print out some details about it


print('seq %s is %i bases long' % (my_seq, len(my_seq)))
print('reverse complement is %s' % my_seq.reverse_complement())
print('protein translation is %s' % my_seq.translate())


print(my_dna)