
:>"^
                 @   s  d  Z  d d l m Z m Z d d l m Z d d l m Z d d l m	 Z	 d d l
 m Z d d l m Z m Z d d l m Z m Z m Z Gd	 d
   d
 e	  Z d d   Z d e d d d d  Z d d   Z d e d d  Z e d d  Z d d   Z d d   Z d d   Z d d   Z e d  d!  Z d" d# d$  Z e d% d&  Z e d' d(  Z  d) d*   Z! e d+ d,  Z" e d- d.  Z# d/ d0   Z$ d1 d2   Z% e& d3 k rd d4 l' m( Z( e(   d S)5zCode for dealing with coding sequence.

CodonSeq class is inherited from Seq class. This is the core class to
deal with sequences in CodonAlignment in biopython.

    )divisionprint_function)permutations)log)Seq)	SeqRecord)generic_dna_ungap)CodonAlphabetdefault_codon_alphabetdefault_codon_tablec               @   s   e  Z d  Z d Z d e d d d d  Z d d   Z d	 d
   Z d d   Z e	 d d d d d  Z
 e d d  Z d d   Z e	 d d d  Z d d d  Z e e d d d   Z d S)CodonSeqaK  CodonSeq is designed to be within the SeqRecords of a CodonAlignment class.

    CodonSeq is useful as it allows the user to specify
    reading frame when translate CodonSeq

    CodonSeq also accepts codon style slice by calling
    get_codon() method.

    **Important:** Ungapped CodonSeq can be any length if you
    specify the rf_table. Gapped CodonSeq should be a
    multiple of three.

    >>> codonseq = CodonSeq("AAATTTGGGCCAAATTT", rf_table=(0,3,6,8,11,14))
    >>> print(codonseq.translate())
    KFGAKF

    test get_full_rf_table method

    >>> p = CodonSeq('AAATTTCCCGG-TGGGTTTAA', rf_table=(0, 3, 6, 9, 11, 14, 17))
    >>> full_rf_table = p.get_full_rf_table()
    >>> print(full_rf_table)
    [0, 3, 6, 9, 12, 15, 18]
    >>> print(p.translate(rf_table=full_rf_table, ungap_seq=False))
    KFPPWV*
    >>> p = CodonSeq('AAATTTCCCGGGAA-TTTTAA', rf_table=(0, 3, 6, 9, 14, 17))
    >>> print(p.get_full_rf_table())
    [0, 3, 6, 9, 12.0, 15, 18]
    >>> p = CodonSeq('AAA------------TAA', rf_table=(0, 3))
    >>> print(p.get_full_rf_table())
    [0, 3.0, 6.0, 9.0, 12.0, 15]

     -Nc             C   s  t  j |  | j   d | | |  _ t | t  s@ t d   | d k r
|  j j | d  } t	 |   d d k r t
 d   t t d d	   t t	 |     |  _ x|  j D]L } |  j | | d  | j k r t
 d
 j |  j | | d     q Wn t | t t f  s+t d   t d d   | D  sPt d   |  j j | d  } xN | D]F } | | | d  | j k rlt
 d j | | | d     qlW| |  _ d S)zInitialize the class.alphabetz0Input alphabet should be a CodonAlphabet object.Nr      r   zJSequence length is not a multiple of three (i.e. a whole number of codons)c             S   s   |  d d k S)Nr   r    )xr   r   </tmp/pip-build-ww9dw3qa/biopython/Bio/codonalign/codonseq.py<lambda>S   s    z#CodonSeq.__init__.<locals>.<lambda>z2Sequence contain codon not in the alphabet ({0})! z)rf_table should be a tuple or list objectc             s   s   |  ] } t  | t  Vq d  S)N)
isinstanceint).0ir   r   r   	<genexpr>b   s    z$CodonSeq.__init__.<locals>.<genexpr>zSElements in rf_table should be int that specify the codon positions of the sequencez7Sequence contain undefined letters from alphabet ({0})!)r   __init__uppergap_charr   r
   	TypeError_datareplacelen
ValueErrorlistfilterrangerf_tablelettersformattupleall)selfdatar   r   r&   Zseq_ungappedr   r   r   r   r   9   s0    	 	$	zCodonSeq.__init__c             C   s   t  |  j | d t S)Nr   )r   r   r   )r+   indexr   r   r   __getitem__n   s    zCodonSeq.__getitem__c                s   t  d d    j D  d k r. t d   t | t  r~ | d
 k rf  j | d | d d  S j | d d  SnJ t t    d       f d d   } | |  } t | d	  j Sd S)z&Get the index codon from the sequence.c             S   s   h  |  ] } | d   q S)r   r   )r   r   r   r   r   	<setcomp>t   s   	 z%CodonSeq.get_codon.<locals>.<setcomp>   z}frameshift detected. CodonSeq object is not able to deal with codon sequence with frameshift. Please use normal slice option.r   Nc                sH     |  } d } x1 | D]) } |  j  | d | d d  7} q W| S)Nr   r   )r   )pZaa_slicecodon_slicer   )aa_indexr+   r   r   cslice   s
    
