### Navigation

- index
- modules |
- Pydna 0.1 documentation »
- Module code »

# Source code for pydna.amplify

    #!/usr/bin/env python
    # -*- coding: utf-8 -*-
    # doctest: +NORMALIZE_WHITESPACE
    '''This module provides functions for PCR.

    Primers with 5' tails as well as inverse PCR on
    circular templates are handled correctly.

    '''

    import datetime
    import itertools
    import math
    import re
    import string
    import sys
    import textwrap
    import collections
    import warnings
    import copy

    from StringIO                       import StringIO
    from Bio.SeqUtils.CheckSum          import seguid
    from math                           import log10
    from math                           import log
    from Bio                            import SeqIO
    from Bio.Seq                        import Seq
    from Bio.Seq                        import reverse_complement
    from Bio.Alphabet.IUPAC             import unambiguous_dna
    from Bio.Alphabet.IUPAC             import ambiguous_dna
    from Bio.SeqRecord                  import SeqRecord
    from Bio.SeqUtils                   import GC
    from Bio.SeqUtils.MeltingTemp       import Tm_staluc
    from Bio.SeqFeature                 import SeqFeature
    from Bio.SeqFeature                 import FeatureLocation
    from Bio.SeqFeature                 import ExactPosition
    from Bio.SeqRecord                  import SeqRecord
    from pydna.dsdna                    import Dseq, Dseqrecord, rc, parse

    def _annealing_positions(primer, template, limit=15):
        '''Finds the annealing position(s) for primer on template where the
        primer anneals with at least limit nucleotides in the 3' part.
        start is a position (integer)
        footprint1 and tail1 are strings.

        ::

            <------------- start ---->
         5'-...gctactacacacgtactgactgcctccaagatagagtcagtaaccacactcgat...3'
               ||||||||||||||||||||||||||||||||||||||||||||||||
                                      3'-gttctatctcagtcattggtgtATAGTG-5'

                                                            <tail>
                                         <---footprint----->
                                         <--------- primer ------>

        Parameters
        ----------
        primer : string
            The primer sequence 5'-3'

        template : string
            The template sequence 5'-3'

        limit : int = 15, optional
            footprint needs to be at least of length limit.

        Returns
        -------
        describe : list of tuples (int, string, string)
            [ (start1, footprint1, tail1), (start2, footprint2, tail2),..., ]

        '''

        if len(primer)<limit:
            return []
        prc = rc(primer)
        head = prc[:limit]
        positions = [m.start() for m in re.finditer('(?={0})'.format(head), template, re.I)]
        if positions:
            tail = prc[limit:]
            length = len(tail)
            revtail = tail[::-1]
            results = []
            for match_start in positions:
                tm = template[match_start+limit:match_start+limit+length]
                footprint = rc(template[match_start:match_start+limit]+"".join([b for a,b in itertools.takewhile(lambda x: x[0].lower()==x[1].lower(),zip(tail, tm))]))
                results.append((match_start, footprint, primer[: len(primer) - len(footprint) ]))
            return results
        return []

    class Amplicon(object):
        '''Holds information about a PCR reaction involving two
        primers and one template. This class is used by the
        Anneal class and is not meant to be instantiated directly.

        Parameters
        ----------
        forward_primer : SeqRecord(Biopython)
            SeqRecord object holding the forward (sense) primer

        reverse_primer : SeqRecord(Biopython)
            SeqRecord object holding the reverse (antisense) primer

        template : Dseqrecord
            Dseqrecord object holding the template (circular or linear)

        saltc : float, optional
            saltc = monovalent cations (mM) (Na,K..)
            default value is 50mM
            This is used for Tm calculations.

        fc : float, optional
            primer concentration (nM)
            default set to 1000nM = 1µM
            This is used for Tm calculations.

        rc : float, optional
            primer concentration (nM)
            default set to 1000nM = 1µM
            This is used for Tm calculations.

        '''

        def __init__(self,
                      forward_primer,
                      reverse_primer,
                      template,
                      saltc=50,
                      fc=1000,
                      rc=1000):

            self.forward_primer             = forward_primer
            self.reverse_primer             = reverse_primer
            self.fc                         = fc
            self.rc                         = rc
            self.saltc                      = saltc
            self.template                   = template
            self.product                    = None

        def __repr__(self):
            '''returns a short string representation of the object
           '''
            return "Amplicon({})".format(self.__len__())

        def __len__(self):
            '''returns the lenght of the PCR product
           '''
            length  = len(self.forward_primer.primer)
            length += len(self.reverse_primer.primer)
            if self.reverse_primer.pos>self.forward_primer.pos:
                length += abs(self.reverse_primer.pos-self.forward_primer.pos)
            else:
                length += len(self.template)-abs(self.reverse_primer.pos-self.forward_primer.pos)
            return length

        def pcr_product(self):
            '''Returns the PCR product as a Dseqrecord object.
           Primers are marked as SeqFeatures(Biopython)
           associated with the sequence.

           '''
            if self.product:
                return self.product

            begin = self.forward_primer.pos
            end   = self.reverse_primer.pos

            if self.template.circular:
                tmpl=Dseqrecord(Dseq(self.template.seq.watson,
                                     self.template.seq.crick,
                                     ovhg=0,
                                     linear = True))
                tmpl.features = copy.deepcopy(self.template.features)
                tmpl=tmpl+tmpl
                for feature in tmpl.features:
                    if feature.location.start > feature.location.end:
                        featureEnd = (feature.location.end+len(self.template))%(2*len(self.template))
                        feature.location._end = ExactPosition(featureEnd)
                if begin>end:
                    end=end+len(self.template)
            else:
                tmpl=self.template

            self.product = ( Dseqrecord(self.forward_primer.primer) +
                             tmpl[begin:end] +
                             Dseqrecord(self.reverse_primer.primer.reverse_complement()) )

            features = tmpl[begin-len(self.forward_primer.footprint):end+len(self.reverse_primer.footprint)].features

            self.product.features = [f._shift(len(self.forward_primer.tail)) for f in features]

            # description = Genbank LOCUS max 16 chars

            length = len(self.product)

            self.product.name = "{0}bp_PCR_prod".format(length)[:16]
            self.product.id = "{0}bp {1}".format( str(length)[:14],self.product.seguid() )
            self.product.description="Primers {0} {1}".format( self.forward_primer.primer.name,
                                                      self.reverse_primer.primer.name)

            self.product.features.append(SeqFeature(FeatureLocation(0,len(self.forward_primer.primer)),
                                           type ="primer_bind",strand = 1,
                                           qualifiers = {"note":[self.forward_primer.primer.name]}))
            self.product.features.append(SeqFeature(FeatureLocation(length - len(self.reverse_primer.primer),length),
                                           type ="primer_bind",
                                           strand = -1,
                                           qualifiers = {"note":[self.reverse_primer.primer.name]}))

            return self.product

        def flankup(self,flankuplength=50):
            '''Returns a Dseqrecord object containing flankuplength bases upstream of the forward primer footprint,
           Truncated if the template is not long enough.

           ::

            <--- flankup --->

                      5TAATAAactactgactatct3
                             ||||||||||||||
            acgcattcagctactgtactactgactatctatcg

           '''
            return self.template.seq[self.forward_primer.pos-flankuplength-len(self.forward_primer.footprint):self.forward_primer.pos-len(self.forward_primer.footprint)]

        def flankdn(self,flankdnlength=50):
            '''Returns a Dseqrecord object containing flankdnlength bases downstream of the reverse primer footprint.
           Truncated if the template is not long enough.

           ::

                                           <---- flankdn ------>

                            3actactgactatctTAATAA5
                             ||||||||||||||
            acgcattcagctactgtactactgactatctatcgtacatgtactatcgtat

           '''
            return self.template.seq[self.reverse_primer.pos+len(self.reverse_primer.footprint):self.reverse_primer.pos+flankdnlength+len(self.reverse_primer.footprint)]

        def _tm(self):
            '''Tm calculations according to SantaLucia 1998
           Four variables are set.

           '''
            self.tmf = Tm_staluc(str(self.forward_primer.footprint),dnac=50, saltc=self.saltc)
            self.tmr = Tm_staluc(str(self.reverse_primer.footprint),dnac=50, saltc=self.saltc)
            '''
           Ta calculation for enzymes with dsDNA binding domains (dbd)
           like Pfu-Sso7d, Phusion or Phire
           https://www.finnzymes.fi/tm_determination.html

           '''
            self.tmf_dbd = tmbresluc(str(self.forward_primer.footprint),primerc=self.fc)
            self.tmr_dbd = tmbresluc(str(self.reverse_primer.footprint),primerc=self.rc)

        def detailed_figure(self):
            '''::

            5gctactacacacgtactgactg3
             |||||||||||||||||||||| tm 52.6 (dbd) 58.3
            5gctactacacacgtactgactg...caagatagagtcagtaaccaca3
            3cgatgatgtgtgcatgactgac...gttctatctcagtcattggtgt5
                                      |||||||||||||||||||||| tm 49.1 (dbd) 57.7
                                     3gttctatctcagtcattggtgt5

           Returns
           -------
           figure:string
                A string containing a text representation of the primers
                annealing on the template.

           Notes
           -----
           tm is the melting temperature (Tm) calculated according to
           SantaLucia 1998 [1] for each primer.

           (dbd) is Tm calculation for enzymes with dsDNA binding domains
           like Pfu-Sso7d [2]. See [3] for more information.

           References
           ----------

           .. [1] J. SantaLucia, “A Unified View of Polymer, Dumbbell, and Oligonucleotide DNA Nearest-neighbor Thermodynamics,” Proceedings of the National Academy of Sciences 95, no. 4 (1998): 1460.
           .. [2] M. Nørholm, “A Mutant Pfu DNA Polymerase Designed for Advanced Uracil-excision DNA Engineering,” BMC Biotechnology 10, no. 1 (2010): 21, doi:10.1186/1472-6750-10-21.
           .. [3] http://www.thermoscientificbio.com/webtools/tmc/

           '''

            self._tm()
            f =   '''
                {sp1}5{faz}...{raz}3
                 {sp3}{rap} tm {tmr} (dbd) {tmr_dbd}
                {sp3}3{rp}5
                5{fp}3
                 {fap:>{fplength}} tm {tmf} (dbd) {tmf_dbd}
                {sp2}3{fzc}...{rzc}5
                '''.format( fp       = self.forward_primer.primer.seq,
                            fap      = "|"*len(self.forward_primer.footprint),
                            fplength = len(self.forward_primer.primer.seq),
                            tmf      = round(self.tmf,1),
                            tmr      = round(self.tmr,1),
                            tmf_dbd  = round(self.tmf_dbd,1),
                            tmr_dbd  = round(self.tmr_dbd,1),
                            rp       = self.reverse_primer.primer.seq[::-1],
                            rap      = "|"*len(self.reverse_primer.footprint),
                            rplength = len(self.reverse_primer.primer.seq),
                            faz      = self.forward_primer.footprint,
                            raz      = self.reverse_primer.footprint.reverse_complement(),
                            fzc      = self.forward_primer.footprint.complement(),
                            rzc      = self.reverse_primer.footprint[::-1],
                            sp1       = " "*(len(self.forward_primer.primer.seq)-len(self.forward_primer.footprint)),
                            sp2       = " "*(len(self.forward_primer.primer.seq)-len(self.forward_primer.footprint)),
                            sp3      = " "*(3+len(self.forward_primer.primer.seq))
                           )
            return textwrap.dedent(f).strip("\n")

        def pcr_program(self):
            '''Returns a string containing a text representation of two proposed
           PCR programs. The first program is adapted for Taq poymerase, while
           the second is adapted for Pfu-Sso7d.

           ::

            Taq (rate 30 nt/s)
            Three-step|         30 cycles     |      |SantaLucia 1998
            94.0°C    |94.0°C                 |      |SaltC 50mM
            __________|_____          72.0°C  |72.0°C|
            04min00s  |30s  \         ________|______|
                      |      \ 46.0°C/ 0min 1s|10min |
                      |       \_____/         |      |
                      |         30s           |      |4-8°C

            Pfu-Sso7d (rate 15s/kb)
            Three-step|          30 cycles   |      |Breslauer1986,SantaLucia1998
            98.0°C    |98.0°C                |      |SaltC 50mM
            __________|_____          72.0°C |72.0°C|Primer1C   1µM
            00min30s  |10s  \ 61.0°C ________|______|Primer2C   1µM
                      |      \______/ 0min 0s|10min |
                      |        10s           |      |4-8°C

           '''
            self._tm()
            if not self.product:
                self.pcr_product()
            # Ta calculation according to
            # Rychlik, Spencer, and Rhoads, 1990, Optimization of the anneal
            # ing temperature for DNA amplification in vitro
            # http://www.ncbi.nlm.nih.gov/pubmed/2003928
            GC_prod=GC(str(self.product.seq))
            tml = min(self.tmf,self.tmr)
            #print GC(str(self.product.seq)), self.saltc/1000.0, len(self.product)
            tmp = 81.5 + 0.41*GC(str(self.product.seq)) + 16.6*log10(self.saltc/1000.0) - 675/len(self.product)
            ta = 0.3*tml+0.7*tmp-14.9
            # Fermentas recombinant taq
            taq_extension_rate = 30  # seconds/kB PCR product length
            extension_time_taq = taq_extension_rate * len(self.product) / 1000 # seconds
            f  = textwrap.dedent(u'''
                                    Taq (rate {rate} nt/s)
                                    Three-step|         30 cycles     |      |SantaLucia 1998
                                    94.0°C    |94.0°C                 |      |SaltC {saltc:2}mM
                                    __________|_____          72.0°C  |72.0°C|
                                    04min00s  |30s  \         ________|______|
                                              |      \ {ta}°C/{0:2}min{1:2}s|10min |
                                              |       \_____/         |      |
                                              |         30s           |      |4-8°C
                                 '''.format(rate    = taq_extension_rate,
                                            ta      = math.ceil(ta),
                                            saltc   = self.saltc,
                                            *divmod(extension_time_taq,60)))

            PfuSso7d_extension_rate = 15 #seconds/kB PCR product
            extension_time_PfuSso7d = PfuSso7d_extension_rate * len(self.product) / 1000  # seconds

            # Ta calculation for enzymes with dsDNA binding domains like Pfu-Sso7d
            # https://www.finnzymes.fi/tm_determination.html

            length_of_f = len(self.forward_primer.footprint)
            length_of_r = len(self.reverse_primer.footprint)

            if (length_of_f>20 and length_of_r>20 and self.tmf_dbd>=69.0 and self.tmr_dbd>=69.0) or (self.tmf_dbd>=72.0 and self.tmr_dbd>=72.0):
                f+=textwrap.dedent( u'''
                                        Pfu-Sso7d (rate {rate}s/kb)
                                        Two-step|    30 cycles |      |Breslauer1986,SantaLucia1998
                                        98.0°C  |98.0C         |      |SaltC {saltc:2}mM
                                        _____ __|_____         |      |Primer1C {fc:3}µM
                                        00min30s|10s  \  72.0°C|72.0°C|Primer2C {rc:3}µM
                                                |      \_______|______|
                                                |      {0:2}min{1:2}s|10min |4-8°C
                                     '''.format(rate = PfuSso7d_extension_rate,
                                                fc = self.fc,
                                                rc = self.rc,
                                                saltc = self.saltc,
                                                *divmod(extension_time_PfuSso7d,60)))
            else:

                if (length_of_f>20 and length_of_r>20):
                    ta = min(self.tmf_dbd,self.tmr_dbd)+3
                else:
                    ta = min(self.tmf_dbd,self.tmr_dbd)

                f+=textwrap.dedent( u'''
                                        Pfu-Sso7d (rate {rate}s/kb)
                                        Three-step|          30 cycles   |      |Breslauer1986,SantaLucia1998
                                        98.0°C    |98.0°C                |      |SaltC {saltc:2}mM
                                        __________|_____          72.0°C |72.0°C|Primer1C {fc:3}µM
                                        00min30s  |10s  \ {ta}°C ________|______|Primer2C {rc:3}µM
                                                  |      \______/{0:2}min{1:2}s|10min |
                                                  |        10s           |      |4-8°C
                                     '''.format(rate = PfuSso7d_extension_rate,
                                                ta   = math.ceil(ta),
                                                fc   = self.fc/1000,
                                                rc   = self.rc/1000,
                                                saltc= self.saltc,
                                                *divmod(extension_time_PfuSso7d,60)))
            return f

    class Anneal(object):
        '''

        Parameters
        ----------
        primers : iterable containing SeqRecord objects
            Primer sequences 5'-3'.

        template : Dseqrecord object
            The template sequence 5'-3'.

        homology_limit : int, optional
            limit length of the annealing part of the primers.

        max_product_size: int, optional
            PCR products longer than max_product_size are not reported

        Attributes
        ----------
        amplicons: list
            A list of Amplicon objects, one for each primer pair that may form a PCR product.

        Examples
        --------
        >>> import pydna
        >>> template = pydna.Dseqrecord("tacactcaccgtctatcattatctactatcgactgtatcatctgatagcac")
        >>> from Bio.SeqRecord import SeqRecord
        >>> p1=pydna.read(">p1\\ntacactcaccgtctatcattatc", ds = False)
        >>> p2 = pydna.read(">p2\\ncgactgtatcatctgatagcac", ds = False).reverse_complement()
        >>> pydna.Anneal((p1,p2), template)
        Anneal(amplicons = 1)
        >>> pydna.Anneal((p1,p2), template).amplicons
        [Amplicon(51)]
        >>> amplicon_list = pydna.Anneal((p1,p2), template).amplicons
        >>> amplicon = amplicon_list.pop()
        >>> amplicon
        Amplicon(51)
        >>> print amplicon.detailed_figure()
        5tacactcaccgtctatcattatc...cgactgtatcatctgatagcac3
                                   |||||||||||||||||||||| tm 50.6 (dbd) 60.5
                                  3gctgacatagtagactatcgtg5
        5tacactcaccgtctatcattatc3
         ||||||||||||||||||||||| tm 49.4 (dbd) 58.8
        3atgtgagtggcagatagtaatag...gctgacatagtagactatcgtg5
        >>> prod = amplicon.pcr_product()
        >>> prod.annotations['date'] = '02-FEB-2013'
        >>> print prod
        Dseqrecord
        circular: False
        size: 51
        ID: 51bp U96+TO06Y6pFs74SQx8M1IVTBiY
        Name: 51bp_PCR_prod
        Description: Primers p1 <unknown name>
        Number of features: 2
        /date=02-FEB-2013
        Dseq(-51)
        taca..gcac
        atgt..cgtg
        >>>

        '''
        def __init__(self,
                      primers,
                      template,
                      homology_limit=13,
                      max_product_size=15000):

            self.homology_limit = homology_limit

            primers=[p for p in primers if p.seq]

            self.template = Dseqrecord(template)

            primer = collections.namedtuple("primer","primer pos footprint tail")

            self.fwd_primers = []
            self.rev_primers = []
            self.amplicons   = []

            twl = len(self.template.seq.watson)
            tcl = len(self.template.seq.crick)

            if self.template.linear:
                tw = self.template.seq.watson
                tc = self.template.seq.crick
            else:
                tw = self.template.seq.watson+template.seq.watson
                tc = self.template.seq.crick+template.seq.crick

            for p in primers:
                self.fwd_primers.extend((primer(p,
                                                tcl-pos - min(self.template.seq.ovhg, 0),
                                                Seq(fp), Seq(tl))
                                        for pos, fp, tl in _annealing_positions(
                                                            str(p.seq),
                                                            tc,
                                                            homology_limit) if pos<tcl))
                self.rev_primers.extend((primer(p,
                                                pos + max(0, self.template.seq.ovhg),
                                                Seq(fp), Seq(tl))
                                         for pos, fp, tl in _annealing_positions(
                                                                         str(p.seq),
                                                                         tw,
                                                                         homology_limit) if pos<twl))
            for fp in self.fwd_primers:
                for rp in self.rev_primers:
                    tail_lengths = len(fp.tail)+len(rp.tail)
                    if (rp.pos >= fp.pos ) and (rp.pos-fp.pos)+len(fp.primer)+len(rp.primer) <= max_product_size:
                        self.amplicons.append(Amplicon(fp, rp, self.template))
                    elif self.template.circular and fp.pos > rp.pos and len(template)-fp.pos-rp.pos < max_product_size:
                        self.amplicons.append(Amplicon(fp, rp, self.template))

        def report(self):
            '''This method is an alias of str(Annealobj).
            Returns a short string representation.
           '''
            return self.__str__()

        def __repr__(self):
            ''' return a short string represenation '''
            return "Anneal(amplicons = {})".format(len(self.amplicons))

        def __str__(self):
            '''return a short report describing if or where primer
           anneal on the template.
           '''
            mystring = "Template {name} {size} nt {top}:\n".format(name=self.template.name,
                                                                        size=len(self.template),
                                                                        top={True:"circular",
                                                                             False:"linear"}[self.template.circular]
                                                                        )
            if self.fwd_primers:
                for p in self.fwd_primers:
                    mystring += "Primer {name} anneals at position {pos}\n".format(name=p.primer.name, pos=p.pos)
            else:
                mystring += "No forward primers anneal...\n"
            if self.rev_primers:
                for p in self.rev_primers:
                    mystring += "Primer {name} anneals reverse at position {pos}\n".format(name=p.primer.name, pos=p.pos)
            else:
                 mystring += "No reverse primers anneal...\n"
            return mystring.strip()

        def products(self):
            ''' returns all pcr product as a list of Dseqrecords '''
            return [a.product for a in self.amplicons]

        def featured_template(self):
            '''returns the template Dseqrecord object with the footprint of all
           annealing primers marked as SeqFeatures.

           '''
            template = self.template
            for primer in self.fwd_primers:
                if primer.pos-len(primer.footprint)>0:
                    start = primer.pos-len(primer.footprint)
                    end = primer.pos
                    template.features.append(SeqFeature(FeatureLocation(start,end+1),
                                                        type ="primer_bind",
                                                        strand = 1,
                                                        qualifiers = {"note":[fp.name],
                                                                      "ApEinfo_fwdcolor":["green"],
                                                                      "ApEinfo_revcolor":["red"]}))
                else:
                    start = len(template)-len(primer.footprint)+primer.pos
                    end = start+len(primer.footprint)-len(template)

                    suba = SeqFeature( FeatureLocation(start,len(template)),
                                       type ="primer_bind",
                                       strand = 1,
                                       qualifiers = {"note":[primer.primer.name]} )

                    subb = SeqFeature( FeatureLocation(0,end),
                                       type ="primer_bind",
                                       strand = 1,
                                       qualifiers = {"note":[primer.primer.name]})

                    template.features.append(SeqFeature(FeatureLocation(start,end),
                                                        type ="primer_bind",
                                                        sub_features = [suba,subb],
                                                        location_operator= "join",
                                                        strand = 1,
                                                        qualifiers = {"note":[primer.primer.name]}))
            for primer in self.rev_primers:
                if primer.pos+len(primer.footprint)<=len(template):
                    start = primer.pos
                    end = primer.pos + len(footprint)
                    template.features.append(SeqFeature(FeatureLocation(start-1,end),
                                                        type ="primer_bind",
                                                        strand = -1,
                                                        qualifiers = {"note":[rp.name],
                                                                      "ApEinfo_fwdcolor":["green"],
                                                                      "ApEinfo_revcolor":["red"]}))
                else:
                    start = primer.pos
                    end = primer.pos+len(primer.footprint)-len(template)
                    suba = SeqFeature(FeatureLocation(start,len(template)),
                                      type ="primer_bind",
                                      strand = -1,
                                      qualifiers = {"note":[primer.primer.name]})
                    subb = SeqFeature(FeatureLocation(0,end),
                                      type ="primer_bind",
                                      strand = -1,
                                      qualifiers = {"note":[primer.primer.name]})
                    template.features.append(SeqFeature(FeatureLocation(start,end),
                                                        type ="primer_bind",
                                                        sub_features = [suba,subb],
                                                        location_operator= "join",
                                                        strand = -1,
                                                        qualifiers = {"note":[primer.primer.name]}))
            return template

    def pcr(*args,  **kwargs):
        '''pcr is a convenience function for Anneal to simplify its usage,
        especially from the command line. If more than one PCR product is
        formed, an exception is raised.

        args is any iterable of sequences or an iterable of iterables of sequences.
        args will be greedily flattened.

        Parameters
        ----------

        args : iterable containing sequence objects
            Several arguments are also accepted.

        limit : int = 13, optional
            limit length of the annealing part of the primers.

        Notes
        -----

        sequences in args could be of type:

        string
        Seq
        SeqRecord
        Dseqrecord

        The last sequence will be interpreted as the template
        all preceeding sequences as primers.

        This is a powerful function, use with care!

        Returns
        -------

        product : Dseqrecord
            a Dseqrecord object representing the PCR product.
            The direction of the PCR product will be the same as
            for the template sequence.

        Examples
        --------

        >>> import pydna
        >>> template = pydna.Dseqrecord("tacactcaccgtctatcattatctactatcgactgtatcatctgatagcac")
        >>> from Bio.SeqRecord import SeqRecord
        >>> p1 = pydna.read(">p1\\ntacactcaccgtctatcattatc", ds = False)
        >>> p2 = pydna.read(">p2\\ncgactgtatcatctgatagcac", ds = False).reverse_complement()
        >>> pydna.pcr(p1, p2, template)
        Dseqrecord(-51)
        >>> pydna.pcr([p1, p2], template)
        Dseqrecord(-51)
        >>> pydna.pcr((p1,p2,), template)
        Dseqrecord(-51)
        >>>

        '''

        import itertools
        from Bio.SeqRecord import SeqRecord
        # flatten args
        output = []
        stack = []
        stack.extend(reversed(args))
        while stack:
            top = stack.pop()
            if hasattr(top, "__iter__") and not isinstance(top, SeqRecord):
                stack.extend(reversed(top))
            else:
                output.append(top)
        new=[]

        for s in output:
            if isinstance(s, Seq):
                s = SeqRecord(s)
            elif isinstance(s, SeqRecord):
                pass
            elif hasattr(s, "watson"):
                s=s.watson
            elif isinstance(s, basestring):
                s = SeqRecord(Seq(s))
            else:
                raise TypeError("the record property needs to be a string, a Seq object or a SeqRecord object")
            new.append(s)

        anneal_primers = Anneal(new[:-1],
                                new[-1],
                                **kwargs)

        if anneal_primers:
            if len(anneal_primers.amplicons) == 1:
                return anneal_primers.amplicons.pop().pcr_product()
            elif len(anneal_primers.amplicons) == 0:
                raise Exception("No PCR products! {}".format(anneal_primers.report()))
            else:
                raise Exception("PCR not specific! {}".format(anneal_primers.report()))
        else:
            raise Exception(anneal_primers.report())
        return

    def basictm(primer, *args, **kwargs):
        '''Returns the melting temperature (Tm) of the primer using
        the basic formula.

        | Tm = (wA+xT)*2 + (yG+zC)*4 assumed 50mM monovalent cations
        |
        | w = number of A in primer
        | x = number of T in primer
        | y = number of G in primer
        | z = number of C in primer

        Parameters
        ----------
        primer : string
            Primer sequence 5'-3'

        Returns
        -------
        tm : int

        Examples
        --------
        >>> from pydna.amplify import basictm
        >>> basictm("ggatcc")
        20
        >>>

        '''
        primer = str(primer).lower()
        return (primer.count("a") + primer.count("t"))*2 + (primer.count("g") + primer.count("c"))*4

    # http://www.promega.com/techserv/tools/biomath/calc11.htm#melt_results

    def tmstaluc98(primer, dnac=50, saltc=50, **kwargs):
        '''Returns the melting temperature (Tm) of the primer using
        the nearest neighbour algorithm. Formula and thermodynamic data
        is taken from SantaLucia 1998. This implementation gives the same
        answer as the one provided by Biopython (See Examples).

        Thermodynamic data used:

        =====  ====  ====
        pair   dH    dS
        =====  ====  ====
        AA/TT  7.9   22.2
        AT/TA  7.2   20.4
        TA/AT  7.2   21.3
        CA/GT  8.5   22.7
        GT/CA  8.4   22.4
        CT/GA  7.8   21.0
        GA/CT  8.2   22.2
        CG/GC  10.6  27.2
        GC/CG  9.8   24.4
        GG/CC  8.0   19.9
        =====  ====  ====

        Parameters
        ----------
        primer : string
            Primer sequence 5'-3' in UPPERCASE

        Returns
        -------
        tm : float
            tm of the primer

        References
        ----------
        .. [1] J. SantaLucia, “A Unified View of Polymer, Dumbbell, and Oligonucleotide DNA Nearest-neighbor Thermodynamics,” Proceedings of the National Academy of Sciences 95, no. 4 (1998): 1460.

        Examples
        --------

        >>> from pydna.amplify import tmstaluc98
        >>> from Bio.SeqUtils.MeltingTemp import Tm_staluc
        >>> tmstaluc98("ACGTCATCGACACTATCATCGAC")
        54.55597724052518
        >>> Tm_staluc("ACGTCATCGACACTATCATCGAC")
        54.555977240525124
        >>>

        '''

        nntermsl={  "AA": (7.9  , 22.2),
                    "TT": (7.9  , 22.2),
                    "AT": (7.2  , 20.4),
                    "TA": (7.2  , 21.3),
                    "CA": (8.5  , 22.7),
                    "TG": (8.5  , 22.7),
                    "GT": (8.4  , 22.4),
                    "AC": (8.4  , 22.4),
                    "CT": (7.8  , 21.0),
                    "AG": (7.8  , 21.0),
                    "GA": (8.2  , 22.2),
                    "TC": (8.2  , 22.2),
                    "CG": (10.6 , 27.2),
                    "GC": (9.8  , 24.4),
                    "GG": (8    , 19.9),
                    "CC": (8    , 19.9),
                    "A" : (0    , 0   ),
                    "C" : (0    , 0   ),
                    "G" : (0    , 0   ),
                    "T" : (0    , 0   )  }

        helixinit = {   "G": (-0.1 ,2.8),
                        "C": (-0.1 ,2.8),
                        "A": (-2.3, -4.1),
                        "T": (-2.3, -4.1) }

        dH, dS = helixinit[primer[0]]
        H ,  S = helixinit[primer[-1]]
        dH = dH+H
        dS = dS+S

        for p in range(len(primer)):
            dn = primer[p:p+2]
            H,S = nntermsl[dn]
            dH+=H
            dS+=S
        R = 1.987 # universal gas constant in Cal/degrees C*Mol
        k = (dnac/4.0)*1e-9
        dS = dS-0.368*(len(primer)-1)*math.log(float(saltc)/1e3)
        tm = ((1000* (-dH))/(-dS+(R * (math.log(k)))))-273.15
        return tm

    def tmbreslauer86(primer, dnac=500.0, saltc=50,thermodynamics=False):
        '''Returns the melting temperature (Tm) of the primer using
        the nearest neighbour algorithm. Formula and thermodynamic data
        is taken from Breslauer 1986. These data are no longer widely used.

        Breslauer 1986, table 2 [1]

        =====  ===== ====   ===
        pair   dH    dS     dG
        =====  ===== ====   ===
        AA/TT  9.1   24.0   1.9
        AT/TA  8.6   23.9   1.5
        TA/AT  6.0   16.9   0.9
        CA/GT  5.8   12.9   1.9
        GT/CA  6.5   17.3   1.3
        CT/GA  7.8   20.8   1.6
        GA/CT  5.6   13.5   1.6
        CG/GC  11.9  27.8   3.6
        GC/CG  11.1  26.7   3.1
        GG/CC  11.0  26.6   3.1
        =====  ===== ====   ===

        Parameters
        ----------
        primer : string
            Primer sequence 5'-3' in UPPERCASE

        Returns
        -------
        tm : float

        References
        ----------
        .. [1] K.J. Breslauer et al., “Predicting DNA Duplex Stability from the Base Sequence,” Proceedings of the National Academy of Sciences 83, no. 11 (1986): 3746.

        Examples
        ---------

        >>> from pydna.amplify import tmbreslauer86
        >>> tmbreslauer86("ACGTCATCGACACTATCATCGAC")
        64.28863985851899

        '''

        nntermbr={  "AA": (9.1   ,24.0   ,1.9),
                    "TT": (9.1   ,24.0   ,1.9),
                    "AT": (8.6   ,23.9   ,1.5),
                    "TA": (6.0   ,16.9   ,0.9),
                    "CA": (5.8   ,12.9   ,1.9),
                    "TG": (5.8   ,12.9   ,1.9),
                    "GT": (6.5   ,17.3   ,1.3),
                    "AC": (6.5   ,17.3   ,1.3),
                    "CT": (7.8   ,20.8   ,1.6),
                    "AG": (7.8   ,20.8   ,1.6),
                    "GA": (5.6   ,13.5   ,1.6),
                    "TC": (5.6   ,13.5   ,1.6),
                    "CG": (11.9  ,27.8   ,3.6),
                    "GC": (11.1  ,26.7   ,3.1),
                    "GG": (11.0  ,26.6   ,3.1),
                    "CC": (11.0  ,26.6   ,3.1),
                    "A" : (0     , 0     ,0),
                    "C" : (0     , 0     ,0),
                    "G" : (0     , 0     ,0),
                    "T" : (0     , 0     ,0),     }
        dH=3.4
        dS=12.4
        dG=0
        for p in range(len(primer)):
            dn = primer[p:p+2]
            H,S,G = nntermbr[dn]
            dG+=G
            dH+=H
            dS+=S

        R = 1.9872          # universal gas constant in Cal/degrees C*Mol
        k = dnac*1E-9/2.0
        dH = dH - 5
        dS = dS-0.368*(len(primer)-1)*math.log(float(saltc)/1E3) # SantaLucia salt correction formula
        tm = 1000 * -dH /(-dS + R * math.log(k) )  - 273.15 # degrees Celsius

        if thermodynamics:
            return tm,dH,dS
        else:
            return tm

    def tmbresluc(primer, primerc=500.0, saltc=50, thermodynamics=False):
        '''Returns the tm for a primer using a formula adapted to polymerases
        with a DNA binding domain.

        Parameters
        ----------

        primer : string
            primer sequence 5'-3'

        primerc : float
           concentration (nM)

        saltc : float, optional
           Monovalent cation concentration (mM)

        thermodynamics : bool, optional
            prints details of the thermodynamic data to stdout. For
            debugging only.

        Returns
        -------
        tm : float
            the tm of the primer

        tm,dH,dS : tuple
            tm and dH and dS used for the calculation

        '''

        import collections

        dHBr = collections.defaultdict(dict)
        dSBr = collections.defaultdict(dict)

        dHBr[0][0] = -9100
        dSBr[0][0] = -24
        dHBr[0][1] = -7633.3000000000002
        dSBr[0][1] = -20.699999999999999
        dHBr[0][2] = -6500
        dSBr[0][2] = -17.300000000000001
        dHBr[0][3] = -8500
        dSBr[0][3] = -22.899999999999999
        dHBr[0][6] = -7800
        dSBr[0][6] = -20.800000000000001
        dHBr[0][7] = -8066.6999999999998
        dSBr[0][7] = -21.699999999999999
        dHBr[0][10] = -8200
        dSBr[0][10] = -22.399999999999999
        dHBr[0][12] = -7800
        dSBr[0][12] = -20.699999999999999
        dHBr[0][13] = -8000
        dSBr[0][13] = -21.5
        dHBr[0][17] = -8450
        dSBr[0][17] = -22.399999999999999
        dHBr[0][18] = -7150
        dSBr[0][18] = -19.100000000000001
        dHBr[0][19] = -8600
        dSBr[0][19] = -23.899999999999999
        dHBr[0][21] = -7800
        dSBr[0][21] = -20.699999999999999
        dHBr[0][22] = -8850
        dSBr[0][22] = -24
        dHBr[0][23] = -8000
        dSBr[0][23] = -21.5
        dHBr[0][24] = -7550
        dSBr[0][24] = -20.600000000000001
        dHBr[1][0] = -5800
        dSBr[1][0] = -14.4
        dHBr[1][1] = -8866.7000000000007
        dSBr[1][1] = -21.800000000000001
        dHBr[1][2] = -9233.2999999999993
        dSBr[1][2] = -22.300000000000001
        dHBr[1][3] = -7722.1999999999998
        dSBr[1][3] = -19.199999999999999
        dHBr[1][6] = -9566.7000000000007
        dSBr[1][6] = -22.399999999999999
        dHBr[1][7] = -7611.1000000000004
        dSBr[1][7] = -19.100000000000001
        dHBr[1][10] = -8683.2999999999993
        dSBr[1][10] = -21.600000000000001
        dHBr[1][12] = -7516.6999999999998
        dSBr[1][12] = -18.399999999999999
        dHBr[1][13] = -8100
        dSBr[1][13] = -20
        dHBr[1][17] = -7683.3000000000002
        dSBr[1][17] = -18.399999999999999
        dHBr[1][18] = -9400
        dSBr[1][18] = -22.399999999999999
        dHBr[1][19] = -7800
        dSBr[1][19] = -20.699999999999999
        dHBr[1][21] = -8200
        dSBr[1][21] = -19.699999999999999
        dHBr[1][22] = -6800
        dSBr[1][22] = -17.600000000000001
        dHBr[1][23] = -8100
        dSBr[1][23] = -20
        dHBr[1][24] = -8516.7000000000007
        dSBr[1][24] = -21.5
        dHBr[2][0] = -5800
        dSBr[2][0] = -12.9
        dHBr[2][1] = -10233.299999999999
        dSBr[2][1] = -25.100000000000001
        dHBr[2][2] = -11000
        dSBr[2][2] = -26.600000000000001
        dHBr[2][3] = -8500
        dSBr[2][3] = -20.5
        dHBr[2][6] = -11900
        dSBr[2][6] = -27.800000000000001
        dHBr[2][7] = -8200
        dSBr[2][7] = -20.100000000000001
        dHBr[2][10] = -9850
        dSBr[2][10] = -24.300000000000001
        dHBr[2][12] = -8400
        dSBr[2][12] = -19.800000000000001
        dHBr[2][13] = -9125
        dSBr[2][13] = -22
        dHBr[2][17] = -8850
        dSBr[2][17] = -20.399999999999999
        dHBr[2][18] = -11450
        dSBr[2][18] = -27.199999999999999
        dHBr[2][19] = -7800
        dSBr[2][19] = -20.800000000000001
        dHBr[2][21] = -9566.7000000000007
        dSBr[2][21] = -22.399999999999999
        dHBr[2][22] = -6800
        dSBr[2][22] = -16.899999999999999
        dHBr[2][23] = -9125
        dSBr[2][23] = -22
        dHBr[2][24] = -9400
        dSBr[2][24] = -23.699999999999999
        dHBr[3][0] = -6900
        dSBr[3][0] = -18.100000000000001
        dHBr[3][1] = -8000
        dSBr[3][1] = -20.300000000000001
        dHBr[3][2] = -7733.3000000000002
        dSBr[3][2] = -19.199999999999999
        dHBr[3][3] = -7722.1999999999998
        dSBr[3][3] = -20
        dHBr[3][6] = -8200
        dSBr[3][6] = -20.100000000000001
        dHBr[3][7] = -7566.6999999999998
        dSBr[3][7] = -19.699999999999999
        dHBr[3][10] = -8133.3000000000002
        dSBr[3][10] = -20.899999999999999
        dHBr[3][12] = -7316.6999999999998
        dSBr[3][12] = -18.699999999999999
        dHBr[3][13] = -7725
        dSBr[3][13] = -19.800000000000001
        dHBr[3][17] = -7550
        dSBr[3][17] = -19.100000000000001
        dHBr[3][18] = -7966.6999999999998
        dSBr[3][18] = -19.600000000000001
        dHBr[3][19] = -8066.6999999999998
        dSBr[3][19] = -21.699999999999999
        dHBr[3][21] = -7611.1000000000004
        dSBr[3][21] = -19.100000000000001
        dHBr[3][22] = -7483.3000000000002
        dSBr[3][22] = -19.899999999999999
        dHBr[3][23] = -7725
        dSBr[3][23] = -19.800000000000001
        dHBr[3][24] = -7900
        dSBr[3][24] = -20.5
        dHBr[6][0] = -5600
        dSBr[6][0] = -13.5
        dHBr[6][1] = -9533.2999999999993
        dSBr[6][1] = -23.5
        dHBr[6][2] = -11100
        dSBr[6][2] = -26.699999999999999
        dHBr[6][3] = -7700
        dSBr[6][3] = -19.100000000000001
        dHBr[6][6] = -11000
        dSBr[6][6] = -26.600000000000001
        dHBr[6][7] = -7733.3000000000002
        dSBr[6][7] = -19.199999999999999
        dHBr[6][10] = -8750
        dSBr[6][10] = -22
        dHBr[6][12] = -8350
        dSBr[6][12] = -20.100000000000001
        dHBr[6][13] = -8550
        dSBr[6][13] = -21
        dHBr[6][17] = -8300
        dSBr[6][17] = -20.100000000000001
        dHBr[6][18] = -11050
        dSBr[6][18] = -26.699999999999999
        dHBr[6][19] = -6500
        dSBr[6][19] = -17.300000000000001
        dHBr[6][21] = -9233.2999999999993
        dSBr[6][21] = -22.300000000000001
        dHBr[6][22] = -6050
        dSBr[6][22] = -15.4
        dHBr[6][23] = -8550
        dSBr[6][23] = -21
        dHBr[6][24] = -8800
        dSBr[6][24] = -22
        dHBr[7][0] = -6966.6999999999998
        dSBr[7][0] = -17.899999999999999
        dHBr[7][1] = -8233.2999999999993
        dSBr[7][1] = -20.800000000000001
        dHBr[7][2] = -7700
        dSBr[7][2] = -19.100000000000001
        dHBr[7][3] = -7988.8999999999996
        dSBr[7][3] = -20.399999999999999
        dHBr[7][6] = -8500
        dSBr[7][6] = -20.5
        dHBr[7][7] = -7722.1999999999998
        dSBr[7][7] = -20
        dHBr[7][10] = -8500
        dSBr[7][10] = -21.699999999999999
        dHBr[7][12] = -7333.3000000000002
        dSBr[7][12] = -18.5
        dHBr[7][13] = -7916.6999999999998
        dSBr[7][13] = -20.100000000000001
        dHBr[7][17] = -7733.3000000000002
        dSBr[7][17] = -19.199999999999999
        dHBr[7][18] = -8100
        dSBr[7][18] = -19.800000000000001
        dHBr[7][19] = -8500
        dSBr[7][19] = -22.899999999999999
        dHBr[7][21] = -7722.1999999999998
        dSBr[7][21] = -19.199999999999999
        dHBr[7][22] = -7733.3000000000002
        dSBr[7][22] = -20.399999999999999
        dHBr[7][23] = -7916.6999999999998
        dSBr[7][23] = -20.100000000000001
        dHBr[7][24] = -8100
        dSBr[7][24] = -21
        dHBr[10][0] = -5800
        dSBr[10][0] = -15.199999999999999
        dHBr[10][1] = -8183.3000000000002
        dSBr[10][1] = -20.199999999999999
        dHBr[10][2] = -8350
        dSBr[10][2] = -20.100000000000001
        dHBr[10][3] = -7333.3000000000002
        dSBr[10][3] = -18.5
        dHBr[10][6] = -8400
        dSBr[10][6] = -19.800000000000001
        dHBr[10][7] = -7316.6999999999998
        dSBr[10][7] = -18.699999999999999
        dHBr[10][10] = -8100
        dSBr[10][10] = -20.199999999999999
        dHBr[10][12] = -7075
        dSBr[10][12] = -17.699999999999999
        dHBr[10][13] = -7587.5
        dSBr[10][13] = -18.899999999999999
        dHBr[10][17] = -7100
        dSBr[10][17] = -17.5
        dHBr[10][18] = -8375
        dSBr[10][18] = -19.899999999999999
        dHBr[10][19] = -7800
        dSBr[10][19] = -20.699999999999999
        dHBr[10][21] = -7516.6999999999998
        dSBr[10][21] = -18.399999999999999
        dHBr[10][22] = -6800
        dSBr[10][22] = -17.899999999999999
        dHBr[10][23] = -7587.5
        dSBr[10][23] = -18.899999999999999
        dHBr[10][24] = -8075
        dSBr[10][24] = -20.399999999999999
        dHBr[12][0] = -7450
        dSBr[12][0] = -18.5
        dHBr[12][1] = -8933.2999999999993
        dSBr[12][1] = -22.899999999999999
        dHBr[12][2] = -8750
        dSBr[12][2] = -22
        dHBr[12][3] = -8500
        dSBr[12][3] = -21.699999999999999
        dHBr[12][6] = -9850
        dSBr[12][6] = -24.300000000000001
        dHBr[12][7] = -8133.3000000000002
        dSBr[12][7] = -20.899999999999999
        dHBr[12][10] = -9025
        dSBr[12][10] = -23.300000000000001
        dHBr[12][12] = -8100
        dSBr[12][12] = -20.199999999999999
        dHBr[12][13] = -8562.5
        dSBr[12][13] = -21.800000000000001
        dHBr[12][17] = -8650
        dSBr[12][17] = -21.399999999999999
        dHBr[12][18] = -9300
        dSBr[12][18] = -23.100000000000001
        dHBr[12][19] = -8200
        dSBr[12][19] = -22.399999999999999
        dHBr[12][21] = -8683.2999999999993
        dSBr[12][21] = -21.600000000000001
        dHBr[12][22] = -7825
        dSBr[12][22] = -20.399999999999999
        dHBr[12][23] = -8562.5
        dSBr[12][23] = -21.800000000000001
        dHBr[12][24] = -8475
        dSBr[12][24] = -22.199999999999999
        dHBr[13][0] = -6625
        dSBr[13][0] = -16.800000000000001
        dHBr[13][1] = -8558.2999999999993
        dSBr[13][1] = -21.5
        dHBr[13][2] = -8550
        dSBr[13][2] = -21
        dHBr[13][3] = -7916.6999999999998
        dSBr[13][3] = -20.100000000000001
        dHBr[13][6] = -9125
        dSBr[13][6] = -22
        dHBr[13][7] = -7725
        dSBr[13][7] = -19.800000000000001
        dHBr[13][10] = -8562.5
        dSBr[13][10] = -21.800000000000001
        dHBr[13][12] = -7587.5
        dSBr[13][12] = -18.899999999999999
        dHBr[13][13] = -8075
        dSBr[13][13] = -20.300000000000001
        dHBr[13][17] = -7875
        dSBr[13][17] = -19.399999999999999
        dHBr[13][18] = -8837.5
        dSBr[13][18] = -21.5
        dHBr[13][19] = -8000
        dSBr[13][19] = -21.5
        dHBr[13][21] = -8100
        dSBr[13][21] = -20
        dHBr[13][22] = -7312.5
        dSBr[13][22] = -19.199999999999999
        dHBr[13][23] = -8075
        dSBr[13][23] = -20.300000000000001
        dHBr[13][24] = -8275
        dSBr[13][24] = -21.300000000000001
        dHBr[17][0] = -7350
        dSBr[17][0] = -18.800000000000001
        dHBr[17][1] = -8583.2999999999993
        dSBr[17][1] = -22.100000000000001
        dHBr[17][2] = -8800
        dSBr[17][2] = -22
        dHBr[17][3] = -8100
        dSBr[17][3] = -21
        dHBr[17][6] = -9400
        dSBr[17][6] = -23.699999999999999
        dHBr[17][7] = -7900
        dSBr[17][7] = -20.5
        dHBr[17][10] = -8475
        dSBr[17][10] = -22.199999999999999
        dHBr[17][12] = -8075
        dSBr[17][12] = -20.399999999999999
        dHBr[17][13] = -8275
        dSBr[17][13] = -21.300000000000001
        dHBr[17][17] = -8375
        dSBr[17][17] = -21.199999999999999
        dHBr[17][18] = -9100
        dSBr[17][18] = -22.899999999999999
        dHBr[17][19] = -7550
        dSBr[17][19] = -20.600000000000001
        dHBr[17][21] = -8516.7000000000007
        dSBr[17][21] = -21.5
        dHBr[17][22] = -7450
        dSBr[17][22] = -19.699999999999999
        dHBr[17][23] = -8275
        dSBr[17][23] = -21.300000000000001
        dHBr[17][24] = -8175
        dSBr[17][24] = -21.300000000000001
        dHBr[18][0] = -5700
        dSBr[18][0] = -13.199999999999999
        dHBr[18][1] = -9883.2999999999993
        dSBr[18][1] = -24.300000000000001
        dHBr[18][2] = -11050
        dSBr[18][2] = -26.699999999999999
        dHBr[18][3] = -8100
        dSBr[18][3] = -19.800000000000001
        dHBr[18][6] = -11450
        dSBr[18][6] = -27.199999999999999
        dHBr[18][7] = -7966.6999999999998
        dSBr[18][7] = -19.600000000000001
        dHBr[18][10] = -9300
        dSBr[18][10] = -23.100000000000001
        dHBr[18][12] = -8375
        dSBr[18][12] = -19.899999999999999
        dHBr[18][13] = -8837.5
        dSBr[18][13] = -21.5
        dHBr[18][17] = -8575
        dSBr[18][17] = -20.199999999999999
        dHBr[18][18] = -11250
        dSBr[18][18] = -26.899999999999999
        dHBr[18][19] = -7150
        dSBr[18][19] = -19.100000000000001
        dHBr[18][21] = -9400
        dSBr[18][21] = -22.399999999999999
        dHBr[18][22] = -6425
        dSBr[18][22] = -16.100000000000001
        dHBr[18][23] = -8837.5
        dSBr[18][23] = -21.5
        dHBr[18][24] = -9100
        dSBr[18][24] = -22.899999999999999
        dHBr[19][0] = -6000
        dSBr[19][0] = -16.899999999999999
        dHBr[19][1] = -6833.3000000000002
        dSBr[19][1] = -16.800000000000001
        dHBr[19][2] = -5600
        dSBr[19][2] = -13.5
        dHBr[19][3] = -6966.6999999999998
        dSBr[19][3] = -17.899999999999999
        dHBr[19][6] = -5800
        dSBr[19][6] = -12.9
        dHBr[19][7] = -6900
        dSBr[19][7] = -18.100000000000001
        dHBr[19][10] = -7450
        dSBr[19][10] = -18.5
        dHBr[19][12] = -5800
        dSBr[19][12] = -15.199999999999999
        dHBr[19][13] = -6625
        dSBr[19][13] = -16.800000000000001
        dHBr[19][17] = -5900
        dSBr[19][17] = -14.9
        dHBr[19][18] = -5700
        dSBr[19][18] = -13.199999999999999
        dHBr[19][19] = -9100
        dSBr[19][19] = -24
        dHBr[19][21] = -5800
        dSBr[19][21] = -14.4
        dHBr[19][22] = -7550
        dSBr[19][22] = -20.5
        dHBr[19][23] = -6625
        dSBr[19][23] = -16.800000000000001
        dHBr[19][24] = -7350
        dSBr[19][24] = -18.800000000000001
        dHBr[21][0] = -6833.3000000000002
        dSBr[21][0] = -16.800000000000001
        dHBr[21][1] = -9133.2999999999993
        dSBr[21][1] = -23.100000000000001
        dHBr[21][2] = -9533.2999999999993
        dSBr[21][2] = -23.5
        dHBr[21][3] = -8233.2999999999993
        dSBr[21][3] = -20.800000000000001
        dHBr[21][6] = -10233.299999999999
        dSBr[21][6] = -25.100000000000001
        dHBr[21][7] = -8000
        dSBr[21][7] = -20.300000000000001
        dHBr[21][10] = -8933.2999999999993
        dSBr[21][10] = -22.899999999999999
        dHBr[21][12] = -8183.3000000000002
        dSBr[21][12] = -20.199999999999999
        dHBr[21][13] = -8558.2999999999993
        dSBr[21][13] = -21.5
        dHBr[21][17] = -8533.2999999999993
        dSBr[21][17] = -20.899999999999999
        dHBr[21][18] = -9883.2999999999993
        dSBr[21][18] = -24.300000000000001
        dHBr[21][19] = -7633.3000000000002
        dSBr[21][19] = -20.699999999999999
        dHBr[21][21] = -8866.7000000000007
        dSBr[21][21] = -21.800000000000001
        dHBr[21][22] = -7233.3000000000002
        dSBr[21][22] = -18.699999999999999
        dHBr[21][23] = -8558.2999999999993
        dSBr[21][23] = -21.5
        dHBr[21][24] = -8583.2999999999993
        dSBr[21][24] = -22.100000000000001
        dHBr[22][0] = -7550
        dSBr[22][0] = -20.5
        dHBr[22][1] = -7233.3000000000002
        dSBr[22][1] = -18.699999999999999
        dHBr[22][2] = -6050
        dSBr[22][2] = -15.4
        dHBr[22][3] = -7733.3000000000002
        dSBr[22][3] = -20.399999999999999
        dHBr[22][6] = -6800
        dSBr[22][6] = -16.899999999999999
        dHBr[22][7] = -7483.3000000000002
        dSBr[22][7] = -19.899999999999999
        dHBr[22][10] = -7825
        dSBr[22][10] = -20.399999999999999
        dHBr[22][12] = -6800
        dSBr[22][12] = -17.899999999999999
        dHBr[22][13] = -7312.5
        dSBr[22][13] = -19.199999999999999
        dHBr[22][17] = -7175
        dSBr[22][17] = -18.699999999999999
        dHBr[22][18] = -6425
        dSBr[22][18] = -16.100000000000001
        dHBr[22][19] = -8850
        dSBr[22][19] = -24
        dHBr[22][21] = -6800
        dSBr[22][21] = -17.600000000000001
        dHBr[22][22] = -8200
        dSBr[22][22] = -22.199999999999999
        dHBr[22][23] = -7312.5
        dSBr[22][23] = -19.199999999999999
        dHBr[22][24] = -7450
        dSBr[22][24] = -19.699999999999999
        dHBr[23][0] = -6625
        dSBr[23][0] = -16.800000000000001
        dHBr[23][1] = -8558.2999999999993
        dSBr[23][1] = -21.5
        dHBr[23][2] = -8550
        dSBr[23][2] = -21
        dHBr[23][3] = -7916.6999999999998
        dSBr[23][3] = -20.100000000000001
        dHBr[23][6] = -9125
        dSBr[23][6] = -22
        dHBr[23][7] = -7725
        dSBr[23][7] = -19.800000000000001
        dHBr[23][10] = -8562.5
        dSBr[23][10] = -21.800000000000001
        dHBr[23][12] = -7587.5
        dSBr[23][12] = -18.899999999999999
        dHBr[23][13] = -8075
        dSBr[23][13] = -20.300000000000001
        dHBr[23][17] = -7875
        dSBr[23][17] = -19.399999999999999
        dHBr[23][18] = -8837.5
        dSBr[23][18] = -21.5
        dHBr[23][19] = -8000
        dSBr[23][19] = -21.5
        dHBr[23][21] = -8100
        dSBr[23][21] = -20
        dHBr[23][22] = -7312.5
        dSBr[23][22] = -19.199999999999999
        dHBr[23][23] = -8075
        dSBr[23][23] = -20.300000000000001
        dHBr[23][24] = -8275
        dSBr[23][24] = -21.300000000000001
        dHBr[24][0] = -5900
        dSBr[24][0] = -14.9
        dHBr[24][1] = -8533.2999999999993
        dSBr[24][1] = -20.899999999999999
        dHBr[24][2] = -8300
        dSBr[24][2] = -20.100000000000001
        dHBr[24][3] = -7733.3000000000002
        dSBr[24][3] = -19.199999999999999
        dHBr[24][6] = -8850
        dSBr[24][6] = -20.399999999999999
        dHBr[24][7] = -7550
        dSBr[24][7] = -19.100000000000001
        dHBr[24][10] = -8650
        dSBr[24][10] = -21.399999999999999
        dHBr[24][12] = -7100
        dSBr[24][12] = -17.5
        dHBr[24][13] = -7875
        dSBr[24][13] = -19.399999999999999
        dHBr[24][17] = -7375
        dSBr[24][17] = -17.600000000000001
        dHBr[24][18] = -8575
        dSBr[24][18] = -20.199999999999999
        dHBr[24][19] = -8450
        dSBr[24][19] = -22.399999999999999
        dHBr[24][21] = -7683.3000000000002
        dSBr[24][21] = -18.399999999999999
        dHBr[24][22] = -7175
        dSBr[24][22] = -18.699999999999999
        dHBr[24][23] = -7875
        dSBr[24][23] = -19.399999999999999
        dHBr[24][24] = -8375
        dSBr[24][24] = -21.199999999999999

        saltc = float(saltc)/1000
        pri  = primerc/10E7
        dS = -12.4
        dH = -3400

        STR = primer.lower();

        for i in range(len(STR)-1):
            n1=ord(STR[i])
            n2=ord(STR[i+1])
            dH += dHBr[n1 - 97][n2 - 97]
            dS += dSBr[n1 - 97][n2 - 97]

        tm = (dH / (1.9872 * math.log(pri / 1600) + dS) + (16.6 * math.log(saltc)) / math.log(10)) - 273.15

        if thermodynamics:
            return tm,dH,dS
        else:
            return tm

    '''

    Examples
    --------

    >>> from pydna import read
    >>> from pydna import Anneal
    >>> template = read(">MyTemplate\\ngctactacacacgtactgactgcctccaagatagagtcagtaaccaca")
    >>> fp = read(">ForwardPrimer\\ngctactacacacgtactgactg", ds=False )
    >>> rp = read(">ReversePrimer\\ntgtggttactgactctatcttg", ds=False )
    >>> primers = (fp,rp)
    >>> p = Anneal(primers, template)
    >>> p
    Anneal(amplicons = 1)
    >>> print p.report()
    Template MyTemplate 48 nt linear:
    Primer ForwardPrimer anneals at position 21
    Primer ReversePrimer anneals reverse at position 27
    >>> p.amplicons
    [Amplicon(48)]
    >>> prod = p.amplicons.pop()
    >>> prod
    Amplicon(48)
    >>> prod.pcr_product()
    Dseqrecord(-48)
    >>> prod.pcr_product().seq
    Dseq(-48)
    gctactacacacgtactgac...agatagagtcagtaaccaca
    cgatgatgtgtgcatgactg...tctatctcagtcattggtgt
    >>> print prod.detailed_figure()
    5gctactacacacgtactgactg3
     |||||||||||||||||||||| tm 52.6 (dbd) 58.3
    5gctactacacacgtactgactg...caagatagagtcagtaaccaca3
    3cgatgatgtgtgcatgactgac...gttctatctcagtcattggtgt5
                              |||||||||||||||||||||| tm 49.1 (dbd) 57.7
                             3gttctatctcagtcattggtgt5
    >>>

    '''

    if __name__=="__main__":
        import doctest
        doctest.testmod()

        t = Dseqrecord("gctactacacacgtactgactgCCnnnTCcaagatagagtcagtaaccaca",circular=False)
        p1 = SeqRecord(Seq("NUNgctactacacacgtactgactg"))
        p2 = SeqRecord(Seq("NCNtgtggttactgactctatct"))
        assert pcr(p1,p2,t).seq.watson == "NUNgctactacacacgtactgactgCCnnnTCcaagatagagtcagtaaccacaNGN"
        t = Dseqrecord("gctactacacacgtactgactgCCnnnTCcaagatagagtcagtaaccaca",circular=True)
        p1 = SeqRecord(Seq("NANgctactacacacgtactgactg"))
        p2 = SeqRecord(Seq("NCNtgtggttactgactctatct"))
        assert pcr(p1,p2,t).seq.watson == "NANgctactacacacgtactgactgCCnnnTCcaagatagagtcagtaaccacaNGN"

        '''

        NANgctactacacacgtactgactg
           cgatgatgtgtgcatgactgacGGnnnAGgttctatctcagtcattggtgt

                                 -------------29--------------

           --------------31---------------
           gctactacacacgtactgactgCCnnnTCcaagatagagtcagtaaccaca
                                          tctatctcagtcattggtgtNCN

        '''

        t = Dseqrecord("caagatagagtcagtaaccacagctactacacacgtactgactgCCnnnTC", circular=True)
        p1 = SeqRecord(Seq("NANgctactacacacgtactgactg"))
        p2 = SeqRecord(Seq("NCNtgtggttactgactctatct"))
        assert pcr(p1,p2,t).seq.watson == "NANgctactacacacgtactgactgCCnnnTCcaagatagagtcagtaaccacaNGN"

        '''
           --
           caagatagagtcagtaaccacagctactacacacgtactgactgCCnnnTC
             tctatctcagtcattggtgtNCN
                              NANgctactacacacgtactgactg
           gttctatctcagtcattggtgtcgatgatgtgtgcatgactgacGGnnnAG
                                                       -------
           ............................................

        '''

        import textwrap

        oligos_table_3='''
        GCGCGC
        CGTCGACG
        GAAGCTTC
        GGAATTCC
        GGTATACC
        GCGAATTCGC
        CAAAAAG
        CAAACAAAG
        CAAAAAAAG
        CAAATAAAG
        CAAAGAAAG
        CGCGTACGCGTACGCG'''

        for oligo in textwrap.dedent(oligos_table_3).splitlines():
            pass #print tmbreslauer86(oligo, thermodynamics = True)[3]

        oligos = '''
                    CAGTCAGTACGTACGTGTACTGCCGTA
                    CAGTCAGTACGTACGTGTACTGCCGTAGG
                    CAGTCAGTACGTACGTGTACTGCCG
                 '''

        for oligo in textwrap.dedent(oligos).splitlines():
            print oligo
            print tmbreslauer86(oligo)#, thermodynamics = True)
            print tmbresluc(oligo)    #, thermodynamics = True)
            print

         #                                  FZ     Primer3Plus
         #  CAGTCAGTACGTACGTGTACTGCCGTA     68.83  68.69
         #  CAGTCAGTACGTACGTGTACTGCCGTAGG   71.5   71.70
         #  CAGTCAGTACGTACGTGTACTGCCG       68.7   68.59

        s = parse('''
            >524_pFA6aF (29-mer)
            cacatacgatttaggtgacactatagaac
            >523_AgTEF1tpR (21-mer)
            ggttgtttatgttcggatgtg
            ''')

        from pydna import Dseq

        f= Dseq("cacatacgatttaggtgacactatagaaccacatccgaacataaacaacc",
                "gctaaatccactgtgatatcttggtgtaggcttgtatttgttgg"[::-1])

        print Anneal(s,Dseqrecord(f)).amplicons.pop().detailed_figure()

    #    from pydna import read

    #    template = Dseqrecord("tacactcaccgtctatcattatctactatcgactgtatcatctgatagcac")
    #    p1 = read(">p1\ntacactcaccgtctatcattatc", ds= False)
    #    p2 = read(">p1\ncgactgtatcatctgatagcac", ds = False).reverse_complement()
    #    amplicon_list = Anneal((p1,p2), template).amplicons
    #    amplicon = amplicon_list.pop()
    #    print amplicon.detailed_figure()

### Navigation

- index
- modules |
- Pydna 0.1 documentation »
- Module code »