.gitignore

@@ -80,3 +80,14 @@ sftp-config.json
 *.svg
 
 !data/
+
+# OS generated files #
+######################
+.DS_Store
+.DS_Store?
+._*
+.Spotlight-V100
+.Trashes
+Icon?
+ehthumbs.db
+Thumbs.db

.travis_test.sh

-#export PYTHONPATH="$PYTHONPATH:./test_install"
 nosetests test/test_treetime.py

__init__.py

@@ -3,8 +3,10 @@ from treetime.clock_tree import ClockTree
 from treetime.treetime import TreeTime
 from treetime.treetime import ttconf as treetime_conf
 from treetime.gtr import GTR
-from treetime.io import *
+from treetime.treetime import plot_vs_years
+from treetime.treetime import treetime_to_newick
 from treetime import seq_utils
 from treetime.merger_models import Coalescent
 import datetime
 import treetime.seq_utils as seq_utils
+from treetime.utils import numeric_date

ancestral_reconstruction.py → bin/ancestral_reconstruction.py

 #!/usr/bin/env python
 from __future__ import print_function, division
 import numpy as np
-from treetime import TreeAnc
+from treetime import TreeAnc, GTR
 from Bio import Phylo, AlignIO
+from Bio import __version__ as bioversion
 
 if __name__=="__main__":
     ###########################################################################
@@ -13,45 +14,110 @@ if __name__=="__main__":
             description='Reconstructs ancestral sequences and maps mutations to the tree.'
                         ' The output consists of a file ending with _ancestral.fasta with ancestral sequences'
                         ' and a tree ending with _mutation.nexus with mutations added as comments'
-                        ' like _A45G_... The inferred GTR model is written to stdout')
+                        ' like _A45G_..., number in SNPs used 1-based index by default.'
+                        ' The inferred GTR model is written to stdout')
     parser.add_argument('--aln', required = True, type = str,  help ="fasta file with input sequences")
-    parser.add_argument('--tree', required = True, type = str,  help ="newick file with tree")
+    parser.add_argument('--tree', type = str,  help ="newick file with tree, "
+                                                     "will attempt to build tree if none given.")
+
+    parser.add_argument('--gtr', type = str, default='infer', help="GTR model to use. "
+        " Type 'infer' to infer the model from the data. Or, specify the model type. "
+        " If the specified model requires additional options, use '--gtr_args' to specify those")
+
+    parser.add_argument('--gtr_params', type=str, nargs='+', help="GTR parameters for the model "
+        "specified by the --gtr argument. The parameters should be feed as 'key=value' list of parameters. "
+        "Example: '--gtr K80 --gtr_params kappa=0.2 pis=0.25,0.25,0.25,0.25'. See the exact definitions of "
+        " the parameters in the GTR creation methods in treetime/nuc_models.py or treetime/aa_models.py")
+
     parser.add_argument('--prot', default = False, action="store_true", help ="protein alignment")
-    parser.add_argument('--marginal', default = False, action='store_true', help='marginal instead of joint ML reconstruction')
-    parser.add_argument('--infer_gtr', default = False, action='store_true', help='infer substitution model')
+    parser.add_argument('--marginal', default = False, action="store_true", help ="marginal reconstruction of ancestral sequences")
+    parser.add_argument('--zero_based', default = False, action='store_true', help='zero based SNP indexing')
     parser.add_argument('--keep_overhangs', default = False, action='store_true', help='do not fill terminal gaps')
+    parser.add_argument('--verbose', default = 1, type=int, help='verbosity of output 0-6')
     params = parser.parse_args()
 
+
+    ###########################################################################
+    ### CHECK FOR TREE, build if not in place
+    ###########################################################################
+    if params.tree is None:
+        from treetime.utils import tree_inference
+        import os,shutil
+        params.tree = os.path.basename(params.aln)+'.nwk'
+        print("No tree given: inferring tree")
+        tmp_dir = 'ancestral_reconstruction_tmp_files'
+        tree_inference(params.aln, params.tree, tmp_dir = tmp_dir)
+        if os.path.isdir(tmp_dir):
+            shutil.rmtree(tmp_dir)
+
+    ###########################################################################
+    ### GTR SET-UP
+    ###########################################################################
+    model = params.gtr
+    gtr_params = params.gtr_params
+    if model == 'infer':
+        gtr = GTR.standard('jc')
+        infer_gtr = True
+    else:
+        try:
+            kwargs = {}
+            if gtr_params is not None:
+                for param in gtr_params:
+                    keyval = param.split('=')
+                    if len(keyval)!=2: continue
+                    if keyval[0] in ['pis', 'pi', 'Pi', 'Pis']:
+                        keyval[1] = map(float, keyval[1].split(','))
+                    elif keyval[0] not in ['alphabet']:
+                        keyval[1] = float(keyval[1])
+                    kwargs[keyval[0]] = keyval[1]
+            else:
+                print ("GTR params are not specified. Creating GTR model with default parameters")
+
+
+            gtr = GTR.standard(model, **kwargs)
+            infer_gtr = False
+        except:
+            print ("Could not create GTR model from input arguments. Using default (Jukes-Cantor 1969)")
+            gtr = GTR.standard('jc')
+            infer_gtr = False
+
+
     ###########################################################################
     ### ANCESTRAL RECONSTRUCTION
     ###########################################################################
