Bioinformatics algorithms in Python
≈ 5 min
Write Python that works with DNA and protein sequences – then use it to see how the algorithms of bioinformatics work, from hash tables to alignments, trees and networks.
strings → lists → loops → functions → hashing → graphs → patterns → alignment → BLAST → multiple alignment → trees → networks
What you will do
Part 1 (chapters 1–4) is Python with DNA and protein sequences: strings, lists, decisions, loops and functions, used for real jobs – GC content, reverse complements, the charge of a protein, melting temperatures. Part 2 (chapters 5–12) uses that Python to look inside the algorithms of bioinformatics: hash tables, graph searches, regular expressions, Smith–Waterman alignment, BLAST, multiple alignments, phylogenetic trees and protein networks.
Each chapter takes 10–20 minutes; the whole practical about three hours. Boxes marked Optional are extensions for when you have time.
How this page works
On the right What it is Notebook Real Python running inside this web page – nothing to install. Each chapter has its own page of cells, already set up with the code you need; your cells are saved in this browser. Files The data files of the practical (sequences, alignments), which Python can open by name. Cheat sheet The Python of this practical on one page.
Exercises are numbered like 1.2 . The cell for each one is ready in the notebook (Go to the cell shows it). When your code is right, the exercise ticks itself as soon as you run the cell. Check tells you what is still wrong; after your first Check, Show a solution puts a working answer under your cell – try first, it is how the Python sticks.
Python counts from 0: the first letter of a sequence is at position 0. Errors are normal – every one comes with a 💡 hint in plain English.
Activity · Start Python
The first time, Python takes 10–30 seconds to download into the page. Wait until the notebook (on the right) says Python … ready .
Run the first cell of the notebook: click ▶ next to it, or click in it and press Shift +Enter . show the cell
Run the next two cells. The second uses the variable the first one made – and it would still be there in a later chapter.
Make a mistake on purpose: change len(my_dna) to len(my_DNA) and run the cell. Read the error and the 💡 hint under it, then put it back.
Getting help
Ask a demonstrator, look at the Cheat sheet , or use an AI tool your university allows – but read and test any code it gives you, and never paste personal or patient data into one.
Your name or student number (goes on your downloaded answers):
DNA as text: strings
≈ 15 min
To Python a DNA sequence is text – a string. Join, repeat, measure and count it; work out its GC content; find a TATA box and the complementary strand.
Sequences as text
DNA is a chain of nucleotides with one of four bases: A denine, C ytosine, G uanine and T hymine. Sequences are always written from the 5′ end to the 3′ end. To Python, a sequence is just text – a string , written between quotes: dna = 'AGTCGAC'. Run the first cell of this chapter’s notebook show the cell to see what you can do with one.
Activity · Strings
1.1 Join and repeat. Store dna_diagram_a followed by dna_diagram_b in dna_joined, and a sequence of 10 T, then 30 G, then 27 A in dna_long – without typing 67 letters.
GC content
G–C base pairs (three hydrogen bonds, and stronger stacking between neighbouring pairs) hold the two strands together more tightly than A–T pairs (two hydrogen bonds), so DNA with more G and C needs more heat to separate its strands. GC-rich stretches of a genome are often rich in genes. The GC content is the percentage of G and C:
GC% = 100 × (NG + NC ) / Ntotal
Activity · Count
1.2 GC content. Work out the percentage of G and C in dna with the formula above, and store it in gc. In Python, × is written *, and the brackets matter. The f-string in the cell prints it to one decimal place: {gc:.1f}.
In a new cell, work out the GC content of 10*'T' + 30*'G' + 27*'A' + 17*'C'. What is it, to one decimal place?
56.0%: 47 of its 84 letters are G or C.
Changing case, replacing letters
.upper(), .lower() and .replace(old, new) each return a new string – the original stays as it was show the cell . Methods like these belong to strings and are written after a dot.
Activity · Find and change
1.3 Highlight a TATA box. In many eukaryotic genes a TATA box, TATAAA, lies about 25–30 bases before the start of transcription: TATA-binding protein binds it and positions RNA polymerase II. Count how many times it occurs in dna (n_tata), and make highlighted: the sequence in lower case with the TATA box in upper case. Hint: dna.lower() first, then .replace() the lower-case box with the upper-case one.1.4 The complementary strand. The two strands of DNA pair A with T and G with C. Finish the code so that comp_strand is the complementary strand of strand. Replacing A with T straight away would not work – the next step would turn those T back into A – which is why the cell changes A to a lower-case t first.
More with strings
Extensions
1.5 RNA. RNA has U (uracil) where DNA has T. Make a variable rna_sequence with at least three U and no letters other than A, U, C and G.1.6 Base composition. For dna = 'AAATTGCCCC', print the percentage of each base, like 30.0% A, 40.0% C, 10.0% G, 20.0% T. Then try it on 10*'A' + 'G' + 4*'C'.1.7 Both strands. Print a sequence and its complement as two strands, with a bar between each pair:5'-AGTCGAC-3'
|||||||
3'-TCAGCTG-5' The complementary strand runs the other way: from 3′ to 5′. Text that contains a ' goes between double quotes: "5'-".
Lists, slices and proteins
≈ 15 min
Lists keep things in order. Strings and lists share the same tools – len(), positions and slices – enough to read a sequence backwards and find its reverse complement; .count() then gives the charge of a protein.
Lists
A list keeps several things in order: organisms = ['fly', 'rice', 'mouse']. Like a string, it has a length, len(organisms), and numbered positions – starting from 0 show the cell .
Organism Scientific name Chromosomes (2n) fruit fly Drosophila melanogaster 8 rice Oryza sativa 24 mouse Mus musculus 40 wheat Triticum aestivum 42 human Homo sapiens 46 potato Solanum tuberosum 48 great white shark Carcharodon carcharias 82 (crucian) carp Carassius carassius 100 red king crab Paralithodes camtschaticus 208
Activity · Lists
2.1 A list of chromosome numbers. Make a list chromosomes with the chromosome numbers of the eight organisms, in order: fly 8, rice 24, mouse 40, wheat 42, human 46, potato 48, great white shark 82, carp 100. Then add the red king crab’s 208 with .append().
Positions and slices
x[i] is the item at position i (0 is the first, −1 the last). A slice x[start:end] takes the items from start up to – but not including – end ; x[::-1] reads the whole thing backwards. Strings and lists work the same way show the cell . One difference: a list can be changed in place (fruits[1] = 'blueberry', del fruits[0]), a string cannot – strings are immutable , so you make a new one instead.
Activity · Slices
2.2 Parts of a sequence. From seq_a, take the first letter, the last letter, the first codon (three letters) and the whole sequence backwards.2.3 The reverse complement. Sequences are always written 5′→3′, so the other strand of a sequence, written the same way, is its complement read backwards : the reverse complement . Make revcomp for seq (it should be GCCAAATTTT).
seq_a = 'ATGGGGCATGCATGC': what is seq_a[6:9]?
CAT: positions 6, 7 and 8 (the 7th to 9th letters).
Protein sequences
Proteins are chains of 20 kinds of amino acid, each with a one-letter code – so a protein sequence is a string too, just with 20 letters instead of 4. The charges of the side chains say a lot about where a protein works. Apamin , a neurotoxin in bee venom that blocks a potassium channel, has the sequence CNCKAPETALCARRCQQH (UniProt P01500).
Activity · Charge
2.4 The charge of a toxin. At neutral pH, R (arginine) and K (lysine) carry a positive charge, and D (aspartic acid) and E (glutamic acid) a negative one. H (histidine) is only partly charged, but simple counts like this one include it with the positives. Count the charged residues of apamin (n_charged) and work out its overall charge (total_charge = positives − negatives). Hint: protein.count('R') counts one letter. Test your code on protein = 'RHKDER' first (6 charged, charge +2), then on apamin.
What percentage of apamin’s residues are charged, to one decimal place?
27.8%: 5 of its 18 residues.
More with lists
Extensions
2.5 Groups of amino acids. Make lists of the one-letter codes: aa_positive (R, H, K), aa_negative (D, E), aa_polar (S, T, N, Q), aa_special (C, G, P) and aa_hydrophobic (A, I, L, M, F, W, Y, V). Then join them with +, as you join strings, into all_aa: it should hold all 20 amino acids – check with len().2.6 Correct a table. Lists can be changed in place: organisms2[0] = 'fruit fly' replaces the first item. The two lists in the cell have five mistakes between them: correct them in place, using the table above.
Decisions and loops
≈ 15 min
if, elif and else let a program choose; a for loop repeats a step for every base of a sequence or every item of a list. Together they give the reverse complement of any sequence – and test it on twenty at once.
Decisions: if, elif, else
An if runs its indented lines only when its condition is true; elif (else if) tries another condition; else catches everything left. Note the colon at the end of each of these lines, the indentation (4 spaces) of the lines that belong to them, and == to compare (= stores a value).
if base == 'A':
comp = 'T'
elif base == 'T':
comp = 'A'
# … an elif for G, and one for C …
else:
raise ValueError('not a DNA base: ' + base)
A program that meets a letter it cannot handle should stop with a clear message instead of carrying on with a wrong answer: that is what raise ValueError(...) does.
Activity · Decide
3.1 One base. Complete the if … elif … else so that comp is the complement of base, and any other letter raises an error: raise ValueError('not a DNA base: ' + base). Try it with each base – and with 'W': for that one, the red ValueError is what should happen.
Loops: for
for base in sequence: runs the indented lines once for each letter, with base set to that letter show the cell . A loop that builds up a result starts from an empty string '' (or an empty list, or 0) and adds to it each time round; n += 1 is short for n = n + 1. for i in range(3): counts 0, 1, 2, and base in 'GC' tests whether a letter is one of G and C.
Activity · Loop
3.2 Backwards with a loop. Use the loop to build reverse_sequence from dna_forward: each base goes in front of what you have so far. No [::-1] this time.3.3 Reverse complement in one loop. Put the if/elif of 3.1 inside the loop – indented one more level – to work out the complement of each base of dna_seq, and add it at the front of rev_comp.3.4 Test on 20 sequences. Programs are tested on examples with known answers. Copy your loop from 3.3 into the for loop over the test cases, going through test_seq and building result, run it, and get 0 failures. The pasted lines must be indented one level more: select them and press Tab (Shift +Tab takes a level away). A loop inside a loop is called a nested loop.
The slice [::-1] reverses a sequence in one step. Why write the reverse complement as a loop at all?
The loop shows how the job is actually done, base by base, and makes room for decisions on each base – unknown letters, RNA, lower case. Most real code would use the short form (or a function from a library such as Biopython), but knowing what it does inside is what lets you check it.
More: tuples
Extensions
3.5 Tuples. The test cases are tuples , written with ( ): like lists, but they cannot be changed. Make test_data, a list of (input, expected output) tuples for 2n + 3 with the inputs 0, 1, -1, 10 and -1.5, and loop over it printing each pair.
Writing your own functions
≈ 20 min
A function packs code under a name, so that it can be used again and again – like len(). Write complement and reverse_complement, then percentage_gc and a melting-temperature predictor, each with a docstring and tests.
Your own functions
You have used Python’s functions: len(), print(), sum(). A function of your own is made with def show the cell :
def double(number):
"""Return twice the number.""" # the docstring: what it does
return 2 * number # the result
number is a parameter : when you call double(2), it is 2 inside the function.return sends the result back; a function without it returns None.Defining a function does nothing visible – run the cell, then call it.
Functions stop you copying and pasting code: if you find yourself pasting the same lines a third time, make them a function.
Activity · Complements
4.1 complement(). Turn your code from chapter 1 into a function: complement(sequence) returns the complementary sequence. Give it a docstring that says what it does – help(complement) shows it.4.2 reverse_complement() reuses complement(). Write reverse_complement(sequence) so that it calls complement() and reverses the result – no copied code. If you later teach complement about RNA, reverse_complement learns it too (try it in 4.7). Inside the function, store the complement under a new name, such as comp: complement = complement(sequence) would hide the function.
Activity · GC and melting temperature
4.3 percentage_gc() and its tests. Write percentage_gc(sequence) (chapter 1’s formula) and run the cell: the loop tests it on sequences whose GC content was worked out by hand and prints PASS or FAIL for each. Decimal numbers are rarely exactly equal, so the test compares them rounded to 3 decimal places with round().4.4 A melting temperature. When DNA is heated, its two strands separate (they ‘melt’). For sequences shorter than 14 bases, the Wallace rule gives a rough melting temperature: tm = 2 * (n_a + n_t) + 4 * (n_g + n_c). For 14 bases or more , use tm = 64.9 + 41 * (n_g + n_c - 16.4) / len(sequence) (Sambrook and Russell, 2001). Put the two formulas under the if and the else. Expected: G 4, GG 8, A 2, GA 6, 14 G 57.9, 14 A 16.9.
Why should reverse_complement call complement instead of containing a copy of its code?
So that there is one place to change and test. If complement is fixed or extended – RNA, lower case, ambiguity codes – reverse_complement gets the change automatically; two copies would drift apart.
More with functions
Extensions
4.5 A default argument. A parameter can have a default value, used when the call leaves it out – the cell has an example. Write complement(sequence, reverse=False) so that complement('AAAG') gives TTTC and complement('AAAG', reverse=True) gives CTTT.4.6 Keep lower case. Sequences often mark regions in lower case, e.g. ctgcTATAAAaggctg. Make complement() keep the case of each letter (gacgATATTTtccgac) with str.maketrans('ACGTacgt', 'TGCAtgca') and .translate() – the cell shows how they work.4.7 RNA. RNA has U where DNA has T, and U pairs with A. Teach complement() that U becomes A, so that complement('UUUCG') gives AAAGC. Then, without touching reverse_complement(), run it on 'UUUCG': it should give CGAAA. Because it calls complement(), it learnt about U too.4.8 The real thing: Biopython. Biopython does all of this, and much more. Run its cell show the cell and compare the answers with your functions: its Wallace rule is your short formula. Real primer design uses the nearest-neighbour method, which looks at each pair of neighbouring bases.
Hash tables: finding things in one step
≈ 12 min
Python’s dictionaries find a key in about one step on average, whatever their size, because they are hash tables. Build a tiny one, make a collision, plot how a filling table slows down, and compare hash functions for codons.
Finding things in one step
Looking a key up in a list means checking item after item. A hash table instead turns the key into a slot number with a hash function and goes straight there – on average one step, however big the table. Python's dictionaries ({'A': 'T'}) and sets work this way, and so do BLAST's word lists, k-mer counters and read mappers.
The simplest hash function for numbers is the remainder: value % table_size show the cell . Two keys that get the same slot collide ; with linear probing the second one goes into the next free slot, so a search has to keep looking until it finds the key or an empty slot. Table sizes are often prime numbers, like 17: then keys that share a pattern – all even, or all multiples of 10 – still spread over all the slots.
Activity · Slots
5.1 A collision. 27 sits in slot 27 % 17 = 10. Find a number between 100 and 200 that the hash function sends to the same slot – a collision – and store it in clash.5.2 Plot the cost of a filling table. With linear probing , a key whose slot is taken goes into the next free slot. The fuller the table, the longer the search: looking for a key that is not there takes about ½ (1 + 1/(1 − α)²) probes, where α (the load factor ) is the fraction of slots in use (Knuth). Plot it for α from 0 to 0.95, with labelled axes and a title.
From your plot (or the formula): about how many probes does looking for a missing key take when the table is 90% full?
About 50 (½ × (1 + 1/0.1²) = 50.5) – against 2.5 when it is half full. That is why hash tables are kept well under full and grown when they fill up.
Hashing text
Keys are usually text, so the letters are first turned into numbers with ord() show the cell . Adding them up is the obvious way – but then GUA and AGU, made of the same letters, always collide. Weighting each letter by its position, or writing the codes in binary one after another, takes the order into account: str_hash_03 writes G, U and A as 1000111, 1010101 and 1000001, reads 100011110101011000001 as one binary number, 1,174,209, and takes the remainder after dividing by 27 – slot 6.
Activity · Hash functions
5.3 Spread the codons. Count how many of the 27 slots each hash function uses for the 64 codons: n_simple is done for you; fill in n_02 and n_03 the same way.
Dictionaries in practice
You rarely write a hash function yourself: counts = {} and counts[word] = counts.get(word, 0) + 1 count words with Python’s own, very good, hash table. Chapter 9’s BLAST keeps the words of its query in a dictionary in the same way.
Graphs: breadth first and depth first
≈ 12 min
A graph is nodes joined by edges – proteins that touch, reads that overlap. Store one as adjacency lists, then visit every node in two orders: breadth first, with a queue, and depth first, with a stack.
Graphs
A graph is a set of nodes joined by edges : proteins that touch, metabolites linked by reactions, overlapping DNA reads in a genome assembly. The simplest way to store one is a list of adjacency lists : graph[i] lists the neighbours of node i show the cell . In an undirected graph each edge is listed at both ends.
a ─── b a: b, c
│ │ b: a, d
c ─── d c: a, d
d: c, b
Activity · Adjacency lists
6.1 A diagonal. In an undirected graph an edge appears twice: d in a’s list and a in d’s. Add the diagonal a–d.
Visiting every node
A traversal visits every node that can be reached from a start. Both kinds below keep a list of nodes still to visit; the main difference is which one they take next draw new_graph :
Breadth first uses a queue : the node waiting longest goes next, so the search spreads out in layers – first the start’s neighbours, then theirs show the cell .Depth first uses a stack : the node added last goes next, so the search runs as far as it can along one path before backing up show the cell .
Activity · Traversals
6.2 Breadth first. Run breadth_first on new_graph from node 0, and store the order in which it visited the nodes in bfs_order. On paper, draw arrows on the picture of new_graph in that order.6.3 Depth first. Now go depth first through new_graph from node 0: dfs_order. Add arrows of another colour to your drawing.6.4 A dead end. Remove the edges d–e and d–f – from the lists at both ends – and go depth first again. Where does the search now have to back up?
Depth first through new_graph goes a, c, e, f, d, b. Why does it reach b last, although b is a neighbour of a?
From a, both b and c go on the stack. A stack gives back the item added last, so c comes next, and its neighbours go on top of b. Each step takes the newest neighbour – c, e, f, d – so the search runs a → c → e → f → d without backing up, and b is reached last, as a neighbour of d. The b that a put on the stack is skipped later, because it has been visited by then.
More: paths
Extensions
6.5 The shortest path. Breadth first finds every node by the fewest edges. Write path_to(parents, node) that follows the parents back from node to the start and returns the path, start first: path_to(breadth_first(new_graph, 0), 5) should be [0, 1, 3, 5].
The two functions follow M. L. Hetland, Python Algorithms (2014), chapter 5.
Patterns: regular expressions
≈ 15 min
Restriction enzymes cut DNA at short patterns, some with ambiguous letters: GANTC, GGWCC. Regular expressions describe such patterns. Use them to find restriction sites – and to map them on the 4,361 base pairs of the plasmid pBR322.
Patterns with variation
Restriction enzymes cut DNA at short recognition sites. Some sites are exact (GAATTC for EcoRI), others allow choices, written with ambiguity codes: N is any base, W is A or T, R is A or G. A regular expression describes such a pattern, and Python’s re module finds it show the cell :
In the pattern Matches GAATTCexactly these letters GG(A|T)CC or GG[AT]CCGGACC or GGTCC GC[ACGT]GCGC, any base, GC .any one character (also letters that are not bases!) A+ A* A?one or more A; zero or more; zero or one (TCTG){3,4}TCTG 3 or 4 times [^ACGT]any character except A, C, G, T
re.search(pattern, text) gives the first match as a match object , or None if there is none; re.finditer gives all of them as match objects show the cell ; re.findall gives all of them as text – but if the pattern has brackets ( ), only the part in the brackets: write GG[AT]CC rather than GG(A|T)CC to count with it. Writing patterns as raw strings, r'...', keeps any backslashes in them as they are.
Activity · Find a pattern
7.1 Where HinfI cuts. The restriction enzyme HinfI cuts at GANTC, where N is any base. Write the pattern as a regular expression and store, for each sequence, whether HinfI cuts it.7.2 A run of A. Use the match object m to find the length of the run of A (polya_length).
The restriction map of a plasmid
Restriction maps show where enzymes cut a piece of DNA – the first thing to know before cloning into it.
Activity · pBR322
7.3 Read pBR322. pBR322 is one of the first cloning vectors: a circular plasmid of 4,361 base pairs see the file . Join its sequence lines – everything but the header – into one string pbr322. Hint: ''.join(lines[1:]).7.4 Restriction sites of pBR322. Count the sites of each enzyme in pbr322 with len(re.findall(pattern, pbr322)). The published map of pBR322 has one site each for HindIII (at position 29, counting from 1), EcoRV (185) and EcoRI (4359); the last lines of the cell print where the search finds the first two.7.5 A circle. The map of pBR322 shows one EcoRI site, at position 4359 – but the search found none. The plasmid is a circle and its sequence file had to start somewhere: the site is split between the end and the start. Search circle, which repeats the first bases after the end, and store where the site starts in ecori_position.
How many times does AvaII (GGWCC) cut pBR322?
Why could a search of the sequence file not find the EcoRI site of pBR322, and what should a program that searches plasmids do?
pBR322 is circular, but the file stores it as a line that has to start somewhere – here in the middle of the EcoRI site (GAA|TTC: GAA is at the end of the file, TTC at its start). A program for circular DNA should also search across the join, e.g. by adding the first few bases to the end (fewer than the length of the longest site), as Biopython does with linear=False.
More: Biopython and overlapping sites
Extensions
7.6 The real thing: Biopython. Biopython knows hundreds of restriction enzymes, both strands, and circles. Run its cell show the cell . It finds 42 BisI sites where re.findall found 35: some BisI sites overlap (GCAGCAGC holds two), and findall carries on after each match. A lookahead, (?=GC[ACGT]GC), counts overlapping matches too. Biopython also gives where each enzyme cuts , not where its site starts.
Local alignment: Smith–Waterman
≈ 12 min
Two sequences that share only a motif align badly from end to end. Smith–Waterman finds the best-matching part: it fills a matrix of scores, never letting a score drop below zero, then traces back from the highest one.
Global and local
A global alignment lines up two sequences from end to end; a local alignment finds the parts that match best and ignores the rest – what you want for two genes that share one domain, or a read against a genome. The two sequences here share a motif but little else.
An alignment is scored column by column. Here: +1 when the letters match, −1 for a mismatch or a gap.
Activity · Scoring
8.1 Score an alignment. Write simple_score_alignment(align_a, align_b): go through the two strings column by column – zip(align_a, align_b) gives the pairs – and add up simple_scoring_scheme for each.
The Smith–Waterman algorithm
Smith–Waterman fills a matrix V : V [i, j] is the best score of any alignment that ends with letter i of one sequence and letter j of the other. Each cell is the best of four choices: the cell diagonally above-left plus the score of aligning the two letters; the cell above or to the left plus a gap; or 0 – start a new alignment here instead of carrying a bad one on. That 0 is what makes it local the matrix .
The alignment is read by the traceback : start at the highest score and step back to whichever neighbour it came from, until a 0 the traceback .
What is the highest score in the Smith–Waterman matrix of seq_a and seq_b?
8: the local alignment TATGCTGGCG / TATGC-GGCG has 9 matches and one gap.
Activity · Align
8.2 One more mutation. Align seq_a with seq_x, which differs from seq_b at one position (C instead of G), and look at the score and the alignment.
Real aligners use the same algorithm with better scores: a table of scores for each pair of letters, and gaps that are expensive to open but cheaper to extend. Biopython’s aligner, scored like the EMBOSS tools needle (global) and water (local), finds the same local alignment run it .
Run the Biopython cell. How do the global and the local alignment differ, and which suits these two sequences?
The global alignment has to place every letter of both sequences, so it contains the poorly matching ends (with mismatches and end gaps); the local one keeps only TATGCTGGCG / TATGC-GGCG, the part the sequences share. For two sequences that only share a motif, the local alignment is the meaningful one; both tools agree with your Smith–Waterman on that part.
More: the price of a gap
Extensions
8.3 Expensive gaps. Write gap4_scoring: like simple_scoring_scheme, but −4 for a gap. Run Smith–Waterman with it. What happens to the local alignment, and why?
Database search: a small BLAST
≈ 12 min
BLAST does not align the query with every database sequence. It first finds short shared words – hits – with a hash table of the query, and extends only those. Run a small BLAST step by step on a database of five sequences.
How BLAST is fast
Smith–Waterman against every sequence of a database would take far too long. BLAST (the Basic Local Alignment Search Tool) uses a shortcut – seeds :
Make a dictionary of every word of w letters in the query, with its positions (a hash table, as in chapter 5).
Go through each database sequence, looking up each of its words in the dictionary: a word found there is a hit .
Extend only the hits, left and right, while the alignment stays good; keep the best extension.
This small version (the MyBlast program of Rocha and Ferreira’s Bioinformatics Algorithms ) searches a database of five DNA sequences, seqBlast.txt , with words of 11 letters – the classic word size for DNA (megablast, BLAST’s default for DNA today, uses 28) – and counts identical letters instead of using a score table show the cell .
Activity · Seeds
9.1 Hits. For each sequence in db, count the hits – words of 11 letters it shares with the query – and collect the numbers in hits_per_sequence.
Extending a hit show the cell : the alignment grows while at least half of the newly added letters match; it never has gaps. Real BLAST scores with a substitution table, stops when the score falls a set amount below its best (the X-drop), allows gaps in a second step, and works out how often a score that good would turn up by chance (the E-value).
Activity · The best alignment
9.2 Show the alignment. Print the aligned parts of the query and of the best sequence, with highlight_line between them. Use the numbers in result: the query part is query[query_start:query_start + size]. It should look like this:cgacgacggcgacgaatgatg
|||||||| ||||||||||||
cgacgacgacgacgaatgatg
Why does BLAST look the database words up in a dictionary of the query’s words, instead of aligning the query with every database sequence?
A dictionary look-up takes one step, so finding the hits takes one pass through the database; only the few places with a hit are then aligned. Full dynamic programming against every sequence would cost query length × database length steps. The price: an alignment without any shared word of w letters is missed.
More: word size
Extensions
9.3 Word size. Run the cell: it searches with words of 8, 11 and 14 letters. Store the result for 14 in result_14. Why does 14 find nothing, and what do smaller words cost?
Multiple alignments: scores, profiles and consensus
≈ 15 min
A multiple alignment lines up a family of sequences. BLOSUM62 scores each pair of amino acids; a profile counts the residues of each column; a consensus takes the majority; the sum of pairs scores the whole alignment.
Scoring amino acids: BLOSUM62
Related proteins differ, but not at random: an isoleucine is often replaced by a leucine, rarely by a glycine. A substitution table gives a score for every pair of amino acids – positive for pairs often seen in related proteins, negative for rare ones. BLOSUM62 , the default of BLAST for proteins, was counted from conserved blocks of related proteins (Henikoff and Henikoff, 1992). Each score is a log-odds: twice the logarithm (base 2) of how much more often the pair is seen in related proteins than by chance – which is why the scores along an alignment can simply be added. Here it is a dictionary of dictionaries show the cell .
Activity · Scores for pairs
10.1 Most similar to W. Look up BLOSUM62['W'][aa] for every amino acid aa except W, and keep the one with the highest score in most_similar_to_w.10.2 Gaps in a pair. Complete pairwise_score_gapped: −1 when one of the two is a gap, 0 when both are gaps (in a multiple alignment two sequences can both have a gap in a column), else the table score.
A family of sequences
A multiple alignment lines up a whole family. Crambin is a small, very stable protein from the seeds of Abyssinian kale; seven related proteins from other plants are aligned with it show the cell . Three ways to summarise such an alignment:
a profile : for each column, the fraction of each residue – used to search for more members of a family (PSI-BLAST, HMMER);
a consensus : the most common residue of each column, or X where no residue is in at least half the sequences;
a conservation line under the alignment, in Clustal’s style: * for (nearly) fully conserved columns, : for similar residues, . for weakly similar ones.
Activity · Summarise
10.3 Consensus. Run the cell that prints the alignment with its conservation line and consensus show the cell . Then store the consensus with a stricter threshold – a residue must be in at least 75% of the sequences – in strict_consensus.
Compare column 33 of the alignment (counting from 1), IKKLIKKQ, with column 18, LLFFFVVI. Why does column 33 get a consensus letter (K) but no conservation mark, while column 18 gets X but a .?
The consensus only counts the most common residue: K is in 4 of the 8 sequences – half, enough for the 0.5 threshold. Conservation looks at all the residues of a column, weighted by how similar BLOSUM62 says they are: column 33 mixes K with I and L, which score badly with it (−3 and −2; only Q is similar, +1), so its conservation is low (0.13), while column 18 holds only L, F, V and I – all hydrophobic, and similar to each other (0.28: weakly similar).
Scoring a whole alignment
Programs that build multiple alignments (Clustal, MUSCLE, MAFFT) need a way to tell a better alignment of the same sequences from a worse one; the sum of pairs is the classic score. It grows with the number of sequences, so it compares alignments of the same set only.
Activity · Sum of pairs
10.4 Sum of pairs. The sum of pairs scores a whole multiple alignment: the gapped pairwise score of every pair of sequences, added up. Use two loops – the second starting after the first (for j in range(i + 1, len(alignment))) – so that each pair counts once. The crambin alignment scores 2964.
Phylogenetic trees with Biopython
≈ 15 min
Distances between aligned sequences give a tree: the closer two sequences, the closer their branches. Build trees of human tRNAs and of myoglobins with Biopython – and identify a mystery animal from its myoglobin.
From distances to a tree
Similar sequences usually share a recent ancestor. A distance method measures how different each pair of aligned sequences is – here simply the fraction of positions that differ – and joins the closest first. UPGMA repeatedly merges the two closest groups, averaging their distances to the rest; neighbour joining corrects for lineages that change at different speeds. Distance methods are quick. The main alternatives – maximum parsimony (the tree that needs the fewest changes) and maximum likelihood or Bayesian methods (the tree that best explains the alignment under a model of evolution) – are slower, but use every position of the alignment directly.
First, 6 of the 20 threonine tRNA genes of the human genome, aligned see the file . Each carries an anticodon that reads threonine codons: AGU (reading ACU), CGU (ACG) or UGU (ACA) – written AGT, CGT and TGT in the gene names the alignment the UPGMA tree . Four codons code for threonine, but there are only three anticodons: no human tRNA gene has GGU, because the AGU tRNA, its A modified to inosine, also reads ACC by wobble pairing.
Activity · tRNAs
11.1 Neighbour joining. Build the tree of the tRNAs again with the neighbour-joining method ('nj' instead of 'upgma') and draw it. UPGMA assumes every lineage changes at the same rate and gives a rooted tree; neighbour joining does not assume that, and its tree has no real root (Biopython draws it from an arbitrary point).
Look at the anticodons in the names. Do the tRNAs group by anticodon?
Not entirely: TGT-2, TGT-3 and TGT-4 form one group, but TGT-1-1 sits with CGT-1-1 and AGT-3-1, from which it differs at only 4 and 8 of the 74 positions (against 22 for the other TGT genes). The anticodon is a few bases; the rest of the molecule, which must work with the same enzymes, shows that TGT-1-1 is more closely related to the AGT and CGT genes – a single change in the anticodon can switch a tRNA from one codon to another.
Myoglobins and a mystery
Activity · Myoglobins
11.2 A tree of myoglobins. Myoglobin stores oxygen in muscle, in nearly all vertebrates. The 8 myoglobins – human, whale, platypus, horse, gummy shark, blue marlin, tuna and carp – were aligned with MUSCLE see the alignment . Build and draw a UPGMA tree of them.11.3 Seal or fish? An animal washed up on a Norfolk beach was too decomposed for DNA testing, but a myoglobin sequence could be read from its muscle protein. It might be a seal. Draw the tree with the unknown myoglobin.11.4 The closest myoglobin. Find the myoglobin with the smallest distance to the unknown one and store its name in closest.
More: search instead of a tree
Extensions
11.5 Search with the small BLAST. If you have run chapter 9’s cells in this session, its best_alignment is still in Python. The last cell searches the unknown myoglobin against the eight with words of 3 letters (BLAST’s word size for proteins) show the cell . Run it: which fish comes top this time, and why could that differ from the tree? (Carp: BLAST scores the best ungapped stretch, the tree the whole alignment.)
Strictly, Biopython’s 'upgma' averages the distances of the two groups it merges without weighting them by their sizes (the method called WPGMA). For these data the trees have the same shape either way.
Networks: the proteins of the ribosome
≈ 12 min
NetworkX stores graphs as dictionaries of dictionaries – chapter 6’s adjacency lists – and adds hundreds of algorithms. Build a graph, then explore the network of contacts between the proteins of the ribosome: its subunits, its hubs, and how it differs from chance.
NetworkX
NetworkX stores a graph as a dictionary of dictionaries – chapter 6’s adjacency lists – and adds hundreds of algorithms: shortest paths, clusters, centrality … Nodes can be anything (names, numbers), and nodes and edges can carry attributes show the cell .
A graph of your own
Extensions
12.1 Build a graph. Add the nodes 1 to 50 (range(1, 51)) with their attribute kind, then the edges between consecutive numbers. integer_graph.add_node(3, kind='odd') adds one node with its attribute.
The proteins of the ribosome
The ribosome builds proteins from the messenger RNA’s code. It is mostly RNA, decorated with some 80 proteins, many of which reach out to touch each other with long extensions. Poirot and Timsit (2016) listed these contacts in the structure of the yeast ribosome; turned into a network, they show how its proteins could pass signals build the network .
STM1 is not one of the ribosome’s own proteins but a factor that binds ribosomes when yeast cells are starved; all its contacts here are with small-subunit proteins, so the S in its name happens to put it in the right group.
Activity · Subunits
12.2 Two subunits. The ribosome is two large particles of RNA and protein, the large and the small subunit, with a cleft between them for the messenger RNA. Make large_node_list – the proteins without an S in their name – and draw them in another colour.
Hubs and chance
Activity · Degree
12.3 The hubs. The degree of a node is its number of edges. Store the names of all the proteins with the highest degree in best_connected (a list or a set – there is a tie).
What is the average degree of the ribosome network, to two decimal places?
3.71: 126 edges among 68 proteins, and 2 × 126 / 68 = 3.71 – every edge counts once for each of its two ends.
In a directed graph (nx.DiGraph(): signals, regulation, metabolic reactions) each node has an in-degree and an out-degree; when the edges carry weights, the weighted degree adds up the weights of a node’s edges (G.degree(weight='weight')).
Is this network special, or would random contacts look the same? Compare it with a random network of the same size and density compare .
How does the ribosome network differ from the random one?
Their degrees are much alike (on average 3.7, at most 8 in both), but the ribosome network has far more triangles – three proteins that all touch each other: 45, against 6 in this random network (about 8 on average over many random networks of this size). Proteins touch each other in groups, as neighbours on the surface of a particle do; random contacts rarely close a triangle. The random network also falls into 3 pieces, two of them single proteins that touch nothing; the ribosome network is one piece – but only because, for this practical, the few small separate pieces in the paper’s table were left out.
Triangles measure clustering . Many real networks are “small worlds”: clustered like this, yet with short paths between any two nodes; nx.watts_strogatz_graph makes random networks of that kind to compare with.
More: across the cleft, cliques
Extensions
12.4 Across the cleft. Which proteins touch the other subunit? Find the edges whose two ends are in different subunits (GBP and STM1 go with the small subunit), and store the protein of highest degree among their ends in cross_hub. Hint: max(ends, key=contacts_graph.degree) picks the end with the highest degree.12.5 Cliques. A clique is a group in which every node touches every other. list(nx.find_cliques(contacts_graph)) lists them; store the size of the largest in largest_clique.
The first part of this chapter follows the NetworkX workshop of Jordi Torrents (PyData Barcelona, 2017). Data: Poirot O. and Timsit Y. (2016) Neuron-like networks between ribosomal proteins within the ribosome. Scientific Reports 6: 26485.
Summary
≈ 5 min
From strings to networks: what each chapter showed, and where the same algorithms work inside the tools of bioinformatics.
What you did, and where it is used
Chapter You used Inside 1–4 strings, lists, if, loops, functions every bioinformatics script; Biopython’s Seq does the same jobs 5 hash tables (dictionaries, sets) BLAST’s word lists, k-mer counters, read mappers, genome assemblers 6 breadth-first and depth-first search genome assembly graphs, shortest paths in networks, ontologies 7 regular expressions motif searches (PROSITE patterns), restriction maps, parsing files 8 Smith–Waterman EMBOSS water, SSEARCH, the last step of BLAST and of read mappers 9 seed and extend BLAST, DIAMOND, read mappers (BWA, minimap2) 10 substitution tables, profiles, sum of pairs Clustal, MUSCLE, MAFFT, PSI-BLAST, HMMER 11 distance trees (UPGMA, neighbour joining) Clustal’s guide trees, quick phylogenies; MEGA, PHYLIP 12 networks STRING, Cytoscape, pathway and interaction analysis
Activity · Keep your work
Download your answers: My answers → Download my answers (top right).
Download your notebook, all chapters in one: the ⤓ button above the notebook → All chapters in one Jupyter notebook . It opens in Jupyter, VS Code or Google Colab.
The Cheat sheet tab has the Python of this practical on one page. To practise more, the extensions marked Optional in each chapter are a good start; then try your own sequences – everything here works on real data.