wget https://molevol.mbl.edu/images/3/3c/MSAlab.zip

*Most of this text is from the MBL tutorial for MSA. I am recording what I find and important points, in case the site gets lost over time.

Alignments with nucleotides: MUSCLE V. mafft

mafft --retree 2 1ped.fasta > mafft_dna.fasta

muscle -verbose -log muscle_dna.log -in 1ped.fasta -out muscle_dna.fasta

Compare the alignments resulting from A steps and B steps. Are they different? How many columns are in each the MAFFT or the MUSCLE alignment? What may be wrong with both? (Hint: these are protein coding genes).

They were not well-aligned, there were indels that were not multiples of three and they are protein coding genes!

Do the topologies and/or branch lengths differ? (Hint: look up some species names to get an idea of the expected topology!) Yes--very much so.

Bottom line: Don't do nucleotide alignments on protein-based sequences!

Comparing what happens with gaps:

muscle -verbose -log muscle_gap-20.log -in 1ped_aa.fasta -out muscle_aa_gap-20.fasta -gapopen -20

muscle -verbose -log muscle_gap-1.log -in 1ped_aa.fasta -out muscle_aa_gap-1.fasta -gapopen -1

Compare the modified gap penalty alignments to the default one. How do they differ? Which has the most gaps? Can you guess the order of magnitude of the default gap penalty? Can you find in the log files of MUSCLE the gap penalty use?

The higher gap penalty makes for a much more compressed alignment and helps to take care of this over "indel" issue. The default for MUSCLE is -2.9.

Alignment extension

   papara -t example.tre -s aln.phy -q sequence1.fasta -n run1

   papara -t example.tre -s aln.phy -q sequence2.fasta -n run2

How does the alignment look? Any concerns? How does the requirement for testing for sequence homology before site alignment apply here?

There are big tails at the ends of the new sequence. The alignment may have been trimmed previously, or the previous group used a different set of primers. The middle part is probably homologous, but a lot of times the ends may not line up perfectly. Always check your alignment.

   pagan --ref-seqfile aln.fasta -t example.tre -q sequence1.fasta -o pagan1

  pagan --ref-seqfile aln.fasta -t example.tre -q sequence2.fasta -o pagan2

Pagan requires input alignments in fasta rather than phylip format, but we have provided that in aln.fasta Try typing only 'pagan' to see the full explanation of parameters Compare the pagan alignments for the two sequences you downloaded.

What happens when you try to align the second sequence? HINT: Look at the output to screen. Why?

Adding the second sequence fails. Pagan is helping us to prevent alignments of completely random sequence! Thanks Pagan!

Both papara and pagan can align short reads as well as full length sequences, but caution about homology is particularly appropriate in these cases.

Codon Alignments

This is not part of the exercises, it's just for your future information.

As you now know, it is not appropriate to align nucleotides of protein coding regions. In the exercises above, you translated the nucleotides to amino acids which you could use to infer trees. But sometimes you want to analyze nucleotides that have been aligned by codon. So how do you go from an amino acid alignment back to codon-aligned nucleotides? We suggest the Pal2Nal server. You will upload your protein alignment and nucleotide sequences, and it will spit out the codon alignment. Please be aware that your nucleotides must be multiples of three (i.e. a full open reading frame).

FASTA tutorial

Identifying homologs and non-homologs; effects of scoring matrices and algorithms; using domain annotations

1. Use the FASTA search page [pgm] to compare Honey bee glutathione transferase D1 NP_001171499/ H9KLY5_APIME [seq] (gi|295842263) to the PIR1 Annotated protein sequence database. Be sure to press, not .

1. Take a look at the output.
   1. How long is the query sequence?
   1. 217 amino acids
   2. How many sequences are in the PIR1 database?
   1. 13,143 sequences
   3. What scoring matrix was used?
   1. BL50 matrix
   4. What were the gap penalties? (what is the penalty for a one-residue gap? two residues?)
   1. one residue: -12 open gap, the -14 for two residues
   5. What are each of the numbers after the description of the library sequence? Which one is best for inferring homology?

   1. sp|P20432|GSTT1_DROME Glutathione S-transferase 1-1; DD ( 209)

   2. 209 is referring to the length of sequence that it matches

   6. How similar is the highest scoring sequence? What is the difference between %_id and %_sim? Why is there no 100% identity match?
   1. 81% similar,  but 57% identity
   2. The sequence is not in the database
   7. Looking at an alignment, where are the boundaries of the alignment (the best local region)? How many gaps are in the best alignment? The second best?
   1. Local alignment, score gets worse if you add more residue so important to consider
   2. There are no gaps in the first, there are 14 in the second