-    model = 'aa' if params.prot else 'Jukes-Cantor'
-    treeanc = TreeAnc(params.tree, aln=params.aln, gtr=model, verbose=4, fill_overhangs=not params.keep_overhangs)
-    treeanc.infer_ancestral_sequences('ml', infer_gtr=params.infer_gtr,
+    treeanc = TreeAnc(params.tree, aln=params.aln, gtr=gtr, verbose=4, fill_overhangs=not params.keep_overhangs)
+    treeanc.infer_ancestral_sequences('ml', infer_gtr=infer_gtr,
                                        marginal=params.marginal)
 
     ###########################################################################
     ### OUTPUT and saving of results
     ###########################################################################
-    if params.infer_gtr:
+
+    model = 'aa' if params.prot else 'Jukes-Cantor'
+    if infer_gtr:
         print('\nInferred GTR model:')
         print(treeanc.gtr)
 
     outaln_name = '.'.join(params.aln.split('/')[-1].split('.')[:-1])+'_ancestral.fasta'
     AlignIO.write(treeanc.get_reconstructed_alignment(), outaln_name, 'fasta')
+    print("--- alignment including ancestral nodes saved as  \n\t %s\n"%outaln_name)
 
     # decorate tree with inferred mutations
     terminal_count = 0
+    offset = 0 if params.zero_based else 1
     for n in treeanc.tree.find_clades():
         if n.up is None:
             continue
         n.confidence=None
-        if n.is_terminal() and len(n.name)>40:
+        # due to a bug in older versions of biopython that truncated filenames in nexus export
+        # we truncate them by hand and make them unique.
+        if n.is_terminal() and len(n.name)>40 and bioversion<"1.69":
             n.name = n.name[:35]+'_%03d'%terminal_count
             terminal_count+=1
         if len(n.mutations):
-            n.comment= '&mutations="' + '_'.join([a+str(pos)+d for (a,pos, d) in n.mutations])+'"'
+            n.comment= '&mutations="' + '_'.join([a+str(pos + offset)+d for (a,pos, d) in n.mutations])+'"'
 
     # write tree to file
     outtree_name = '.'.join(params.tree.split('/')[-1].split('.')[:-1])+'_mutation.nexus'
     Phylo.write(treeanc.tree, outtree_name, 'nexus')
+    print("--- tree saved in nexus format as  \n\t %s\n"%outtree_name)

bin/homoplasy_scanner.py

+#!/usr/bin/env python
+from __future__ import print_function, division
+import numpy as np
+from treetime import TreeAnc, GTR
+from Bio import Phylo, AlignIO
+from Bio import __version__ as bioversion
+import os,shutil
+
+if __name__=="__main__":
+    ###########################################################################
+    ### parameter parsing
+    ###########################################################################
+    import argparse
+    parser = argparse.ArgumentParser(
+            description='Reconstructs ancestral sequences and maps mutations to the tree.'
+                        ' The tree is then scanned for homoplasies. An excess number of homoplasies'
+                        ' might suggest contamination, recombination, culture adaptation or similar. ')
+    parser.add_argument('--aln', required = True, type = str,  help ="fasta file with input nucleotide sequences")
+    parser.add_argument('--tree', type = str,  help ="newick file with tree (optional if tree builders installed)")
+    parser.add_argument('--const', type = int, default=0, help ="number of constant sites not included in alignment")
+    parser.add_argument('--rescale', type = float, default=1.0, help ="rescale branch lengths")
+    parser.add_argument('--detailed', required = False, action="store_true",  help ="generate a more detailed report")
+    parser.add_argument('--gtr', required=False, type = str, default='infer', help="GTR model to use. "
+        " Type 'infer' to infer the model from the data. Or, specify the model type. "
+        " If the specified model requires additional options, use '--gtr_args' to specify those")
+
+    parser.add_argument('--gtr_params', type=str, nargs='+', help="GTR parameters for the model "
+        "specified by the --gtr argument. The parameters should be feed as 'key=value' list of parameters. "
+        "Example: '--gtr K80 --gtr_params kappa=0.2 pis=0.25,0.25,0.25,0.25'. See the exact definitions of "
+        " the parameters in the GTR creation methods in treetime/nuc_models.py. Only nucleotide models supported at present")
+
+    parser.add_argument('--prot', default = False, action="store_true", help ="protein alignment")
+    parser.add_argument('--zero_based', default = False, action='store_true', help='zero based SNP indexing')
+    parser.add_argument('-n', default = 10, type=int, help='number of mutations/nodes that are printed to screen')
+    parser.add_argument('--verbose', default = 1, type=int, help='verbosity of output 0-6')
+    params = parser.parse_args()
+
+
+    ###########################################################################
+    ### CHECK FOR TREE, build if not in place
+    ###########################################################################
+    if params.tree is None:
+        from treetime.utils import tree_inference
+        params.tree = os.path.basename(params.aln)+'.nwk'
+        print("No tree given: inferring tree")
+        tmp_dir = 'homoplasy_scanner_tmp_files'
+        tree_inference(params.aln, params.tree, tmp_dir = tmp_dir)
+        if os.path.isdir(tmp_dir):
+            shutil.rmtree(tmp_dir)
+
+    ###########################################################################
+    ### GTR SET-UP
+    ###########################################################################
+    model = params.gtr
+    gtr_params = params.gtr_params
+    if model == 'infer':
+        gtr = GTR.standard('jc')
+        infer_gtr = True
+    else:
+        infer_gtr = False
+        try:
+            kwargs = {}
+            if gtr_params is not None:
+                for param in gtr_params:
+                    keyval = param.split('=')
+                    if len(keyval)!=2: continue
+                    if keyval[0] in ['pis', 'pi', 'Pi', 'Pis']:
+                        keyval[1] = map(float, keyval[1].split(','))
+                    elif keyval[0] not in ['alphabet']:
+                        keyval[1] = float(keyval[1])
+                    kwargs[keyval[0]] = keyval[1]
+            else:
+                print ("GTR params are not specified. Creating GTR model with default parameters")
+
+
+            gtr = GTR.standard(model, **kwargs)
+        except:
+            print ("Could not create GTR model from input arguments. Using default (Jukes-Cantor 1969)")
+            gtr = GTR.standard('jc')
+
+
+    ###########################################################################
+    ### ANCESTRAL RECONSTRUCTION
+    ###########################################################################
+    treeanc = TreeAnc(params.tree, aln=params.aln, gtr=gtr, verbose=1,
+                      fill_overhangs=True)
+    L = treeanc.aln.get_alignment_length() + params.const
+    N_seq = len(treeanc.aln)
+    N_tree = treeanc.tree.count_terminals()
+    if params.rescale!=1.0:
+        for n in treeanc.tree.find_clades():
+            n.branch_length *= params.rescale
+            n.mutation_length = n.branch_length
+
+    print("read alignment from file %s with %d sequences of length %d"%(params.aln,N_seq,L))
+    print("read tree from file %s with %d leaves"%(params.tree,N_tree))
+    print("\ninferring ancestral sequences...")
+
+    treeanc.infer_ancestral_sequences('ml', infer_gtr=infer_gtr, marginal=False)
+    print("...done.")
+
+    ###########################################################################
+    ### analysis of reconstruction
+    ###########################################################################
+    from collections import defaultdict
+    from scipy.stats import poisson
+    offset = 0 if params.zero_based else 1
+
+    # construct dictionaries gathering mutations and positions
+    mutations = defaultdict(list)
+    positions = defaultdict(list)
+    terminal_mutations = defaultdict(list)
+    for n in treeanc.tree.find_clades():
+        if n.up is None:
+            continue
+
+        if len(n.mutations):
+            for (a,pos, d) in n.mutations:
+                if '-' not in [a,d]:
+                    mutations[(a,pos+offset,d)].append(n)
+                    positions[pos+offset].append(n)
+            if n.is_terminal():
+                for (a,pos, d) in n.mutations:
+                    if '-' not in [a,d]:
+                        terminal_mutations[(a,pos+offset,d)].append(n)
+
+    # gather homoplasic mutations by strain
+    mutation_by_strain = defaultdict(list)
+    for n in treeanc.tree.get_terminals():
+        for a,pos,d in n.mutations:
+            if pos in positions and len(positions[pos])>1:
+                mutation_by_strain[n.name].append([(a,pos+offset,d), len(positions[pos])])
+
+
+    # total_branch_length is the expected number of substitutions
+    # corrected_branch_length is the expected number of observable substitutions
+    # (probability of an odd number of substitutions at a particular site)
+    total_branch_length = treeanc.tree.total_branch_length()
+    corrected_branch_length = np.sum([np.exp(-x.branch_length)*np.sinh(x.branch_length)
+                                      for x in treeanc.tree.find_clades()])
+    corrected_terminal_branch_length = np.sum([np.exp(-x.branch_length)*np.sinh(x.branch_length)
+                                      for x in treeanc.tree.get_terminals()])
+    expected_mutations = L*corrected_branch_length
+    expected_terminal_mutations = L*corrected_terminal_branch_length
+
+    # make histograms and sum mutations in different categories
+    multiplicities = np.bincount([len(x) for x in mutations.values()])
+    total_mutations = np.sum([len(x) for x in mutations.values()])
+
+    multiplicities_terminal = np.bincount([len(x) for x in terminal_mutations.values()])
+    terminal_mutation_count = np.sum([len(x) for x in terminal_mutations.values()])
+
+    multiplicities_positions = np.bincount([len(x) for x in positions.values()])
+    multiplicities_positions[0] = L - np.sum(multiplicities_positions)
+
+    ###########################################################################
+    ### Output the distribution of times particular mutations are observed
+    ###########################################################################
+    print("\nThe TOTAL tree length is %1.3e, expecting %1.1f mutations vs an observed %d"
+          %(total_branch_length,expected_mutations,total_mutations))
+    print("Of these %d mutations,"%total_mutations
+            +"".join(['\n\t - %d occur %d times'%(n,mi)
+                      for mi,n in enumerate(multiplicities) if n]))
+    # additional optional output this for terminal mutations only
+    if params.detailed:
+        print("\nThe TERMINAL branch length is %1.3e, expecting %1.1f mutations vs an observed %d"
+              %(corrected_terminal_branch_length,expected_terminal_mutations,terminal_mutation_count))
+        print("Of these %d mutations,"%terminal_mutation_count
+                +"".join(['\n\t - %d occur %d times'%(n,mi)
+                          for mi,n in enumerate(multiplicities_terminal) if n]))
+
+
+    ###########################################################################
+    ### Output the distribution of times mutations at particular positions are observed
+    ###########################################################################
+    print("\nOf the %d positions in the genome,"%L
+            +"".join(['\n\t - %d were hit %d times (expected %1.2f)'%(n,mi,L*poisson.pmf(mi,1.0*total_mutations/L))
+                      for mi,n in enumerate(multiplicities_positions) if n]))
+
+
+    # compare that distribution to a Poisson distribution with the same mean
+    p = poisson.pmf(np.arange(10*multiplicities_positions.max()),1.0*total_mutations/L)
+    print("\nlog-likelihood difference to Poisson distribution with same mean: %1.3e"%(
+            - L*np.sum(p*np.log(p+1e-100))
+            + np.sum(multiplicities_positions*np.log(p[:len(multiplicities_positions)]+1e-100))))
+
+
+    ###########################################################################
+    ### Output the mutations that are observed most often
+    ###########################################################################
+    print("\n\nThe ten most homoplasic mutations are:\n\tmut\tmultiplicity")
+    mutations_sorted = sorted(mutations.items(), key=lambda x:len(x[1])-0.1*x[0][1]/L, reverse=True)
+    for mut, val in mutations_sorted[:params.n]:
+        if len(val)>1:
+            print("\t%s%d%s\t%d"%(mut[0], mut[1], mut[2], len(val)))
+        else:
+            break
+
+    # optional output specifically for mutations on terminal branches
+    if params.detailed:
+        print("\n\nThe ten most homoplasic mutation on terminal branches are:\n\tmut\tmultiplicity")
+        terminal_mutations_sorted = sorted(terminal_mutations.items(), key=lambda x:len(x[1])-0.1*x[0][1]/L, reverse=True)
+        for mut, val in terminal_mutations_sorted[:params.n]:
+            if len(val)>1:
+                print("\t%s%d%s\t%d"%(mut[0], mut[1], mut[2], len(val)))
+            else:
+                break
+
+    ###########################################################################
+    ### Output strains that have many homoplasic mutations
+    ###########################################################################
+    # TODO: add statistical criterion
+    if params.detailed:
+        print("\n\nTaxons that carry positions that mutated elsewhere in the tree:\n\ttaxon name\t#of homoplasic mutations")
+        mutation_by_strain_sorted = sorted(mutation_by_strain.items(), key=lambda x:len(x[1]), reverse=True)
+        for name, val in mutation_by_strain_sorted[:params.n]:
+            if len(val):
+                print("\t%s\t%d"%(name, len(val)))

bin/mugration.py

+#!/usr/bin/env python
+from __future__ import print_function, division
+import numpy as np
+import pandas as pd
+from treetime import TreeAnc, GTR
+from Bio.SeqRecord import SeqRecord
+from Bio.Seq import Seq
+from Bio.Align import MultipleSeqAlignment
+from Bio import Phylo, AlignIO
+from Bio import __version__ as bioversion
+import os
+
+if __name__=="__main__":
+    ###########################################################################
+    ### parameter parsing
+    ###########################################################################
+    import argparse
+    parser = argparse.ArgumentParser(
+            description='Reconstructs discrete ancestral states, for example '
+                        'geographic location, host, or similar.')
+    parser.add_argument('--tree', required = True, type=str, help ="newick file with tree")
+    parser.add_argument('--attribute', type=str, help ="attribute to reconstruct, e.g. country")
+    parser.add_argument('--states', required = True, type=str, help ="csv or tsv file with discrete characters."
+                                    "\n#name,country,continent\ntaxon1,micronesia,oceania\n...")
+    parser.add_argument('--weights', type=str, help="csv or tsv file with probabilities of that a randomly sampled "
+                        "sequence has a particular state. E.g. population of different continents or countries. E.g.:"
+                        "\n#country,weight\nmicronesia,0.1\n...")
+    # parser.add_argument('--migration', type=str, help="csv or tsv file with symmetric migration/transition rates "
+    #                     "between states. For example passenger flow.")
+    # parser.add_argument('--infer_gtr', action="store_true", help="infer GTR model from tree. "
+    #                                 "Ignored when prop or migration is specified.")
+    parser.add_argument('--confidence', action="store_true", help="output confidence of mugration inference")
+    parser.add_argument('--pc', type=float, default=1.0, help ="pseudo-counts higher numbers will results in 'flatter' models")
+    parser.add_argument('--verbose', default = 1, type=int, help='verbosity of output 0-6')
+    params = parser.parse_args()
+    missing = "?"
+
+    ###########################################################################
+    ### Parse states
+    ###########################################################################
+    if os.path.isfile(params.states):
+        states = pd.read_csv(params.states, sep='\t' if params.states[-3:]=='tsv' else ',',
+                             skipinitialspace=True)
+    else:
+        print("file with states does not exist")
+        exit(1)
+
+    taxon_name = 'name' if 'name' in states.columns else states.columns[0]
+    if params.attribute and params.attribute in states.columns:
+        attr = params.attribute
+    else:
+        attr = states.columns[1]
+
+    leaf_to_attr = {x[taxon_name]:x[attr] for xi, x in states.iterrows()
+                    if x[attr]!=missing}
+    unique_states = sorted(set(leaf_to_attr.values()))
+    nc = len(unique_states)
+    if nc>180:
+        print("mugration: can't have more than 180 states!")
+        exit(1)
+    elif nc<2:
+        print("mugration: only one or zero states found -- this doesn't make any sense")
+        exit(1)
+
+    ###########################################################################
+    ### make a single character alphabet that maps to discrete states
+    ###########################################################################
+    alphabet = [chr(65+i) for i,state in enumerate(unique_states)]
+    missing_char = chr(65+nc)
+    letter_to_state = {a:unique_states[i] for i,a in enumerate(alphabet)}
+    letter_to_state[missing_char]=missing
+    reverse_alphabet = {v:k for k,v in letter_to_state.items()}
+
+    ###########################################################################
+    ### construct gtr model
+    ###########################################################################
+    if params.weights:
+        params.infer_gtr = True
+        tmp_weights = pd.read_csv(params.weights, sep='\t' if params.states[-3:]=='tsv' else ',',
+                             skipinitialspace=True)
+        weights = {row[0]:row[1] for ri,row in tmp_weights.iterrows()}
+        mean_weight = np.mean(weights.values())
+        weights = np.array([weights[c] if c in weights else mean_weight for c in unique_states], dtype=float)
+        weights/=weights.sum()
+    else:
+        weights = np.ones(nc, dtype=float)/nc
+
+    # set up dummy matrix
+    W = np.ones((nc,nc), dtype=float)
+
+    mugration_GTR = GTR.custom(pi = weights, W=W, alphabet = np.array(alphabet))
+    mugration_GTR.profile_map[missing_char] = np.ones(nc)
+    mugration_GTR.ambiguous=missing_char
+
+    ###########################################################################
+    ### set up treeanc
+    ###########################################################################
+    treeanc = TreeAnc(params.tree, gtr=mugration_GTR, verbose=params.verbose)
+    pseudo_seqs = [SeqRecord(id=n.name,name=n.name,
+                   seq=Seq(reverse_alphabet[leaf_to_attr[n.name]] if n.name in leaf_to_attr else missing))
+                   for n in treeanc.tree.get_terminals()]
+    treeanc.aln = MultipleSeqAlignment(pseudo_seqs)
+
+    treeanc.infer_ancestral_sequences(method='ml', infer_gtr=True,
+            store_compressed=False, pc=params.pc, marginal=True, normalized_rate=False,
+            fixed_pi=weights if params.weights else None)
+
+
+    ###########################################################################
+    ### output
+    ###########################################################################
+    print("\nCompleted mugration model inference of attribute '%s' for"%attr,params.tree)
+
+    bname = './'+os.path.basename(params.tree)
+    gtr_name = bname + '.GTR.txt'
+    with open(gtr_name, 'w') as ofile:
+        ofile.write('Character to attribute mapping:\n')
+        for state in unique_states:
+            ofile.write('  %s: %s\n'%(reverse_alphabet[state], state))
+        ofile.write('\n\n'+str(treeanc.gtr)+'\n')
+        print("\nSaved inferred mugration model as:", gtr_name)
+
+    terminal_count = 0
+    for n in treeanc.tree.find_clades():
+        if n.up is None:
+            continue
+        n.confidence=None
+        # due to a bug in older versions of biopython that truncated filenames in nexus export
+        # we truncate them by hand and make them unique.
+        if n.is_terminal() and len(n.name)>40 and bioversion<"1.69":
+            n.name = n.name[:35]+'_%03d'%terminal_count
+            terminal_count+=1
+        n.comment= '&%s="'%attr + letter_to_state[n.sequence[0]] +'"'
+
+    if params.confidence:
+        conf_name = bname+'.confidence.csv'
+        with open(conf_name, 'w') as ofile:
+            ofile.write('#name, '+', '.join(unique_states)+'\n')
+            for n in treeanc.tree.find_clades():
+                ofile.write(n.name + ', '+', '.join([str(x) for x in n.marginal_profile[0]])+'\n')
+        print("Saved table with ancestral state confidences as:", conf_name)
+
+    # write tree to file
+    outtree_name = bname+'.mugration.nexus'
+    Phylo.write(treeanc.tree, outtree_name, 'nexus')
+    print("Saved annotated tree as:",outtree_name)
\ No newline at end of file

temporal_signal.py → bin/temporal_signal.py

@@ -12,15 +12,17 @@ if __name__=="__main__":
     import argparse
     parser = argparse.ArgumentParser(
             description="Calculates the root-to-tip regression and quantifies the 'clock-i-ness' of the tree. "
-                        "It will optionally reroot the tree to maximize the clock-like signal and recalculate branch length.")
+                        "It will reroot the tree to maximize the clock-like "
+                        "signal and recalculate branch length unless run with --keep_root.")
     parser.add_argument('--tree', required = True, type = str,  help ="newick file with tree")
     parser.add_argument('--dates', required = True, type = str,
                         help ="csv with dates for nodes with 'node_name, date' where date is float (as in 2012.15)")
