from Bio.Seq import Seq
my_seq = Seq("GATCG")
for index, letter in enumerate(my_seq):
    print("%i %s" % (index, letter))
print(len(my_seq))

"""
Explanation: Source of the materials: Biopython Tutorial and Cookbook (adapted)
Sequence Objects
Biological sequences are arguably the central object in Bioinformatics, and in this chapter we’ll introduce the Biopython mechanism for dealing with sequences, the Seq object. Chapter 4 will introduce the related SeqRecord object, which combines the sequence information with any annotation, used again in Chapter 5 for Sequence Input/Output.
Sequences are essentially strings of letters like AGTACACTGGT, which seems very natural since this is the
most common way that sequences are seen in biological file formats.
The most important difference between Seq objects and standard Python strings is they have different methods. Although the Seq object supports many of the same methods as a plain string, its translate() method differs by doing biological translation, and there are also additional biologically relevant methods like reverse_complement().
Sequences act like strings
In many ways, we can deal with Seq objects as if they were normal Python strings, for example getting the length, or iterating over the elements:
End of explanation
"""


print(my_seq[0]) #first letter
print(my_seq[2]) #third letter
print(my_seq[-1]) #last letter

"""
Explanation: You can access elements of the sequence in the same way as for strings (but remember, Python counts from zero!):
End of explanation
"""


print("AAAA".count("AA"))
print(Seq("AAAA").count("AA"))

"""
Explanation: The Seq object has a .count() method, just like a string.
Note that this means that like a Python string, this gives a
non-overlapping count:
End of explanation
"""


my_seq = Seq('GATCGATGGGCCTATATAGGATCGAAAATCGC')
print(len(my_seq))
print(my_seq.count("G"))
print(100 * float(my_seq.count("G") + my_seq.count("C")) / len(my_seq))

"""
Explanation: For some biological uses, you may actually want an overlapping count
(i.e. $3$ in this trivial example). When searching for single letters, this
makes no difference:
End of explanation
"""


 from Bio.SeqUtils import GC

my_seq = Seq('GATCGATGGGCCTATATAGGATCGAAAATCGC')
GC(my_seq)

"""
Explanation: While you could use the above snippet of code to calculate a GC\%, note that  the Bio.SeqUtils module has several GC functions already built.  For example:
End of explanation
"""


my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC")
my_seq[4:12]

"""
Explanation: Note that using the Bio.SeqUtils.GC() function should automatically cope with mixed case sequences and the ambiguous nucleotide S which means G or C.
Also note that just like a normal Python string, the Seq object is in some ways ``read-only''.  If you need to edit your sequence, for example simulating a point mutation, look at the Section sec:mutable-seq below which talks about the MutableSeq object.
Slicing a sequence
A more complicated example, let's get a slice of the sequence:
End of explanation
"""


print(my_seq[0::3])
print(my_seq[1::3])
print(my_seq[2::3])

"""
Explanation: Note that ‘Seq‘ objects follow the usual indexing conventions for Python strings, with the first element
of the sequence numbered 0. When you do a slice the first item is included (i.e. 4 in this case) and the last
is excluded (12 in this case).
Also like a Python string, you can do slices with a start, stop and stride (the step size, which defaults to one).  For example, we can get the first, second and third codon positions of this DNA sequence:
End of explanation
"""


my_seq[::-1]

"""
Explanation: Another stride trick you might have seen with a Python string is the use of a -1 stride to reverse the string.  You can do this with a Seq object too:
End of explanation
"""


str(my_seq)

"""
Explanation: Turning Seq objects into strings
If you really do just need a plain string, for example to write to a file, or insert into a database, then this is very easy to get:
End of explanation
"""


print(my_seq)

"""
Explanation: Since calling str() on a Seq object returns the full sequence as a string,
you often don't actually have to do this conversion explicitly.
Python does this automatically in the print function:
End of explanation
"""


fasta_format_string = ">Name\n%s\n" % my_seq
print(fasta_format_string)

"""
Explanation: You can also use the Seq object directly with a %s placeholder when using the Python string formatting or interpolation operator (%):
End of explanation
"""


protein_seq = Seq("EVRNAK")
dna_seq = Seq("ACGT")
protein_seq + dna_seq


"""
Explanation: This line of code constructs a simple FASTA format record (without worrying about line wrapping). Section 4.6 describes a neat way to get a FASTA formatted string from a SeqRecord object, while the more
general topic of reading and writing FASTA format sequence files is covered in Chapter 5.
Concatenating or adding sequences
As of Biopython 1.78, you can add any two Seq objects together.
End of explanation
"""


