17
\$\begingroup\$

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.

Graphical presentation of transition vs transversion with probabilities a and b

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!

\$\endgroup\$
16
  • 2
    \$\begingroup\$ For the purpose of questions to self-contain the necessary informations, a brief explanation of what is a transition and a transversion would be welcome \$\endgroup\$
    – Kaddath
    Feb 10 at 12:33
  • 17
    \$\begingroup\$ It's a shame that P isn't a base, otherwise we could calculate the chance of PPCG mutating to CGCC... \$\endgroup\$
    – Neil
    Feb 10 at 12:40
  • 1
    \$\begingroup\$ tatcactaag | tgtcaatgag | 0.0625 should be corrected to tatcactaag | tgtcaatgag | 0.00625 and the sixth case should also be 0.00625 and the nineth should be 0.000001 I think since it's 6 bs and no as. \$\endgroup\$
    – Kjetil S.
    Feb 10 at 13:52
  • 1
    \$\begingroup\$ Can I output as a fraction e.g. 1/40? \$\endgroup\$
    – Neil
    Feb 10 at 18:54
  • 1
    \$\begingroup\$ I‘ve corrected the fifth and sixth test cases to be more inline with what all the answers are getting, to avoid more confusion \$\endgroup\$ Feb 11 at 15:24

17 Answers 17

11
\$\begingroup\$

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\$
\$\endgroup\$
5
  • \$\begingroup\$ Save a byte by using different digits, e.g. TIO \$\endgroup\$ Feb 11 at 20:20
  • \$\begingroup\$ @JonathanAllan Nice catch, I thought I‘d checked for combinations that didn’t require duplicating the 10 \$\endgroup\$ Feb 11 at 20:41
  • \$\begingroup\$ So, it appears that there is a tie! Should I accept your answer, as @cairdcoinheringaahing acknowledges that most of their solution is ported from yours? \$\endgroup\$
    – Whitehot
    Feb 17 at 15:59
  • 1
    \$\begingroup\$ @Whitehot If you want to accept an answer, we typically use “first to get that that score” as a tiebreaker metric (which I believe was mine). However, we generally discourage accepting answers at all, because code-golf is best thought of as a competition within each language, so I’m competing against other Jelly answers, rather than every other answer, in order to avoid discouraging people who like golfing in longer languages, such as C# or Java. That said, it’s entirely up to you if you want to accept an answer or not \$\endgroup\$ Feb 17 at 16:02
  • \$\begingroup\$ thanks for the info, I wasn't aware! I'll leave it open then :) \$\endgroup\$
    – Whitehot
    Feb 18 at 14:04
9
\$\begingroup\$

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\$
\$\endgroup\$
10
  • \$\begingroup\$ -1 byte taking input as a list of code points \$\endgroup\$ Feb 10 at 14:32
  • 1
    \$\begingroup\$ I assumed that it was ok because of Arnauld's answer, I'm not sure if there are any meta defaults \$\endgroup\$ Feb 10 at 14:40
  • 2
    \$\begingroup\$ 12 bytes (porting my Jelly answer) \$\endgroup\$ Feb 10 at 17:16
  • 1
    \$\begingroup\$ I've got to stop helping you outgolfing me :) \$\endgroup\$ Feb 10 at 17:50
  • 1
    \$\begingroup\$ Would it save to use different digits (like 2, 3, 6, and 7)? \$\endgroup\$ Feb 11 at 20:23
7
\$\begingroup\$

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
\$\endgroup\$
3
  • \$\begingroup\$ Unsure as to why test("attcac", "taagcg") gives you such an output of 0.0000010000000000000002. In theory, the lowest decimal point that can reach is 10^-6, seeing as the sequence is only 6 bases long and only b mutations. Could there be an issue in comparing two strings with no characters in common? \$\endgroup\$
    – Whitehot
    Feb 10 at 13:31
  • 4
    \$\begingroup\$ @Whitehot That‘ll be a floating point issue \$\endgroup\$ Feb 10 at 13:34
  • 4
    \$\begingroup\$ @Whitehot Caird is right. Negative powers of ten cannot be represented as exact values in IEEE-754 format. For instance, \$10^{-1}\$ is stored as \$1.600000023841858\times 2^{-4}\$ (the exponent is exact, but the mantissa is approximated). For this challenge, we divide by 10 several times, so there is an accumulation of small consecutive errors which may or may not be visible in the displayed number depending on the final rounding. (We would always get exact results if we had b=0.125, which is a power of 2.) \$\endgroup\$
    – Arnauld
    Feb 10 at 15:39
6
\$\begingroup\$

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!

\$\endgroup\$
5
\$\begingroup\$

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
\$\endgroup\$
3
  • 1
    \$\begingroup\$ 26 bytes using integer encoding: [:*/%@4 10{~1<~:|@-/"1@#,. Try it online! \$\endgroup\$
    – Jonah
    Feb 11 at 6:05
  • \$\begingroup\$ 28 bytes for current encoding: [:*/%@10 4{~0=5|1#.3 u:~:#,. Try it online! \$\endgroup\$
    – Jonah
    Feb 11 at 14:39
  • \$\begingroup\$ 22 bytes with a different integer encoding: [:*/%@10 4{~0=1#.~:#,. Try it online! \$\endgroup\$
    – Jonah
    Feb 11 at 14:47
5
\$\begingroup\$

Wolfram Language (Mathematica), 27 bytes

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

Try it online!

Uses 0123 for acgt.

\$\endgroup\$
5
\$\begingroup\$

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!

\$\endgroup\$
2
  • 1
    \$\begingroup\$ x==[]or ... saves two bytes. Try it online! (This returns True for the empty input, which is allowed by default) \$\endgroup\$
    – ovs
    Feb 11 at 9:25
  • \$\begingroup\$ @ovs Very nice, thanks! \$\endgroup\$ Feb 11 at 13:56
3
\$\begingroup\$

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).

\$\endgroup\$
1
  • 1
    \$\begingroup\$ @ceilingcat thank you! :D \$\endgroup\$ Feb 12 at 11:07
3
\$\begingroup\$

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
\$\endgroup\$
2
\$\begingroup\$

DESMOS, 125 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

\$\endgroup\$
2
\$\begingroup\$

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(^)
\$\endgroup\$
1
\$\begingroup\$

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!

\$\endgroup\$
1
\$\begingroup\$

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
\$\endgroup\$
1
\$\begingroup\$

Charcoal, 23 bytes

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

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
     E                  Map over characters
         ι              Current character
        ℅               Ordinal
            η           Second input
           §            Character at index
             κ          Current index
          ℅             Ordinal
       ⁻                Subtract
    E                   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
I                       Cast to string
                        Implicitly print
\$\endgroup\$
1
\$\begingroup\$

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
}
\$\endgroup\$
1
\$\begingroup\$

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.

\$\endgroup\$
1
\$\begingroup\$

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
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.