2. What is the highest scoring non-homolog? (The non-homolog with the highest alignment score, or the lowest E()-value.) How would you confirm that your candidate non-homolog was truly unrelated? (Hint - compare your candidate non-homolog with SwissProt for a more comprehensive test.)
1. Highest scoring unrelated sequence should have an E() value of 1 (or around 1). Can do "general re-search of this sequence and make sure it does not find your originally queried value. 
3. Homology (and non-homology) can also be inferred by domain relationships. Try the same Honey bee GSTD1 search [pgm] search using the Annotate Query and Database Annotations: set to show InterPro domains/no features.
   1. Does the domain display change your mind about the highest scoring non-homolog?
   1. Results with E() of 0.11 show completely different shared domains and some with an even higher E() show a shared domain. Often say NODOM Q-values <30 have no domain.
   2. There are three parts to the domain display, the domain structure of the query (top) sequence (if available), the domain structure of the library (bottom) sequence, and the domain alignment boundaries in the middle (inside the alignment box). The boundaries and color of the alignment domain coloring match the Region: sub-alignment scores.
   
   1. Based on the score, the second part of the domain is not significant. Sometimes the homolog is "too far away" (evolutionary distance) to be given a high score.  
   2. Domains have characteristic length (about 100 residues). Look to see if there is room for a domain, if it is not contributing to the score, it could just be too far away.
   3. Use transitivity from other sequences to make decision it is there.
   3. Note that the alignment of Honey bee GSTD1 and SSPA_ECO57 includes portions of both the N-terminal and C-terminal domains, but neither domain is completely aligned. Why do you think the alignments do not include the complete domains?
   4. Is your explanation for the partial domain alignment consistent the the argument that domains have a characteristic length? How might you test whether a complete domain is present?
   5. The FASTA programs can partition an alignment score based domain boundaries, and report the amount of score associated with a domain. In the SSPA_ECO57 alignment, how much of the score comes from the GST-N terminal domain? The GST-C terminal domain? Does this alignment provide strong support for the presence of a GST-C terminal domain? How might you test for the presence of the domain on GSTD1? On SSPA_ECO57?In the subalignment scores, the Q value is -10 * log(p) for the sub-alignment score, so Q=30.0 means p < 0.001.
4. Repeat the GSTD1 search [pgm] using the BLASTP62/-11/-1 scoring matrix BlastP62 (30%) that BLAST uses. Re-examine the GSTD1:SSPA_ECO57 alignment. Are both Glutathione transferase domains present? Look at the alignments to the homologs above and below SSPA_ECO57. Based on those aligments, do you think the Glutathione-S-Trfase C-like domain is really missing? Why did the alignment become so much shorter?
1. SKIP
5. One of the candidate non-homologs is sp|Q9SI20|EF1D2_ARATH, with an E()-value of 0.11.
1. Does the domain structure of EF1D2_ARATH suggest that it could be a glutatione transferase homolog?
1. No, it looks very different
   
2. Only residues 8-60 overlapped
2. Use the General Research to explore the domains contained in EF1D2_ARATH homologs found in SwissProt.
3. Does this secondary search support homology or non-homology?


1. Focusing on the NODOM, it actually just isn't annotated! Scorlling down you can see that a glutathion homolog eventaully does show up in the annotation!!
2. 
   

Trust annotations when they say something, but don't trust them if they don't say something. 

2. Exploring domains and alignment over-extension -- cortactin (SRC8_HUMAN)

Compare SRC8_HUMAN [pgm] (human cortactin) to the SwissProt protein sequence database.

1. Looking at the top five alignments, how many cortactin orthologs do you see? (ortholog, same protein, different species).
1. Four orthologs
2. In the SRC8 HUMAN:CHICK alignment, both the query and the subject (library) sequences align seven cortactin domains and an SH3 domain. In addition, two regions (one before the cortactin domain cluster and one after) are well conserved, but do not have annotated domains (NODOM). Are these non-domain (NODOM) regions as well conserved as the annotated domains?

1. 
   
3. Look at the SRC8_HUMAN:HCLS1_MOUSE alignment. How many cortactin domains does HCLS1_MOUSE contain? How much score does the NODOM spanning the region between cortactin domains and the SH3 domain contribute? Why is it included in the alignment? Is it likely to be homologous?
1. It is a homolog, but not an ortholog--it is a paralog! We already found mouse domain; only 4 shared domain too. 
4. Is the NODOM between the cortactin domains and the SH3 domain likely to be homologous in the SRC8_HUMAN:DBNLB_XENLA alignment?

1. No. The domain score of the second, non annotated domain has a poor score (even negative!). Still needs to be included to include the SH3 domain. No SH3 domain, the alignment would stop.
   
5. In the SRC8_HUMAN:LASP1_HUMAN alignment, the alignment extends to include several Nebulin_repeat domains. Do you think there is a Nebulin_repeat domain in SRC8_HUMAN? Why do you think those domains are aligning?
6. What scoring matrix should be used to reduce over-extension from the SH3 domain?

3. Exploring domains and over-extension with local alignments -- death associated protein kinase (DAPK1_HUMAN)