from Bio.Seq import Seq
list_of_seqs = [Seq("ACGT"), Seq("AACC"), Seq("GGTT")]
concatenated = Seq("")
for s in list_of_seqs:
    concatenated += s
concatenated


"""
Explanation: Deliberately mixing DNA and protein like this is likely a mistake though...
You may often have many sequences to add together, which can be done with a for loop like this:
End of explanation
"""


from Bio.Seq import Seq
contigs = [Seq("ATG"), Seq("ATCCCG"), Seq("TTGCA")]
spacer = Seq("N"*10)
spacer.join(contigs)


"""
Explanation: Like Python strings, Biopython Seq also has a .join method:
End of explanation
"""


dna_seq = Seq("acgtACGT")
print(dna_seq)
print(dna_seq.upper())
print(dna_seq.lower())

"""
Explanation: Changing case
Python strings have very useful upper and lower methods for changing the case. For example,
End of explanation
"""


print("GTAC" in dna_seq)
print("GTAC" in dna_seq.upper())

"""
Explanation: These are useful for doing case insensitive matching:
End of explanation
"""


from Bio.Seq import Seq
my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC")
print(my_seq)
print(my_seq.complement())
print(my_seq.reverse_complement())


"""
Explanation: Nucleotide sequences and (reverse) complements
For nucleotide sequences, you can easily obtain the complement or reverse
complement of a Seq object using its built-in methods:
End of explanation
"""


my_seq[::-1]

"""
Explanation: As mentioned earlier, an easy way to just reverse a Seq object (or a
Python string) is slice it with -1 step:
End of explanation
"""


from Bio.Seq import Seq
protein_seq = Seq("EVRNAK")
protein_seq.complement()


"""
Explanation: If you do accidentally end up trying to do something weird like taking the (reverse)complement of a
protein sequence, the results are biologically meaningless:
End of explanation
"""


coding_dna = Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG")
print(coding_dna)
template_dna = coding_dna.reverse_complement()
print(template_dna)

"""
Explanation: Here the letter “E” is not a valid IUPAC ambiguity code for nucleotides, so was not complemented.
However, “V” means “A”, “C” or “G” and has complement “B“, and so on.
The example in Section 5.5.3 combines the Seq object’s reverse complement method with Bio.SeqIO for sequence input/output.
Transcription
Before talking about transcription, I want to try and clarify the strand issue.
Consider the following (made up) stretch of double stranded DNA which
encodes a short peptide:
|    | DNA coding strand (Crick strand, strand +1)    |    |
| -- | ---------------------------------------------- | -- |
| 5' | ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG        | 3' |
| 3' | TACCGGTAACATTACCCGGCGACTTTCCCACGGGCTATC        | 5' |
|    | DNA template strand (Watson strand, strand -1) |    |
|    |           Transcription                    |    |
| 5' | AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG}       | 3' |
|    | Single stranded messenger RNA                  |    |
The actual biological transcription process works from the template strand, doing a reverse complement (TCAG -> CUGA) to give the mRNA.  However, in Biopython and bioinformatics in general, we typically work directly with the coding strand because this means we can get the mRNA sequence just by switching T -> U.
Now let's actually get down to doing a transcription in Biopython.  First, let's create Seq objects for the coding and template DNA strands:
End of explanation
"""


messenger_rna = coding_dna.transcribe()
messenger_rna

"""
Explanation: These should match the figure above - remember by convention nucleotide sequences are normally read from the 5' to 3' direction, while in the figure the template strand is shown reversed.
Now let's transcribe the coding strand into the corresponding mRNA, using the Seq object's built in transcribe method:
End of explanation
"""


template_dna.reverse_complement().transcribe()

"""
Explanation: As you can see, all this does is switch T $\rightarrow$ U.
If you do want to do a true biological transcription starting with the template strand, then this becomes a two-step process:
End of explanation
"""


messenger_rna.back_transcribe()

"""
Explanation: The Seq object also includes a back-transcription method for going from the mRNA to the coding strand of the DNA.  Again, this is a simple U $\rightarrow$ T substitution:
End of explanation
"""


messenger_rna.translate()

"""
Explanation: Translation
Sticking with the same example discussed in the transcription section above,
now let's translate this mRNA into the corresponding protein sequence - again taking
advantage of one of the Seq object's biological methods:
End of explanation
"""


coding_dna.translate()

"""
Explanation: You can also translate directly from the coding strand DNA sequence:
End of explanation
"""


coding_dna.translate(table="Vertebrate Mitochondrial")

"""
Explanation: You should notice in the above protein sequences that in addition to the end stop character, there is an internal stop as well.  This was a deliberate choice of example, as it gives an excuse to talk about some optional arguments, including different translation tables (Genetic Codes).
The translation tables available in Biopython are based on those from the NCBI (see the next section of this tutorial). By default, translation will use the standard genetic code (NCBI table id 1). Suppose we are dealing with a mitochondrial sequence. We need to tell the translation function to use the relevant genetic code instead:
End of explanation
"""


