Reproduce codeml analysis from Yang Z (MBE, 2005)¶

This cookbook aims at reproducing the results about selection on sites in a given lineage presented in the work Yang Z (MBE, 2005) using ete-evol.

The data files adapted from the examples in the PAML package, large dataset will be used.

This files are:

data/lysozyme/lysozymeLarge.fasta
data/lysozyme/lysozymeLarge.nw

In [11]:

%%bash
ete3 view -t data/lysozyme/lysozymeLarge.nw  --alg data/lysozyme/lysozymeLarge.fasta --alg_type compactseq \
-i data/lysozyme/lysozymeLarge.png

In [12]:

from IPython.display import Image
Image(filename='data/lysozyme/lysozymeLarge.png')

Out[12]:

In particular we aim at reproducing the highlighted part the Table2:

In [9]:

from IPython.display import Image
Image(filename='data/lysozyme/table5_YangZ_2005.png')

Out[9]:

We are thus going to run 6 models.

First we run models M1 and M2:

In [6]:

%%bash
ete3 evol -t data/lysozyme/lysozymeLarge.nw --alg data/lysozyme/lysozymeLarge.fasta --models M1 bsA bsA1  -o /tmp/ete3-tmp/ \
--mark "colobus_Cgu&Can,,probiscis_Nla" --cpu 3 --resume -i data/lysozyme/lysozymeLarge_bsA_bsA1_M1.png

Using: /Users/fransua/.miniconda2/bin/ete3_apps/bin/Slr
Using: /Users/fransua/.miniconda2/bin/ete3_apps/bin/codeml

Running CodeML/Slr (3 CPUs)
  - processing model M1
Model M1 already executed... SKIPPING
  - processing model bsA
       marking branches 22

Model bsA.22-1 already executed... SKIPPING
  - processing model bsA1
       marking branches 22

Model bsA1.22-1 already executed... SKIPPING

LRT

               Null model |        Alternative model | p-value
     ------------------------------------------------------------
                bsA1.22-1 |                 bsA.22-1 | 0.224738
                       M1 |                bsA1.22-1 | 0.170264
                       M1 |                 bsA.22-1 | 0.186889