1. Look up the domain structure of DAPK1_HUMAN at Pfam [pgm].
   1. What are the major (PfamA) domain regions on the protein?
   2. Which of the domains is repeated?
   3. In a local (LALIGN) alignment, where would you expect to see overlapping domains like those in Calmodulin (CALM_HUMAN) and Cortactin (DAPK1_HUMAN)?
2. Use lalign/plalign [pgm] to examine local similarities between DAPK1_HUMAN and itself. Check the options to "annotate sequence 1 domains" and "annotate sequence 2 domains". Do you see the domains you expected from Pfam? Do they map in the same places?

After running GARLI, I wanted to compare
OR137:
OR213: lots of disagreements; redo with 100
OR4: lots of disagreements; redo with 100
OR589: lot
OR6: all agree after 8 runs
OR10: all agree after 8 runs
OR11: tree2 disagrees, 7 others agree, use this as main gene tree
OR51: lots of disagreements; redo with 100
OR52: lots of disagreements; redo with 100

100 is a lot of tree to compare, so I wrote a small script doing so:

setwd("/Volumes/Yango/Hayden_batORs/garli/OR51/trees")

#tell where to open up all of your downloaded files from CIPRES
my.files<-list.files("/Volumes/Yango/Hayden_batORs/garli/OR51/trees/")
#read in all of these trees to a list using read.nexus function
my.trees<-lapply(my.files, read.nexus)

#create a matrix to store all of the tree comparisons
tree.matrix<-matrix(nrow = length(my.trees), ncol = length(my.trees))

for(i in 1:length(my.trees)){ #for every tree in the file
  this.tree<-as.phylo(my.trees[[i]]) #look at that tree
  #compare the topology to all the other trees
  tree.matrix[i,]<-sapply(my.trees, function(my.trees) dist.topo(this.tree,my.trees))
}
tree.matrix

#print out the ones that agree with each other
for(i in 1:length(tree.matrix[1,])){
  this.list<-as.numeric(which(tree.matrix[i,]==0))
  print(this.list)
}
write.csv(tree.matrix, "tree.matrix.csv")

For OR51, are things looking better?

[1] 1

[1] 2

[1] 3

[1]  4 84

[1] 5

[1] 6

[1] 7

[1]  8 15 41 79

[1]  9 14

[1] 10 54

[1] 11

[1] 12

[1] 13 73

[1]  9 14

[1]  8 15 41 79

[1] 16

[1] 17 20 75

[1] 18

[1] 19

[1] 17 20 75

[1] 21 51 68

[1] 22 37

[1] 23 50 71

[1] 24

[1] 25

[1] 26

[1] 27

[1] 28 69 88

[1] 29

[1] 30

[1] 31

[1] 32 46

[1] 33

[1] 34 70

[1] 35 53

[1] 36

[1] 22 37

[1] 38 59

[1] 39 47 78

[1] 40

[1]  8 15 41 79

[1] 42

[1] 43

[1] 44

[1] 45

[1] 32 46

[1] 39 47 78

[1] 48

[1] 49

[1] 23 50 71

[1] 21 51 68

[1] 52

[1] 35 53

[1] 10 54

[1] 55

[1] 56

[1] 57

[1] 58

[1] 38 59

[1] 60

[1] 61

[1] 62 86 87

[1] 63

[1] 64

[1] 65

[1] 66 95

[1] 67

[1] 21 51 68

[1] 28 69 88

[1] 34 70

[1] 23 50 71

[1] 72

[1] 13 73

[1] 74

[1] 17 20 75

[1] 76

[1] 77

[1] 39 47 78

[1]  8 15 41 79

[1] 80

[1] 81

[1] 82

[1] 83

[1]  4 84

[1] 85

[1] 62 86 87

[1] 62 86 87

[1] 28 69 88

[1] 89

[1] 90

[1] 91

[1] 92

[1] 93

[1] 94

[1] 66 95

[1] 96

[1] 97

[1] 98

Womp, womp...not really--only a maximum of 4 trees out of 100 are agreeing.

OR4s looking much better:
[1]  1 22 28 42 75
[1]  2 99
[1]  3 14 18 83 89
 [1]  4 12 17 20 23 24 27 51 58 59 72 77 79 88 96
[1] 5
[1] 6
[1]  7 15 98
[1] 8
[1]  9 19 38 80
[1] 10 50
[1] 11 67
 [1]  4 12 17 20 23 24 27 51 58 59 72 77 79 88 96
[1] 13 86
[1]  3 14 18 83 89
[1]  7 15 98
[1] 16
 [1]  4 12 17 20 23 24 27 51 58 59 72 77 79 88 96
[1]  3 14 18 83 89
[1]  9 19 38 80
 [1]  4 12 17 20 23 24 27 51 58 59 72 77 79 88 96
[1] 21 31 44 62 84
[1]  1 22 28 42 75
 [1]  4 12 17 20 23 24 27 51 58 59 72 77 79 88 96
 [1]  4 12 17 20 23 24 27 51 58 59 72 77 79 88 96
