!pip install -q biopythonSession 5 — BLAST and Profile Methods in code
This notebook reproduces every worked example on the Session 5 page with real code and real data: the BLOSUM62 neighbor-word scores from BLAST’s seeding stage, the X-drop extension rule, a BLASTP search of HBB (beta-globin) against its globin relatives together with the optimal alignment that BLAST’s heuristic missed, a PSSM built from scratch (with and without pseudocounts), and a live fetch of HBB’s real Pfam profile-HMM domain hit.
Setup
Stage 1: seeding — BLOSUM62 neighbor words
For the query word LWG, which other 3-letter words score above a threshold T=12 against it, using the real BLOSUM62 matrix?
from Bio.Align import substitution_matrices
from itertools import product
blosum62 = substitution_matrices.load("BLOSUM62")
AA = sorted(set(blosum62.alphabet) - {"*", "X", "B", "Z"})
def word_score(w1, w2):
return sum(blosum62[a, b] for a, b in zip(w1, w2))
query_word = "LWG"
T = 12
scores = {"".join(w): word_score(query_word, "".join(w)) for w in product(AA, repeat=3)}
neighbors = sorted(scores.items(), key=lambda kv: -kv[1])
above_threshold = [(w, s) for w, s in neighbors if s >= T]
print(f"{len(above_threshold)} words score >= {T} against {query_word!r}; top 10:")
for w, s in above_threshold[:10]:
print(f" {w} {s}")49 words score >= 12 against 'LWG'; top 10:
LWG 21.0
IWG 19.0
MWG 19.0
VWG 18.0
FWG 17.0
AWG 16.0
CWG 16.0
TWG 16.0
YWG 16.0
KWG 15.0Stage 2: extension — the X-drop rule
Track a running score and a running maximum; stop as soon as the running score falls more than X below the running maximum, then trim back to where that maximum actually occurred.
def xdrop_extend(query, target, X):
q = query.replace(" ", "")
t = target.replace(" ", "")
n = min(len(q), len(t))
score = 0
running_max = 0
max_pos = 0
trace = []
for i in range(n):
score += 1 if q[i] == t[i] else -1
if score > running_max:
running_max = score
max_pos = i
trace.append(score)
if running_max - score >= X:
return {
"stopped_at": i,
"peak_score": running_max,
"peak_position": max_pos,
"trimmed_query": q[: max_pos + 1],
"trimmed_target": t[: max_pos + 1],
"trace": trace,
}
return {"stopped_at": None, "peak_score": running_max, "trace": trace}
result = xdrop_extend(
"The quick brown fox jumps over the lazy dog.",
"The quiet brown cat purrs when she sees him.",
X=5,
)
print("Extension stopped at character index:", result["stopped_at"])
print("Peak running score:", result["peak_score"], "at position", result["peak_position"])
print("Trimmed alignment:")
print(" query:", result["trimmed_query"])
print(" target:", result["trimmed_target"])
assert result["trimmed_query"] == "Thequickbrown"
assert result["trimmed_target"] == "Thequietbrown"Extension stopped at character index: 19
Peak running score: 9 at position 12
Trimmed alignment:
query: Thequickbrown
target: ThequietbrownStage 3: a real BLAST search, and what its heuristic missed
The book page’s table comes from a BLASTP search of HBB (beta-globin) against five other human globin-family proteins. It was run with BLAST+, the command-line version of BLAST, when this notebook was written. Its output is reproduced below, so that the notebook needs no command-line tools. The commands are given in a comment for reference.
The sequences themselves are fetched live from UniProt, because the next step aligns two of them again in Python.
import requests
GLOBINS = {
"P68871": "HBB_HUMAN (beta-globin) -- the query",
"P69905": "HBA_HUMAN (alpha-globin)",
"P02144": "MYG_HUMAN (myoglobin)",
"P02008": "HBAZ_HUMAN (zeta-globin, embryonic alpha-like)",
"P02042": "HBD_HUMAN (delta-globin)",
"P02100": "HBE_HUMAN (epsilon-globin, embryonic)",
}
fasta = {}
for acc, label in GLOBINS.items():
r = requests.get(f"https://rest.uniprot.org/uniprotkb/{acc}.fasta", timeout=10)
r.raise_for_status()
fasta[acc] = r.text
print(label)HBB_HUMAN (beta-globin) -- the query
HBA_HUMAN (alpha-globin)
MYG_HUMAN (myoglobin)
HBAZ_HUMAN (zeta-globin, embryonic alpha-like)
HBD_HUMAN (delta-globin)
HBE_HUMAN (epsilon-globin, embryonic)# Run with BLAST+ 2.16.0, on the five sequences fetched above:
# makeblastdb -in globin_db.fasta -dbtype prot -out globin_db
# blastp -query hbb_query.fasta -db globin_db -outfmt "6 sseqid pident length evalue bitscore" -evalue 1000
# blastp ... -window_size 0 -max_hsps 1 (the same search, with the two-hit rule switched off)
BLAST_DEFAULT = """sp|P02042|HBD_HUMAN 93.197 147 9.28e-105 284
sp|P02100|HBE_HUMAN 75.510 147 2.05e-87 240
sp|P69905|HBA_HUMAN 43.448 145 8.95e-38 114
sp|P02008|HBAZ_HUMAN 35.862 145 3.81e-32 100
sp|P02144|MYG_HUMAN 16.000 125 63 10.8"""
BLAST_NO_TWO_HIT = """sp|P02144|MYG_HUMAN 25.517 145 5.38e-11 45.4"""
def show_blast(table):
print(f"{'Hit':<20}{'%identity':>10}{'length':>8}{'E-value':>14}{'bit score':>11}")
for line in table.splitlines():
sseqid, pident, length, evalue, bitscore = line.split("\t")
print(f"{sseqid.split('|')[2]:<20}{float(pident):>10.1f}{length:>8}{float(evalue):>14.3g}{float(bitscore):>11.1f}")
show_blast(BLAST_DEFAULT)Hit %identity length E-value bit score
HBD_HUMAN 93.2 147 9.28e-105 284.0
HBE_HUMAN 75.5 147 2.05e-87 240.0
HBA_HUMAN 43.4 145 8.95e-38 114.0
HBAZ_HUMAN 35.9 145 3.81e-32 100.0
MYG_HUMAN 16.0 125 63 10.8Four confident hits, and myoglobin with an E-value of about 63. Taken at face value, that is no statistical evidence at all, although myoglobin is a known remote homolog. But an E-value only evaluates the alignment BLAST found. BLAST’s default two-hit rule only extends a seed when a second seed lies nearby on the same diagonal. Is there a better alignment that BLAST never extended? Smith-Waterman (Session 4) finds the best local alignment by full dynamic programming, with no shortcuts. Here it is with BLAST’s own scoring (BLOSUM62, a gap costing 11 to open and 1 per residue):
import math
from Bio import Align
from Bio.Align import substitution_matrices
def seq_of(acc):
return "".join(fasta[acc].splitlines()[1:])
aligner = Align.PairwiseAligner(mode="local", substitution_matrix=substitution_matrices.load("BLOSUM62"),
open_gap_score=-12, extend_gap_score=-1) # a gap of length k costs 11 + k, as in BLAST
best = aligner.align(seq_of("P68871"), seq_of("P02144"))[0]
identity = sum(a == b and a != "-" for a, b in zip(best[0], best[1]))
length = best.shape[1]
lam, K = 0.267, 0.041 # BLAST's Karlin-Altschul parameters for BLOSUM62 with these gap costs
bits = (lam * best.score - math.log(K)) / math.log(2)
m = len(seq_of("P68871"))
n = sum(len(seq_of(acc)) for acc in GLOBINS if acc != "P68871")
E = m * n * 2 ** -bits # E-value: the search space times 2 to the minus bit score
print(f"Smith-Waterman, HBB against MYG: raw score {best.score:.0f}, {bits:.1f} bits, "
f"{length} columns, {100 * identity / length:.1f}% identity, E = {E:.1g}")
print()
print("BLAST with the two-hit rule switched off (-window_size 0) finds the same alignment:")
show_blast(BLAST_NO_TWO_HIT)Smith-Waterman, HBB against MYG: raw score 110, 47.0 bits, 145 columns, 25.5% identity, E = 8e-10
BLAST with the two-hit rule switched off (-window_size 0) finds the same alignment:
Hit %identity length E-value bit score
MYG_HUMAN 25.5 145 5.38e-11 45.4The optimal alignment covers almost the whole of both proteins (145 residues, about 25% identity) and is highly significant. BLAST with the two-hit rule switched off finds it too, with almost the same bit score; its E-value is a little smaller because BLAST also corrects for sequence composition and for the ends of the sequences. With the default settings no two seeds happened to land close together on the right diagonal, so the good alignment was never extended: the price of a fast heuristic. More distant relatives than myoglobin give no pairwise signal even with an optimal alignment, which is the case profile methods are built for.
Building a PSSM from scratch
The same 5-sequence toy alignment from the book page, with and without pseudocounts.
import math
msa = ["AGLSP", "AGLTP", "RGISP", "AALSQ", "AALSP"]
# Typical background frequencies for these residues in proteins
background = {
"A": 0.074, "R": 0.042, "G": 0.074, "L": 0.076, "S": 0.081,
"T": 0.062, "P": 0.050, "Q": 0.037, "I": 0.038, "W": 0.013,
}
residues = sorted(r for r in background if r != "W") # residues that actually appear in this MSA
n_pos = len(msa[0])
def counts_at(pos):
col = [seq[pos] for seq in msa]
return {r: col.count(r) for r in residues}
def pssm(pseudocount):
matrix = []
for pos in range(n_pos):
c = counts_at(pos)
total = sum(c.values()) + pseudocount * len(residues)
freqs = {r: (c[r] + pseudocount) / total for r in residues}
matrix.append({r: math.log2(freqs[r] / background[r]) if freqs[r] > 0 else float("-inf") for r in residues})
return matrix
def best_sequence(matrix):
best = [max(col, key=col.get) for col in matrix]
score = sum(matrix[i][best[i]] for i in range(len(best)))
return "".join(best), score
pssm_pseudo = pssm(pseudocount=1)
best_seq, best_score = best_sequence(pssm_pseudo)
print("With pseudocounts, highest-scoring sequence:", best_seq, f"(score {best_score:.2f})")
assert best_seq == "AGLSP", "should recover the alignment's own consensus"With pseudocounts, highest-scoring sequence: AGLSP (score 11.43)# Why pseudocounts matter: a rare residue occurring once can outscore
# a common residue occurring four times, purely because it's rarer overall.
a_count, w_count, n = 4, 1, 5
score_A = math.log2((a_count / n) / background["A"])
score_W = math.log2((w_count / n) / background["W"])
print(f"A (appears in 4/5 sequences): log-odds = {score_A:.2f}")
print(f"W (appears in 1/5 sequences): log-odds = {score_W:.2f}")
assert score_W > score_A, "the rare residue should outscore the common one without pseudocounts"A (appears in 4/5 sequences): log-odds = 3.43
W (appears in 1/5 sequences): log-odds = 3.94A real profile-HMM hit on HBB
Pfam’s profile HMMs, queried live via EBI’s InterPro API, for the same HBB UniProt entry used throughout this course.
r = requests.get("https://www.ebi.ac.uk/interpro/api/entry/pfam/protein/uniprot/P68871/", timeout=10)
r.raise_for_status()
entry = r.json()["results"][0]
meta = entry["metadata"]
loc = entry["proteins"][0]["entry_protein_locations"][0]
frag = loc["fragments"][0]
print("Pfam accession: ", meta["accession"])
print("Domain name: ", meta["name"])
print("Matched region: ", f"residues {frag['start']}-{frag['end']} of {entry['proteins'][0]['protein_length']}")
print("E-value: ", loc["score"]) # InterPro calls the domain E-value "score"Pfam accession: PF00042
Domain name: Globin
Matched region: residues 27-142 of 147
E-value: 1.3e-29