-    parser.add_argument('--aln', required = False, type = str,  help ="fasta file with input sequences")
     parser.add_argument('--infer_gtr', default = False, action='store_true', help='infer substitution model')
-    parser.add_argument('--reroot', required = False, action="store_true", default=False,
-                        help ="reroot the tree to maximize the correlation of root-to-tip distance with sampling time")
-    parser.add_argument('--plot', required = False, action="store_true", default=False,)
+    parser.add_argument('--keep_root', required = False, action="store_true", default=False,
+                        help ="don't reroot the tree. Otherwise, reroot to minimize the "
+                              "the residual of the regression of root-to-tip distance and sampling time")
+    parser.add_argument('--plot', required = False, action="store_true", default=False,
+                        help = "save the root-to-tip regression as a pdf")
     parser.add_argument('--verbose', default = 0, type=int,
                         help='verbosity of output 0-6')
     params = parser.parse_args()
@@ -37,18 +39,18 @@ if __name__=="__main__":
                 dates[name] = float(date)
             except:
                 failed_dates+=1
-                print("couldn't parse date from:",line.strip(),"\n\texpecting float in second column")
+
         if len(dates)<failed_dates:
             print("\n\nDATE PARSING FAILED, ABORTING...")
             import sys
             sys.exit(1)
 