[1] 25
[1] 26
 [1]  4 12 17 20 23 24 27 51 58 59 72 77 79 88 96
[1]  1 22 28 42 75
[1] 29
[1] 30 33
[1] 21 31 44 62 84
[1] 32 37 61 68 81
[1] 30 33
[1]  34 100
[1] 35
[1] 36
[1] 32 37 61 68 81
[1]  9 19 38 80
[1] 39 71
[1] 40 49
[1] 41
[1]  1 22 28 42 75
[1] 43
[1] 21 31 44 62 84
[1] 45
[1] 46
[1] 47
[1] 48
[1] 40 49
[1] 10 50
 [1]  4 12 17 20 23 24 27 51 58 59 72 77 79 88 96
[1] 52
[1] 53
[1] 54
[1] 55
[1] 56 90
[1] 57
 [1]  4 12 17 20 23 24 27 51 58 59 72 77 79 88 96
 [1]  4 12 17 20 23 24 27 51 58 59 72 77 79 88 96
[1] 60
[1] 32 37 61 68 81
[1] 21 31 44 62 84
[1] 63
[1] 64
[1] 65
[1] 66
[1] 11 67
[1] 32 37 61 68 81
[1] 69
[1] 70
[1] 39 71
 [1]  4 12 17 20 23 24 27 51 58 59 72 77 79 88 96
[1] 73
[1] 74
[1]  1 22 28 42 75
[1] 76
 [1]  4 12 17 20 23 24 27 51 58 59 72 77 79 88 96
[1] 78
 [1]  4 12 17 20 23 24 27 51 58 59 72 77 79 88 96
[1]  9 19 38 80
[1] 32 37 61 68 81
[1] 82
[1]  3 14 18 83 89
[1] 21 31 44 62 84
[1] 85
[1] 13 86
[1] 87
 [1]  4 12 17 20 23 24 27 51 58 59 72 77 79 88 96
[1]  3 14 18 83 89
[1] 56 90
[1] 91
[1] 92
[1] 93
[1] 94
[1] 95
 [1]  4 12 17 20 23 24 27 51 58 59 72 77 79 88 96
[1] 97
[1]  7 15 98
[1]  2 99
[1]  34 100

PAML! Oh, good fun.

The worst part about PAML is getting ready to run PAML.

Be sure to keep organized. To run PAML, navigate to /Applications/paml4.8/bin/. There should be all the different PAML binaries listed.


For the null model:

bash-3.2$ sudo ./codeml codeml.ctl 

Tree file with two groups: group 0 is taxa with 