coding_dna.translate(table=2)

"""
Explanation: You can also specify the table using the NCBI table number which is shorter, and often included in the feature annotation of GenBank files:
End of explanation
"""


coding_dna.translate()

coding_dna.translate(to_stop=True)

coding_dna.translate(table=2)

coding_dna.translate(table=2, to_stop=True)

"""
Explanation: Now, you may want to translate the nucleotides up to the first in frame stop codon,
and then stop (as happens in nature):
End of explanation
"""


coding_dna.translate(table=2, stop_symbol="@")

"""
Explanation: Notice that when you use the to_stop argument, the stop codon itself
is not translated - and the stop symbol is not included at the end of your protein
sequence.
You can even specify the stop symbol if you don't like the default asterisk:
End of explanation
"""


gene = Seq("GTGAAAAAGATGCAATCTATCGTACTCGCACTTTCCCTGGTTCTGGTCGCTCCCATGGCA" +
           "GCACAGGCTGCGGAAATTACGTTAGTCCCGTCAGTAAAATTACAGATAGGCGATCGTGAT" + 
           "AATCGTGGCTATTACTGGGATGGAGGTCACTGGCGCGACCACGGCTGGTGGAAACAACAT" + 
           "TATGAATGGCGAGGCAATCGCTGGCACCTACACGGACCGCCGCCACCGCCGCGCCACCAT" + 
           "AAGAAAGCTCCTCATGATCATCACGGCGGTCATGGTCCAGGCAAACATCACCGCTAA")
gene.translate(table="Bacterial")

gene.translate(table="Bacterial", to_stop=True)

"""
Explanation: Now, suppose you have a complete coding sequence CDS, which is to say a
nucleotide sequence (e.g. mRNA -- after any splicing) which is a whole number
of codons (i.e. the length is a multiple of three), commences with a start
codon, ends with a stop codon, and has no internal in-frame stop codons.
In general, given a complete CDS, the default translate method will do what
you want (perhaps with the to_stop option). However, what if your
sequence uses a non-standard start codon? This happens a lot in bacteria --
for example the gene yaaX in E. coli K12:
End of explanation
"""


gene.translate(table="Bacterial", cds=True)

"""
Explanation: In the bacterial genetic code GTG is a valid start codon,
and while it does normally encode Valine, if used as a start codon it
should be translated as methionine. This happens if you tell Biopython your
sequence is a complete CDS:
End of explanation
"""


from Bio.Data import CodonTable
standard_table = CodonTable.unambiguous_dna_by_name["Standard"]
mito_table = CodonTable.unambiguous_dna_by_name["Vertebrate Mitochondrial"]

"""
Explanation: In addition to telling Biopython to translate an alternative start codon as
methionine, using this option also makes sure your sequence really is a valid
CDS (you'll get an exception if not).
Translation Tables
In the previous sections we talked about the Seq object translation method (and mentioned the equivalent function in the Bio.Seq module – see Section 3.13). Internally these use codon table objects derived from the NCBI information at ftp://ftp.ncbi.nlm.nih.gov/entrez/misc/data/gc.prt, also shown on https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi in a much more readable layout.
As before, let’s just focus on two choices: the Standard translation table, and the translation table for Vertebrate Mitochondrial DNA.
End of explanation
"""


standard_table = CodonTable.unambiguous_dna_by_id[1]
mito_table = CodonTable.unambiguous_dna_by_id[2]

"""
Explanation: Alternatively, these tables are labeled with ID numbers 1 and 2, respectively:
End of explanation
"""


print(standard_table)

print(mito_table)

"""
Explanation: You can compare the actual tables visually by printing them:
End of explanation
"""


print(mito_table.stop_codons)
print(mito_table.start_codons)
print(mito_table.forward_table["ACG"])

"""
Explanation: You may find these following properties useful -- for example if you are trying
to do your own gene finding:
End of explanation
"""


from Bio.Seq import Seq
seq1 = Seq("ACGT")
seq2 = Seq("ACGT")



"""
Explanation: Comparing Seq objects
Sequence comparison is actually a very complicated topic, and there is no easy way to decide if two sequences
are equal. The basic problem is the meaning of the letters in a sequence are context dependent - the letter
“A” could be part of a DNA, RNA or protein sequence. Biopython can track the molecule type, so comparing
two Seq objects could mean considering this too.
Should a DNA fragment “ACG” and an RNA fragment “ACG” be equal? What about the peptide
“ACG“? Or the Python string “ACG“? In everyday use, your sequences will generally all be the same type
of (all DNA, all RNA, or all protein). Well, as of Biopython 1.65, sequence comparison only looks at the
sequence and compares like the Python string.
End of explanation
"""


