<jupyter_start><jupyter_text>DNA complement and reverse This notebook was used to benchmark the different approach of computingthe complement of a sequence.Takes about 15ms to analyse a 1e7 long sequence.<jupyter_code># First, let us create a sequence.
from biokit.sequence.benchmark import SequenceBenchmark

def create_sequence(expectedLength=1e6):
    s = SequenceBenchmark()
    return s.create_sequence(expectedLength)
sequence = create_sequence(1e7)<jupyter_output>/home/cokelaer/miniconda3/envs/py3/lib/python3.5/site-packages/matplotlib/font_manager.py:280: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
  'Matplotlib is building the font cache using fc-list. '<jupyter_text>BioPython<jupyter_code>from Bio.Seq import Seq
from Bio.Alphabet import IUPAC
from Bio.SeqUtils import GC
seq1 = Seq(sequence, IUPAC.unambiguous_dna)
GC(seq1)
%%timeit -n 3
seq1 = Seq(sequence, IUPAC.unambiguous_dna)
res = seq1.complement()

len(seq1)<jupyter_output><empty_output><jupyter_text>Python with lists<jupyter_code>class DNAList(object):
    bases = ["A", "C", "T", "G"]
    cbases = {"T":"A", "G":"C", "A":"T", "C":"G"}
    def __init__(self, data):
        self.data = data[:]
    def get_complement(self):
        res = "".join((self.cbases[x] for x in self))
        return res 
    def __getitem__(self, index):
        if type(index) == slice:
            return "".join(self.data[index])
        return self.data[index]
seq2 = DNAList(sequence[0:1000000])  # slow
%%timeit -n 3
res = seq2.get_complement()<jupyter_output>395 ms ± 8.97 ms per loop (mean ± std. dev. of 7 runs, 3 loops each)<jupyter_text>Python with strings This is a simplified version of what is done in BioPython and BioKit<jupyter_code>import string
# with strings as ACGT / TGCA, there is no improvment
trans = bytes.maketrans(b'ACGTagct', b'TGCAtgca')
class DNA(object):
    def __init__(self, sequence):
        self.sequence = sequence
    def complement(self):
        return self.sequence.translate(trans)
d = DNA(sequence)

%%timeit -n 3
res = d.complement()
#fh = open("test.fasta", "w")
#fh.write(res)
#fh.close()<jupyter_output>8.71 ms ± 454 µs per loop (mean ± std. dev. of 7 runs, 3 loops each)<jupyter_text>Another way with Pandas?<jupyter_code>import pandas as pd
d = pd.Series(sequence)
# here we just store the sequence so the timeseries is 
# just a container. Instanciation is 50 times longer 
# than a simple class storing the sequence as a string though
trans = bytes.maketrans(b'ACGTacgt', b'TGCAtgca')
%%timeit -n 3
res = "".join(d).translate(trans)<jupyter_output>8.63 ms ± 1.04 ms per loop (mean ± std. dev. of 7 runs, 3 loops each)<jupyter_text>BioKit<jupyter_code>from biokit.sequence import dna

d = dna.DNA(sequence)
%%timeit -n 10
# maketrans
d.get_complement()<jupyter_output>8.34 ms ± 1.28 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)<jupyter_text>cython test - Cythonise does not make it faster. probably the reason being that maketrans  is alredy optimised- THe code below provides a cython version of the code implemented earlier (CyTransSeq) and a cythonised version of a dictionary implementation.- Neither are faster than the maketrans- lookup and lut table do not help either- references: http://nbviewer.ipython.org/urls/gist.github.com/Jorge-C/d51b48d3e18897c46ea2/raw/73d7e11e4b72d6ba90e0021931afa230e63031e9/cython+sequences.ipynb?create=1<jupyter_code>%load_ext cython

%%cython --annotate
# source http://nbviewer.ipython.org/urls/gist.github.com/Jorge-C/d51b48d3e18897c46ea2/raw/73d7e11e4b72d6ba90e0021931afa230e63031e9/cython+sequences.ipynb?create=1
cimport cython
import numpy as np
cimport numpy as cnp

table = bytes.maketrans(b'ACGTacgt',
                        b'TGCAtgca')

cdef class CyTransSeq(object):
    """Simply defining the class as a cython one, uses translation table"""
    cdef public str seq

    def __cinit__(self, seq):
        self.seq = seq

    cpdef rc(self):
        return list(reversed(self.seq.translate(table)))
    
    cpdef complement(self):
        return self.seq.translate(table)

cdef class CySeq(object):
    """Translation using a dict, cythonized"""
    _complement_map = {
        'A': 'T', 'C':'G', 'G':'C', 'T':'A',
        'a': 't', 'c':'g', 'g':'c', 't':'a'}
    cdef public str seq
    
    def __cinit__(self, seq):
        self.seq = seq

    cdef _rc(self):
        result = []
        for base in reversed(self.seq):
            result.append(self._complement_map[base])
        return result
    
    cdef _complement1(self):
        result = []
        for base in self.seq:
            result.append(self._complement_map[base])
        return result
    def complement1(self):
        return self._complement1()
    def complement2(self):
        return self._complement2()
        
    cdef _complement2(self):
        return [self._complement_map[base] for base in self.seq]

    def rc(self):
        return self._rc()


seq = CyTransSeq(sequence)
%%timeit -n 3
res = seq.complement()<jupyter_output>8.13 ms ± 588 µs per loop (mean ± std. dev. of 7 runs, 3 loops each)<jupyter_text>Plot<jupyter_code>import matplotlib as plt
%matplotlib
import pylab
N = [1e5, 1e6, 5e6, 1e7, 5e7, 1e8]
import time
timesBioKit = []
timesC = []
timesBio = []
timesCyt = []
for el in N:
    # BioKit
    factor = int(el / 50.) #100 is the length of the small string
    subseq = "AGCTTTTCATTCTGACTGCAACGGGCAATATGTCAGTGTCTCGTTGCAAA"
    sequence_in = "".join([subseq]*factor)
    seq = dna.DNA(sequence_in)
    t1 = time.time()
    res = seq.get_complement()
    t2 = time.time()
    timesBioKit.append(t2-t1)
    
    #t1 = time.time()
    #res = seq.get_complement_c()
    #t2 = time.time()
    #timesC.append(t2-t1)
    
    
    # biopython
    seqbio = Seq(sequence_in, IUPAC.unambiguous_dna)
    t1 = time.time()
    seqbio.complement()
    t2 = time.time()
    timesBio.append(t2-t1)
    
    # cython
    seqcyt = CyTransSeq(sequence_in)
    t1 = time.time()
    seqcyt.complement()
    t2 = time.time()
    timesCyt.append(t2-t1)
    
    
    print(el)
    
%pylab inline
pylab.clf()
pylab.loglog(N, timesBioKit, 'o-', N, timesBio, 'gx-', 
             N, timesCyt, 
             'ko-')
pylab.legend(['timesBioKit', 'timesBio', 'times Cython'], 
             loc='best')<jupyter_output>Populating the interactive namespace from numpy and matplotlib
{'fraction_non_alphanumeric': 0.11826605936861945, 'fraction_numerical': 0.03156902779959164, 'mean_word_length': 3.542082738944365, 'pattern_counts': {'":': 4, '<': 26, '<?xml version=': 0, '>': 26, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '24295', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 2438, 'n_tokens_neox': 2272, 'n_words': 600, 'path': 'notebooks/benchmark/DNA sequence complement.ipynb', 'repo_name': 'HabibMrad/biokit', 'size': '6367'}