'z"CodonSeq.get_codon.<locals>.cslicer   )	r!   r&   RuntimeErrorr   r   r   r%   r   r   )r+   r-   r4   r2   r   )r3   r+   r   	get_codonr   s    "zCodonSeq.get_codonc             C   s   t  |  j  S)z,Return the number of codons in the CodonSeq.)r!   r&   )r+   r   r   r   get_codon_num   s    zCodonSeq.get_codon_num*Tc       
      C   s  g  } | r' |  j  j |  j d  } n	 |  j  } | d k rE |  j } d } xA| D]9} t | t  rz | j d  qR n d | | | d  k r| d	 k s | | d k r | } | | | d  j d d  d d  }	 q*| | d k r*| | | d  }	 | } n | | | d  }	 | } |	 | j k rI| j |  qR y | j | j |	  WqR t	 k
 rt
 d j |	    YqR XqR Wd j |  S)
a1  Translate the CodonSeq based on the reading frame in rf_table.

        It is possible for the user to specify
        a rf_table at this point. If you want to include
        gaps in the translated sequence, this is the only
        way. ungap_seq should be set to true for this
        purpose.
        r   Nr0   r   r      zOUnknown codon detected ({0}). Did you forget to specify the ungap_seq argument?r5   r5   )r   r    r   r&   r   floatappendstop_codonsforward_tableKeyErrorr6   r(   join)
r+   codon_tablestop_symbolr&   	ungap_seqZamino_acidsZtr_seqr1   r   codonr   r   r   	translate   s:    
		-		zCodonSeq.translatec             C   s   t  |  j t  S)zConvert DNA to seq object.)r   r   r   )r+   r   r   r   r   toSeq   s    zCodonSeq.toSeqc                s^  |  j  j d d      f d d   |  j D } |  j d g } xQ t d t |  j d d   d  D]) } | j |  j | |  j | d  qh Wg  } d } xt d d	   t t |  j     D]} |  j  | | d
  |  j d
 k r| j | d  n| | d k r.| j |  | d 7} n | | d k rt |  j  | d
 |  j d d   } | d
 k r| j | | |  nM | d k r| j | d | |  n% | d k r| j | d | |  | d 7} n! | | d k r| j | d  y; t |  j  | | d
  j d d   } | | | 8<Wq t k
 rUYq Xq W| S)zReturn full rf_table of the CodonSeq records.

        A full rf_table is different from a normal rf_table in that
        it translate gaps in CodonSeq. It is helpful to construct
        alignment containing frameshift.
        r   r   c                s$   g  |  ] }   | | d    q S)r   r   )r   r   )rC   r   r   
<listcomp>   s   	 z.CodonSeq.get_full_rf_table.<locals>.<listcomp>r   r0   Nc             S   s   |  d d k S)Nr   r   r   )r   r   r   r   r      s    z,CodonSeq.get_full_rf_table.<locals>.<lambda>r   g           r5   )r5   rI   )	r   r    r&   r%   r!   r<   r$   r   	Exception)r+   	codon_lstZrelative_posr   full_rf_table	codon_numZgap_statZthis_lenr   )rC   r   get_full_rf_table   s<    -'+$))	zCodonSeq.get_full_rf_tablec          	   C   s.   |  j    } |  j d | d | d | d d  S)z,Apply full translation with gaps considered.rA   rB   r&   rC   F)rN   rE   )r+   rA   rB   rL   r   r   r   full_translate   s    zCodonSeq.full_translatec             C   s   t  |  j d  rv | s' |  j j } n= | |  j j k rd t d t |  t |  j j j  f   t |  j  } n | s t d   n	 |  j } t |  d k s t | t  r t d t |    t	 t |  j
  j | d  | d |  j S)z;Return a copy of the sequence without the gap character(s).r   z&Gap %s does not match %s from alphabetz3Gap character not given and not defined in alphabetr0   zUnexpected gap character, %sr   r&   )hasattrr   r   r"   reprr	   r!   r   strr   r   r    r&   )r+   Zgapalphar   r   r   ungap   s    %	"!zCodonSeq.ungapc             C   s<   | d k r |  | j  d | S|  | j  d | d | Sd S)z&Get codon sequence from sequence data.Nr   r&   )r   )clsseqr   r&   r   r   r   from_seq
  s    zCodonSeq.from_seq)__name__
__module____qualname____doc__r   r   r.   r7   r8   r   rE   r   rF   rN   rO   rT   classmethodrW   r   r   r   r   r      s    4/(r   c             C   s6  |  j    } g  } xt |  D]\ } } t | t  r | } y t | | d  } Wn t k
 rv | d } Yn Xt |  | |   } t |  d k r | j |  q.| j t | j     q t |  t |  t |  d   d k r| j d  q | j |  t |  t |  d   q W| S)zAList of codons according to full_rf_table for counting (PRIVATE).r0   r   z---)	rN   	enumerater   r   
