Question: Code Golf - bisulfite conversion
7
gravatar for Chris Miller
13 days ago by
Chris Miller20k
Washington University in St. Louis, MO
Chris Miller20k wrote:

Been a while since we had a code golf challenge (shortest script wins, but expect to receive upvotes for cleverness)!

Simplifying slightly, bisulfite conversion means that Cs get converted to Ts, but some bases may be protected by methylgroups and thus remain unconverted. (This happens on the opposite strand too, so some complimentary Gs could get converted to As)

So, given an input sequence read (assume a [ACTGN] string of no more than 150 bp), output all reads that might result from partial or complete bisulfite conversion of that sequence.

EDIT: To be clear, C>T happens on one read, G>A happens on the other. Reads will not exist with both conversions! See this screenshot:

Small example: Input:

AACGCGAA

Output:

AATGTGAA
AATGCGAA
AACGTGAA
AACACAAA
AACACGAA
AACGCAAA
AACGCGAA

(the original string could be in there as well, since all bases could be protected!)

Happy Golfing!

ADD COMMENTlink modified 12 days ago by zx87545.3k • written 13 days ago by Chris Miller20k

So that we can test our outputs, what is the output for ACTGN?

ADD REPLYlink written 13 days ago by zx87545.3k

after I saw the other's answers, i'm not sure if we only need to change the 'C' in the sequence or the 'C' AND the 'G'.

ADD REPLYlink modified 13 days ago • written 13 days ago by Pierre Lindenbaum112k

Both need to change since the methylation and read you are looking at could be from either strand.

ADD REPLYlink modified 13 days ago • written 13 days ago by genomax57k
1

If I am understanding the problem correctly, you can't change C->T and G->A for the same read.

5'- ATGCCCATG -3'
3'- TACGGGTAC -5'

The only conversion chemically occurring is C->T, so there are 3 conversions possible for the upper strand, and two for the lower strand. If my interpretation is correct, Pierre Lindenbaum and jrj.healey solutions are incomplete, and zx8754 solution is incorrect.

edit: my interpretation was right, but all answers have been corrected and updated.

ADD REPLYlink modified 12 days ago • written 13 days ago by h.mon20k

This is correct, it's either C->T or G->A, but not both!

ADD REPLYlink written 13 days ago by Chris Miller20k

OK, R code corrected, gives same output as expected.

ADD REPLYlink written 12 days ago by zx87545.3k

Updated this post with more details and example input/output. Sorry for the confusion!

ADD REPLYlink written 13 days ago by Chris Miller20k
5
gravatar for jrj.healey
13 days ago by
jrj.healey7.5k
United Kingdom
jrj.healey7.5k wrote:

A python solution

Edit II: gained 8 bytes (now 108):

NB This doesn’t solve the problem completely correctly as present (was written before OP clarified).

from itertools import product as p
print(["".join(s) for s in p(*[i+{'C':'t','G':'a'}.get(i,'') for i in ""])])

Edit: a little enhancement at the cost of readability:

As a pseudo-oneliner, without counting the input string as bytes:

from itertools import product
print(["".join(s) for s in product(*[i + {'C':'t'}.get(i, '') for i in ""])])

112 bytes, or 116, if the dictionary should read {'C':'t','G':'a'}.


Not 100% sure if this is correct, but I get the same output as Pierre:

from itertools import product
s = [i + {'C':'t'}.get(i, '') for i in "ATGCCCATG"]
print(["".join(sub) for sub in product(*s)])

127 bytes.

['ATGCCCATG', 'ATGCCtATG', 'ATGCtCATG', 'ATGCttATG', 'ATGtCCATG', 'ATGtCtATG', 'ATGttCATG', 'ATGtttATG']

If the G substitution is also needed (OP was unclear):

from itertools import product
s = [i + {'C':'t','G':'a'}.get(i, '') for i in "ATGCCCATG"]
print(["".join(sub) for sub in product(*s)])

I think I get the same as zx8754:

['ATGCCCATG', 'ATGCCCATa', 'ATGCCtATG', 'ATGCCtATa', 'ATGCtCATG', 'ATGCtCATa', 'ATGCttATG', 'ATGCttATa', 'ATGtCCATG', 'ATGtCCATa', 'ATGtCtATG', 'ATGtCtATa', 'ATGttCATG', 'ATGttCATa', 'ATGtttATG', 'ATGtttATa', 'ATaCCCATG', 'ATaCCCATa', 'ATaCCtATG', 'ATaCCtATa', 'ATaCtCATG', 'ATaCtCATa', 'ATaCttATG', 'ATaCttATa', 'ATatCCATG', 'ATatCCATa', 'ATatCtATG', 'ATatCtATa', 'ATattCATG', 'ATattCATa', 'ATatttATG', 'ATatttATa']

@Chris, it would be useful if you had some 100% correct example output?

ADD COMMENTlink modified 12 days ago • written 13 days ago by jrj.healey7.5k

Very interesting ! Could you detailed please :)

ADD REPLYlink written 13 days ago by Bastien Hervé2.2k
2

Each iterable variation on the input string is stored, with a dictionary guiding the substitutions for each character in the string (for i in "ATG...").

