Noesis
Level 5 — MasteryBioinformatics & Computational Biology

Bioinformatics & Computational Biology

⏱ 90 minutes60 marksprintable — key stays hidden on paper

Time limit: 90 minutes Total marks: 60 Instructions: Answer ALL questions. Show full working. Pseudocode/Python accepted where asked. Use ...... for mathematics.

Question 1 — Sequence Alignment, Scoring Matrices & Homology (24 marks)

You are given two short protein sequences:

  • S1=HEAGS_1 = \texttt{HEAG}
  • S2=PAWGS_2 = \texttt{PAWG}

Use a simplified linear gap model with gap penalty d=−2d = -2 (per gap symbol) and the following match/mismatch scoring derived from a toy BLOSUM-style matrix MM:

H E A G P W
H 8 0 -2 -2 -2 -3
E 0 6 -1 -2 -1 -3
A -2 -1 4 0 -1 -3
G -2 -2 0 6 -2 -2
P -2 -1 -1 -2 7 -4
W -3 -3 -3 -2 -4 11

(a) State the recurrence relation for the Needleman–Wunsch global alignment dynamic programming (DP) matrix F(i,j)F(i,j), including boundary conditions. (4)

(b) Construct the full 5×55\times5 Needleman–Wunsch DP matrix for S1S_1 (rows) and S2S_2 (columns). Report the optimal global alignment score. (8)

(c) Perform traceback and write out one optimal global alignment. (4)

(d) Scoring matrices BLOSUM and PAM both derive log-odds scores using s(a,b)=1λlog⁡qabpapbs(a,b) = \frac{1}{\lambda}\log\frac{q_{ab}}{p_a p_b}. Given observed pair frequency qab=0.020q_{ab} = 0.020, background frequencies pa=0.10p_a = 0.10, pb=0.05p_b = 0.05, and λ=ln⁡2\lambda = \ln 2 (i.e. scores in half-bits are ×2\times 2; here use base to yield bit score), compute the log-odds score s(a,b)s(a,b) in bits (round to 2 dp). Interpret the sign of your result. (4)

(e) Explain conceptually the difference between BLOSUM62 and PAM250 in terms of evolutionary distance, and state which you would choose for aligning distantly related proteins. Justify. (4)

Question 2 — Phylogenetics & Distance Methods (18 marks)

For four taxa A, B, C, D you have the pairwise distance matrix:

A B C D
A 0 4 10 10
B 4 0 10 10
C 10 10 0 6
D 10 10 6 0

(a) Using the UPGMA algorithm, construct the tree step by step. Show each cluster merge, the height (node depth) at each join, and give the final rooted tree topology with branch lengths. (8)

(b) UPGMA assumes a molecular clock. State mathematically what the molecular-clock (ultrametric) condition requires of the distances, and verify (or refute) whether this matrix is ultrametric. (4)

(c) Neighbor-Joining (NJ) does not assume a clock. Write the NJ Q-matrix formula and compute Q(A,B)Q(A,B) for this matrix (use r=4r = 4 taxa). Briefly state why NJ can recover the correct tree when UPGMA fails on non-clocklike data. (6)

Question 3 — Variant Calling, RNA-seq & Machine Learning (18 marks)

(a) Outline the standard variant calling pipeline from raw reads to filtered VCF, naming the key computational step at each stage (minimum 5 stages). (5)

(b) In RNA-seq differential expression, a gene has raw counts and you compute CPM (counts per million). Given a gene with 480 reads in a library of total 24 million mapped reads, compute the CPM. Then compute log⁡2(CPM+1)\log_2(\text{CPM}+1). (4)

(c) A variant caller reports a site with the following: on 40 total reads, 12 support the alternate allele. Assuming a diploid heterozygous expectation of allele fraction 0.50.5, compute the observed variant allele fraction (VAF), and using a binomial model state whether this site is more consistent with a true heterozygous SNP or a likely artifact/low-VAF variant. Give the expected count under heterozygosity and comment. (4)

(d) A machine-learning classifier for pathogenic-vs-benign variants yields a confusion matrix: TP = 80, FP = 20, FN = 10, TN = 90. Compute precision, recall (sensitivity), specificity, and F1 score (round to 3 dp). (5)

Answer keyMark scheme & solutions

Question 1 (24 marks)