print(seq1 == seq2)
print(seq1 == seq1)

"""
Explanation: So, what does Biopython do? Well, the equality test is the default for Python
objects -- it tests to see if they are the same object in memory. This is a
very strict test:
End of explanation
"""


my_seq = Seq("GCCATTGTAATGGGCCGCTGAAAGGGTGCCCGA")

"""
Explanation: As an extension to this, while you can use a Python dictionary with
Seq objects as keys, it is generally more useful to use the sequence a
string for the key.
MutableSeq Objects
Just like the normal Python string, the Seq object is ``read only'', or in Python terminology, immutable.  Apart from wanting the Seq object to act like a string, this is also a useful default since in many biological applications you want to ensure you are not changing your sequence data:
End of explanation
"""


my_seq[5] = "G"

"""
Explanation: Observe what happens if you try to edit the sequence:
End of explanation
"""


mutable_seq = my_seq.tomutable()
mutable_seq

"""
Explanation: However, you can convert it into a mutable sequence (a MutableSeq object) and do pretty much anything you want with it
End of explanation
"""


from Bio.Seq import MutableSeq
mutable_seq = MutableSeq("GCCATTGTAATGGGCCGCTGAAAGGGTGCCCGA")

"""
Explanation: Alternatively, you can create a MutableSeq object directly from a string:
End of explanation
"""


mutable_seq[5] = "C"
print(mutable_seq)
mutable_seq.remove("T")
print(mutable_seq)
mutable_seq.reverse()
print(mutable_seq)

"""
Explanation: Either way will give you a sequence object which can be changed:
End of explanation
"""


new_seq = mutable_seq.toseq()
new_seq

"""
Explanation: Do note that unlike the Seq object, the MutableSeq object's methods like reverse_complement() and reverse() act in-situ!
An important technical difference between mutable and immutable objects in Python means that you can't use a MutableSeq object as a dictionary key, but you can use a Python string or 
a Seq object in this way.
Once you have finished editing your a MutableSeq object, it's easy to get back to a read-only Seq object should you need to:
End of explanation
"""


from Bio.Seq import UnknownSeq
unk = UnknownSeq(20)
unk

print(unk)
print(len(unk))

"""
Explanation: You can also get a string from a MutableSeq object just like from a Seq object.
UnknownSeq Objects
The UnknownSeq object is a subclass of the basic Seq object
and its purpose is to represent a
sequence where we know the length, but not the actual letters making it up.
You could of course use a normal Seq object in this situation, but it wastes
rather a lot of memory to hold a string of a million N'' characters when you could
just store a single letterN'' and the desired length as an integer.
End of explanation
"""


unk_dna = UnknownSeq(20, character= "N")
unk_dna

print(unk_dna)

"""
Explanation: For DNA or RNA sequences, unknown nucleotides are commonly denoted by the letter “N”, while for
proteins “X” is commonly used for unknown amino acids. When creating an ‘UnknownSeq‘, you can specify
the character to be used instead of “?” to represent unknown letters. For example
End of explanation
"""


unk_dna

unk_dna.complement()

unk_dna.reverse_complement()

unk_dna.transcribe()

unk_protein = unk_dna.translate()
unk_protein

print(unk_protein)
print(len(unk_protein))

"""
Explanation: You can use all the usual Seq object methods too, note these give back
memory saving UnknownSeq objects where appropriate as you might expect:
End of explanation
"""


from Bio.Seq import reverse_complement, transcribe, back_transcribe, translate
my_string = "GCTGTTATGGGTCGTTGGAAGGGTGGTCGTGCTGCTGGTTAG"
print(reverse_complement(my_string))
print(transcribe(my_string))
print(back_transcribe(my_string))
print(translate(my_string))

"""
Explanation: You may be able to find a use for the UnknownSeq object in your own
code, but it is more likely that you will first come across them in a
SeqRecord object created by Bio.SeqIO.
Some sequence file formats don't always include the actual sequence, for
example GenBank and EMBL files may include a list of features but for the
sequence just present the contig information.  Alternatively, the QUAL files
used in sequencing work hold quality scores but they never contain a
sequence -- instead there is a partner FASTA file which does have the
sequence.
Working with strings directly
To close this chapter, for those you who really don't want to use the sequence
objects (or who prefer a functional programming style to an object orientated one),
there are module level functions in Bio.Seq will accept plain Python strings,
Seq objects (including UnknownSeq objects) or MutableSeq objects:
End of explanation
"""