+
     ###########################################################################
-    ### FAKING ALIGMENT IF NONE GIVEN
+    ### FAKING ALIGMENT TO APPEASE TREETIME
     ###########################################################################
-    if params.aln is None:
     from Bio import Seq, SeqRecord, Align
-        aln = Align.MultipleSeqAlignment([SeqRecord.SeqRecord(Seq.Seq("AAA"), id=name, name=name)
+    aln = Align.MultipleSeqAlignment([SeqRecord.SeqRecord(Seq.Seq("AAA"), id=node, name=node)
                                     for node in dates])
 
 
@@ -57,18 +59,40 @@ if __name__=="__main__":
     ###########################################################################
     base_name = '.'.join(params.tree.split('/')[-1].split('.')[:-1])
     myTree = TreeTime(dates=dates, tree=params.tree,
-                      aln=params.aln, gtr='JC69', verbose=params.verbose)
-    myTree.one_mutation=0.001
+                      aln=aln, gtr='JC69', verbose=params.verbose)
 
-    if params.reroot:
+    if not params.keep_root:
         myTree.reroot('best')
+
+    d2d = DateConversion.from_tree(myTree.tree)
+    print('\n',d2d)
+    print('The R^2 value indicates the fraction of variation in'
+          '\nroot-to-tip distance explained by the sampling times.'
+          '\nHigher values corresponds more clock-like behavior (max 1.0).')
+
+    print('\nThe rate is the slope of the best fit of the date to'
+          '\nthe root-to-tip distance and provides an estimate of'
+          '\nthe substitution rate. The rate needs to be positive!'
+          '\nNegative rates suggest an inappropriate root.\n\n')
+
+    if not params.keep_root:
         # write rerooted tree to file
         outtree_name = base_name+'_rerooted.newick'
         Phylo.write(myTree.tree, outtree_name, 'newick')
+        print("--- re-rooted tree written to \n\t %s\n"%outtree_name)
 
-
-    d2d = DateConversion.from_tree(myTree.tree)
-    print('\n',d2d)
+    table_fname = base_name+'_rtt.csv'
+    with open(table_fname, 'w') as ofile:
+        ofile.write("#name, date, root-to-tip distance\n")
+        ofile.write("#Dates of nodes that didn't have a specified date are inferred from the root-to-tip regression.\n")
+        for n in myTree.tree.get_terminals():
+            if hasattr(n, "numdate_given"):
+                ofile.write("%s, %f, %f\n"%(n.name, n.numdate_given, n.dist2root))
+            else:
+                ofile.write("%s, %f, %f\n"%(n.name, d2d.numdate_from_dist2root(n.dist2root), n.dist2root))
+        for n in myTree.tree.get_nonterminals(order='preorder'):
+            ofile.write("%s, %f, %f\n"%(n.name, d2d.numdate_from_dist2root(n.dist2root), n.dist2root))
+        print("--- wrote dates and root-to-tip distances to \n\t %s\n"%table_fname)
 
 
     ###########################################################################
@@ -78,10 +102,11 @@ if __name__=="__main__":
         import matplotlib.pyplot as plt
         myTree.plot_root_to_tip(label=False)
         t = np.array([np.min(dates.values()), np.max(dates.values())])
-        plt.plot(t, t*d2d.slope+d2d.intercept,
-                 label='y=%1.4f+%1.5ft, r^2=%1.2f'%(d2d.intercept, d2d.slope, d2d.r_val**2))
+        plt.plot(t, t*d2d.clock_rate+d2d.intercept,
+                 label='y=%1.4f+%1.5ft, r^2=%1.2f'%(d2d.intercept, d2d.clock_rate, d2d.r_val**2))
         plt.legend(loc=2)
         plt.savefig(base_name+'_root_to_tip_regression.pdf')
+        print("--- root-to-tip plot saved to  \n\t %s_root_to_tip_regression.pdf\n"%base_name)
 
 
 

timetree_inference.py → bin/timetree_inference.py

 #!/usr/bin/env python
 from __future__ import print_function, division
 import numpy as np
-from treetime import TreeTime
+from treetime import TreeTime, GTR
 from Bio import Phylo, AlignIO
+from Bio import __version__ as bioversion
 
 if __name__=="__main__":
     ###########################################################################
@@ -10,37 +11,48 @@ if __name__=="__main__":
     ###########################################################################
     import argparse
     parser = argparse.ArgumentParser(
-            description="Reconstructs ancestral sequences and infers a molecular clock tree."
-                        ' The output consists of a file ending with _ancestral.fasta with ancestral sequences'
-                        ' and a tree ending with _mutation.newick with mutations appended to node names'
-                        ' as _A45G_.... The branches of this tree are scaled such that branch length'
-                        ' correspond to times in units of the molecular clock. The molecular clock,'
-                        ' along with the inferred GTR model, is written to stdout')
+            description=\
+"Reconstructs ancestral sequences and infers a molecular clock tree. The"\
+" script produces an alignment file ending on _ancestral.fasta which contains"\
+" the inferred ancestral sequences and a tree file ending on _timetree.nexus."\
+" Inferred mutations are included as comments. The molecular clock, along with the inferred"\
+" GTR model, is written to stdout)")
     parser.add_argument('--aln', required = True, type = str,  help ="fasta file with input sequences")