IndexErrorrR   r!   r<   rT   )ZcodonseqrL   rK   r   kstartend
this_codonr   r   r   _get_codon_list  s"    ,+rc   NG86r0   Nc             C   s  t  |  t  r! t  | t  r! n? t  |  t  rT t  | t  rT |  j }  | j } n t d   t |  j    t | j    k r t d j t |  j    t | j       | d k r d } n$ | d k	 r | d k r t d   | d k r!d	 d l	 } | j
 d
 j |   d } t |   } t |  } g  }	 g  }
 xO t | |  D]> \ } } d | k rUd | k rU|	 j |  |
 j |  qUWd t d t d t d t i } | d k r| | |	 |
 | |  S| | |	 |
 | |  Sd S)a  Calculate dN and dS of the given two sequences.

    Available methods:
        - NG86  - `Nei and Gojobori (1986)`_ (PMID 3444411).
        - LWL85 - `Li et al. (1985)`_ (PMID 3916709).
        - ML    - `Goldman and Yang (1994)`_ (PMID 7968486).
        - YN00  - `Yang and Nielsen (2000)`_ (PMID 10666704).

    .. _`Nei and Gojobori (1986)`: http://www.ncbi.nlm.nih.gov/pubmed/3444411
    .. _`Li et al. (1985)`: http://www.ncbi.nlm.nih.gov/pubmed/3916709
    .. _`Goldman and Yang (1994)`: http://mbe.oxfordjournals.org/content/11/5/725
    .. _`Yang and Nielsen (2000)`: https://doi.org/10.1093/oxfordjournals.molbev.a026236

    Arguments:
     - codon_seq1 - CodonSeq or or SeqRecord that contains a CodonSeq
     - codon_seq2 - CodonSeq or or SeqRecord that contains a CodonSeq
     - w  - transition/transversion ratio
     - cfreq - Current codon frequency vector can only be specified
       when you are using ML method. Possible ways of
       getting cfreq are: F1x4, F3x4 and F61.

    zVcal_dn_ds accepts two CodonSeq objects or SeqRecord that contains CodonSeq as its seq!zBfull_rf_table length of seq1 ({0}) and seq2 ({1}) are not the sameNF3x4ZMLz8cfreq can only be specified when you are using ML methodF1x4F61r   zWUnknown cfreq ({0}). Only F1x4, F3x4 and F61 are acceptable. Use F3x4 in the following.r   rd   ZLWL85ZYN00)rf   re   rg   )r   r   r   rV   r   r!   rN   r6   r(   warningswarnrc   zipr<   _ml_ng86_lwl85_yn00)Z
codon_seq1Z
codon_seq2methodrA   r_   Zcfreqrh   Zseq1_codon_lstZseq2_codon_lstseq1seq2r   jZ	dnds_funcr   r   r   	cal_dn_ds/  s@    			
rs   c          	   C   s;  t  |  d | d | \ } } t  | d | d | \ } } | | d } | | d }	 d d g }
 xH t |  |  D]7 \ } } d d   t |
 t | | d |  D }
 qt W|
 d | } |
 d |	 } | d k  r t d t d d |   } n d } | d k  r+t d t d d |   } n d } | | f S)z$NG86 method main function (PRIVATE).rA   r_   g       @r   c             S   s    g  |  ] \ } } | |  q Sr   r   )r   mnr   r   r   rG   ~  s   	 z_ng86.<locals>.<listcomp>r0   r      g      @g      @g      ?g      g      gUUUUUU?r5   g      ?g      g      gUUUUUU?r5   )_count_site_NG86rj   _count_diff_NG86absr   )rp   rq   r_   rA   S_sites1N_sites1S_sites2N_sites2S_sitesN_sitesSNr   rr   pspndSdNr   r   r   rl   t  s&    		!!rl   c             C   sx  d } d } d } d } d } xM|  D]E} d g  d g  i }	 | j  d d  } | d	 k r^ q% xt |  D]\ }
 } x | D] } | | k r q~ | | k r | | k r t |  } | | |
 <d
 j |  } |	 d j |  q~ | | k r5| | k r5t |  } | | |
 <d
 j |  } |	 d j |  q~ t |  } | | |
 <d
 j |  } |	 d j |  q~ Wqk W| j | } d } } xX |	 d D]L } | | j k r| d 7} q| j | | k r| d 7} q| d 7} qWxX |	 d D]L } | | j k r| | 7} q| j | | k r2| | 7} q| | 7} qW| | d } | | | 7} | | | 7} q% W| | f S)a  Count synonymous and non-synonymous sites of a list of codons (PRIVATE).

    Arguments:
     - codon_lst - A three letter codon list from a CodonSeq object.
       This can be returned from _get_codon_list method.
     - k - transition/transversion rate ratio.

    r   AGTCZ
transitionZtransversionUz---r   r0   r   )r   r   )r   r   )r   r   r   r   )r    r]   r#   r@   r<   r>   r=   )rK   r_   rA   ZS_siteZN_sitepurine
pyrimidine
base_tuplerD   neighbor_codonru   r   rr   Zcodon_charsrb   aaZthis_codon_N_siteZthis_codon_S_siteZneighbor
norm_constr   r   r   rw     s\    	