According to the docs, the imported itertools module is functionally equivalent to nested for-loops (https://docs.python.org/3/library/itertools.html#itertools.product), allowing you to yield permutations of the input iterables, as 'cartesian products'. Its kind of similar to itertools.permutations, but that would give all the shuffled possibilities (though its following that train of thought that got me to this).

I confess I'm black-boxing this a little as it's heavily based upon https://stackoverflow.com/a/29184387/3691040 and I can't take all the credit. You'll be much better guided reading the docs than listening to me :P

ADD REPLYlink modified 13 days ago • written 13 days ago by jrj.healey7.5k

Yeah the product(*s) is quite murky to me. Anyhow it's really good coding !

ADD REPLYlink written 13 days ago by Bastien Hervé2.2k

The (*s) syntax is just telling the product function that it will receive a variable number of args, as the number of iterables returned will depend on the sequence and its length. The inner workings of product I don't have a very strong feeling for, as I suspect it required better linear algebra than I have haha..

ADD REPLYlink written 13 days ago by jrj.healey7.5k
4
gravatar for Pierre Lindenbaum
13 days ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum112k wrote:

In C :

compilation:

gcc -O3 -o biostar341051 biostar341051.c

usage:

$ ./biostar341051 ATGCCCATG
AACGCGAA
AACGCAAA
AACACGAA
AACACAAA
AACGCGAA
AACGTGAA
AATGCGAA
AATGTGAA

Edit: fixed after I saw your comments about 'G'->A

Edit2 : fixed after Chris's comment.

ADD COMMENTlink modified 13 days ago • written 13 days ago by Pierre Lindenbaum112k
4
gravatar for zx8754
13 days ago by
zx87545.3k
London
zx87545.3k wrote:

Using R, 212 bytes:

Note: I am converting C > t and G > a so it is easier to see the changes in the output

The function with whitespaces:

foo <- function(s){
  x = strsplit(s, "")[[1]]
  c(s, unlist(
    lapply(c("C", "G"), function(i){
      y = which(x == i)    
      sapply(seq(y), function(j) combn(y, j, function(k){
        z = x
        z[ k ] = setNames(c("t", "a"), c("C", "G"))[ z[ k ] ]
        paste(z, collapse = "")}))
    })))
  }

Usage:

foo("AACGCGAA")
# [1] "AACGCGAA" "AAtGCGAA" "AACGtGAA" "AAtGtGAA" "AACaCGAA" "AACGCaAA" "AACaCaAA"
ADD COMMENTlink modified 12 days ago • written 13 days ago by zx87545.3k
2
gravatar for Alex Reynolds
12 days ago by
Alex Reynolds26k
Seattle, WA USA
Alex Reynolds26k wrote:

EDIT I may have a correct (Python-based) answer. Please see edited code below.

In addition to being shorter, jrj.healey's answer is also a bit more than twice as fast as building a powerset of substitution indices and reducing:

from functools import reduce
from itertools import product as p
import timeit

e = 'AACGCGAA'
m = {'C':'t','G':'a'}

def f(k,s,t):
    for i in s:
        t[i] = ord(m[k])
    return t.decode('ascii')

def a1():
    return [i for s in [[f(k,x,bytearray(e,'ascii')) for x in [reduce(lambda z,x: z + [y + [x] for y in z], s, [[]]) for s in [[i for i,v in enumerate(e) if v == k] for k in m.keys()]][i] if x != []] for i,k in enumerate(m.keys())] + [[e]] for i in s]

def a2():
    return ["".join(s) for s in p(*[i+{'C':'t','G':'a'}.get(i,'') for i in e])]

print("a1 ->", timeit.timeit(stmt=a1, number=1000))
print(a1)

print("a2 ->", timeit.timeit(stmt=a2, number=1000))
print(a2)

"""
$ python test341051.py
a1 -> 0.015046401007566601
['AAtGCGAA', 'AACGtGAA', 'AAtGtGAA', 'AACaCGAA', 'AACGCaAA', 'AACaCaAA', 'AACGCGAA']
a2 -> 0.006538485002238303
['AACGCGAA', 'AACGCaAA', 'AACGtGAA', 'AACGtaAA', 'AACaCGAA', 'AACaCaAA', 'AACatGAA', 'AACataAA', 'AAtGCGAA', 'AAtGCaAA', 'AAtGtGAA', 'AAtGtaAA', 'AAtaCGAA', 'AAtaCaAA', 'AAtatGAA', 'AAtataAA']
"""

For code golf purposes, the minified version uses 421 characters 310 bytes. Might be worth comparing with R solution on time or memory, to see how it scales.

from functools import reduce
e='AACGCGAA'
m={'C':'t','G':'a'}
def f(k,s,t):
 for i in s:
  t[i]=ord(m[k])
 return t.decode('ascii')
print([i for s in[[f(k,x,bytearray(e,'ascii'))for x in[reduce(lambda z,x:z+[y+[x]for y in z],s,[[]])for s in[[i for i,v in enumerate(e)if v==k]for k in m.keys()]][i]if x!=[]]for i,k in enumerate(m.keys())]+[[e]]for i in s])

More positive controls may be useful.

ADD COMMENTlink modified 12 days ago • written 12 days ago by Alex Reynolds26k
2

Forget efficiency :)

shortest script wins

Use http://bytesizematters.com/ and tick "Ignore whitespace".

ADD REPLYlink modified 12 days ago • written 12 days ago by zx87545.3k

haha okay 310 bytes, then!

ADD REPLYlink written 12 days ago by Alex Reynolds26k

mine still needs some work now Chris has clarified, it actually substitutes too many things :/

ADD REPLYlink written 12 days ago by jrj.healey7.5k

Hmm, mine will work on each key of the substitution map m on its own, one at a time — so either C or G is substituted, but not both — so I think mine should give correct output. It would be useful to have some positive examples though, just to test this out.

ADD REPLYlink modified 12 days ago • written 12 days ago by Alex Reynolds26k

When I edited the post, I added a small example!

ADD REPLYlink written 12 days ago by Chris Miller20k

Okay, I made some edits and re-ran things with output. I think mine matches your example. Some other input or edge cases may be useful to test.

ADD REPLYlink written 12 days ago by Alex Reynolds26k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1895 users visited in the last hour