SUMMARY BY MODEL

 - Model bsA1.22-1
   * Marked branches
      (((((langur_Sen&Sve,langur_Tob&Tfr,Douc_langur_Pne),probiscis_Nla),colobus_Cgu&Can) #1,(((baboon_Pcy,mangabey_Cat),rhesus_Mmu),Allen_Ani,(talapoin_Mta,(patas_Epa,vervet_Cae)))),(((human,chimp_bonobo_gorilla),orangutan_Ppy),gibbon_Ggo),(squirrel_m,(tamarin_Soe,Marmoset_Cja)));

 - Model M1
   * Average omega for all tree: 0.495
   * Sites significantly caracterized
      codon position |   category
     -----------------------------------------------------------
                 2   |   Relaxed (probability > 0.95)
                14   |   Relaxed (probability > 0.99)
                15   |   Relaxed (probability > 0.99)
                17   |   Relaxed (probability > 0.95)
                21   |   Relaxed (probability > 0.95)
                23   |   Relaxed (probability > 0.95)
                29   |   Relaxed (probability > 0.99)
                37   |   Relaxed (probability > 0.99)
                41   |   Relaxed (probability > 0.95)
                47   |   Relaxed (probability > 0.99)
                50   |   Relaxed (probability > 0.99)
                62   |   Relaxed (probability > 0.95)
                67   |   Relaxed (probability > 0.95)
                69   |   Relaxed (probability > 0.95)
                75   |   Relaxed (probability > 0.95)
           78-  79   |   Relaxed (probability > 0.99)
                82   |   Relaxed (probability > 0.95)
           87-  88   |   Relaxed (probability > 0.95)
                90   |   Relaxed (probability > 0.99)
                91   |   Relaxed (probability > 0.95)
                94   |   Relaxed (probability > 0.99)
               101   |   Relaxed (probability > 0.95)
               106   |   Relaxed (probability > 0.99)
          113- 114   |   Relaxed (probability > 0.95)
               115   |   Relaxed (probability > 0.99)
               117   |   Relaxed (probability > 0.95)
               119   |   Relaxed (probability > 0.95)
          121- 122   |   Relaxed (probability > 0.99)
          125- 126   |   Relaxed (probability > 0.99)

 - Model bsA.22-1
   * Marked branches
      (((((langur_Sen&Sve,langur_Tob&Tfr,Douc_langur_Pne),probiscis_Nla),colobus_Cgu&Can) #1,(((baboon_Pcy,mangabey_Cat),rhesus_Mmu),Allen_Ani,(talapoin_Mta,(patas_Epa,vervet_Cae)))),(((human,chimp_bonobo_gorilla),orangutan_Ppy),gibbon_Ggo),(squirrel_m,(tamarin_Soe,Marmoset_Cja)));

In the output above we can see in the LRT the results of the so called "test 1" (bsA vs M1) and "test 2" (bsA vs bsA1), and in both cases, the results are in favor of the null hypothesis just as persented in the article.

The image generated by the ete_evol tool:

In [7]:

from IPython.display import Image
Image(filename='data/lysozyme/lysozymeLarge_bsA_bsA1_M1.png')

Out[7]:

To search for lineages with specific sites evolving under positive selection we could explore the selective pressure on branches using the free-ration (or free-branch) model:

In [12]:

%%bash
ete3 evol -t data/lysozyme/lysozymeLarge.nw --alg data/lysozyme/lysozymeLarge.fasta --models fb  -o /tmp/ete3-tmp/ \
--cpu 1 --resume -i data/lysozyme/lysozymeLarge_fb.png

Using: /Users/fransua/.miniconda2/bin/ete3_apps/bin/Slr
Using: /Users/fransua/.miniconda2/bin/ete3_apps/bin/codeml

Running CodeML/Slr (1 CPUs)
  - processing model fb
Model fb already executed... SKIPPING
SUMMARY BY MODEL

 - Model fb
   * Marked branches
      (((((langur_Sen&Sve #27,langur_Tob&Tfr #28,Douc_langur_Pne #29) #19,probiscis_Nla #20) #10,colobus_Cgu&Can #11) #4,(((baboon_Pcy #30,mangabey_Cat #31) #21,rhesus_Mmu #22) #12,Allen_Ani #13,(talapoin_Mta #23,(patas_Epa #32,vervet_Cae #33) #24) #14) #5) #1,(((human #25,chimp_bonobo_gorilla #26) #15,orangutan_Ppy #16) #6,gibbon_Ggo #7) #2,(squirrel_m #8,(tamarin_Soe #17,Marmoset_Cja #18) #9) #3);

        Branches  =>   omega
              #4  =>   3.516
             #29  =>   0.185
             #12  =>   0.000
             #19  => 999.000
             #27  =>   0.000
             #30  =>   0.000
             #14  => 999.000
             #16  =>   0.000
             #20  =>   0.000
             #18  =>   0.078
             #24  => 999.000
             #28  =>   0.125
             #15  => 999.000
             #33  =>   0.000
             #11  =>   1.370
             #21  =>   0.000
             #25  =>   0.000
             #10  =>   0.856
              #8  =>   0.632
              #5  =>   0.432
             #32  =>   0.000
             #17  =>   0.137
              #1  =>   1.514
              #6  => 999.000
             #22  => 999.000
              #3  =>   0.386
             #13  => 999.000
              #2  => 999.000
             #26  =>   0.000
              #7  =>   0.379
              #9  =>   2.133
             #23  =>   0.366
             #31  =>   0.000

In [11]:

from IPython.display import Image
Image(filename='data/lysozyme/lysozymeLarge_fb.png')

Out[11]:

Many lineage seems to have evolved at high $\omega$ rates. We will use the hominide ancestor.

In [13]:

%%bash
ete3 evol -t data/lysozyme/lysozymeLarge.nw --alg data/lysozyme/lysozymeLarge.fasta --models M1 bsA bsA1  -o /tmp/ete3-tmp/ \
--mark orangutan_Ppy,,gibbon_Ggo --cpu 3 --resume -i data/lysozyme/lysozymeLarge_hominidae_bsA_bsA1_M1.png

Using: /Users/fransua/.miniconda2/bin/ete3_apps/bin/Slr
Using: /Users/fransua/.miniconda2/bin/ete3_apps/bin/codeml

Running CodeML/Slr (3 CPUs)
  - processing model M1
Model M1 already executed... SKIPPING
  - processing model bsA
       marking branches 30

Model bsA.30-1 already executed... SKIPPING
  - processing model bsA1
       marking branches 30

Model bsA1.30-1 already executed... SKIPPING

LRT

               Null model |        Alternative model | p-value
     ------------------------------------------------------------
                bsA1.30-1 |                 bsA.30-1 | 0.008960**
                       M1 |                bsA1.30-1 | 0.092070
                       M1 |                 bsA.30-1 | 0.007952**

SUMMARY BY MODEL

 - Model bsA1.30-1
   * Marked branches
      (((((langur_Sen&Sve,langur_Tob&Tfr,Douc_langur_Pne),probiscis_Nla),colobus_Cgu&Can),(((baboon_Pcy,mangabey_Cat),rhesus_Mmu),Allen_Ani,(talapoin_Mta,(patas_Epa,vervet_Cae)))),(((human,chimp_bonobo_gorilla),orangutan_Ppy),gibbon_Ggo) #1,(squirrel_m,(tamarin_Soe,Marmoset_Cja)));
   * Sites significantly caracterized
      codon position |   category
     -----------------------------------------------------------
                67   |   Relaxed (probability > 0.95)
                79   |   Relaxed (probability > 0.95)
               115   |   Relaxed (probability > 0.95)
               122   |   Relaxed (probability > 0.95)

 - Model M1
   * Average omega for all tree: 0.495
   * Sites significantly caracterized
      codon position |   category
     -----------------------------------------------------------
                 2   |   Relaxed (probability > 0.95)
                14   |   Relaxed (probability > 0.99)
                15   |   Relaxed (probability > 0.99)
                17   |   Relaxed (probability > 0.95)
                21   |   Relaxed (probability > 0.95)
                23   |   Relaxed (probability > 0.95)
                29   |   Relaxed (probability > 0.99)
                37   |   Relaxed (probability > 0.99)
                41   |   Relaxed (probability > 0.95)
                47   |   Relaxed (probability > 0.99)
                50   |   Relaxed (probability > 0.99)
                62   |   Relaxed (probability > 0.95)
                67   |   Relaxed (probability > 0.95)
                69   |   Relaxed (probability > 0.95)
                75   |   Relaxed (probability > 0.95)
           78-  79   |   Relaxed (probability > 0.99)
                82   |   Relaxed (probability > 0.95)
           87-  88   |   Relaxed (probability > 0.95)
                90   |   Relaxed (probability > 0.99)
                91   |   Relaxed (probability > 0.95)
                94   |   Relaxed (probability > 0.99)
               101   |   Relaxed (probability > 0.95)
               106   |   Relaxed (probability > 0.99)
          113- 114   |   Relaxed (probability > 0.95)
               115   |   Relaxed (probability > 0.99)
               117   |   Relaxed (probability > 0.95)
               119   |   Relaxed (probability > 0.95)
          121- 122   |   Relaxed (probability > 0.99)
          125- 126   |   Relaxed (probability > 0.99)

 - Model bsA.30-1
   * Marked branches
      (((((langur_Sen&Sve,langur_Tob&Tfr,Douc_langur_Pne),probiscis_Nla),colobus_Cgu&Can),(((baboon_Pcy,mangabey_Cat),rhesus_Mmu),Allen_Ani,(talapoin_Mta,(patas_Epa,vervet_Cae)))),(((human,chimp_bonobo_gorilla),orangutan_Ppy),gibbon_Ggo) #1,(squirrel_m,(tamarin_Soe,Marmoset_Cja)));
   * Sites significantly caracterized
      codon position |   category
     -----------------------------------------------------------
                79   |   Positively-selected (probability > 0.95)
               122   |   Positively-selected (probability > 0.95)

In [14]:

from IPython.display import Image
Image(filename='data/lysozyme/lysozymeLarge_hominidae_bsA_bsA1_M1.png')

Out[14]:

This time we do find a couple of sites under positive section in the hominide ancestors.