Use a GPT as a programming assistant. Your job is to understand the optimization problem, define the state and control, derive the DP recurrence, and interpret the result.

In a simplified story, a lab takes biological material from a crime scene and from several candidate references, i.e., suspects. From each sample it reads the same short region of DNA, called a locus. A sequencing machine then returns a string of letters such as ACGTCGATAC.

Your task is to determine which candidate sequence is most similar to the crime-scene sequence.

It is tempting to picture a gene that starts with a start codon, contains protein-coding DNA, and ends with a stop codon. That is not the best model for this exercise. Real forensic DNA profiling often uses regions such as short tandem repeats, where a short pattern is repeated and the number of repeats can differ between people.¹¹See Wikipedia: DNA profiling for the forensic setting, and Wikipedia: short tandem repeat for STRs. In this tutorial, we ignore the laboratory and statistical details and keep only the computational object: a DNA string for the crime-scene sample and a DNA string for each candidate.

In this exercise, a DNA sequence is a string over the alphabet \(\{ =A=, =C=, =G=, =T= \}\). For example: ACGTCGATAC.

Two related sequences need not have exactly the same nucleotide at every position. Differences can include:

1. a substitution: one nucleotide is replaced by another;
2. an insertion: one sequence contains an extra nucleotide;
3. a deletion: one sequence lacks a nucleotide.

Insertions and deletions are often called indels. Suppose we compare these two sequences:

Sequence 1: A C G T C
Sequence 2: A C T C

If we compare position 1 with position 1, position 2 with position 2, and so on, then after the missing G every later position is shifted. This can create many apparent mismatches. Instead, we can align the sequences by inserting a gap:

Sequence 1: A C G T C
Sequence 2: A C - T C

The gap means that the G in the first sequence is aligned with no nucleotide in the second sequence. After that, the T and C line up again.

Spend at most the first tutorial hour on this.

We now develop the DP and include a few simple exercises for you to check by hand. The next step is to ask a GPT to build everything we need. The example of this section can then serve as a test case for the code the GPT builds.

Let \(X\) and \(Y\) be the two sequences, with lengths \(n\) and \(m\). The state is \(x=(i,j)\). Here \(i\) points at the first nucleotide of \(X\) that has not yet been aligned, and \(j\) at the first nucleotide of \(Y\) that has not yet been aligned, so \(1\leq i\leq n+1\) and \(1\leq j\leq m+1\). The value \(i=n+1\) means that \(X\) is used up. The start state is \((1,1)\) and \((n+1,m+1)\) is the terminal state.

A control \(u\) tells the algorithm what to align next. There are three possible controls:

1. \(u=(1,1)\): align \(X_i\) with \(Y_j\);
2. \(u=(1,0)\): align \(X_i\) with a gap;
3. \(u=(0,1)\): align a gap with \(Y_j\).

Controls that run past the end of a sequence are not feasible. For example, if \(i=n+1\), then \(u=(1,1)\) and \(u=(1,0)\) are not allowed.

Written out, the control set is \[ U(i,j)=\{(a,b)\in\{(1,1),(1,0),(0,1)\}: i+a\leq n+1,\ j+b\leq m+1\}. \]

The transition function is \[ T((i,j),(a,b))=(i+a,j+b). \]

Thus \(T\) updates the state after applying a control.

Every control increases \(i+j\) by at least one, so no state can be visited twice, and any sequence of controls reaches \((n+1,m+1)\) in at most \(n+m\) steps. The position \((i,j)\) therefore plays the role of the time index, and we write \(g\) and \(J^{*}\) without a stage index \(k\).

Use the following simple cost system: \[ \text{match}=0,\qquad \text{mismatch}=1,\qquad \text{gap}=2. \] For two nucleotides \(a\) and \(b\), define the nucleotide cost

A gap has cost \(2\). These numbers are modeling choices. They are not biological constants.

The terminal cost is \[ g_N(n+1,m+1)=0. \]

Let \(J^*(i,j)\) be the minimum remaining cost from state \((i,j)\) to the terminal state \((n+1,m+1)\).

After filling the table, the value \(J^*(1,1)\) gives the minimum cost. It does not yet show the alignment.

To recover an optimal alignment, start at state \((1,1)\). At each state, choose a control \(u\) that achieves the minimum in the recurrence. At state \((i,j)\):

