This tutorial was generated from a Jupyter notebook. You can download the notebook here.
There are two fundamental approaches to coding, a tactical and a strategic: The tactical focuses on getting things done; it's quick and dirty, whereas the strategic is more concerned about code maintance and reusability. Both concepts have their benefits and often there is no clear-cut boundary between them. No matter which line of attack the following coding philosophies will avert some pain:
While functions help to avoid duplicated code object-oriented programming (OOP) goes further:
Classes form the core of Python's OOP. They are devices that create new objects. This is called instantiation. I.e., an object is an instance of a class. These objects have three characteristics which makes the use of classes particularly attractive:
+ operator has a different meaning for floats vs. strings.This is all a bit theoretical, so let's look at an example.
In this example, we will make a class for biological sequences. The goal here is just to lift the covers off how classes are built, so you have a basic understanding what you are doing when you use the dot notation on objects.
class Biosequence(object):
"""Biological sequence class with methods."""
# The special method __init__() is run when the class is instantiated
def __init__(self, seq='', material='dna'):
"""
Instantiate Biosequence object with sequence and material.
Material is either 'dna', 'rna', or 'protein'.
"""
self.seq = seq
self.material = material.lower()
# Define a method
def seq_len(self):
"""The length of the sequence."""
return len(self.seq)
Methods are functions contained within a class. In contrast to functions they need to include the term self within the parentheses. Because when Python calls a method the current object becomes the first argument of the method. Using the word self is a convention; theoretically any other word could be used.
The method we have written for the Biosequence class is trivial, but it demonstrates the concept.
# Instantial Biosequence
s = Biosequence(seq='ATGAAGGGTCC')
# Call the methods
print('The sequence has', s.seq_len(), 'bases.')
In addition to methods, an object also has attributes, which are values associated with it. This class has two attributes, seq and material, and has two methods, seq_len() and the special method __init__().
Let's create a subclass Nucleotides. To make it inherit all of the defined methods and attributes of Biosequence (its parent class), we refer to Biosequence in the class definition. We will write a single method that uses the Marmur rule of thumb for computing the melting temperature of a stretch of double-stranded DNA,
where $N_i$ denotes the number of times nucleotide $i$ appears in the sequence. Remember, the class inherits all of the methods from its parent, so we only need to write methods specific to nucleic acids. This is part of the DRY principle.
class Nucleotides(Biosequence):
"""Nucleotide sequences and related methods."""
# Add a nucleic acid-specific method
def T_marmur(self):
"""Melting temperature by Marmur rule of thumb."""
# Convert to upper case
seq_up = self.seq.upper()
return 2 * (seq_up.count('A') + seq_up.count('T')) \
+ 4 * (seq_up.count('G') + seq_up.count('C'))
So, we now have a new class for nucleotides that inherited all of the methods and attributes of its parent, Biosequence. It also has the new melting temperature method. Let's put it to use.
# Instantiate object
s = Nucleotides('AAAGGTTTTTTTTTTTC', material='dna')
# Compute result
s.T_marmur()
Let's try it with another sequence.
s = Nucleotides('aaattg')
s.T_marmur()
Let's add the Wallace rule of thumb for melting temperature in degrees Celsius.
\begin{align} T_m = 64.9 + 41\,\frac{N_G + N_C - 16.4}{N_A + N_T + N_G + N_C}. \end{align}This, of course, works only for sequences that have some reasonable GC content. Clearly, a poly-AT sequence should not melt below freezing. Again, we want to return the results so we can actually use this class.
class Nucleotides(Biosequence):
"""Nucleotide sequences and related methods."""
def T_marmur(self):
"""Melting temperature by Marmur rule of thumb."""
# Convert to upper case
seq_up = self.seq.upper()
return 2 * (seq_up.count('A') + seq_up.count('T')) \
+ 4 * (seq_up.count('G') + seq_up.count('C'))
def T_wallace(self):
"""Melting temperature by Wallace rule of thumb."""
# Convert to upper case
seq_up = self.seq.upper()
return 64.9 + 41 * (seq_up.count('G') + seq_up.count('C') - 16.4) \
/ len(self.seq)
Now let's compare the two melting temperatures.
# Instatiate class
s = Nucleotides('GGGTTTCCCTACA')
# Compute melting temperatures
print('T_m via Marmur rule:', s.T_marmur(), '°C')
print('T_m via Wallace rule:', s.T_wallace(), '°C')
Now that we can calculate melting temperatures, we can introduce a class that deals with primer design. It inherits the Nucleotides class, and then has methods to choose a primer.
class Primer(Nucleotides):
"""Primer design."""
def marmur_primer(self, Tm):
"""
Primer design based on Marmur rule of thumb for melting
temperature.
`Tm`: the desired melting temperature in deg. C of the
primer-template duplex.
"""
# Start with a short primer of 10 nucleotides
n = 10
# Start at 5' end with forward primer
forward_p = Nucleotides(self.seq[:n])
# Keep adding based to the primer until we hit the desired Tm
while (forward_p.T_marmur() < Tm and n < len(self.seq)):
n += 1
forward_p.seq = self.seq[:n]
return forward_p.seq
Let's try to design a primer for a sequence. We'll take our melting temperature to be 42°C, or course.
# Instantiate object with template sequence
s = Primer('AACCCCCCAAATTTTTTTTTTGAAAAAAAAAACATATTCTTCTCTCGGGGGGGGG')
s.marmur_primer(42)
We would also like to compute a reverse primer. To do this, we need to be able to compute the reverse complement of a sequence. We have already done that in the last exercise, so we add that functionality to the Nucleotides class. Remember that since the method is aware of self, and the object has a seq attribute, we do not need to pass the seq argument into the function that computes the reverse complement.
class Nucleotides(Biosequence):
"""Nucleotide sequences and related methods."""
def T_marmur(self):
"""Melting temperature by Marmur rule of thumb."""
# Convert to upper case
seq_up = self.seq.upper()
return 2 * (seq_up.count('A') + seq_up.count('T')) + \
4 * (seq_up.count('G') + seq_up.count('C'))
def T_wallace(self):
"""Melting temperature by Wallace rule of thumb."""
# Convert to upper case
seq_up = self.seq.upper()
return 64.9 + 41 * (seq_up.count('G') + seq_up.count('C') - 16.4) / \
len(self.seq)
def reverse_complement(self):
"""Compute reverse complement of the sequence."""
# Initialize reverse complement as reverse of sequence
rev_comp = self.seq.upper()[::-1]
# Replace bases with complement
rev_comp = rev_comp.replace('A', 't')
rev_comp = rev_comp.replace('T', 'a')
rev_comp = rev_comp.replace('C', 'g')
rev_comp = rev_comp.replace('G', 'c')
return rev_comp.upper()
We now need to update the Primer class to compute the reverse primer. We already wrote most of the code to do this, so we just need to tweak the class a bit.
class Primer(Nucleotides):
"""Primer design."""
def forward_and_reverse_marmur_primers(self, Tm):
"""
Forward and reverse primer design using Marmur rule
of thum for melting temperature.
`Tm`: the desired melting temperature in deg. C of the
primer-template duplex.
"""
# Compute forward primer
forward_p = self.marmur_primer(Tm)
# Compute reverse complement of strand
rev_comp = Primer(self.reverse_complement())
# Compute reverse primer
reverse_p = rev_comp.marmur_primer(Tm)
return forward_p, reverse_p
def marmur_primer(self, Tm):
"""
Primer design based on Marmur rule of thumb for melting
temperature.
`Tm`: the desired melting temperature in deg. C of the
primer-template duplex.
"""
# Start with a short primer of 10 nucleotides
n = 10
# Start at 5' end with forward primer
forward_p = Nucleotides(self.seq[:n])
# Keep adding based to the primer until we hit the desired Tm
while (forward_p.T_marmur() < Tm and n < len(self.seq)):
n += 1
forward_p.seq = self.seq[:n]
return forward_p.seq
Let's go through each line of the method forward_and_reverse_marmur_primers(). First, we just compute the forward primer using the method we already wrote. Note that because we are in a method within the class definition, we need to use self, i.e., self.marmum_primer(Tm). Next, we instantiate a new Primer instance that has the reverse complement of the strand of interest as its sequence. We then compute the forward primer of that (which is the reverse complement of our input sequence). We then return both primers.
Let's give it a whirl!
# Instantiate sequence
s = Primer('AACCCCCCAAATTTTTTTTTTGAAAAAAAAAACATATTCTTCTCTCGGGGGGGGG')
# Compute primers
s.forward_and_reverse_marmur_primers(Tm=42)
This works, that's great! Notice that we tried to reuse code we already wrote as much as possible. This is part of the DRY principle. It also is a great help in debugging!
In keeping with the DRY principle, we can recycle the Biosequence class to handle protein sequences. We'll make a new class with that computes the net charge of residues in a protein (assuming histidine is not charged).
class Protein(Biosequence):
"""Protein sequences."""
def netcharge(self):
"""Compute the net charge of a protein"""
# Conver to upper case
sequence = self.seq.upper()
return sequence.count('K') + sequence.count('R') - \
sequence.count('E') - sequence.count('D')
Again, we inherit all of the methods from the Biosequence class. We can instantiate and compute net charges!
# Specify the sequence
seq = 'MKKILLSVLTAFVAVVLAACGGNSDSKTLNSLDKIKQNGVVRIGVFGDKPPFGYVDEKGNNQGYDIALAK\
RIAKELFGDENKVQFVLVEAANRVEFLKSNKVDIILANFTQTPQRAEQVDFCSPYMKVALGVAVPKDSNI\
TSVEDLKDKTLLLNKGTTADAYFTQNYPNIKTLKYDQNTETFAALMDKRGDALSHDNTLLFAWVKDHPDF\
KMGIKELGNKDVIAPAVKKGDKELKEFIDNLIIKLGQEQFFHKAYDETLKAHFGDDVKADDVVIEGGKI'
# Instantiate
p = Protein(seq)
# Compute net charge
p.netcharge()