How likely is this mutation?

Question

In my previous bioinformatics challenge, I asked you to mutate a DNA sequence. This time, I'd like you to evaluate how likely a mutation, or a series of mutations, is.

The two types of substitutions are transitions and transversions, and due to the chemical structure of DNA bases, transitions are more likely to occur than transversions. A transition is when a base is turned into one of the same size (purine -> purine or pyrimidine -> pyrimidine), and a transversion involves two bases of different sizes (purine <-> pyrimidine). Following Kimura's model, postulated in 1980, we can define a as the probability of a transition occurring for each unit of time, and b as the probability of a transversion occurring.

See diagram below. A and G are purines, C and T are pyrimidines.

Although the exact values of a and b change from organism to organism (and even between different areas in an organism's genome), we can set a=0.25 and b=0.1 for this challenge.

Given two DNA strings of same length as input, I would like you to calculate how likely it would be for string B to be a mutated version of string A. This is code golf, so fewest bytes wins!

Test cases:

Input String A | Input string B | Output probability
tgcctatc       | tgcctata       |  0.1
aggttcctt      | gggttcctt      |  0.25
aactgg         | aaccgg         |  0.25
atgccct        | atcgcct        |  0.01
tatcactaag     | tgtcaatgag     |  0.00625
ctcctgggca     | cttttgcgca     |  0.00625
ctgtgtgct      | cagagagca      |  0.0001
atcgctttca     | ttggctttca     |  0.01
attcac         | taagcg         |  0.000001
attttcattg     | attttttacg     |  0.000625

Specifications:

• The input DNA strings can be as strings of characters, binary or numerical.
• We will assume that all of the mutations are independent events.
• The output can be in plain notation (0.00015) or scientific (\$1.5\times10^{-4}\$).

If you have any questions or require more specifications, send them my way in the comments!

Answer 1

Jelly, 9 bytes

^ị“½¥¢‘İP

Try it online!

Takes input as 2376 instead of acgt respectively. The Footer does this for you

-1 byte thanks to Jonathan Allan!

How it works

^ị“½¥¢‘İP - Main link. Takes the first input on the left and the second on the right
^         - XOR the corresponding elements
  “½¥¢‘   - Yield [10, 4, 1]
 ị        - Index into this array, 1-indexed and modularly
       İ  - Take the inverse of each
        P - Take the product

Why this works

Characters	As integers	XOR	1-index into [10,4,1]	Inverse
aa	2, 2	\$0\$	\$1\$	\$1\$
cc	3, 3	\$0\$	\$1\$	\$1\$
gg	7, 7	\$0\$	\$1\$	\$1\$
tt	6, 6	\$0\$	\$1\$	\$1\$
ac	2, 3	\$1\$	\$10\$	\$0.1\$
ag	2, 7	\$5\$	\$4\$	\$0.25\$
at	2, 6	\$4\$	\$10\$	\$0.1\$
cg	3, 7	\$4\$	\$10\$	\$0.1\$
ct	3, 6	\$5\$	\$4\$	\$0.25\$
gt	7, 6	\$1\$	\$10\$	\$0.1\$

Answer 2

05AB1E, 22 18 16 13 12 10 9 bytes

^Tη4ªsèzP

-4 bytes by porting @Arnauld's JavaScript answer (his first version)
-2 bytes by taking the inputs as lists of codepoints (thanks @cairdCoinheringaahing)
-3 bytes by taking the inputs as [0,1,2,3] for agct respectively
-1 byte thanks to @cairdCoinheringaahing by porting his Jelly answer (and now taking [0,1,2,3] for acgt)
-1 byte thanks to a tip of @JonathanAllan (and now taking [2,3,7,6] for acgt)

Inputs as two lists of integers [2,3,7,6] for acgt respectively.

Try it online or try it online with strings input or verify all test cases.

Explanation:

^          # Bitwise-XOR the (implicit) input-lists at the same positions together
 T         # Push 10
  η        # Pop and push its prefixes: [1,10]
   4ª      # Append a 4: [1,10,4]
     sè    # Index (0-based modular) the XOR-ed values into this list
       z   # Take 1/v for each
        P  # And take the product of the list
           # (after which the result is output implicitly as result)

Here the transformations per pair:

Character-pair	values	XOR-ed	0-based index (modulo-3)	indexed into [1,10,4]	1/v
aa	\$2,2\$	\$0\$	\$0\$	\$1\$	\$1\$
cc	\$3,3\$	\$0\$	\$0\$	\$1\$	\$1\$
gg	\$7,7\$	\$0\$	\$0\$	\$1\$	\$1\$
tt	\$6,6\$	\$0\$	\$0\$	\$1\$	\$1\$
ac	\$2,3\$	\$1\$	\$1\$	\$10\$	\$0.1\$
ag	\$2,7\$	\$5\$	\$2\$	\$4\$	\$0.25\$
at	\$2,6\$	\$4\$	\$1\$	\$10\$	\$0.1\$
cg	\$3,7\$	\$4\$	\$1\$	\$10\$	\$0.1\$
ct	\$3,6\$	\$5\$	\$2\$	\$4\$	\$0.25\$
gt	\$7,6\$	\$1\$	\$1\$	\$10\$	\$0.1\$

Answer 3

JavaScript (ES6),  52  51 bytes

Expects two lists of codepoints as (a)(b).

a=>b=>a.map((c,i)=>p/=!(c^=b[i])|22%c&4||10,p=1)&&p

Try it online!

How?

Given the code points p and q of the nucleobase characters, we use the formula:

(22 MOD (p XOR q)) AND 4

whose result is 4 for a transition pair (a <-> g or c <-> t) or 0 otherwise. Conveniently, 4 is the inverse of the probability of said transition.

 c0  | c1  |  p  |  q  | x = p^q | 22%x | &4
-----+-----+-----+-----+---------+------+----
 'a' | 'a' |  97 |  97 |     0   |  NaN |  0
 'a' | 'c' |  97 |  99 |     2   |    0 |  0
 'a' | 'g' |  97 | 103 |     6   |    4 |  4
 'a' | 't' |  97 | 116 |    21   |    1 |  0
 'c' | 'c' |  99 |  99 |     0   |  NaN |  0
 'c' | 'g' |  99 | 103 |     4   |    2 |  0
 'c' | 't' |  99 | 116 |    23   |   22 |  4
 'g' | 'g' | 103 | 103 |     0   |  NaN |  0
 'g' | 't' | 103 | 116 |    19   |    3 |  0
 't' | 't' | 116 | 116 |     0   |  NaN |  0

Answer 4

R, 77 73 bytes

(input directly as DNA sequence)

function(t,f,`+`=utf8ToInt,d=abs(+t-+f)).1^sum(d>0)*2.5^sum(d%in%c(6,17))

Try it online!

---

R, 53 50 48 bytes

(input as vector of 1 for A, 2 for C, 3 for G, 4 for T)

function(t,f,d=(t-f)^2).1^sum(d>0)*2.5^sum(d==4)

Try it online!

Answer 5

J, 34 bytes

%@10 4*/@({~(>ag`ct)e.~/:~"1)~:#,.

Try it online!

How it works

%@10 4*/@({~(>ag`ct)e.~/:~"1)~:#,.
                                ,. join the characters to pairs
                               #   take the ones
                             ~:    where both are different
                       /:~"1       sort the pairs
                    e.~            is it in …
            (>ag`ct)               'ag' or 'ct'?
          {~                       map 0 and 1 to
%@10 4                             0.1 and 0.25
      */@                          fold with multiplication

Answer 6

Wolfram Language (Mathematica), 27 bytes

1##&@@1[.1,.25,.1][[#-#2]]&

Try it online!

Uses 0123 for acgt.

Answer 7

Python 3, 57 bytes

f=lambda x,y:x==[]or[1,.1,.25,.1][x.pop()^y.pop()]*f(x,y)

Try it online!

Port of my Jelly answer

-2 bytes thanks to ovs!

Answer 8

C (gcc), 79 71 bytes

-8 bytes thanks to ceilingcat

float r;f(a,b,l)int*a,*b;{for(r=1;l--;++a,++b)*a-*b?r/=*a**b-2?10:4:0;}

Try it online!

The function receives two arrays of integers. The bases are mapped as follows

 A  G
-2 -1

 C  T
 1  2

This mapping allows the following decision making

*a-*b?r/=*a**b-2?10:4:0;

In a nutshell: the product of two bases is 2 only if they have the same size (A<->G or C<->T).

Answer 9

Jelly,  8  6 bytes

-2 bytes since we may take input mapped from actg to arbitrary numbers.

^«9‘İP

A dyadic Link accepting two lists of integers, a:4 c:9 t:10 g:7 which yields the probability.

Try it online! Or see the test-suite.

How?

^«9‘İP - Link: a, b               e.g.  [4, 9,10, 7], [4, 4, 4, 4]
^      - bitwise XOR (vectorises)       [0,13,14, 3]
  9    - nine                           9
 «     - minimum                        [0, 9, 9, 3]
   ‘   - increment                      [1,10,10, 4]
    İ  - inverse                        [1, 0.1, 0.25, 0.1]
     P - product                        0.0025

Answer 10

DESMOS, ¹²⁵ 109 bytes

$$f(k,l)=\prod_{n=1}^{\operatorname{length}\left(k\right)}\left\{k-l=0,\left|k-l\right|=2:.25,.1\right\}\left[n\right]$$

\prod_{n=1}^{\operatorname{length}\left(k\right)}\left\{k-l=0,\left|k-l\right|=2:.25,.1\right\}\left[n\right]

try it online:
https://www.desmos.com/calculator/6nxdztewy3

Explanation:
\$\prod_{n=1}^{\operatorname{length}\left(k\right)}{}\left[n\right]\$ takes the product of the list starting with n=1
\$k-l=0\$ if k-l = 0 returns 1
\$\left|k-l\right|=2:.25,.1\$ else if abs(k-l) = 2 returns 0.25 else return 0.1

Answer 11

Vyxal, 10 7 bytes

꘍9v∵›ĖΠ

Try it Online!

-3 by porting Jonathan Allan's Jelly answer

Explained

꘍9v∵›ĖΠ
꘍         # Bitwise XOR of the two input lists
 9v∵      # vectorise(min, ^, 9)
    ›     # ^ + 1 (vectorises)
     Ė    # 1 / ^ (vectorises)
      Π   # product(^)

Answer 12

Perl 5, 81 bytes

sub{$i="@_";$r=1;$_=ord($1^$3),$r/=10-/0/*9-/3|6/*6while$i=~s/(.)(.* )(.)/$2/;$r}

Try it online!

Answer 13

Scala, 52 51 bytes

_.lazyZip(_).map(_.^(_)*7%8/3*7/4*3+1)./:(1.0)(_/_)

Try it online!

Uses 0123 for actg. Like many of the other answers, it xors the values and then hashes them to get the right probabilities.

Here is a table for the XOR, borrowed from Kevin Cruijssen's answer (modified for actg instead of acgt).

Character pair	values	XOR
aa	0,0	0
cc	1,1	0
gg	2,2	0
tt	3,3	0
ac	0,1	1
at	0,2	2
ag	0,3	3
ct	1,2	3
cg	1,3	2
gt	3,2	1

Basically, we get 0 when there isn't a mutation, 1 or 2 for purine -> pyrimidine or vice versa, or 3 for purine -> purine or pyrimidine -> pyrimidine. The goal is to get 0 to 1, 1 and 2 to 10, and 3 to 4, and then later reduce by dividing these inverses of the probabilities.

Original value	*7	%8	/3	*7	/4	*3	+1
0	0	0	0	0	0	0	1
1	7	7	2	14	3	9	10
2	14	6	2	14	3	9	10
3	21	5	1	7	1	3	4

Answer 14

Charcoal, 23 bytes

Ｉ∕¹ΠＥＥθ⁻℅ι℅§ηκ∨¬ι⎇﹪ι⁶χ⁴

Try it online! Link is to verbose version of code. Takes the RNA translations as input, so expects characters from the set acgu. Works with both lower and upper case but not mixed case. Explanation:

      θ                 First input
     Ｅ                  Map over characters
         ι              Current character
        ℅               Ordinal
            η           Second input
           §            Character at index
             κ          Current index
          ℅             Ordinal
       ⁻                Subtract
    Ｅ                   Map over ordinal differences
                ι       Current difference
               ¬        Is zero
              ∨         Logical Or
                   ι    Current difference
                  ﹪     Modulo
                    ⁶   Literal `6`
                 ⎇      If nonzero then
                     χ  Predefined variable `10`
                      ⁴ Else literal `4`
   Π                    Take the product
  ¹                     Literal `1`
 ∕                      Divide
Ｉ                       Cast to string
                        Implicitly print

Answer 15

PHP, 78 bytes

Accepts input as arrays of integers [0,1,2,3], where a->0, c->1, g->2, t->3

function($o,$p){$s=1;foreach($o as$k=>$v){$s/='104'[$p[$k]^$v]?:10;}return$s;}

Try it online!

Explanation

function($o, $p) {        // function accepting two arrays of integers
  $s=1;                   // initial probability
  foreach($o as $k=>$v) { // loop through first input as key, value pairs
    $s /=                 // divide by... and assign result to $s
    '104'                 // string of possible division integers
                          // 1 => no change in probability
                          // 0 => special case
                          // 4 => inverse of 0.25
    [$p[$k] ^ $v]         // take an integer from string offset based on XOR value
                          // of loop's 'value' and corresponding value in second input
    ?:10;                 // evaluate boolean value of integer and return it when "truthy"
                          // else return 10 (inverse of 0.1) when "falsy"
                          // 1) integer is 0 (special case) 
                          // 2) integer is undefined (XOR was 3, string offset doesn't exist)
  }
  return $s;              // return probability
}

Answer 16

Retina, 124 72 bytes

O$`.
$.%`
¶

L`..
A`(.)\1
%O`.
ag|ct
4
..
10
%`$
$*
¶

~`.+
.+¶1/$$.($&_

Try it online! Takes input on separate lines but link includes test suite that converts from comma-separated inputs for convenience. Outputs as a fraction 1/n. Explanation:

O$`.
$.%`
¶

L`..

Transpose the inputs.

A`(.)\1
%O`.
ag|ct
4
..
10

Get the probability of each mutation.

%`$
$*
¶

~`.+
.+¶1/$$.($&_

Take the product and prefix 1/ to the result.

Previous 124-byte solution used decimal output:

O$`.
$.%`
¶

L`..
%O`.
(.)\1
1/1
ag|ct
25/100
[a-z].
1/10
+`\d+(.+)¶(\d+)/1
$.(*$2*)$1
{`/1$

^\B
0
(.)(\.(\d+))?/10
.$1$3/1

Try it online! Takes input on separate lines but link includes test suite that converts from comma-separated inputs for convenience. Explanation:

O$`.
$.%`
¶

L`..

Transpose the inputs.

%O`.
(.)\1
1/1
ag|ct
25/100
[a-z].
1/10

Get the probability of each mutation.

+`\d+(.+)¶(\d+)/1
$.(*$2*)$1

Multiply all of the probabilities together.

{`/1$

^\B
0
(.)(\.(\d+))?/10
.$1$3/1

Convert the fraction to a decimal.

Answer 17

Haskell, 53 bytes

a?b|a==b=1|a/b<0=10|1>0=4
a!b=foldl1(/)$zipWith(?)a b

Try it online!

• input : a-> 1, g-> 2, c-> -1, t -> -2
• ? compares bases to get the inverse of probabilities
• ! zips sequences using ? to compare each pair and folds by dividing