rw   c                s  t  |  t  s  t  | t  rD t d j t |   t |     t |   d k sh t |  d k r t d j t |   t |     d d g } |  d k s | d k r | Sd   t   f d
 d   |  D  s t d j |     t   f d d   | D  s"t d j |    |  | k r2| Sg  } xD t t	 |  |   D]- \ } } | d | d k rN| j
 |  qNWt d d d  } t |  d k rd d   t	 | | |  | d |  D } nt |  d k r| j | } x| D] } |  d |  | | |  | d d  }	 d d   t	 | | |  |	 d | d d  D } d d   t	 | | |	 | d | d d  D } qWn]t |  d k r| j | } t t d d d g d   }
 g  } x|
 D]} |  d | d  | | d |  | d d d  } | d | d  | | d | | d d d  } | j
 | | f  d d   t	 | | |  | | d d  D } d d   t	 | | | | | d d   D } d d   t	 | | | | | d d!  D } qW| S)"zCount differences between two codons, three-letter string (PRIVATE).

    The function will take multiple pathways from codon1 to codon2
    into account.
    zM_count_diff_NG86 accepts string object to represent codon ({0}, {1} detected)r   z7codon should be three letter string ({0}, {1} detected)r   z---r   r   r   r   c             3   s   |  ] } |   k Vq d  S)Nr   )r   r   )r   r   r   r     s    z#_count_diff_NG86.<locals>.<genexpr>zNUnrecognized character detected in codon1 {0} (Codons consist of A, T, C or G)c             3   s   |  ] } |   k Vq d  S)Nr   )r   r   )r   r   r   r     s    zNUnrecognized character detected in codon2 {0} (Codons consist of A, T, C or G)r0   c             S   sX   d } } t  t t | j j |  | g    d k rD | | 7} n
 | | 7} | | f S)z4Compare two codon accounting for different pathways.r   r0   )r!   setmapr>   get)codon1codon2rA   weightsdZndr   r   r   compare_codon  s    

z'_count_diff_NG86.<locals>.compare_codonc             S   s    g  |  ] \ } } | |  q Sr   r   )r   r   rr   r   r   r   rG     s   	 z$_count_diff_NG86.<locals>.<listcomp>rA   rH   Nc             S   s    g  |  ] \ } } | |  q Sr   r   )r   r   rr   r   r   r   rG     s   	 r   g      ?c             S   s    g  |  ] \ } } | |  q Sr   r   )r   r   rr   r   r   r   rG     s   	 c             S   s    g  |  ] \ } } | |  q Sr   r   )r   r   rr   r   r   r   rG     s   	 c             S   s    g  |  ] \ } } | |  q Sr   r   )r   r   rr   r   r   r   rG     s   	 c             S   s    g  |  ] \ } } | |  q Sr   r   )r   r   rr   r   r   r   rG     s   	 )r   r   r   r   gUUUUUU?gUUUUUU?gUUUUUU?)r   rR   r   r(   typer!   r6   r*   r]   rj   r<   r   r>   r#   r   )r   r   rA   r   diff_posr   r_   r   	codon2_aa
temp_codonpaths	tmp_codonr1   tmp1tmp2r   )r   r   rx     sr     	$			"
	%*		66rx   c          	   C   sd  t  |  } d d g } d d g } d d g } x |  | D]u } | | }	 xb |	 D]Z }
 |
 d k ru | d d 7<qR |
 d k r | d d 7<qR |
 d k rR | d d 7<qR Wq; Wt |  d t |  d t |  d g } d g d } xr t |  |  D]a \ } } | d k s | d k s | | k r4q q d	 d
   t | t | | d |  D } q Wd d
   t | | d  D } | d d  } | d d  } d d
   t | |  D } d d
   | D } d | d | d | d | d | d | d d | d } d | d | d | d | d | d d | d d | d } | | f S)zlLWL85 method main function (PRIVATE).

    Nomenclature is according to Li et al. (1985), PMID 3916709.
    r   0r0   24g       @r:   z---c             S   s    g  |  ] \ } } | |  q Sr   r   )r   r   rr   r   r   r   rG   >  s   	 z_lwl85.<locals>.<listcomp>	fold_dictc             S   s    g  |  ] \ } } | |  q Sr   r   )r   r   rr   r   r   r   rG   C  s   	 rH   Nr   c          	   S   sP   g  |  ]F \ } } d t  d  d d | |  d t  d  d d |   q S)g      ?rH   r0   rv   g      ?g      ?)r   )r   r   rr   r   r   r   rG   F  s   	c             S   s,   g  |  ]" } d t  d  d d |   q S)g      ?rH   r0   g      ?)r   )r   r   r   r   r   rG   H  s   	 )_get_codon_foldsumrj   _diff_codon)rp   rq   r_   rA   codon_fold_dictZfold0Zfold2Zfold4rD   fold_numfLZPQr   r   PQr   Br   r   r   r   r   rm   $  s>    