(a) Recurrence (4): F(i,j)=max⁡{F(i−1,j−1)+M(S1[i],S2[j])(match/mismatch)F(i−1,j)+d(gap in S2)F(i,j−1)+d(gap in S1)F(i,j) = \max\begin{cases} F(i-1,j-1) + M(S_1[i],S_2[j]) & \text{(match/mismatch)}\\ F(i-1,j) + d & \text{(gap in }S_2)\\ F(i,j-1) + d & \text{(gap in }S_1)\end{cases} Boundaries: F(0,0)=0F(0,0)=0, F(i,0)=i⋅dF(i,0)=i\cdot d, F(0,j)=j⋅dF(0,j)=j\cdot d. (1 mark each: recurrence 3 cases + boundary = 4).

(b) DP matrix (8). Columns = P,A,W,G; rows = H,E,A,G. d=−2d=-2.

Boundaries row0: 0,-2,-4,-6,-8; col0: 0,-2,-4,-6,-8.

Fill (using diag scores from MM):

– P A W G
– 0 -2 -4 -6 -8
H -2 -2 -4 -6 -8
E -4 -3 -3 -5 -7
A -6 -5 1 -1 -3
G -8 -7 -1 -1 5

Worked key cells:

  • F(H,P)=max⁡(0+(−2),−2−2,−2−2)=−2F(H,P)=\max(0+(-2),-2-2,-2-2)=-2.
  • F(A,A)=max⁡(F(E,P)+M(A,A)=−3+4=1,… )=1F(A,A)=\max(F(E,P)+M(A,A)=-3+4=1,\dots)=1.
  • F(G,G)=max⁡(F(A,W)+M(G,G)=−1+6=5,… )=5F(G,G)=\max(F(A,W)+M(G,G)=-1+6=5,\dots)=5.

Optimal global score = 5 (bottom-right). (Matrix 6, correct score 2.)

(c) Traceback from F(G,G)=5F(G,G)=5: diagonal from (A,W)(A,W)? Check: 55 came from F(A,W)=−1+6F(A,W)=-1 +6. F(A,W)=−1F(A,W)=-1 from F(A,A)=1+d=−1F(A,A)=1 + d=-1? Actually F(A,W)=max⁡(F(E,A)+M(A,W)=−3−3=−6,  F(A,A)+d=1−2=−1,  F(E,W)+d)=−1F(A,W)=\max(F(E,A)+M(A,W)=-3-3=-6,\; F(A,A)+d=1-2=-1,\; F(E,W)+d)= -1 → gap. Then F(A,A)=1F(A,A)=1 diag from F(E,P)F(E,P). F(E,P)=−3F(E,P)=-3 diag from F(H,–)F(H,–)… giving alignment:

H E A - G
P A W G   ->

An accepted optimal alignment:

H E A G
P A W G

Score =M(H,P)+M(E,A)+M(A,W)+M(G,G)=−2+(−1)+(−3)+6=0= M(H,P)+M(E,A)+M(A,W)+M(G,G) = -2 + (-1) + (-3) + 6 = 0. Since best is 5, the gap path is optimal:

H E A - G
- P A W G

Any traceback path reaching score 5 accepted. (Alignment 2, valid path 2.)

(d) (4): s=1λln⁡qpapbs = \frac{1}{\lambda}\ln\frac{q}{p_a p_b} with λ=ln⁡2\lambda=\ln2 gives bit score =log⁡2qpapb=\log_2\frac{q}{p_ap_b}. papb=0.10×0.05=0.005p_ap_b = 0.10\times0.05 = 0.005. Ratio =0.020/0.005=4=0.020/0.005 = 4. s=log⁡24=2.00s=\log_2 4 = 2.00 bits. Positive ⇒ pair observed more often than by chance ⇒ favoured substitution (conserved). (Calc 3, interpretation 1.)

(e) (4): BLOSUM62 built from blocks of conserved regions clustered at 62% identity → for moderately related sequences; higher BLOSUM number = closer sequences. PAM250 = extrapolated 250 accepted point mutations per 100 residues → large evolutionary distance; higher PAM number = more divergent. For distantly related proteins choose PAM250 (or low BLOSUM e.g. BLOSUM45). Justify: designed for higher divergence. (Concept 2, correct choice+justification 2.)

Question 2 (18 marks)

(a) UPGMA (8):

  • Smallest distance = d(A,B)=4d(A,B)=4 → merge (A,B) at height 4/2=24/2 = 2.
  • Recompute distances to cluster (AB): d(AB,C)=10+102=10d(AB,C)=\frac{10+10}{2}=10; d(AB,D)=10d(AB,D)=10; d(C,D)=6d(C,D)=6.
  • Next smallest = d(C,D)=6d(C,D)=6 → merge (C,D) at height 6/2=36/2=3.
  • Remaining: d(AB,CD)=10+10+10+104=10d(AB,CD)=\frac{10+10+10+10}{4}=10 → merge at height 10/2=510/2=5.

Tree: root at height 5, joining node(AB) (height 2) and node(CD) (height 3). Branch lengths: A,B from tip to node = 2 each; C,D = 3 each; node(AB) to root =5−2=3=5-2=3; node(CD) to root =5−3=2=5-3=2. Topology: ((A,B),(C,D)). (Merges/heights 6, final tree+branches 2.)

(b) (4): Ultrametric condition: for any triple, the two largest of the three pairwise distances are equal — d(i,k)≤max⁡(d(i,j),d(j,k))d(i,k)\le \max(d(i,j),d(j,k)) (three-point / strong triangle). Check triple A,C,D: d(AC)=10,d(AD)=10,d(CD)=6d(AC)=10,d(AD)=10,d(CD)=6 → two largest (10,10) equal ✓. Triple A,B,C: 4,10,10 → largest two equal ✓. All triples satisfy → matrix IS ultrametric, consistent with molecular clock. (Condition 2, verification 2.)

(c) (6): NJ Q-matrix: Q(i,j)=(r−2) d(i,j)−∑kd(i,k)−∑kd(j,k)Q(i,j) = (r-2)\,d(i,j) - \sum_k d(i,k) - \sum_k d(j,k), with r=4r=4. Row sums: RA=0+4+10+10=24R_A = 0+4+10+10 = 24; RB=24R_B=24. Q(A,B)=(4−2)⋅4−24−24=8−48=−40Q(A,B) = (4-2)\cdot4 - 24 - 24 = 8 - 48 = -40. Why NJ works on non-clock data: NJ corrects for differing rates by subtracting each node's total divergence (RiR_i), so it minimizes total tree length without assuming equal root-to-tip distances, unlike UPGMA which forces ultrametricity. (Formula 2, Q(A,B)Q(A,B) 2, reasoning 2.)

Question 3 (18 marks)

(a) (5), 1 mark each (any 5):

  1. QC/trimming of raw reads (FastQC, Trimmomatic).
  2. Alignment/mapping to reference (BWA/Bowtie2) → BAM.
  3. Sort + mark duplicates (samtools/Picard).
  4. Base quality score recalibration (BQSR, GATK).
  5. Variant calling (GATK HaplotypeCaller / bcftools) → VCF.
  6. Variant filtering (hard filters / VQSR) and annotation.

(b) (4): CPM=48024,000,000×106=48024=20\text{CPM} = \frac{480}{24{,}000{,}000}\times10^6 = \frac{480}{24} = 20 CPM. log⁡2(20+1)=log⁡221=4.392\log_2(20+1)=\log_2 21 = 4.392 (2 marks + 2 marks).

(c) (4): VAF =12/40=0.30=12/40 = 0.30. Expected count under het (AF 0.5) =40×0.5=20= 40\times0.5 = 20 reads. Observed 12 < 20; VAF 0.30 is below the 0.5 heterozygous expectation. Binomially, 12/40 is on the low side but not extreme — consistent with a real but low-VAF variant (e.g. allelic imbalance/subclonal) rather than clean germline het; a strict het filter expecting ~0.5 might flag it. (VAF 1, expected 1, comment 2.)

(d) (5):

  • Precision =TPTP+FP=80100=0.800= \frac{TP}{TP+FP}=\frac{80}{100}=0.800.
  • Recall =TPTP+FN=8090=0.889=\frac{TP}{TP+FN}=\frac{80}{90}=0.889.
  • Specificity =TNTN+FP=90110=0.818=\frac{TN}{TN+FP}=\frac{90}{110}=0.818.
  • F1 =2⋅0.8×0.8890.8+0.889=2⋅0.71111.6889=0.842=2\cdot\frac{0.8\times0.889}{0.8+0.889}=2\cdot\frac{0.7111}{1.6889}=0.842. (1 each + F1 2.)
[
  {"claim":"BLOSUM log-odds bit score = 2.00", "code":"import math\nq=0.020; pa=0.10; pb=0.05\ns=math.log(q/(pa*pb),2)\nresult = abs(s-2.0)<1e-9"},
  {"claim":"NJ Q(A,B) = -40", "code":"r=4; dAB=4; RA=0+4+10+10; RB=4+0+10+10\nQ=(r-2)*dAB - RA - RB\nresult = (Q==-40)"},
  {"claim":"CPM=20 and log2(CPM+1)=4.392", "code":"import math\ncpm=480/24_000_000*1e6\nl=math.log(cpm+1,2)\nresult = abs(cpm-20)<1e-9 and abs(l-4.392)<0.001"},
  {"claim":"F1 score = 0.842", "code":"TP=80;FP=20;FN=10\nprec=TP/(TP+FP); rec=TP/(TP+FN)\nf1=2*prec*rec/(prec+rec)\nresult = abs(f1-0.842)<0.001"}
]