• if the chosen control is \(u=(1,1)\), align \(X_i\) with \(Y_j\);
• if the chosen control is \(u=(1,0)\), align \(X_i\) with a gap;
• if the chosen control is \(u=(0,1)\), align a gap with \(Y_j\).

There can be more than one optimal alignment if there are ties. In that case, any optimal traceback is acceptable, as long as its cost equals \(J^*(1,1)\).

Start with this as soon as you are sure you can run the DP algorithm during the exam.

Have a GPT write a program²²Remember: you are the architect, hence responsible for the final product. that does the following:

1. read the crime-scene sequence; see the data files below.
2. read all candidate sequences;
3. implement global sequence alignment using the DP recurrence above;
4. develop small tests for your alignment code before using the candidate data;
5. compute the optimal alignment cost between the crime-scene sequence and every candidate;
6. rank the candidates by optimal alignment cost, with lower costs ranked first;
7. identify the best matching candidate;
8. recover and display an optimal alignment for the best candidate;
9. make a visual inspection of the recovered alignment; use green to mark gaps and red to mark nucleotide differences.
10. explain where substitutions and gaps occur in that alignment.

Use these files:³³FASTA is a plain-text format for sequence data. Each record starts with a header line beginning with >, followed by one or more lines containing the sequence. See Wikipedia: FASTA format.

data/crime_scene.fasta
data/candidates.fasta

The first file contains one sequence:

crime_scene

The second file contains the candidate reference sequences:

candidate_01
candidate_02
...

You should open them in an editor and check the contents.⁴⁴Checking the primary data is fundamental in data science.

Make a bunch of tests yourself. Testing here is extremely important: we do not want to put the wrong person behind bars because we were too lazy to thoroughly test the code.⁵⁵Please realize: in your job your boss pays you not to be sloppy. Developing testing skills really matters. As a hint, start with very small sequences where you can compute the answer by hand, such as \(AAA\) aligned with \(AAA\), \(AA\), or \(AAT\). Also try empty sequences.⁶⁶Don’t shy away from using super dumb tests. You will be amazed how many simple errors you can catch by this strategy. Add your tests to a FASTA file of your own, for example test_sequences.fasta. Use it to check that your FASTA parser and your alignment code work together. FASTA headers can contain extra text after the identifier, so for your own test file you may store the expected cost in the header:

>test_01_x pair=test_01 expected_cost=0
AAA
>test_01_y pair=test_01 expected_cost=0
AAA

The file tests/test_alignment.py contains a very small test suite for the provided code. From the dna_alignment directory, run it with:

python -m unittest discover -s tests

Read the tests before writing your own. They are intentionally simple.

After completing the main exercise, try changing the gap cost. Compare:

match = 0, mismatch = 1, gap = 2

with:

match = 0, mismatch = 1, gap = 4

Before running the program, predict two things: what happens to the costs, and what happens to the ranking of the candidates.

Here are some interesting questions that we do not address here⁷⁷This is a toy example after all..

1. How much overlap is normal for such DNA sequences from one human?
2. What would be a good threshold for saying that two sequences are sufficiently similar or different?
3. Can we use sequences that do not really overlap—so we are sure the villain is not among the candidates—while still using partial overlap to get hints about whether a family member of a candidate could be the villain?

If you find the above too simple, here are some interesting optional challenges.

• Find a highly similar subsequence inside a longer sequence; the classic DP algorithm is Smith-Waterman;
• Dynamic time warping can align two time series even if one is stretched or compressed in time;⁸⁸This can be used to detect gravitational waves.
• Speech-to-text alignment, where parts of an audio signal or model output are aligned with words, phonemes, or characters;
• The translation-alignment problem discussed in the book is applied to short text. How would you deal with entire books?
• Read Wikipedia: DNA profiling and figure out how real forensic DNA profiling works.

Last ten minutes.

1. Discuss which prompts you used to get the results.
2. How did you design your tests? Which tests did you make?
3. What was simple?
4. What was difficult?
5. What should you look out for?
6. What mistakes did the GPT make? How did you discover them?
7. What mistakes did you make? How did you discover them?
8. How confident are you in your product? Any ideas for increasing your confidence?

Please understand that we (= Nicky and Gijs) also don’t know the answers to these questions. GPTs are new for all of us.

Hopefully you have seen that with a good list of instructions, a GPT can build the code without too many problems. In the following tutorials, we expect you to use prompting more and more to get the right instructions in the right sequence. A second important point is learning how to instruct a GPT to make good test suites.