-$ 	BFrm   c             C   sV   d d   } i  } x3 |  j  D]( } d | k r | | |  j   | | <q Wd | d <| S)zFClassify different position in a codon into different folds (PRIVATE).c       
      S   s>  d d d d h } d } t  |   } xt |  D]\ } } | t |  } g  } xX | D]P }	 |	 | | <y | j | d j |   WqZ t k
 r | j d  YqZ XqZ W| j | |   d k r | d 7} nX | j | |   d k r | d 7} n2 | j | |   d k r | d 7} n t d   | | | <q1 W| S)Nr   r   r   r   r   stopr   r   r0   rH   r   r   r   z3Unknown Error, cannot assign the position to a fold)r0   rH   )r#   r]   r   r<   r@   r?   countr6   )
rD   r>   basefoldZcodon_base_lstr   bZ
other_baser   rr   r   r   r   find_fold_classP  s*    
z(_get_codon_fold.<locals>.find_fold_classr   z---)r>   )rA   r   Z
fold_tablerD   r   r   r   r   N  s    
r   c             C   s7  d } } } } } } | |  }	 d }
 d } xt  t |  |   D]\ } \ } } | | k r | |
 k r | |
 k r |	 | d k r | d 7} nN |	 | d k r | d 7} n1 |	 | d	 k r | d 7} n t d
 |	 |   | | k rv| | k rv| | k rv|	 | d k r(| d 7} nN |	 | d k rE| d 7} n1 |	 | d	 k rb| d 7} n t d
 |	 |   | | k rF | |
 k r| | k s| | k rF | |
 k rF |	 | d k r| d 7} qF |	 | d k r| d 7} qF |	 | d	 k r	| d 7} qF t d
 |	 |   qF W| | | | | | f S)zCount number of different substitution types between two codons (PRIVATE).

    returns tuple (P0, P2, P4, Q0, Q2, Q4)

    Nomenclature is according to Li et al. (1958), PMID 3916709.
    r   r   r   r   r   r   r0   r   r   zUnexpected fold_num %d)r   r   )r   r   )r]   rj   r6   )r   r   r   ZP0ZP2ZP4ZQ0ZQ2ZQ4r   r   r   ru   r   rr   r   r   r   r   q  s>    
($$$r   c       ,   
      s  d d l  m } d d l m } d d d d d d d d i d d d d d d d d i d d d d d d d d i g } t |  } | t  } | t  }	 x |  | D] }
 |
 d k r | d |
 d d	 7<| d	 |
 d	 d	 7<| d
 |
 d
 d	 7<| |
 } xY t |  D]K \  } | d k r>| |
  d	 7<q| d k r|	 |
  d	 7<qWq Wt | j    } t |	 j    } xA t	 | |	  D]0 \  } |  | |  <|	  | |	  <qWt
 |  | d | } t | |  t |	 |  f } | | d | | d	 | | } xQ t d  D]C  t |  j      f d d   |  j   D |  <q.W| t  } x: t | j j    | j D]  d  k rd |  <qWx  |  | D]  |  d	 7<qWt |  | | d | d | \ } } } t | |  | d | d | \ } } } | | d
 } | | d
 } d d d d d d d d i d d d d d d d d i g } xK t d
  D]=  x4 d! D], } |  | |  | d
 |  | <qWqWd d g } xH t	 |  |  D]7 \  } d d   t	 | t  | d |  D } qW| d | } | d	 | }  t |  | | }! t d	 d" |   t d	 d# |  }" d% t d	 d& |!  }# d  d d g }$ xNt d  D]@}% d d   t | j j    | j D }& t | | |" |& |  }' | |' |#  }( d d d d g } d d g }) i    x_ t	 |  |  D]N \  }  d k rD| d k rD  j  | f d     | f d	 7<qDWxS   D]K  t  d  d	 |( |& |  }*    f d d   t	 | |*  D } qW| d | | d	 | f | d
 | | d | f f } g  }+ x9 t	 | |  D]( \ } }* |+ j t | |* d d  q:W|+ d d | | | |+ d	 d | | | }# |+ d	 |+ d }" t t  f d d   d d   t	 |+ |$  D   r|+ d	 |+ d f S|+ }$ qWd  S)'zsYN00 method main function (PRIVATE).

    Nomenclature is according to Yang and Nielsen (2000), PMID 10666704.
    r   )defaultdict)expmr   r   r   r   z---r0   rH   r   r   rA   r   c                s#   i  |  ] \ } } |   |  q Sr   r   )r   rr   r_   )totr   r   