-    parser.add_argument('--tree', required = True, type = str,  help ="newick file with tree")
     parser.add_argument('--dates', required = True, type = str,
                         help ="csv with dates for nodes with 'node_name, date' where date is float (as in 2012.15)")
-    parser.add_argument('--infer_gtr', default = True, action='store_true', help='infer substitution model')
+    # parser.add_argument('--infer_gtr', default = True, action='store_true', help='infer substitution model')
+    parser.add_argument('--tree', type = str,  help ="newick file with tree")
+    parser.add_argument('--gtr', default='infer', type = str, help="GTR model to use. "
+        " Type 'infer' to infer the model from the data. Or, specify the model type. "
+        "Optionally, feed the arguments with the '--gtr_args' option")
+    parser.add_argument('--gtr_params', type=str, nargs='+', help="GTR parameters for the model "
+        "specified by the --gtr argument. The parameters should be feed as 'key=value' list of parameters. "
+        "Example: '--gtr K80 --gtr_params kappa=0.2 pis=0.25,0.25,0.25,0.25'. See the exact definitions of "
+        " the parameters in the GTR creation methods.")
+
     parser.add_argument('--reroot', required = False, type = str, default='best',
                         help ="reroot the tree. Valid arguments are 'best', 'midpoint', or a node name")
-    parser.add_argument('--resolve_polytomies', default = True, action='store_true',
-                        help='resolve polytomies using temporal information.')
+    parser.add_argument('--optimize_branch_length', default = False, action='store_true',
+                        help="Reoptimize branch length. Note that branch length optimized by treetime are only accurate at short evolutionary distances.")
+    parser.add_argument('--keep_polytomies', default = False, action='store_true',
+                        help="Don't resolve polytomies using temporal information.")
     parser.add_argument('--relax',nargs='*', default = False,
                         help='use an autocorrelated molecular clock. Prior strength and coupling of parent '
                              'and offspring rates can be specified e.g. as --relax 1.0 0.5')
     parser.add_argument('--max_iter', default = 2, type=int,
-                        help='maximal number of iterations the inference cycle is run')
-    parser.add_argument('--verbose', default = 3, type=int,
+                        help='maximal number of iterations the inference cycle is run. Note that for polytomy resolution and coalescence models max_iter should be at least 2')
+    parser.add_argument('--verbose', default = 1, type=int,
                         help='verbosity of output 0-6')
-    parser.add_argument('--Tc', default = 0.0, type=float,
+    parser.add_argument('--Tc', default = "0.0", type=str,
                         help='coalescent time scale -- sensible values are on the order of the average '
-                             'hamming distance of contemporaneous sequences')
+                             'hamming distance of contemporaneous sequences. In addition, "opt" '
+                             '"skyline" are valid options and estimate a constant coalescent rate'
+                             'or a piecewise linear coalescent rate history')
     parser.add_argument('--plot', default = False, action='store_true',
-                        help='plot the tree on a time axis')
+                        help='plot the tree on a time axis and save as pdf')
     params = parser.parse_args()
     if params.relax==[]:
         params.relax=True
-    if params.Tc<1e-5: params.Tc=False
 
     ###########################################################################
     ### PARSING DATES
@@ -54,28 +66,97 @@ if __name__=="__main__":
             except:
                 continue
 
+
+    ###########################################################################
+    ### CHECK FOR TREE, build if not in place
+    ###########################################################################
+    if params.tree is None:
+        from treetime.utils import tree_inference
+        import os,shutil
+        params.tree = os.path.basename(params.aln)+'.nwk'
+        print("No tree given: inferring tree")
+        tmp_dir = 'timetree_inference_tmp_files'
+        tree_inference(params.aln, params.tree, tmp_dir = tmp_dir)
+        if os.path.isdir(tmp_dir):
+            shutil.rmtree(tmp_dir)
+
+
+    ###########################################################################
+    ### GTR SET-UP
+    ###########################################################################
+    model = params.gtr
+    gtr_params = params.gtr_params
+    if model == 'infer':
+        gtr = GTR.standard('jc')
+        infer_gtr = True
+    else:
+        try:
+            kwargs = {}
+            if gtr_params is not None:
+                for param in gtr_params:
+                    keyval = param.split('=')
+                    if len(keyval)!=2: continue
+                    if keyval[0] in ['pis', 'pi', 'Pi', 'Pis']:
+                        keyval[1] = map(float, keyval[1].split(','))
+                    elif keyval[0] not in ['alphabet']:
+                        keyval[1] = float(keyval[1])
+                    kwargs[keyval[0]] = keyval[1]
+            else:
+                print ("GTR params are not specified. Creating GTR model with default parameters")
+
+
+            gtr = GTR.standard(model, **kwargs)
+            infer_gtr = False
+        except:
+            print ("Could not create GTR model from input arguments. Using default (Jukes-Cantor 1969)")
+            gtr = GTR.standard('jc')
+            infer_gtr = False
+
+
+
+    ###########################################################################
+    # PARSING OPTIONS
+    ###########################################################################
+    try:
+        Tc = float(params.Tc)
+        if Tc<1e-5:
+            Tc = None
+    except:
+        if params.Tc in ['opt', 'skyline']:
+            Tc = params.Tc
+        else:
+            Tc = None
+
     ###########################################################################
-    ### ANCESTRAL RECONSTRUCTION AND SET-UP
+    ### SET-UP and RUN
     ###########################################################################
     myTree = TreeTime(dates=dates, tree=params.tree,
-                       aln=params.aln, gtr='JC69', verbose=params.verbose)
+                       aln=params.aln, gtr=gtr, verbose=params.verbose)
     myTree.run(root=params.reroot, relaxed_clock=params.relax,
-               resolve_polytomies=params.resolve_polytomies,
-               Tc=params.Tc, max_iter=params.max_iter)
+               resolve_polytomies=(not params.keep_polytomies),
+               Tc=Tc, max_iter=params.max_iter,
+               use_input_branch_length = (not params.optimize_branch_length))
 
     ###########################################################################
     ### OUTPUT and saving of results
     ###########################################################################
-    if params.infer_gtr:
+    if infer_gtr:
         print('\nInferred GTR model:')
         print(myTree.gtr)
 
     print(myTree.date2dist)
-    base_name = '.'.join(params.aln.split('/')[-1].split('.')[:-1])
 