(((((((((((((((((((((Artibeus_fimbriatus,(Artibeus_jamaicensis,Artibeus_obscurus)),Artibeus_inopinatus),Artibeus_concolor),((((Artibeus_watsoni,Artibeus_incomitatus),Artibeus_phaeotis),Artibeus_cinereus),(Artibeus_glaucus,Artibeus_gnomus))),((Phyllops_falcatus,(Ariteus_flavescens,Ardops_nichollsi)),((Pygoderma_bilabiatum,Sphaeronycteris_toxophyllum),Centurio_senex))),Ectophylla_alba),Enchisthenes_hartii),((((((((Platyrrhinus_nigellus,Platyrrhinus_dorsalis),(Platyrrhinus_infuscus,Platyrrhinus_aurarius)),Platyrrhinus_albericoi),Platyrrhinus_lineatus),Vampyrodes_caraccioli),(Vampyressa_thyone,Mesophylla_macconnelli)),(((Chiroderma_villosum,Chiroderma_improvisum),Chiroderma_doriae),(Vampyressa_brocki,Vampyressa_bidens))),Uroderma_magnirostrum)),(((Sturnira_bogotensis,Sturnira_tildae),Sturnira_magna),(Sturnira_oporaphilum,Sturnira_lilium))),Rhinophylla_fischerae),(((((Carollia_perspicillata,Carollia_brevicauda),Carollia_sowelli),Carollia_subrufa),Carollia_castanea),(Glyphonycteris_sylvestris,Trinycteris_nicefori))),((((Lonchophylla_hesperia,Lionycteris_spurrelli),Platalina_genovensium),Lonchophylla_chocoana),Lonchophylla_mordax)),((((((Glossophaga_longirostris,Glossophaga_commissarisi),Glossophaga_soricina),Leptonycteris_curasoae),((Phyllonycteris_poeyi,(Erophylla_sezekorni,Erophylla_bombifrons)),(Brachyphylla_pumila,Brachyphylla_cavernarum))$1),((Musonycteris_harrisoni,Choeroniscus_godmani),((Anoura_latidens,Anoura_geoffroyi),Anoura_caudifer))),(Lonchorhina_aurita,Lonchorhina_inusitata))),(((Phyllostomus_discolor,(Lophostoma_brasiliense,(Lophostoma_evotis,Lophostoma_silvicola))),(Tonatia_saurophila,Tonatia_bidens)),Mimon_cozumelae)),((Diaemus_youngi,Desmodus_rotundus),Diphylla_ecaudata)),((Micronycteris_hirsuta,(Micronycteris_minuta,Micronycteris_schmidtorum)),Lampronycteris_brachyotis)),(((Pteronotus_davyi,Pteronotus_quadridens),Pteronotus_parnellii),(Mormoops_megalophylla,Mormoops_blainvillei)$1)),(Noctilio_leporinus,Furipterus_horrens)$1),(Thyroptera_lavali,Thyroptera_tricolor)$1),Tadarida_brasilense#1,Chilonatalus_micropus#1),Cormura_brevirostris#1,((Megaderma_lyra,(Rhinolophus_ferrumequinum,Hipposideros_pratti)),((Acerodon_celebensis,Eonycteris_spelaea),Cynopterus_sphinx))$1);

codeml.CladeC.ctl
      seqfile = bat_data_match.phylip * sequence data filename
     treefile = bat_tree_match_cladeModelC_2groups.tre      * tree structure file name
      outfile = trpc2cladeModelC_2groups_0point1.txt           * main result file name

        noisy = 9  * 0,1,2,3,9: how much rubbish on the screen
      verbose = 0  * 0: concise; 1: detailed, 2: too much
      runmode = 0  * 0: user tree

      seqtype = 1  * 1:codons; 2:AAs; 3:codons-->AAs
    CodonFreq = 2  * 0:1/61 each, 1:F1X4, 2:F3X4, 3:codon table
        clock = 0  * 0:no clock, 1:clock; 2:local clock; 3:CombinedAnalysis

  * Model C: model = 3  NSsites = 2 
    * Model D: model = 3  NSsites = 3

        model = 3  * 3 = clade models
      NSsites = 2  * choose "2" or "3"
 
        icode = 0  * 0:universal code; 1:mammalian mt; 2-10:see below

    fix_kappa = 0  * 1: kappa fixed, 0: kappa to be estimated
        kappa = 2.5  * initial or fixed kappa

    fix_omega = 0  * 1: omega or omega_1 fixed, 0: estimate 
        omega = 0.1 * initial or fixed omega, for codons or codon-based AAs

    fix_alpha = 1  * 0: estimate gamma shape parameter; 1: fix it at alpha
        alpha = 0. * initial or fixed alpha, 0:infinity (constant rate)

        ncatG = 3  * # of categories in dG of NSsites models

        getSE = 0  * 0: don't want them, 1: want S.E.s of estimates
 RateAncestor = 0  * (0,1,2): rates (alpha>0) or ancestral states (1 or 2)

   Small_Diff = .2e-6
    cleandata = 0  * remove sites with ambiguity data (1:yes, 0:no)?

We are ready to go. Run the control file.

bash-3.2$ sudo ./codeml codeml.CladeC.ctl 

Running GARLI from our Model-O-Matic code.

If we remember from our model-o-matic output and sort by the lowest AIC, we can see the best model for OR6 is CodonM0.0.EQU+4dG. It is interesting how all codon models are preferred to the nucleotide models! This turned out to be the case in every set of alignments.

Things are looking good. I'll keep you posted on the results.

Here is the following best-fit models for each of the OR gene families (see summary files with AIC and likelihoods in /Volumes/Yango/Hayden_batORs/model_o_matic/excel_summaries):

OR137      CodonM0.0.EQU+4dG

OR213      CodonM0.0F64+4dG

OR4          CodonM0.0.F64+4dG

OR589      CodonM0.0.F64+4dG

OR6          CodonM0.0.EQU + 4dG

OR10        CodonM0.0.F64+4dG

OR11        CodonM0.0F64+4dG

OR51        CodonM0.0F64+4dG

OR52        CodonM0.0F64+4dG

Those highlighted in pink were run in model-O-matic as the -100normal parameter, as the program would crash if it was run as -normal. My thoughts are there are too many taxa, but other than than that I am at a loss.

Now we can apply this model to GARLI using our tree output. There are two options for running GARLI (at least that I use). One is on the CIPRES server, the other can be from command line. I am currently having CPU usage problems on CIPRES so I am going to run it from command.

Installing GARLI is easy if you have homebrew.

105-238:~ loloyohe$ brew install garli

Should be ready to go.

Let's take a second to interpret the output from

105-238:OR6 loloyohe$ garli garli.conf

Running GARLI Version 2.01.1067 (18 May 2012)

->Single processor version for 64-bit OS<-

##############################################################

 This is GARLI 2.0, the first "official" release including 

          partitioned models.  It is a merging of

   official release 1.0 and beta version GARLI-PART 0.97

  Briefly, it includes models for nucleotides, amino acids,

 codons, and morphology-like characters, any of which can be 

  mixed together and applied to different subsets of data.

    General program usage is extensively documented here:

            http://www.nescent.org/wg/garli/

      see this page for details on partitioned usage:

  http://www.nescent.org/wg/garli/Partition_testing_version

   and this page for details on Mkv mophology model usage:

    http://www.nescent.org/wg/garli/Mkv_morphology_model

         PLEASE LET ME KNOW OF ANY PROBLEMS AT:

                garli.support@gmail.com

##############################################################

Compiled Apr 27 2015 16:06:23 using GNU gcc compiler version 4.9.2

Using NCL version 2.1.17

#######################################################

Reading config file garli.conf

###################################################

READING OF DATA

Attempting to read data file in Nexus format (using NCL):

Phyllos_OR6_Funct_tAlign_noStops.nex ...

Reading DATA block...storing implied block: TAXA

storing read block: DATA

 successful

###################################################

PARTITIONING OF DATA AND MODELS

GARLI data subset 1

CHARACTERS block #1 ("Untitled DATA Block 1")

Data read as Nucleotide data,

modeled as Codon data

Gaps or ambiguity codes found within codon for taxon Lsi_29_OR6.

Codons coded as missing for that taxon: 6 120 235 236 

Gaps or ambiguity codes found within codon for taxon Phe_43_OR6.

Codons coded as missing for that taxon: 6 

Gaps or ambiguity codes found within codon for taxon Sti_17_OR6.

Codons coded as missing for that taxon: 185 

Summary of data:

  35 sequences.

  6 constant characters.

  232 parsimony-informative characters.

  5 uninformative variable characters.

  243 total characters.

  243 unique patterns in compressed data matrix.

Pattern processing required < 1 second

###################################################

NOTE: Unlike many programs, the amount of system memory that Garli will

use can be controlled by the user.

(This comes from the availablememory setting in the configuration file.

Availablememory should NOT be set to more than the actual amount of 

physical memory that your computer has installed)

For this dataset:

 Mem level availablememory setting

  great     >= 150 MB

  good approx 149 MB to 99 MB

  low approx 98 MB to 42 MB

  very low approx 42 MB to 31 MB

the minimum required availablememory is 31 MB

You specified that Garli should use at most 512.0 MB of memory.

Garli will actually use approx. 224.4 MB of memory

**Your memory level is: great (you don't need to change anything)**

#######################################################

Found outgroup specification:  1

#######################################################

STARTING RUN

>>>Search rep 1 (of 2)<<<

MODEL REPORT - Parameters are at their INITIAL values (not yet optimized)

Model 1

  Number of states = 61 (codon data, standard code)

  Nucleotide Relative Rate Matrix Assumed by Codon Model:     2 rates (transition and transversion) K param = 4.0000

  Equilibrium State Frequencies: equal (1/61 = 0.01639, fixed)

  Rate Heterogeneity Model:

    4 nonsynonymous rate categories, rate and proportion of each estimated

     (this is effectively the M3 model of PAML)

      dN/dS Proportion

      0.2500 0.2500

      0.5000 0.2500

      0.7500 0.2500

      1.0000 0.2500

Starting with seed=287687

creating likelihood stepwise addition starting tree...

number of taxa added:

 4  5  6  7  8  9  10  11  12  13  14  15  16  17

We can set the exact same thing up in CIPRES. 

I am using modelOmatic to determine which models best fit the sequences, and it include comparisons of codon models. Installing and running modelOmatic. For awhile I thought I was having Yosemite issues but because they are all binaries, everything seems to checking out okay. Lots of "duh" moments were had.

To install on Mac OSX Yosemite:
[1] Download the binaries and move to where you store your programs.

mv ~/Downloads/modelomatic_OSX /usr/local/bin/

[2] Set modelOmatic into your path.

vi .bash_profile  #in your home directory, this opens up your bash profile

#press "i" to edit and then paste this into your .bash_profile

export PATH=/usr/local/bin/modelomatic_OSX:$PATH 

#press ESC, then type "wq!"; restart the terminal

[3] ModelOmatic requires a tree file to run. It can make a NJ tree if you have bioNJ set up. Download the binaries for bioNJ here.

[4] Move the bioNJ into where you store your programs. Set your path. You have to rename your path in lower case so that modelOmatic can read it.

mv ~/Downloads/BIONJ/ /usr/local/bin/ 

mv /usr/local/bin/BIONJ/BIONJ /usr/local/bin/BIONJ/bionj

#add to your path variable
vi .bash_profile

export PATH=/usr/local/bin/BIONJ/BIONJ:$PATH

export PATH=/usr/local/bin/BIONJ:$PATH

The format for input 

./ModelOMatic <data> bionj <output> <genetic_code> fast

I am running it as modelomatic_OSX because it is in my path. 

Using the bionj option tell modelOmatic to make its own tree using bioNJ.

Specify in <output> slot where to put your output folder.

The <genetic_code> has been tricky. It can choke and tell you that you have stop codons if you don't use the correct number:

     0:  Universal code
     1:  Vertebrate mt
     2:  Yeast mt
     3:  Mould mt

     4:  Invertebrate mt
     5:  Ciliate nuclear
     6:  Echinoderm mt
     7:  Euplotid mt

     8:  Alternative yeast nuclear
     9:  Ascidian mt

     10: Blepharisma nuclear

     11: Everything codes (64 character state-space)

My adviser was having trouble but then realized she was using cytb so it actually should have been "1" instead of "0". 

ModelOmatic uses .phylip file. Allegedly can use other formats but I have not tried. However, it actually seems ridiculously flexible in reading a .phylip file.

OR_alignments loloyohe$ modelomatic_OSX Phyllos_OR6_Funct_tAlign.phylip bionj test.txt 1

---------------------------------------------------------------------------------------------

  ModelOMatic (v1.01 (release)).

A program for choosing substitutions models for phylogenetic inference.

Written by Simon Whelan.

Contributions from James Allen, Ben Blackburne and David Talavera.

---------------------------------------------------------------------------------------------

Data: <Phyllos_OR6_Funct_tAlign.phylip>: 35 sequences of length 729 (DataMatrix: 35 x 570)

Checking for sparse (>85% gaps) sequences

Starting with 35 ... after removal there are 35 sequences

Creating start tree ... 

Assertion failed: (0), function Error, file tools.cxx, line 235.

Unrecognised codon: AGA for data codon position 241 in sequence 1...Abort trap: 6

We had stop codons. Okay, so I translated the sequences, looked for where the stop codons were, and then I replaced them in the original .phylip file using "---".  Note in the future that all .phylip sequence files in /Volumes/Yango/Hayden_ORs/OR_alignments have the stop codon replaced, while the .nex files have the original sequence alignments. Lets try to get it running...

modelomatic_OSX  /Volumes/Yango/Hayden_batORs/OR_alignments/Phyllos_OR4_Funct_tAlign.phylip bionj /Volumes/Yango/Hayden_batORs/Phyllos_OR4_Funct_tAlign_mOm_results.txt 0 fast

---------------------------------------------------------------------------------------------

  ModelOMatic (v1.01 (release)).

A program for choosing substitutions models for phylogenetic inference.

Written by Simon Whelan.

Contributions from James Allen, Ben Blackburne and David Talavera.

---------------------------------------------------------------------------------------------

Data: </Volumes/Yango/Hayden_batORs/OR_alignments/Phyllos_OR4_Funct_tAlign.phylip>: 127 sequences of length 753 (DataMatrix: 127 x 666)

Checking for sparse (>85% gaps) sequences

Starting with 127 ... after removal there are 127 sequences

Creating start tree ...  estimated using bionj (1.15646s)

Optimisation settings: fast 

Scanning for modelomatic.ini file ... done

Working with genetic code: Universal

>>> Doing model analysis <<< 

RY Done  (3.56514s)

NT Done  (10.6517s)

AA Done  (71.2184s).........

Codon Done  (232.733s)

Outputting results to <Phyllos_OR4_Funct_tAlign_mOm_results.txt>

Successful exit

####now for the others

modelomatic_OSX /Volumes/Yango/Hayden_batORs/OR_alignments/Phyllos_OR6_Funct_tAlign.phylip bionj /Volumes/Yango/Hayden_batORs/Phyllos_OR6_Funct_tAlign_mOm_results.txt 0 fast

modelomatic_OSX  /Volumes/Yango/Hayden_batORs/OR_alignments/Phyllos_OR10_Funct_tAlign.phylip bionj /Volumes/Yango/Hayden_batORs/Phyllos_OR10_Funct_tAlign_mOm_results.txt 0 fast

modelomatic_OSX  /Volumes/Yango/Hayden_batORs/OR_alignments/Phyllos_OR51_Funct_tAlign.phylip bionj /Volumes/Yango/Hayden_batORs/Phyllos_OR51_Funct_tAlign_mOm_results.txt 0 fast

Had a hiccup with this data:

modelomatic_OSX  /Volumes/Yango/Hayden_batORs/OR_alignments/Phyllos_OR2_13_Funct_tAlign.phylip bionj /Volumes/Yango/Hayden_batORs/model_o_matic/Phyllos_OR2_13_Funct_tAlign_mOm_results.txt 0 normal

---------------------------------------------------------------------------------------------

  ModelOMatic (v1.01 (release)).

A program for choosing substitutions models for phylogenetic inference.

Written by Simon Whelan.

Contributions from James Allen, Ben Blackburne and David Talavera.

---------------------------------------------------------------------------------------------

Data: </Volumes/Yango/Hayden_batORs/OR_alignments/Phyllos_OR2_13_Funct_tAlign.phylip>: 264 sequences of length 1530 (DataMatrix: 264 x 943)

Checking for sparse (>85% gaps) sequences

Starting with 264 ... after removal there are 264 sequences

Creating start tree ...  estimated using bionj (4.68751s)

Optimisation settings: normal 

Scanning for modelomatic.ini file ... done

Working with genetic code: Universal

>>> Doing model analysis <<< 

Broken Prob(): -0.000104812

modelomatic_OSX  /Volumes/Yango/Hayden_batORs/OR_alignments/Phyllos_OR137_Funct_tAlign.phylip bionj /Volumes/Yango/Hayden_batORs/Phyllos_OR137_Funct_tAlign_mOm_results.txt 0 fast

---------------------------------------------------------------------------------------------

  ModelOMatic (v1.01 (release)).

A program for choosing substitutions models for phylogenetic inference.

Written by Simon Whelan.

Contributions from James Allen, Ben Blackburne and David Talavera.

---------------------------------------------------------------------------------------------

Data: </Volumes/Yango/Hayden_batORs/OR_alignments/Phyllos_OR137_Funct_tAlign.phylip>: 452 sequences of length 804 (DataMatrix: 452 x 766)

Checking for sparse (>85% gaps) sequences

Starting with 452 ... after removal there are 452 sequences

Creating start tree ...  estimated using bionj (4.7425s)

Optimisation settings: fast 

Scanning for modelomatic.ini file ... done

Working with genetic code: Universal

>>> Doing model analysis <<< 

RY Done  (8.05337s)

NT Done  (44.7082s)

Broken Prob(): -40.8111

modelomatic_OSX  /Volumes/Yango/Hayden_batORs/OR_alignments/Phyllos_OR52_Funct_tAlign.phylip bionj /Volumes/Yango/Hayden_batORs/Phyllos_OR52_Funct_tAlign_mOm_results.txt 0 fast

---------------------------------------------------------------------------------------------

  ModelOMatic (v1.01 (release)).

A program for choosing substitutions models for phylogenetic inference.

Written by Simon Whelan.

Contributions from James Allen, Ben Blackburne and David Talavera.

---------------------------------------------------------------------------------------------

Data: </Volumes/Yango/Hayden_batORs/OR_alignments/Phyllos_OR51_Funct_tAlign.phylip>: 152 sequences of length 750 (DataMatrix: 152 x 712)

Checking for sparse (>85% gaps) sequences

Starting with 152 ... after removal there are 152 sequences

Creating start tree ...  estimated using bionj (1.32059s)

Optimisation settings: fast 

Scanning for modelomatic.ini file ... done

Working with genetic code: Universal

>>> Doing model analysis <<< 

RY Done  (4.69988s)

NT Done  (13.8241s)

AA Done  (85.6409s).........

...

Broken Prob(): -2.64535e+41

I think the problem is in running the "fast" option with large datasets. When I use "trimfast" and "normal", they both work. 

modelomatic_OSX  /Volumes/Yango/Hayden_batORs/OR_alignments/Phyllos_OR51_Funct_tAlign.phylip bionj /Volumes/Yango/Hayden_batORs/model_o_matic/Phyllos_OR51_Funct_tAlign_mOm_results.txt 0 trimfast

---------------------------------------------------------------------------------------------

  ModelOMatic (v1.01 (release)).

A program for choosing substitutions models for phylogenetic inference.

Written by Simon Whelan.

Contributions from James Allen, Ben Blackburne and David Talavera.

---------------------------------------------------------------------------------------------

Data: </Volumes/Yango/Hayden_batORs/OR_alignments/Phyllos_OR51_Funct_tAlign.phylip>: 152 sequences of length 750 (DataMatrix: 152 x 712)

Checking for sparse (>85% gaps) sequences

Starting with 152 ... after removal there are 152 sequences

Creating start tree ...  estimated using bionj (1.58694s)

Optimisation settings: fast trim=10

TRIMMING: Tree contains 152 sequences (>TrimTree=10) ...

          Loose optimisation of start tree under JC ... done

          Obtaining greedy start tree with 10 sequences ... done

          Reinitialising objects for trimmed data ... done

          New files available: data = </Volumes/Yango/Hayden_batORs/model_o_matic/Phyllos_OR51_Funct_tAlign_mOm_results.txt.trim.data>; tree = </Volumes/Yango/Hayden_batORs/model_o_matic/Phyllos_OR51_Funct_tAlign_mOm_results.txt.trim.tree>

Scanning for modelomatic.ini file ... done

Working with genetic code: Universal

>>> Doing model analysis <<< 

RY Done  (1.08313s)

NT Done  (1.18703s)

AA Done  (6.64464s).........

Codon Done  (24.4936s)

Outputting results to </Volumes/Yango/Hayden_batORs/model_o_matic/Phyllos_OR51_Funct_tAlign_mOm_results.txt>

Successful exit

However, word of caution. The "trimfast" option gives a vastly different result:

Normal it is! With fewer parameters, its best just to wait it out.