We present a probabilistic logic program to generate an educational puzzle that introduces the basic principles of next generation sequencing, gene nding and the translation of genes to proteins following the central dogma in biology. In the puzzle, a secret "protein word" must be found by assembling DNA from fragments (reads), locating a gene in this sequence and translating the gene to a protein. Sampling using this program generates random instance of the puzzle, but it is possible constrain the di culty and to customize the secret protein word. Because of these constraints and the randomness of the generation process, sampling may fail to generate a satisfactory puzzle. To avoid failure we employ a strategy using adaptive probabilities which change in response to previous steps of generative process, thus minimizing the risk of failure.
Assemble-yourself 3 is a paper-printable \bioinformatics\ puzzle which is
intended to be fun as well as educational { it introduces concepts in genetics and
Next Generation Sequencing (NGS) [
The game is generated by a probabilistic logic program and many aspects of the game can be con gured prior to generation. For instance, the secret protein 3 https://github.com/cth/assemble-yourself word can be changed and the minimal read depth | the number of reads covering any position | can be con gured. In order to ensure satisfaction of constraints in a computationally e cient way, we introduce a heuristic scheme to adapt random switch probabilities based previous states of the program. 2
In the style of logic programming, we begin with the goal. The goal of the game is to nd a particular protein word. For a given game instance, a mutated version of the speci ed protein word | where some letters (amino acids) are randomly substituted | is generated and reported in the background story, c.f., gure 1. Recently an interesting protein with the amino acid sequence ILP was found in the bacteria S. Equencia. It is now to be determined if a homologue exists in the species B. Ionformatica. To determine this a lab ampli cied a relevant part of the DNA of B. Ionformatica using PCR primers anking the gene in S. Equencia which are believed to be highly conserved also in B. Ionformatica, although the sequence of B. Ionformatica is currently not known. The ampli ed DNA was sequenced using Ullamini LoSeq next generation sequencing tech. The quality of the reads are not perfect { read errors resulting in random \mutations\ are expected in one out of twenty bases. As a bioinformatician you are given the task to nd out if B. Ionformatica has a homologue of the protein ILP and determine how its amino acid sequence di ers in B. Ionformatica. However, the high performance moon grid engine supercluster is currently down (as it sometimes is) and you have to do it all by hand. Fortunately, you have printed all the reads. You task is as follows: 1) Perform de-novo assembly of all the reads, 2) Find open reading frames that may contain a gene, 3) Find the amino acid sequence of any such gene to determine if it could be a homologue to ILP, 4) Report your nding and claim eternal fame.
The amino acids are re ected in a corresponding the DNA sequence of nucleotide bases forming a gene, where each triplet of DNA bases (called a codon) correspond to an amino acid. Only part of the DNA sequence is used to encode the protein word. It is anked by random DNA on both sides, but the o set of the encoded protein is randomly determined. A given amino acid may be represented by one of multiple coden triplets, a particular codon may indicate the start of the gene and another is used to indicate the end of it (see gure 2, and e.g., https://en.wikipedia.org/wiki/Genetic code for details). The DNA sequence can be read either left-to-right or right-to-left, but in the latter case the DNA sequence should be complimented, i.e., the bases should be replaced by the one they pair with in the double-helical DNA molecule (A $ T and G $ C). Yet another complication is that codon triplets may start in any position.
Part of the puzzle is to nd possible genes in DNA, translate the nucleotide bases into an amino sequence to nd the secret protein word. Once found the game is completed. However, before doing this, the player is faced with solving the problem of assembling the sequencing from a number of overlapping subsequences of the sequence (reads). This mimics the situation that arise when using NGS technologies, where the assembly process is handled by computationally intensive algorithms. The user however, is tasked with doing this by hand for a small DNA sequence. A scaled down instance of the game is shown in gure 3.
T C A G TTT Phe F TCT Ser S TAT Tyr Y TGT Cys Y T TTC Phe F TCC Ser S TAC Tyr Y TGC Cys Y C T
T CTT Leu L CCT Pro P CAT His H CGT Arg R T h CTC Leu L CCC Pro P CAC His H CGC Arg R C ird CTA Leu L CCA Pro P CAA Gln Q CGA Arg R A ab CTG Leu L CCG Pro P CAG Gln Q CGG Arg R G se ATT Ile I ACT Thr T AAT Asn N AGT Ser S T in ATC Ile I ACC Thr T AAC Asn N AGC Ser S C co d ATA Ile I ACA Thr T AAA Lys K AGA Arg R A on ATG Met M ACG Thr T AAG Lys K AGG Arg R G GTT Val V GCT Ala A GAT Asp D GGT Gly G T GTC Val V GCC Ala A GAC Asp D GGC Gly G C G
GTA Val V GCA Ala A GAA Glu E GGA Gly G A GTG Val V GCG Ala A GAG Glu E GGG Gly G G
We have written the program to generate puzzles in the PRISM language [
The part of the puzzle which involves the protein word corresponding amino acid sequence and its underlying DNA sequence can be generated by a series of probabilistic choices on where to place the protein word, on generation of the
T A T G G A A A T G T T A C C T A A G A A A A T G G A A T C T A C C T T T A C C A C C T T T A C C T G G A A A T G G A A
C G T T G A C C T T
L
R protein word with probability of point mutations as well as corresponding and additional point mutations in the underlying DNA sequence.
We must, however, control the number of mutations that are introduced in
the protein word. If no mutations occur, the word will be the same as the the
homologue protein in the background story, whereas if too many mutations
occur, it may be di cult to determine if a discovered protein is correct. Hence, we
constrain the number of mutations in the protein word to a speci ed range. To
address this constraint we use PRISMs soft_msw construct [
Note that it is not strictly necessary to rely on backtracking to solve this particular problem. An alternative failure-free approach would be to randomly choose the number of mutations from the speci ed range, randomly determine at which positions these mutations should occur and then introduce these mutations in the protein word.
In the generation of reads another global constraint come into play { each position of the DNA sequence must be covered by a speci ed minimum number of reads { the depth. At the same time, the total number of reads is constrained and generally kept to a minimum in order to balance the level fun and di culty in the puzzle.
Reads are generated by a recursive predicate in which the termination case of the predicate speci es the condition that all positions must have the required minimum depth and the recursive case generates a random read, aligns it to the DNA sequence and updates a depth vector, d1 : : : dn, for each position 1 : : : n in the DNA sequence. It is not possible for this predicate to fail, but it may produce too many reads causing an assertion to fail in the calling predicate.
Initially we used soft_msw in our implementation, which would backtrack upon failure, but this approach exhibited poor performance for larger board sizes with a low cap on the total number of reads due to the potentially extensive thrashing behaviour associated with backtracking, where partial solutions leading to failures are repeatedly revisited.
The predicate responsible for placing reads is shown below, placeread(Seq, Part1, Part2,ReadSize,Depths) :length(Seq,L), LMax is L - ReadSize, findall(X,between(0,LMax,X),AllLen), findall(D2,( between(0,LMax,X), length(C1,X), length(C2,ReadSize), append(C1,C2,C3), append(C3,_,Depths),
D2 is min(C2))), MinDepths), inverse_depths_probs(MinDepths,Probs), random_select(AllLen,Probs,L1), L #= L1+L2, length(Part2,L2), append(Part1,Part2,Seq).
The placeread/4 predicate attempts to construct a list of depths from the part of the sequence before the read, the part covering the read, and the part after the read for all combinations of possible read placements. For each of possible read placement, it then nds the minimum depth of any position covered by the read. The call to inverse_depth_probs/2 constructs a probability distribution over the possible read positions from the minimum depths. The probability of placing a new read of length r starting at position i, is given by, P (pos = i) = ( 1 ; n Phn=w1ri wh if Pin=1 di = 0: otherwise: , where wi =
Pjn=1r min dj : : : dj+r min di : : : di+r
The probabilistic selection of read positions is not implemented using the msw construct, but instead using the random_select construct, which given a list of probabilities which serves as an probability distribution over a corresponding list of outcomes, randomly selects one of the outcomes according to the probability distribution. Instead of having a xed probability distribution, we adapt the probabilities in response to previous placements. This ensures that the probability of placing a new read at a given position is inversely proportional to the minimum depth of the positions covered by the read. This approach is not as declarative, since the possible outcomes and the probability distribution can only be determined during program execution. A disadvantage that comes with this lack of declarativity, i.e., that random choices occur outside the de ned PRISM model, is that other forms of inference which rely on a xed model are no longer possible.
The approach is very fast (almost instantaneous), but may occasionally fail. It rarely happens, but when it does, the procedure can just be restarted. It results in fewer reads and a distribution of depths/reads which appears uniformly random, except in the ends near the edges of the board which are characterized by lower depths. 4
Our game has many similarities to Gigsaw [
Using PRISM programs to sample biological sequence data has been done
before, e.g., in [
A heuristic approach to generate programs with low probability of failure,
as embodied by adaptive probability distributions is more e cient than the
soft_msw approach in our case. However, the scheme is speci c to our
application and cannot easily be generalized. As a point of future research, it would
be useful to develop generic, but heuristically informed methods for PLP-based
random sampling which are less prone to thrashing than the pure soft_msw
backtracking approach. Perhaps inspiration can be drawn from the methods of
constraint programming [