+    if Tc=='skyline':
+        skyline = myTree.merger_model.skyline_inferred(gen=50)
+        print("inferred skyline assuming 50 generations per year:")
+        for (x,y) in zip(skyline.x, skyline.y):
+            print("%1.3f\t%1.3f"%(x,y))
+
+
+    base_name = '.'.join(params.aln.split('/')[-1].split('.')[:-1])
     # plot
     if params.plot:
-        from treetime.io import plot_vs_years
+        from treetime.treetime import plot_vs_years
         import matplotlib.pyplot as plt
         plt.ion()
         leaf_count = myTree.tree.count_terminals()
@@ -84,6 +165,7 @@ if __name__=="__main__":
                                        +('...' if  len(x.mutations)>10 else '')) if leaf_count<30 else ''
         plot_vs_years(myTree, show_confidence=False, label_func = label_func) #, branch_labels=branch_label_func)
         plt.savefig(base_name+'_tree.pdf')
+        print("--- saved tree as pdf in \n\t %s\n"%(base_name+'_tree.pdf'))
     else:
         # convert branch length to years (this is implicit in the above plot)
         myTree.branch_length_to_years()
@@ -91,14 +173,23 @@ if __name__=="__main__":
     # decorate tree with inferred mutations
     outaln_name = base_name+'_ancestral.fasta'
     AlignIO.write(myTree.get_reconstructed_alignment(), outaln_name, 'fasta')
+    print("--- alignment including ancestral nodes saved as  \n\t %s\n"%outaln_name)
+
+    terminal_count = 0
     for n in myTree.tree.find_clades():
         if n.up is None:
             continue
+        n.confidence=None
+        # due to a bug in older versions of biopython that truncated filenames in nexus export
+        # we truncate them by hand and make them unique.
+        if n.is_terminal() and len(n.name)>40 and bioversion<"1.69":
+            n.name = n.name[:35]+'_%03d'%terminal_count
+            terminal_count+=1
         if len(n.mutations):
-            n.name+='_'+'_'.join([a+str(pos)+d for (a,pos, d) in n.mutations])
+            n.comment= '&mutations="' + '_'.join([a+str(pos)+d for (a,pos, d) in n.mutations])+'"'
 
-    # write tree to file. Branch length will now be scaled such that node
-    # positions correspond to sampling times.
-    outtree_name = base_name+'_timetree.newick'
-    Phylo.write(myTree.tree, outtree_name, 'newick')
+    # write tree to file
+    outtree_name = '.'.join(params.tree.split('/')[-1].split('.')[:-1])+'_timetree.nexus'
+    Phylo.write(myTree.tree, outtree_name, 'nexus')
 
+    print("--- tree saved in nexus format as  \n\t %s\n"%outtree_name)

data/H3N2_NA_allyears_NA.20.metadata.csv

-name,numdate_given
+name, date
 A/Hawaii/02/2013|KF789866|05/28/2013|USA|12_13|H3N2/1-1409,2013.40520192
 A/Boston/DOA2_107/2012|CY148382|11/01/2012|USA|12_13|H3N2/1-1409,2012.83778234
 A/Oregon/15/2009|GQ895004|06/25/2009|USA|08_09|H3N2/1-1409,2009.48186174

data/H3N2_NA_allyears_NA.200.metadata.csv

-name,numdate_given
+name, date
 A/Texas/UR06_0358/2007|CY027565|02/20/2007|USA|06_07|H3N2/5-1413, 2007.13963039
 A/Western_Australia/8/2000|CY015662|10/17/2000|Australia||H3N2/1-1409, 2000.79671458
 A/Australia/NHRC0008/2003|CY091279|07/08/2003|Australia||H3N2/4-1412, 2003.5174538

data/Zika_country_weights.csv

+#country,weight
+american_samoa,55
+brazil,200000
+china,1000000
+colombia,50000
+cuba,12000
+dominican_republic,10000
+ecuador,16000
+french_guiana,280
+french_polynesia,276
+guadeloupe,400
+guatemala,17000
+haiti,10000
+honduras,9000
+martinique,385
+mexico,124000
+panama,4000
+puerto_rico,3400
+suriname,558
\ No newline at end of file

data/Zika_metadata.csv

