import sandbox_widget48 Biopython
Somebody already wrote the classes you need for sequences, sequence files, alignments, BLAST searches and protein structures. This chapter shows you how to find your way around them.
Classes someone else wrote
In Chapter 38 you wrote a Gene class that kept a name and a DNA sequence together, and you gave it __str__ and __len__ so that print and len knew what to do with it. Biopython is what you get when people who need that kind of class every day keep at it for more than twenty-five years. It is a package (Chapter 43), which means a folder of ordinary .py files full of class statements just like yours.
The class you will use most is Seq, which you import with from Bio.Seq import Seq. Calling Seq("ATGGCC") runs its __init__ and hands you a Seq object, just as Point(4, 3) handed you a Point. Peek behind the curtain and you find the same special methods you wrote yourself: __str__ and __repr__ decide how it prints, __len__ makes len work, __contains__ makes in work, __add__ makes + work, and __getitem__ makes indexing and slicing work (which is exactly what your Gene could not do when you tried gene[0]). That is why a Seq behaves almost like a string. On top of that, it has methods no string has, like reverse_complement and translate. “Hey sequence, translate yourself!”
A SeqRecord is Biopython’s grown-up version of your Gene. It keeps a Seq in its seq attribute, along with an id, a name and a description, and you reach all of them with the same dot you used for p.x. You will rarely make one yourself, because reading a FASTA or GenBank file hands them to you. Whenever a Biopython object puzzles you, ask it what you would ask any object: what type are you (type), what data do you carry (your attributes), and what can you do (your methods)? Running help on it answers the last two.
Biopython is big, and this chapter only shows you the parts you are most likely to need. Everything here is described in much more detail in the Biopython Tutorial, and that tutorial, not the assistant, is what you check against when you are unsure what a Biopython function does.
Sequences
Let us make a Seq and poke at it. Type the cell below yourself (no copy-paste) and run it.
from Bio.Seq import Seq
dna = Seq("ATGGCCATTGTAATGGGCCGCTGA")
dnaThe last line of the cell displayed Seq('ATGGCCATTGTAATGGGCCGCTGA'). That is the __repr__ method at work, and it tells you both the type and the content. print calls __str__ instead, which gives you just the letters:
print(dna)Exercise 48-1
Everything you know about strings works on a Seq. Decide what each line prints before you run the cell. Pay special attention to the last line: when you slice a Seq, what type of object do you get back?
%%sandbox
from Bio.Seq import Seq
dna = Seq("ATGGCCATTGTAATGGGCCGCTGA")
print(len(dna))
print(dna[:3])
print(dna[-3:])
print("GCC" in dna)
print(dna.count("G"))
print(dna[:3] == "ATG")
print(type(dna[:3]))Exercise 48-2
A Seq is like a string in one more way. Remember from Section 18.0.0.15 that you cannot change a string in place? Predict what happens when you try with a Seq, then run it and read the last line of the error message.
%%sandbox
from Bio.Seq import Seq
dna = Seq("ATGGCC")
dna[0] = "C"Just like with strings, the way around it is to make a new sequence: Seq("C") + dna[1:]. The + works because Seq has an __add__ method, the same special method you gave your Point class.
Biology a string does not know
So far a Seq has behaved like a string. Here is where it earns its keep. Predict what each line prints before you run it. You know what a complement, a reverse complement, and a transcript are, so you can work out the answers on paper first.
print(dna.complement())
print(dna.reverse_complement())
print(dna.transcribe())
print(dna.translate())“Hey sequence, reverse complement yourself!” Each method returns a new Seq and leaves dna itself unchanged. The translation ends in *, which is how Biopython writes a stop codon, just as you did in the translation project (Chapter 33).
The translate method takes keyword arguments that change how it works. to_stop=True stops at the first stop codon and leaves the * out. table picks the genetic code by its NCBI number. Table 1 is the standard code you used in the translation project, and table 2 is the code used in vertebrate mitochondria, where TGA codes for tryptophan (W) and only TAA, TAG, AGA and AGG mean stop.
Exercise 48-3
Decide what each of the three lines prints, then run the cell to check. The sequence has two stop codons in the standard code. How many does it have in the mitochondrial code?
%%sandbox
from Bio.Seq import Seq
orf = Seq("ATGTGATAA")
print(orf.translate())
print(orf.translate(to_stop=True))
print(orf.translate(table=2))Exercise 48-4
In the translation project you wrote your own translate_orf function, and the project’s tests are the contract it had to satisfy. Here are the three test inputs from that test file run through Biopython instead. Predict what Biopython prints for each, and then compare with what your own translate_orf returns for the same three inputs (run it in your project folder).
The two agree on two of the inputs and disagree on one. Which one, and who is right? Look at the test file again before you decide. Each function does what its own contract says, and the contracts are different.
%%sandbox
from Bio.Seq import Seq
for orf in ["ATGCCCATGTGA", "ATGATNATGTGA", "ATGTGA"]:
print(orf, Seq(orf).translate())GC content
Biopython keeps a box of small sequence utilities in the module Bio.SeqUtils. One of them computes GC content, which you also wrote yourself as a method on your Gene class. This one is a function rather than a method, so you hand it the sequence: “Hey function, compute the GC fraction of this sequence!”
from Bio.SeqUtils import gc_fraction
gc_fraction(dna)Exercise 48-5
In Chapter 38 you tested your own gc_content method on 'GGCC', where the answer must be 1.0, and on 'ATAT', where it must be 0.0. Run the same two tests on gc_fraction. Then predict what gc_fraction returns for the lowercase sequence 'ggcc', and check. Would your own method have passed that third test?
%%sandbox
from Bio.SeqUtils import gc_fraction
print(gc_fraction("GGCC"))Records
A sequence on its own is rarely enough. You want to know what it is called and where it came from. A SeqRecord packages a Seq together with that information, the same way your Gene packaged a name with a sequence:
from Bio.SeqRecord import SeqRecord
record = SeqRecord(dna, id="orf1", description="a small made-up ORF")
print(record)Printing the record called its __str__ method, which lists what it holds. Each item is an attribute that you can reach with the dot:
print(record.id)
print(record.description)
print(record.seq)
print(len(record))len(record) works because SeqRecord has a __len__ method that hands the question on to the sequence it holds. That is exactly what the __len__ of your Gene class did with its string.
Exercise 48-6
Objects can hold other objects, and methods can return new objects. Predict the type of each of the three expressions below before you run the cell. The last one is the surprising one: what do you think Biopython does to the id and the seq when you slice a whole record?
%%sandbox
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
record = SeqRecord(Seq("ATGGCCATTGTAATGGGCCGCTGA"), id="orf1", description="a small made-up ORF")
print(type(record.seq))
print(type(record.seq.translate()))
print(type(record[3:9]))
print(record[3:9].id, record[3:9].seq)Reading and writing sequence files
In the projects you read sequences from files by opening them and looping over the lines. Biopython does that work for you in the module SeqIO, short for sequence input/output. It reads and writes many file formats, and the two you will meet most are FASTA and GenBank.
Let us first write a few records to a FASTA file, so we have something to read. SeqIO.write takes the records, a file name, and the name of the format:
from Bio import SeqIO
my_records = [
SeqRecord(Seq("ATGGCCATTGTAATGGGCCGCTGA"), id="orf1", description="first"),
SeqRecord(Seq("ATGAAACCCGGGTTTTAA"), id="orf2", description="second"),
SeqRecord(Seq("ATGTGATAA"), id="orf3", description="third"),
]
SeqIO.write(my_records, "my_orfs.fasta", "fasta")SeqIO.write returned 3, the number of records it wrote. Open my_orfs.fasta in VS Code and have a look. It is an ordinary text file, and you could have written it with print(..., file=f) as in Section 30.0.0.4. Now read it back with SeqIO.parse:
for rec in SeqIO.parse("my_orfs.fasta", "fasta"):
print(rec.id, len(rec), rec.seq.translate(to_stop=True))SeqIO.parse hands you one record at a time as the for-loop asks for it, just like a file object serves one line at a time (Section 30.0.0.5). That is why it can read a file with millions of sequences without filling up your memory. If you do want all the records at once, turn them into a list with list(SeqIO.parse("my_orfs.fasta", "fasta")).
Exercise 48-7
SeqIO has a sibling of parse called read, which you use when the file holds exactly one record. It then returns that record directly, so you do not need a loop. What do you think happens if you use read on a file with three records? Decide before you run it.
%%sandbox
from Bio import SeqIO
rec = SeqIO.read("my_orfs.fasta", "fasta")
print(rec.id)Fetching records from NCBI
Where did that GenBank file come from? From the NCBI databases, which you can search in the browser at ncbi.nlm.nih.gov. The module Bio.Entrez lets your program do the same searches and downloads. NCBI asks for three things in return. You must tell them your email address, so they can contact you if your program misbehaves. You may send at most three requests per second (Biopython makes sure of that for you). And if you need to send more than a hundred requests in a row, you should do it on a weekend or outside US office hours.
The cell below downloads the record we have been using. Put your own email address in it before you run it. This cell, and the other cells in the chapter that talk to NCBI, need an internet connection.
from Bio import Entrez
Entrez.email = "your.name@post.au.dk" # use your own address
with Entrez.efetch(db="nucleotide", id="NM_148900.4", rettype="gb", retmode="text") as handle:
fetched = SeqIO.read(handle, "genbank")
print(fetched.description)
SeqIO.write(fetched, "data/NM_148900.4.gb", "genbank")Entrez.efetch gives you a handle, an object that behaves like an open file, except that the content comes from NCBI rather than from your disk. That is why you can hand it to SeqIO.read just as you would a file name, and why you use it in a with statement that closes it for you afterward.
If you do not know the accession number, you search for it first with Entrez.esearch. The search term uses the same syntax as the search box on the NCBI website, and Entrez.read turns NCBI’s answer into a dictionary:
with Entrez.esearch(db="nucleotide", term="FOXP2[Gene] AND Homo sapiens[Organism] AND srcdb_refseq[prop]", retmax=10) as handle:
result = Entrez.read(handle)
print(result["Count"])
print(result["IdList"])Exercise 48-8
Change the search so it finds the FOXP2 sequences of the chimpanzee (Pan troglodytes) instead. Then fetch one of them with Entrez.efetch, read it with SeqIO.read, and print its description and length. Predict first: do you expect it to be about as long as the human one?
Exercise 48-9
The search term above contains srcdb_refseq[prop]. Ask the assistant what that part of the term means and what would change if you left it out. Then settle it by running the search with and without it and comparing the two values of result["Count"]. Write in your logbook whether the numbers matched the explanation.
GenBank files
A FASTA file holds little more than names and sequences. A GenBank file holds a sequence together with everything the database knows about it. The data folder of your course folder has the GenBank file for the messenger RNA of the human gene FOXP2. FOXP2 is a transcription factor, and mutations in its gene cause a rare inherited disorder of speech and language (Lai et al., 2001). Open NM_148900.4.gb in VS Code and scroll through it before you go on. Then read it:
record = SeqIO.read("data/NM_148900.4.gb", "genbank")
print(record.id)
print(record.name)
print(record.description)
print(len(record))The rest of the header of the file ends up in record.annotations, which is an ordinary dictionary:
print(record.annotations["organism"])
print(record.annotations["molecule_type"])Exercise 48-10
record.annotations is a dictionary, and you know how to see all the keys of a dictionary. Print them all. Then print the value for the key that tells you when the record was last updated.
SeqIO.convert reads a file in one format and writes it in another, so converting the GenBank file to FASTA is one line:
SeqIO.convert("data/NM_148900.4.gb", "genbank", "NM_148900.4.fasta", "fasta")Exercise 48-11
Read NM_148900.4.fasta back in with SeqIO.read and print its id, name and description. Predict first: which of the three are the same as for the GenBank record above, and which changed? Open the FASTA file in VS Code to see why. The FASTA format has only a single header line to keep that information in.
Features
Scroll down to the FEATURES section of the GenBank file. It lists regions of the sequence and says what each one is: the gene, its exons, the coding sequence (CDS), and more. Biopython reads each one into a SeqFeature object and puts them all in the list record.features:
for feature in record.features:
print(feature.type, feature.location)Exercise 48-12
Find the line in the GenBank file where the CDS starts. The file says 557..2755, but Biopython printed [556:2755](+). Is one of them wrong? Before you read on, work out how the two can describe the same bases. Hint: how does Python count positions, and what does a slice [556:2755] include?
Both are right. GenBank counts from 1 and includes both ends. Biopython translates that into Python’s way of counting, so that you can use the numbers directly as a slice. The (+) says the feature is on the forward strand.
Let us find the coding sequence. We loop over the features and keep the one whose type is CDS:
for feature in record.features:
if feature.type == "CDS":
cds = feature
print(cds.qualifiers["gene"])
print(cds.qualifiers["product"])The qualifiers attribute is a dictionary with the extra information listed under each feature in the file. Notice that each value is a list of strings, even when there is only one string in it. The extract method cuts the feature’s part out of the record’s sequence. For a feature on the reverse strand, it also reverse complements it, which is why you should use extract rather than slicing yourself:
cds_seq = cds.extract(record.seq)
protein = cds_seq.translate(to_stop=True)
print(protein)Exercise 48-13
The GenBank file carries the database’s own translation of the coding sequence, under cds.qualifiers["translation"]. That means you can check your protein against it. Decide what each of the two lines below prints before you run it. If you predicted the same answer for both, look again at what you just learned about the values in qualifiers.
print(protein == cds.qualifiers["translation"])
print(protein == cds.qualifiers["translation"][0])Exercise 48-14
Count how many of the features are exons, and add up the lengths of all the exons (len works on a feature). Before you run your code, predict how the sum compares to len(record). Then explain the answer, remembering that this is the sequence of a messenger RNA.
Pairwise alignment
To compare two sequences, you align them. In Biopython that is the job of a PairwiseAligner object. You make one, and then you ask it to align two sequences:
from Bio import Align
aligner = Align.PairwiseAligner()
alignments = aligner.align("GAACT", "GAT")
print(alignments.score)
print(len(alignments))
for alignment in alignments:
print(alignment)There are two different alignments with the same, best score, and the aligner gives you both. The numbers in the margins are positions in each sequence, counted Python’s way from 0.
The score comes from the aligner’s scoring settings, which are attributes on the aligner object. By default a match scores 1, a mismatch scores 0, and every gap position scores -1:
print(aligner.match_score)
print(aligner.mismatch_score)
print(aligner.open_gap_score)
print(aligner.extend_gap_score)Exercise 48-15
Use the four scores above to compute the score of each of the two alignments by hand, position by position, the way you would do substitution and reduction steps: replace each column by its score, then add them up. Do you get the score the aligner printed?
Change the attributes and you change the alignment. Setting mode to "local" asks for the best-scoring part of the two sequences rather than the whole of both. For proteins, you usually score with a substitution matrix like BLOSUM62, which gives similar amino acids a positive score even if they are not identical:
from Bio.Align import substitution_matrices
protein_aligner = Align.PairwiseAligner()
protein_aligner.mode = "local"
protein_aligner.substitution_matrix = substitution_matrices.load("BLOSUM62")
protein_aligner.open_gap_score = -10
protein_aligner.extend_gap_score = -0.5
best = protein_aligner.align("HEAGAWGHEE", "PAWHEAE")[0]
print(best.score)
print(best)Exercise 48-16
Predict what the score of aligning "GAACT" with "GAT" becomes if you set aligner.mode = "local", keeping the default scores. Write down the local alignment you expect to win before you run your code.
Exercise 48-17
Look at the line in the middle of a printed alignment, the one made of |, . and -. Ask the assistant what each of the three symbols means. Then check its answer against one of the alignments you printed yourself, column by column. In the protein alignment, is every . really a mismatch between two different amino acids?
BLAST
An alignment compares two sequences you already have. BLAST finds the sequences in a database that are most similar to yours, and then aligns your sequence to each of them. The function Blast.qblast sends your sequence to NCBI’s BLAST servers. You give it the BLAST program to use (blastp searches a protein database with a protein), the database (swissprot is a carefully curated protein database), and your sequence. Here we search with the FOXP2 protein you translated above.
A BLAST search at NCBI takes at least a minute, and when many people are using the servers, your search can wait in their queue for a lot longer. NCBI also asks that you do not send more than one search every ten seconds. So save the result to a file and work with the file from then on. The cell below shows how. You do not need to run it, because we have run it for you and put the result in data/foxp2_blastp.xml.
proteinfrom Bio import Blast
Blast.email = "your.name@post.au.dk" # use your own address
stream = Blast.qblast("blastp", "swissprot", protein, hitlist_size=20)
with open("data/foxp2_blastp.xml", "wb") as f:
f.write(stream.read())
stream.close()The "wb" opens the file for writing bytes rather than text, because that is what the stream from NCBI delivers. The result is an XML file. Open it in VS Code if you are curious, but you will be happy that you do not have to read it yourself. Blast.read does that:
from Bio import Blast
blast_record = Blast.read("data/foxp2_blastp.xml")
print(len(blast_record))The result is built from objects inside objects, just like the records and features above. The blast_record holds one Hit for each database sequence that BLAST found, and you loop over it just like a list. Each Hit holds one or more HSPs (high-scoring segment pairs), which are the alignments between your sequence and that database sequence. Each HSP has its statistics in a dictionary called annotations:
for hit in blast_record:
hsp = hit[0]
print(hit.target.id, hsp.annotations["evalue"], hit.target.description[:60])Exercise 48-18
Pick out the first hit and its first HSP, and print the HSP itself. Before you run it, decide what printing an HSP shows you. You have seen printed alignments above, so what would you expect it to look like?
Exercise 48-19
hsp.annotations["identity"] is the number of identical positions in the alignment, and hsp.length is the length of the alignment. Write a loop that prints the percent identity of the first HSP of every hit, with the description of the hit. Which hits are FOXP2 from other species, and which are other genes of the same family (FOXP1, FOXP3, FOXP4)?
The first hit is human FOXP2 itself, yet it is not 100% identical to our protein. Predict why before you print its HSP and look. Hint: compare len(protein) with the length of the database sequence, which the printed HSP shows you, and remember that our record is transcript variant 4.
Exercise 48-20
Ask the assistant how to print the E-value of every hit in a BLAST result with Biopython, and give it nothing else. Many answers still use an older interface, Bio.Blast.NCBIXML, where the E-value lives in an attribute called hsp.expect. Run the code it gives you on data/foxp2_blastp.xml. If it fails, read the error and say which attribute or function does not exist in the Biopython you have installed. Settle it against the BLAST chapter of the Biopython Tutorial, and write in your logbook what the assistant claimed and what actually ran.
Protein structures
The Protein Data Bank (PDB) holds the measured three-dimensional structures of proteins, each identified by a four-character code. Structure 1HMP is the human enzyme hypoxanthine-guanine phosphoribosyltransferase (HGPRT), which recycles purine bases. Children born without a working copy of this enzyme have Lesch-Nyhan syndrome. The class PDBList downloads structures from the PDB:
from Bio.PDB import PDBList
pdblist = PDBList()
pdblist.retrieve_pdb_file("1HMP", pdir=".", file_format="pdb")That saves the file as pdb1hmp.ent. We have put a copy in your data folder as data/1hmp.pdb, so you do not need the download. A PDB file is a text file with one line per atom. Have a look at the first lines:
with open("data/1hmp.pdb") as f:
for i in range(10):
print(f.readline(), end="")The PDBParser class reads such a file into a Structure object (QUIET=True keeps it from complaining about small irregularities in the file):
from Bio.PDB import PDBParser
parser = PDBParser(QUIET=True)
structure = parser.get_structure("1HMP", "data/1hmp.pdb")
print(structure.header["name"])A structure is objects inside objects inside objects. A Structure holds one or more models, a Model holds chains (one per protein molecule in the crystal), a Chain holds residues, and a Residue holds atoms. Each level is a container that you index with square brackets and loop over with a for-loop, which tells you that each class has a __getitem__ and an __iter__ method. You index a model by number, a chain by its letter, a residue by its number in the protein sequence, and an atom by its name:
model = structure[0]
for chain in model:
print(chain.id, len(chain))
residue = model["A"][4]
print(residue.resname)
for atom in residue:
print(atom.get_id(), atom.coord)Each atom carries its coordinates in space. Now for a last peek behind the curtain. Atom objects have an __sub__ method, which is the one Python calls for the - operator. It does not subtract anything. It returns the distance between the two atoms in ångström, just like the distance method you wrote for your Point class:
ca4 = model["A"][4]["CA"]
ca5 = model["A"][5]["CA"]
print(ca4 - ca5)Exercise 48-21
Predict what ca4 - ca4 and ca5 - ca4 give, then run them. If you remember the exercise in Chapter 38 where you asked whether p1.distance(p2) and p2.distance(p1) give the same answer, you already know the answer to the second one.
Exercise 48-22
CA is the alpha carbon of each amino acid, and neighboring alpha carbons in a chain are held at about the same distance by the peptide bond. Loop over the residues of chain A and print the distance from each alpha carbon to the next one. Some residues in the chain are not amino acids (water molecules and the bound GMP), and they have no alpha carbon. Use in to skip those: "CA" in residue works because Residue has a __contains__ method. Predict what most of the distances will be before you run your loop. Then find the smallest and the largest distance. How much do they vary?
The class PPBuilder finds the stretches of the chain that are connected polypeptides. Each one knows its own amino acid sequence, and get_sequence gives it to you as a Seq object, so everything you learned at the start of the chapter works on it too:
from Bio.PDB import PPBuilder
polypeptides = PPBuilder().build_peptides(model["A"])
peptide = polypeptides[0]
print(len(polypeptides))
print(peptide.get_sequence())A Ramachandran plot
The shape of a protein backbone is set by two rotation angles at each amino acid, called phi (\(\phi\)) and psi (\(\psi\)). A polypeptide computes them for you with get_phi_psi_list, which gives one pair per residue, in radians. The first residue has no phi and the last has no psi, so those come back as None. Plotting psi against phi for every residue gives a Ramachandran plot, where residues in alpha-helices and beta-sheets each gather in their own region. The math module from ?sec-standard-library converts radians to degrees:
import math
import matplotlib.pyplot as plt
import seaborn as sns
phi_list = []
psi_list = []
for phi, psi in peptide.get_phi_psi_list():
if phi is not None and psi is not None:
phi_list.append(math.degrees(phi))
psi_list.append(math.degrees(psi))
sns.kdeplot(x=phi_list, y=psi_list, cmap="Reds", fill=True)
plt.scatter(phi_list, psi_list, marker="x", color="black")
plt.xlim(-180, 180)
plt.ylim(-180, 180)
plt.title("Ramachandran plot of 1HMP chain A")
plt.xlabel(r"$\phi$")
plt.ylabel(r"$\psi$")
plt.show()Exercise 48-23
Ask the assistant in which region of a Ramachandran plot residues in alpha-helices fall, and in which region residues in beta-sheets fall, as ranges of phi and psi in degrees. Do not take its word for it. Settle it against the Ramachandran plot in your biochemistry textbook, and then find the two regions in your own plot. Write in your logbook whether the assistant’s ranges matched the textbook.
General exercises
Exercise 48-24
Write a function first_protein(filename) that reads a GenBank file with a single record, finds its first CDS feature, and returns the translated protein as a string, stopping at the stop codon. Test it on NM_148900.4.gb by comparing its return value to the translation stored in the file itself.
Exercise 48-25
Write down on paper what a function translate_fasta(filename) should return for my_orfs.fasta: a dictionary that maps the id of every record to its protein, as a string, stopping at the first stop codon. Then ask the assistant to draft the function using SeqIO, and give it your description of what it should return, but not your answer. Read every line of the draft before you run it, and settle it by comparing what it returns to the dictionary you wrote on paper.
Exercise 48-26
Ask the assistant for a different way to find the CDS feature in record.features than the for-loop with an if-statement used in this chapter. Run both on the FOXP2 record and check that they find the same feature. Which of the two can you read more easily, and what does each one do if the record has two CDS features, or none at all? Try both cases to find out, rather than taking either your own or the assistant’s word for it.