Two's complement in Python (shift left on many bits with rounding)

Question

How we could code the reverse complete of a DNA sequence from its code?

A DNA sequence can contain 4 different characters A, C, G, T; where A is the complement of T and C is the complement of G.

A reverse complement of A DNA sequence is the complement of a sequence but in an inverse way (we compute the complement of each character from right to left).

Example: the complement of (AA) is: TT, the complement of (AC) is GT and so on...

In general, using python we code a sequence by mapping each character to a number going from 0 to 3,

    {A:0, C:1, G:2, T:3}

then the coding of AA is: 0, the coding of AC is:

    AC = 0*4^0+1*4^1 = 4

the coding of GT is:

    GT = 2*4^0+3*4^1 = 14 

How could I transform the code of each sequence to its reverse complement in python without creating a dictionary? For the above example: convert 4 to 14? and 0 to 15 ...

Answer 1

the reverse of a list in python

>>> xs = [1,2,3]
>>> reversed(xs)
<listreverseiterator object at 0x10089c9d0>
>>> list(reversed(xs))
[3, 2, 1]
>>>

def complement(x):
   return ~x & 15 # as 15 == int('1111', 2) 

the 15 is a bitmask. It represents the binary 1111. We then use the binary and operator.

>>> "{0:b}".format(complement(int('1111',2)))
'0'
>>> "{0:b}".format(complement(int('0001',2)))
'1110'
>>> "{0:b}".format(complement(int('1001',2)))
'110'

>>> xs = [int('1111',2), int('1001',2), int('0110',2), int('1011',2)]
>>> map(complement, xs)
[0, 6, 9, 4]
>>> list(reversed(map(complement, xs)))
[4, 9, 6, 0]

Basing your example where

given a sequence of 6 characters: ACCGTT, the complement of A is: T, and the complement of C is G; so the reverse complement of ACCGTT is: AACGGT.

assume that you have c complemnt function complement and a reverse function reverse.

we have reverse(ACCGTT) = TTGCCA and complement(ACCGTT) = TGGCAA . Reversing a list after calling a function on each element is the same as calling a function on each element on a list.

complement(reverse(ACCGTT)) = reverse(complement(ACCGTT))

So the other part of the question is that you want to map

{A:0, C:1, G:2, T:3}
A -> T | 0 -> 3
T -> A | 3 -> 0
C -> G | 1 -> 2
G -> C | 2 -> 1

which in binary would be

a = int('00', 2) # 0
c = int('01', 2) # 1
g = int('10', 2) # 2
t = int('11', 2) # 3

def complement(x):
    return ~x & 3 # this 3 is the same as int('11', 2)

def reverse_complement(list_of_ints):
    return list(reversed(map(complement, list_of_ints)))

Answer 2

Your symbol set is too small for a hash map to actually be efficient. And mixing two's complement into your problem has just caused confusion.

symbols = 'ACGT'
complements = symbols[::-1]   # reverse order
import string
table = string.maketrans(symbols, complements)
sample = 'ACCGTT'
print(sample[::-1].translate(table))
# output: AACGGT

Converting to some bitpacked format would take less space but require a lot more special handling, as you'd need to track sizes separately, perform arbitrarily wide shifts and so on. Python can certainly do it, in particular with int() accepting many bases and creating arbitrary width results, but it's likely a counterproductive detour.

digits = string.digits[:len(symbols)]
length = len(sample)
digitmap = string.maketrans(symbols, digits)
number = int(sample.translate(digitmap), len(digits))

def reversemapnumber(function=id, number=0, radix=0b100, length=0):
    result = 0
    for i in range(length):
        number,digit = divmod(number, radix)
        result = result*radix + function(digit)
    return result
revcomplemented = reversemapnumber(function=lambda x: 3-x,
        number=number, length=length)
# binary form
print('{:0{}b}'.format(revcomplemented, length*2))
# back to text form
print(''.join(symbols[(revcomplemented>>i)&0b11]
    for i in range(2*length-2, -2, -2)))

In that jumble of code I've used division rather than shifts to be somewhat more generic (supporting radix not a power of two), but the printing examples rely on the width exactly. In the end it's just tricky and unclear.