+#name,date,country
+PuertoRico/ZF8/2016|PuertoRico/ZF8/2016|2016-06-13|puerto_rico,2016-06-13,puerto_rico
+PRVABC59|KX601168|2015-12-01|puerto_rico,2015-12-01,puerto_rico
+Martinique/ZF1/2016|Martinique/ZF1/2016|2016-03-07|martinique,2016-03-07,martinique
+NL00013|KU937936|2016-02-11|suriname,2016-02-11,suriname
+V17271|KU758877|2015-12-15|french_guiana,2015-12-15,french_guiana
+BeH819015|KU365778|2015-07-01|brazil,2015-07-01,brazil
+GZ01|KU820898|2016-02-14|china,2016-02-14,china
+GD01|KU740184|2016-02-15|china,2016-02-15,china
+GDZ16001|KU761564|2016-02-12|china,2016-02-12,china
+Z16019|KU955590|2016-02-26|china,2016-02-26,china
+GZ02/2016|KX056898|2016-02-25|china,2016-02-25,china
+Z1106033|KU312312|2015-10-02|suriname,2015-10-02,suriname
+BeH815744|KU365780|2015-07-01|brazil,2015-07-01,brazil
+BeH818995|KU365777|2015-07-01|brazil,2015-07-01,brazil
+BeH819966|KU365779|2015-07-01|brazil,2015-07-01,brazil
+SSABR1|KU707826|2015-07-01|brazil,2015-07-01,brazil
+Aedessp/MEX/MEX_2_81/2016|KX446950|2016-01-01|mexico,2016-01-01,mexico
+Aedessp/MEX/MEX_I_7/2016|KX446951|2016-01-01|mexico,2016-01-01,mexico
+MEX_I_7|KX247632|2015-11-15|mexico,2015-11-15,mexico
+FB_GWUH_2016|KU870645|2016-02-02|guatemala,2016-02-02,guatemala
+V103344|KU501216|2015-12-01|guatemala,2015-12-01,guatemala
+V8375|KU501217|2015-11-01|guatemala,2015-11-01,guatemala
+HND/R103451/2015|KX694534|2016-07-01|honduras,2016-07-01,honduras
+V103451|KX262887|2016-01-06|honduras,2016-01-06,honduras
+Rio_U1|KU926309|2016-01-14|brazil,2016-01-14,brazil
+PAN/CDC_259249_V1_V3/2015|KX156775|2015-12-11|panama,2015-12-11,panama
+PAN/CDC_259359_V1_V3/2015|KX156774|2015-12-18|panama,2015-12-18,panama
+PAN/CDC_259364_V1_V2/2015|KX156776|2015-12-18|panama,2015-12-18,panama
+FLR|KU820897|2015-12-15|colombia,2015-12-15,colombia
+COL/FLR/2015|KX087102|2015-12-15|colombia,2015-12-15,colombia
+COL/UF_1/2016|KX247646|2016-02-09|colombia,2016-02-09,colombia
+PAN/BEI_259634_V4/2016|KX198135|2016-07-01|panama,2016-07-01,panama
+COL/FCC00093/2015|KX548902|2015-10-07|colombia,2015-10-07,colombia
+MEX/InDRE/Sm/2016|KU922960|2016-02-25|mexico,2016-02-25,mexico
+MRS_OPY_Martinique_PaRi_2015|KU647676|2015-12-15|martinique,2015-12-15,martinique
+Brazil/PE243/2015|KX197192|2015-07-01|brazil,2015-07-01,brazil
+SPH2015|KU321639|2015-03-15|brazil,2015-03-15,brazil
+Haiti/1225/2014|KU509998|2014-12-12|haiti,2014-12-12,haiti
+Haiti/1/2016|KX051563|2016-02-05|haiti,2016-02-05,haiti
+Brazil_ZKV2015|KU497555|2015-11-30|brazil,2015-11-30,brazil
+Brazil/2016/INMI1|KU991811|2016-03-06|brazil,2016-03-06,brazil
+PHE_semen_Guadeloupe|KX673530|2016-04-21|guadeloupe,2016-04-21,guadeloupe
+Cuba/ZF10/2016|Cuba/ZF10/2016|2016-06-02|cuba,2016-06-02,cuba
+Dominican_Republic/2016/PD1|KU853012|2016-02-01|dominican_republic,2016-02-01,dominican_republic
+BeH823339|KU729217|2015-07-01|brazil,2015-07-01,brazil
+Rio_S1|KU926310|2016-01-29|brazil,2016-01-29,brazil
+Natal_RGN|KU527068|2015-07-01|brazil,2015-07-01,brazil
+Ecuador/EC062/2016|Ecuador/EC062/2016|2016-04-15|ecuador,2016-04-15,ecuador
+Ecuador/EC089/2016|Ecuador/EC089/2016|2016-04-15|ecuador,2016-04-15,ecuador
+Paraiba_01|KX280026|2015-07-01|brazil,2015-07-01,brazil
+Bahia04|KX101062|2015-06-15|brazil,2015-06-15,brazil
+Bahia08|KU940227|2015-07-15|brazil,2015-07-15,brazil
+Bahia05|KX101063|2015-12-15|brazil,2015-12-15,brazil
+HS_2015_BA_01|KX520666|2015-07-01|brazil,2015-07-01,brazil
+Bahia09|KU940224|2015-08-01|brazil,2015-08-01,brazil
+Bahia12|KX101067|2015-05-15|brazil,2015-05-15,brazil
+Bahia07|KU940228|2015-07-01|brazil,2015-07-01,brazil
+Bahia03|KX101061|2015-05-15|brazil,2015-05-15,brazil
+Bahia11|KX101064|2015-04-15|brazil,2015-04-15,brazil
+Bahia15|KX101065|2016-01-15|brazil,2016-01-15,brazil
+Bahia01|KX101066|2015-05-15|brazil,2015-05-15,brazil
+Bahia02|KX101060|2015-05-15|brazil,2015-05-15,brazil
+BeH828305|KU729218|2015-07-01|brazil,2015-07-01,brazil
+1_0049_PF|KX447510|2013-12-15|french_polynesia,2013-12-15,french_polynesia
+1_0117_PF|KX447518|2013-12-15|french_polynesia,2013-12-15,french_polynesia
+1_0016_PF|KX447520|2014-01-15|french_polynesia,2014-01-15,french_polynesia
+PF13/251013_18|KX369547|2013-10-25|french_polynesia,2013-10-25,french_polynesia
+1_0181_PF|KX447512|2013-12-15|french_polynesia,2013-12-15,french_polynesia
+1_0134_PF|KX447513|2013-12-15|french_polynesia,2013-12-15,french_polynesia
+1_0199_PF|KX447519|2013-11-15|french_polynesia,2013-11-15,french_polynesia
+1_0015_PF|KX447511|2014-01-15|french_polynesia,2014-01-15,french_polynesia
+SZ01/2016/China|KU866423|2016-07-01|american_samoa,2016-07-01,american_samoa
+Zhejiang04|KX117076|2016-02-17|china,2016-02-17,china
+SZ_WIV01|KU963796|2016-07-01|american_samoa,2016-07-01,american_samoa
+ZJ03|KU820899|2016-02-17|china,2016-02-17,china
+CN/SZ02/2016|KX185891|2016-02-17|american_samoa,2016-02-17,american_samoa
+ZKC2/2016|KX253996|2016-02-16|china,2016-02-16,china
+Z16006|KU955589|2016-02-16|china,2016-02-16,china
+SMGC_1|KX266255|2016-02-14|china,2016-02-14,china
+1_0035_PF|KX447514|2014-01-15|french_polynesia,2014-01-15,french_polynesia
+1_0080_PF|KX447521|2014-02-15|french_polynesia,2014-02-15,french_polynesia
+H/PF/2013|KJ776791|2013-11-28|french_polynesia,2013-11-28,french_polynesia
+1_0087_PF|KX447509|2013-12-15|french_polynesia,2013-12-15,french_polynesia
+1_0111_PF|KX447516|2014-01-15|french_polynesia,2014-01-15,french_polynesia
+1_0030_PF|KX447515|2013-11-15|french_polynesia,2013-11-15,french_polynesia
+1_0038_PF|KX447517|2014-01-15|french_polynesia,2014-01-15,french_polynesia

data/Zika_tree.newick