<dictcomp>  s   	 z_yn00.<locals>.<dictcomp>r   r_   c             S   s    g  |  ] \ } } | |  q Sr   r   )r   rt   ru   r   r   r   rG     s   	 z_yn00.<locals>.<listcomp>g      @rv   gh㈵>   c             S   s"   g  |  ] } d  | k r |  q S)r   r   )r   r   r   r   r   rG     s   	 c                s(   g  |  ] \ } } | |     q Sr   r   )r   rt   ru   )codon_npathr   r   r   rG     s   	 tTc                s
   |    k  S)Nr   )r   )	tolerancer   r   r     s    z_yn00.<locals>.<lambda>c             S   s&   g  |  ] \ } } t  | |   q Sr   )ry   )r   r   rr   r   r   r   rG     s   	 N)r   r   r   r   gUUUUUU?gUUUUUU?g      gUUUUUU?)collectionsr   scipy.linalgr   r   r   r]   r   valuesrj   _get_TV_get_kappa_tr%   itemsr#   r>   keysr=   _count_site_YN00rx   r   _get_Q
setdefault_count_diff_YN00r<   r*   r   ),rp   rq   r_   rA   r   r   fcodonr   Z	fold0_cntZ	fold4_cntrD   r   r   Zf0_totalZf4_totalrr   TVZk04kappapirz   r{   ZbfreqSN1r|   r}   ZbfreqSN2r   r~   ZbfreqSNr   r   r   r   r1   wr   ZdSdN_pretemprK   r   r   sitestvZdSdNr   )r   r   r   r   r   rn     s    !
"+#!.	&	 )8 7rn   c             C   s  d } d	 } d d g } d } x t  |  |  D] \ } } d | | f k r. x t  | |  D] \ }	 }
 |	 |
 k rw n` |	 | k r |
 | k r | d d 7<n7 |	 | k r |
 | k r | d d 7<n | d d 7<| d 7} q\ Wq. W| d | | d | f S)
zGet TV (PRIVATE).

    Arguments:
     - T - proportions of transitional differences
     - V - proportions of transversional differences

    r   r   r   r   r   z---r0   )r   r   )r   r   )rj   )
codon_lst1
codon_lst2rA   r   r   r   r   r   r   r   rr   r   r   r   r     s     r   Fc       	      C   s  |  d |  d |  d <|  d |  d |  d <d |  d |  d |  d |  d d |  d |  d |  d |  d |  d |  d |  d |  d d | d d |  d |  d | d	 d |  d |  d |  d |  d |  d |  d } d | d d |  d |  d } d t  |  } d t  |  } | | d } | d k rd |  d |  d |  d |  d |  d |  d | |  d |  d |  d |  d } | Sd |  d |  d d | |  d d |  d |  d d | |  d d |  d |  d | } | Sd S)zmCalculate kappa (PRIVATE).

    The following formula and variable names are according to PMID: 10666704
    r   r   Yr   r   RrH   r0   r   g      ?Frv   Ng      g      )r   )	r   r   r   r   r   ar   ZkappaF84Z
kappaHKY85r   r   r   r   '  s    7"Wbr   c          	   C   s  t  |   t  |  k r= t d t  |   t  |  f   n t  |   } d } d } d } | j }	 | j }
 i  } x_ t |  |  D]N \ } } | d k r | d k r | j | | f d  | | | f d 7<q Wd } } d d d d d d d d i d d d d d d d d i g } xj| j   D]\\ } } | d } d } } xt d	  D]} x| D] } | | | k r{qb| d
 |  | | | d d
  } | |
 k rqb| | } | | | k r| | k r| | 9} n& | | | k r	| | k r	| | 9} |	 | |	 | k r@| | 7} | d | | | 7<qb| | 7} | d | | | 7<qbWqUW| | | 7} | | | 7} q(Wd	 | | | } | | 9} | | 9} x? | D]7 } t | j	    } x | D] } | | | <qWqW| | | f S)a  Site counting method from Ina / Yang and Nielsen (PRIVATE).

    Method from `Ina (1995)`_ as modified by `Yang and Nielsen (2000)`_.
    This will return the total number of synonymous and nonsynonymous sites
    and base frequencies in each category. The function is equivalent to
    the ``CountSites()`` function in ``yn00.c`` of PAML.

    .. _`Ina (1995)`: https://doi.org/10.1007/BF00167113
    .. _`Yang and Nielsen (2000)`: https://doi.org/10.1093/oxfordjournals.molbev.a026236

    z?Length of two codon_lst should be the same (%d and %d detected)r   r   r   r   z---r   r0   r   N)r   r   )r   r   )r   r   r   r   )
r!   r6   r>   r=   rj   r   r   r%   r   r   )r   r   r   r_   rA   lengthr   r   r   Z
codon_dictr   r   r   rr   r~   r   ZfreqSNZ
codon_pairnpathrD   SNposr   r   r   r   r   r   r   r   r   @  s^    		
!

&





r   c                s8  t   t  s  t   t  rD t d j t   t      t   d k sh t   d k r t d j t   t      d d d d g } d }  d k s  d k r | Sd   t   f d
 d    D  s t d j     t   f d d    D  s.t d j       k r>| Sg  } xD t t	     D]- \ } }	 |	 d |	 d k rZ| j
 |  qZWd d d  }
 t |  d k rd } d d   t	 | |
   | d |   D } nt |  d k ru| j  }   f d d   | D } g   xx | D]p } t t | j  |  g   } | | d | d f | | d | d f f }  j
 | d | d  q,W f d d    D  xt |  D] \ } }  d |   |  | d d  } d d   t	 | |
  | | | d  | d  D } d d   t	 | |
  | | | d  | d  D } qWn]t |  d k r| j  } t t d d d g d   } g   g  } x| D]}  d | d   | d  | d d d  } | d | d   | d | | d d d  } | j
 | | f  t t | j  | |  g   } | | d | d f | | d | d f | | d | d f f }  j
 | d | d | d  qW f d d    D  x t	 |  |  D] \ } } }	 d d   t	 | |
  | d |	 d | d | d  D } d d   t	 | |
 | d | d |	 d | d | d  D } d d   t	 | |
 | d  |	 d | d | d  D } qW | j k s | j k rd d g } n5 | j  | j  k r(d d g } n d d g } | S) a&  Count differences between two codons (three-letter string; PRIVATE).

    The function will weighted multiple pathways from codon1 to codon2
    according to P matrix of codon substitution. The proportion
    of transition and transversion (TV) will also be calculated in
    the function.
    zM_count_diff_YN00 accepts string object to represent codon ({0}, {1} detected)r   z7codon should be three letter string ({0}, {1} detected)r   z---r   r   r   r   c             3   s   |  ] } |   k Vq d  S)Nr   )r   r   )r   r   r   r     s    z#_count_diff_YN00.<locals>.<genexpr>zNUnrecognized character detected in codon1 {0} (Codons consist of A, T, C or G)c             3   s   |  ] } |   k Vq d  S)Nr   )r   r   )r   r   r   r     s    zNUnrecognized character detected in codon2 {0} (Codons consist of A, T, C or G)r0   c       	      S   s  d } d } | j  } | j } |  | k s6 | | k r |  | | k rf | | | k rf d d | d g S|  | | k r | | | k r d d | d g Sd d d | g Sn | |  | | k r0|  | | k r | | | k r | d d d g S|  | | k r| | | k r| d d d g Sd | d d g Snp |  | | k r`| | | k r`d d | d g S|  | | k r| | | k rd d | d g Sd d d | g Sd  S)Nr   r   r   r   r   )r   r   )r   r   )r>   r=   )	r   r   diffrA   r   r   r   Zdicr   r   r   r   count_TV  s*    		      z"_count_diff_YN00.<locals>.count_TVc             S   s    g  |  ] \ } } | |  q Sr   r   )r   r1   qr   r   r   rG     s   	 z$_count_diff_YN00.<locals>.<listcomp>rH   c                s:   g  |  ]0 }   d  |   |   | d d    q S)Nr0   r   )r   r   )r   r   r   r   rG     s   	c                s$   g  |  ] } d  | t      q S)rH   )r   )r   r   )	path_probr   r   rG     s   	 Nc             S   s    g  |  ] \ } } | |  q Sr   r   )r   r1   r   r   r   r   rG     s   	 r   c             S   s    g  |  ] \ } } | |  q Sr   r   )r   r1   r   r   r   r   rG     s   	 c                s$   g  |  ] } d  | t      q S)r   )r   )r   r   )r   r   r   rG     s   	 c             S   s    g  |  ] \ } } | |  q Sr   r   )r   r1   r   r   r   r   rG     s   	 c             S   s    g  |  ] \ } } | |  q Sr   r   )r   r1   r   r   r   r   rG     s   	 c             S   s    g  |  ] \ } } | |  q Sr   r   )r   r1   r   r   r   r   rG     s   	 )r   r   r   r   )r   rR   r   r(   r   r!   r6   r*   r]   rj   r<   r>   r#   r   r-   r   r=   )r   r   r   rK   rA   r   siter   r   r_   r   Zprobr   r   Z	codon_idxru   r   r   r1   r   r   rr   r   )r   r   r   r   r   r     s    	 	$			"2
!*66$%"#'#r   c             C   s2  d d l  m } d d l m } |   } t |  | | d | } xC t |  |  D]2 \ } }	 d | |	 f k rQ | | |	 f d 7<qQ Wd d   t | j j    | j	 D }
 | | |
 | d	 d
  } | | d d d g d d d d d d } | j
 \ } } } t | | | |
 |  } d } } x t |
  D] \ } } x t |
  D] \ }	 } | |	 k rLyY | j | | j | k r| | | | | |	 f 7} n | | | | | |	 f 7} WqLt k
 rYqLXqLWq3W| | 9} | | 9} | | d d d g d d d d d d } | j
 \ } } } t | | | |
 |  } d } } x t |
  D] \ } } x t |
  D] \ }	 } | |	 k rryY | j | | j | k r| | | | | |	 f 7} n | | | | | |	 f 7} Wqrt k
 rYqrXqrWqYW| d 9} | d 9} | | } | | } | | f S)z"ML method main function (PRIVATE).r   )Counter)minimizerA   z---r0   c             S   s"   g  |  ] } d  | k r |  q S)r   r   )r   r   r   r   r   rG   
  s   	 z_ml.<locals>.<listcomp>c          
   S   s/   t  |  d |  d |  d | | d | d | S)z'Temporary function, params = [t, k, w].r   r0   rH   rK   rA   )_likelihood_func)paramsr   	codon_cntrK   rA   r   r   r   func  s    z_ml.<locals>.funcg?rH   ro   zL-BFGS-BZbounds绽|=r   
   Ztolgh㈵>r   r   r   r   r   r   r   )r   r   r   r   r   r   r   r0   r0   )r   r   r   )r   r   Zscipy.optimizer   _get_pirj   r#   r>   r   r=   r   r   r]   r?   )rp   rq   cmethodrA   r   r   r   r   r   rr   rK   r   Zopt_resr   r_   r   r   ZSdZNdc1c2ZrhoSZrhoNr   r   r   r   r   rk      sd    		 		

 

	

 



rk   c          
      s(  i  } | d k r d d d d d d d d i } x= |  | D]1 } | d k r; x | D] } | | d 7<qT Wq; Wt  | j        f d	 d
   | j   D } x| j j   | j D]< } d | k r | | d | | d | | d | | <q Wn)| d k rd d d d d d d d i d d d d d d d d i d d d d d d d d i g } x` |  | D]T } | d k ri| d | d d 7<| d | d d 7<| d | d d 7<qiWxQ t d  D]C } t  | | j        f d d
   | | j   D | | <qWxt | j j    | j D]H } d | k r2| d | d | d | d | d | d | | <q2Wn | d k r$x4 | j j   | j D] } d | k rd | | <qWx, |  | D]  } | d k r| | d 7<qWt  | j        f d d
   | j   D } | S)zObtain codon frequency dict (pi) from two codon list (PRIVATE).

    This function is designed for ML method. Available counting methods
    (cfreq) are F1x4, F3x4 and F64.
    rf   r   r   r   r   r   z---r0   c                s#   i  |  ] \ } } |   |  q Sr   r   )r   rr   r_   )r   r   r   r   [  s   	 z_get_pi.<locals>.<dictcomp>r   rH   re   r   c                s#   i  |  ] \ } } |   |  q Sr   r   )r   rr   r_   )r   r   r   r   k  s   	 rg   g?c                s#   i  |  ] \ } } |   |  q Sr   r   )r   rr   r_   )r   r   r   r   x  s   	 )r   r   r   r>   r   r=   r%   r#   )rp   rq   r   rA   r   r   r   cr   )r   r   r   J  sL    	1!+=r   c             C   s  |  | k r d S|  | j  k s. | | j  k r2 d S|  | k sJ | | k rN d Sd	 } d
 } g  } xK t t |  |   D]4 \ }	 \ }
 } |
 | k rv | j |	 |
 | f  qv Wt |  d k r d S| j |  | j | k rQ| d d | k r| d d | k r| | | S| d d | k rF| d d | k rF| | | S| | Sn| | d d | k r| d d | k r| | | | S| d d | k r| d d | k r| | | | S| | | Sd S)a   Q matrix for codon substitution (PRIVATE).

    Arguments:
     - i, j  : three letter codon string
     - pi    : expected codon frequency
     - k     : transition/transversion ratio
     - w     : nonsynonymous/synonymous rate ratio
     - codon_table: Bio.Data.CodonTable object

    r   r   r   r   r   rH   r0   N)r   r   )r   r   )r=   r]   rj   r<   r!   r>   )r   rr   r   r_   r   rA   r   r   r   ru   r   r   r   r   r   _q|  s2    (((((r   c             C   s#  d d l  } t |  } | j | | f  } xg t |  D]Y } xP t |  D]B }	 | |	 k rM t | | | |	 |  | | d | | | |	 f <qM Wq: Wd }
 xu t |  D]g } t | | d d  f  | | | f <y% |
 |  | | | | | f 7}
 Wq t k
 rYq Xq W| |
 } | S)z*Q matrix for codon substitution (PRIVATE).r   NrA   )Znumpyr!   Zzerosr%   r   r   r?   )r   r_   r   rK   rA   nprM   r   r   rr   Znucl_substitutionsr   r   r   r     s"    '%	
r   c          	   C   s   d d l  m } t | | | | |  } | | |   }	 d }
 x t |  D] \ } } x t |  D] \ } } | | f | k rd |	 | | f | | d k r |
 | | | f d 7}
 qd |
 | | | f t | | |	 | | f  7}
 qd WqK W|
 S)z,Likelihood function for ML method (PRIVATE).r   )r   )r   r   r   r]   r   )r   r_   r   r   r   rK   rA   r   r   r   Z
likelihoodr   r   rr   r   r   r   r   r     s    8r   __main__)run_doctest))r[   
__future__r   r   	itertoolsr   mathr   ZBio.Seqr   ZBio.SeqRecordr   ZBio.Alphabetr   r	   ZBio.codonalign.codonalphabetr
   r   r   r   rc   rs   rl   rw   rx   rm   r   r   rn   r   r   r   r   rk   r   r   r   r   rX   Z
Bio._utilsr   r   r   r   r   <module>   s<   D>X*#/kA~J22