+(((((((((((('PuertoRico/ZF8/2016|PuertoRico/ZF8/2016|2016-06-13|puerto_rico':0.00087,'PRVABC59|KX601168|2015-12-01|puerto_rico':0.0):0.00028,'Martinique/ZF1/2016|Martinique/ZF1/2016|2016-03-07|martinique':0.00068):0.00095,('NL00013|KU937936|2016-02-11|suriname':0.00096,'V17271|KU758877|2015-12-15|french_guiana':0.00058):0.00019,'BeH819015|KU365778|2015-07-01|brazil':0.00113):0.00019,(('GZ01|KU820898|2016-02-14|china':0.00023,('GD01|KU740184|2016-02-15|china':0.0,'GDZ16001|KU761564|2016-02-12|china':0.0):0.00015):0.00067,('Z16019|KU955590|2016-02-26|china':0.00019,'GZ02/2016|KX056898|2016-02-25|china':0.0):0.00073):0.00109):9e-05,'Z1106033|KU312312|2015-10-02|suriname':0.00127):0.00019,(('BeH815744|KU365780|2015-07-01|brazil':0.00019,'BeH818995|KU365777|2015-07-01|brazil':0.0):0.00038,'BeH819966|KU365779|2015-07-01|brazil':0.0,'SSABR1|KU707826|2015-07-01|brazil':0.0):0.00085):0.00047,((('Aedessp/MEX/MEX_2_81/2016|KX446950|2016-01-01|mexico':0.00056,('Aedessp/MEX/MEX_I_7/2016|KX446951|2016-01-01|mexico':0.00019,'MEX_I_7|KX247632|2015-11-15|mexico':0.0):0.00037):0.00075,'FB_GWUH_2016|KU870645|2016-02-02|guatemala':0.00178,('V103344|KU501216|2015-12-01|guatemala':0.0001,'V8375|KU501217|2015-11-01|guatemala':0.0):0.00088,('HND/R103451/2015|KX694534|2016-07-01|honduras':0.00028,'V103451|KX262887|2016-01-06|honduras':9e-05):0.00047):0.00075,'Rio_U1|KU926309|2016-01-14|brazil':0.00141):9e-05):9e-05,(((((('PAN/CDC_259249_V1_V3/2015|KX156775|2015-12-11|panama':0.00028,'PAN/CDC_259359_V1_V3/2015|KX156774|2015-12-18|panama':0.00024,'PAN/CDC_259364_V1_V2/2015|KX156776|2015-12-18|panama':9e-05):0.00061,(('FLR|KU820897|2015-12-15|colombia':0.0,'COL/FLR/2015|KX087102|2015-12-15|colombia':0.0):0.00022,'COL/UF_1/2016|KX247646|2016-02-09|colombia':0.00015):0.00056):0.00018,'PAN/BEI_259634_V4/2016|KX198135|2016-07-01|panama':0.00057):0.00019,('COL/FCC00093/2015|KX548902|2015-10-07|colombia':0.00124,('MEX/InDRE/Sm/2016|KU922960|2016-02-25|mexico':0.00038,'MRS_OPY_Martinique_PaRi_2015|KU647676|2015-12-15|martinique':0.0):0.00067):0.0001):0.0008,((('Brazil/PE243/2015|KX197192|2015-07-01|brazil':0.00093,'SPH2015|KU321639|2015-03-15|brazil':0.0):0.00046,'Haiti/1225/2014|KU509998|2014-12-12|haiti':0.0):0.00028,'Haiti/1/2016|KX051563|2016-02-05|haiti':0.00075):0.00075):0.00027,('Brazil_ZKV2015|KU497555|2015-11-30|brazil':0.00212,'Brazil/2016/INMI1|KU991811|2016-03-06|brazil':0.00095):0.00017):0.00049,((((('PHE_semen_Guadeloupe|KX673530|2016-04-21|guadeloupe':0.00114,'Cuba/ZF10/2016|Cuba/ZF10/2016|2016-06-02|cuba':0.00087):9e-05,'Dominican_Republic/2016/PD1|KU853012|2016-02-01|dominican_republic':0.00048):0.00133,'BeH823339|KU729217|2015-07-01|brazil':0.00144):0.00018,('Rio_S1|KU926310|2016-01-29|brazil':0.00188,'Natal_RGN|KU527068|2015-07-01|brazil':0.00159):9e-05):9e-05,(('Ecuador/EC062/2016|Ecuador/EC062/2016|2016-04-15|ecuador':0.0002,'Ecuador/EC089/2016|Ecuador/EC089/2016|2016-04-15|ecuador':0.0001):0.00134,'Paraiba_01|KX280026|2015-07-01|brazil':0.00084):9e-05):0.00019,(((((('Bahia04|KX101062|2015-06-15|brazil':0.00068,'Bahia08|KU940227|2015-07-15|brazil':0.0,'Bahia05|KX101063|2015-12-15|brazil':0.0):0.00033,(('HS_2015_BA_01|KX520666|2015-07-01|brazil':0.00048,'Bahia09|KU940224|2015-08-01|brazil':0.0):0.00016,'Bahia12|KX101067|2015-05-15|brazil':0.0,'Bahia07|KU940228|2015-07-01|brazil':0.0):0.00016,'Bahia03|KX101061|2015-05-15|brazil':0.00067):0.00021,'Bahia11|KX101064|2015-04-15|brazil':0.00029,'Bahia15|KX101065|2016-01-15|brazil':0.0):0.00021,'Bahia01|KX101066|2015-05-15|brazil':0.00115):0.00026,'Bahia02|KX101060|2015-05-15|brazil':0.00046):0.00115,'BeH828305|KU729218|2015-07-01|brazil':0.00123):9e-05):0.00037,'1_0049_PF|KX447510|2013-12-15|french_polynesia':0.00038,'1_0117_PF|KX447518|2013-12-15|french_polynesia':0.00033):0.00028,'1_0016_PF|KX447520|2014-01-15|french_polynesia':0.00102,('PF13/251013_18|KX369547|2013-10-25|french_polynesia':0.00056,'1_0181_PF|KX447512|2013-12-15|french_polynesia':0.00019):9e-05,('1_0134_PF|KX447513|2013-12-15|french_polynesia':0.00038,'1_0199_PF|KX447519|2013-11-15|french_polynesia':0.00011):9e-05,'1_0015_PF|KX447511|2014-01-15|french_polynesia':0.00057):9e-05,(('SZ01/2016/China|KU866423|2016-07-01|american_samoa':0.0002,'Zhejiang04|KX117076|2016-02-17|china':0.00019,'SZ_WIV01|KU963796|2016-07-01|american_samoa':0.00019,'ZJ03|KU820899|2016-02-17|china':9e-05,'CN/SZ02/2016|KX185891|2016-02-17|american_samoa':0.0):0.00047,('ZKC2/2016|KX253996|2016-02-16|china':0.0,'Z16006|KU955589|2016-02-16|china':0.0,'SMGC_1|KX266255|2016-02-14|china':0.0):0.00018):0.00235,('1_0035_PF|KX447514|2014-01-15|french_polynesia':0.00048,'1_0080_PF|KX447521|2014-02-15|french_polynesia':0.00047,'H/PF/2013|KJ776791|2013-11-28|french_polynesia':0.00028,'1_0087_PF|KX447509|2013-12-15|french_polynesia':0.00019):9e-05,'1_0111_PF|KX447516|2014-01-15|french_polynesia':0.00066,'1_0030_PF|KX447515|2013-11-15|french_polynesia':0.00048):0.00015,'1_0038_PF|KX447517|2014-01-15|french_polynesia':0.00118):0.001;

data/ebola.csv → data/ebola.metadata.csv

+name, date
 EM_COY_2015_015982, 2015.30
 G3676, 2014.40
 EM_COY_2015_015980, 2015.30

debian/TODO

+Replace
+
+    Testsuite: autopkgtest-pkg-python
+
+by a real test suite running build time tests

debian/changelog

+python-treetime (0.2.4-1) unstable; urgency=medium
+
+  * d/watch: Point to github since upstream has started using release tags
+  * New upstream version
+  * Point Vcs fields to salsa.debian.org
+  * Standards-Version: 4.1.4
+  * debhelper 11
+  * Testsuite: autopkgtest-pkg-python
+  * Delete get-orig-source target
+
+ -- Andreas Tille <tille@debian.org>  Mon, 16 Apr 2018 14:47:01 +0200
+
 python-treetime (0.0+20170607-1) unstable; urgency=medium
 
   * Initial release (Closes: #864823)

debian/compat

-10
+11

debian/control

@@ -2,14 +2,15 @@ Source: python-treetime
 Maintainer: Debian Med Packaging Team <debian-med-packaging@lists.alioth.debian.org>
 Uploaders: Andreas Tille <tille@debian.org>
 Section: science
+Testsuite: autopkgtest-pkg-python
 Priority: optional
-Build-Depends: debhelper (>= 10),
+Build-Depends: debhelper (>= 11~),
                dh-python,
                python-all,
                python-setuptools
-Standards-Version: 3.9.8
-Vcs-Browser: https://anonscm.debian.org/cgit/debian-med/python-treetime.git
-Vcs-Git: https://anonscm.debian.org/git/debian-med/python-treetime.git
+Standards-Version: 4.1.4
+Vcs-Browser: https://salsa.debian.org/med-team/python-treetime
+Vcs-Git: https://salsa.debian.org/med-team/python-treetime.git
 Homepage: https://github.com/neherlab/treetime
 X-Python-Version: >= 2.6
 

debian/get-orig-source

-#!/bin/sh -e
-
-COMPRESS=xz
-
-NAME=`dpkg-parsechangelog | awk '/^Source/ { print $2 }'`
-MVERSION=`dpkg-parsechangelog | awk '/^Version:/ { print $2 }' | sed 's/^\([0-9\.]\+\)[+~][-0-9]\+$/\1/'`
-
-mkdir -p ../tarballs
-cd ../tarballs
-# need to clean up the tarballs dir first because upstream tarball might
-# contain a directory with unpredictable name
-rm -rf *
-git clone --quiet https://github.com/neherlab/treetime $NAME
-cd $NAME
-VERSION=${MVERSION}+`date -d @$(git show --format="%at" | head -n1) +%Y%m%d`
-# for esthetical reasons set file timestamps (if git-restore-mtime is installed)
-git restore-mtime || true
-cd ..
-TARDIR=${NAME}-${VERSION}
-mv ${NAME} ${TARDIR}
-rm -rf ${TARDIR}/.git
-
-GZIP="--best --no-name" tar --owner=root --group=root --mode=a+rX -caf "$NAME"_"$VERSION".orig.tar.${COMPRESS} "${TARDIR}"
-rm -rf ${TARDIR}