Use prokka to annotate each strain
332298
def annotatethreads(self):
        
        # Move the files to subfolders and create objects
        self.runmetadata = createobject.ObjectCreation(self)
        # Fix headers
        self.headers()
        printtime('Performing prokka analyses', self.start)
        # Create and start threads
        for i in range(self.cpus):
            # Send the threads to the appropriate destination function
            threads = Thread(target=self.annotate, args=())
            # Set the daemon to true - something to do with thread management
            threads.setDaemon(True)
            # Start the threading
            threads.start()
        for sample in self.metadata.samples:
            # Create the prokka attribute in the metadata object
            setattr(sample, 'prokka', GenObject())
            sample.prokka.outputdir = os.path.join(sample.general.outputdirectory, 'prokka')
            if not os.path.isdir(sample.prokka.outputdir):
                os.makedirs(sample.prokka.outputdir)
            # TODO Incorporate MASH/rMLST/user inputted genus, species results in the system call
            # Create the system call
            # prokka 2014-SEQ-0275.fasta --force --genus Escherichia --species coli --usegenus --addgenes
            # --prefix 2014-SEQ-0275 --locustag EC0275 --outputdir /path/to/sequences/2014-SEQ-0275/prokka
            sample.prokka.command = 'prokka {} ' \
                                    '--force ' \
                                    '--genus {} ' \
                                    '--species {} ' \
                                    '--usegenus ' \
                                    '--addgenes ' \
                                    '--prefix {} ' \
                                    '--locustag {} ' \
                                    '--outdir {}' \
                .format(sample.general.fixedheaders,
                        self.genus, self.species, sample.name, sample.name, sample.prokka.outputdir)
            self.queue.put(sample)
        self.queue.join()
458327
def map_all(self, prot_alignment, nucl_sequences):
        
        zipped = itertools.zip_longest(prot_alignment, nucl_sequences)
        for p, n in zipped:
            if p is None:
                raise ValueError("Exhausted protein sequences")
            elif n is None:
                raise ValueError("Exhausted nucleotide sequences")
            yield self.map_alignment(p, n)
415252
def create_column_index(annotations):
    
    _column_index = OrderedDict({'Column Name' : annotations['Column Name']})
    categorical_rows = annotation_rows('C:', annotations)
    _column_index.update(categorical_rows)
    numerical_rows = {name: [float(x) if x != '' else float('NaN') for x in values]
            for name, values in annotation_rows('N:', annotations).items()} # to floats
    _column_index.update(numerical_rows)
    column_index = pd.MultiIndex.from_tuples(list(zip(*_column_index.values())), names=list(_column_index.keys()))
    if len(column_index.names) == 1:
        # flatten single-level index
        name = column_index.names[0]
        column_index = column_index.get_level_values(name)
    return column_index
709498
def annotatedcore(self):
        
        logging.info('Calculating annotated core')
        # Determine the total number of core genes
        self.total_core()
        # Iterate through all the samples, and process all Escherichia
        for sample in self.metadata:
            if sample.general.bestassemblyfile != 'NA':
                # Create a set to store the names of all the core genes in this strain
                sample[self.analysistype].coreset = set()
                if sample.general.referencegenus == 'Escherichia':
                    # Add the Escherichia sample to the runmetadata
                    self.runmetadata.samples.append(sample)
                    # Parse the BLAST report
                    try:
                        report = sample[self.analysistype].report
                        self.blastparser(report=report,
                                         sample=sample,
                                         fieldnames=self.fieldnames)
                    except KeyError:
                        sample[self.analysistype].coreset = list()
        # Create the report
        self.reporter()
2112
def parse(self, xml_file):
        ""

        recipes = []

        with open(xml_file, "rt") as f:
            tree = ElementTree.parse(f)

        for recipeNode in tree.iter():
            if self.to_lower(recipeNode.tag) != "recipe":
                continue

            recipe = Recipe()
            recipes.append(recipe)

            for recipeProperty in list(recipeNode):
                tag_name = self.to_lower(recipeProperty.tag)

                if tag_name == "fermentables":
                    for fermentable_node in list(recipeProperty):
                        fermentable = Fermentable()
                        self.nodes_to_object(fermentable_node, fermentable)
                        recipe.fermentables.append(fermentable)

                elif tag_name == "yeasts":
                    for yeast_node in list(recipeProperty):
                        yeast = Yeast()
                        self.nodes_to_object(yeast_node, yeast)
                        recipe.yeasts.append(yeast)

                elif tag_name == "hops":
                    for hop_node in list(recipeProperty):
                        hop = Hop()
                        self.nodes_to_object(hop_node, hop)
                        recipe.hops.append(hop)

                elif tag_name == "miscs":
                    for misc_node in list(recipeProperty):
                        misc = Misc()
                        self.nodes_to_object(misc_node, misc)
                        recipe.miscs.append(misc)

                elif tag_name == "style":
                    style = Style()
                    recipe.style = style
                    self.nodes_to_object(recipeProperty, style)

                elif tag_name == "mash":

                    for mash_node in list(recipeProperty):
                        mash = Mash()
                        recipe.mash = mash

                        if self.to_lower(mash_node.tag) == "mash_steps":
                            for mash_step_node in list(mash_node):
                                mash_step = MashStep()
                                self.nodes_to_object(mash_step_node, mash_step)
                                mash.steps.append(mash_step)
                        else:
                            self.nodes_to_object(mash_node, mash)

                else:
                    self.node_to_object(recipeProperty, recipe)

        return recipes
1265393
def write_strain_functional_genes(self, force_rerun=False):
        """"""
        if len(self.df_orthology_matrix) == 0:
            raise RuntimeError('Empty orthology matrix, please calculate first!')
        ref_functional_genes = [g.id for g in self.reference_gempro.functional_genes]
        log.info('Building strain specific models...')
        result = []
        for s in tqdm(self.strain_ids):
            result.append(self._write_strain_functional_genes(s, ref_functional_genes, self.df_orthology_matrix, force_rerun=force_rerun))

        for strain_id, functional_genes in result:
            self.strain_infodict[strain_id]['functional_genes'] = functional_genes
1057932
def print_table(seqs, id2name, name):
    
    itable = open('%s.itable' % (name.rsplit('.', 1)[0]), 'w')
    print('\t'.join(['#sequence', 'gene', 'model', 'insertion', 'gene position', 'model position', 'length', 'orf?', 'intron?', 'orf?intron?', 'insertion', 'orf', 'intron']), file=itable)
    for seq, info in list(seqs.items()):
        gene, model, insertions = info
        name = id2name[seq]
        for i, ins in enumerate(insertions, 1):
            gene_pos, model_pos, length, iseq, \
                    orfs, introns, orfs_b, introns_b, orf_annotations = ins
            # append annotation to orf header
            for orf in orfs:
                parts = orf[0].split()
                annotation = orf_annotations[parts[0].split('>')[1]]
                orf[0] = '%s %s %s' % (parts[0], annotation, ' '.join(parts[1:]))
            # get orf position
            gene_pos = '-'.join([str(j) for j in gene_pos])
            # check if orf, intron is present
            if orfs_b is True or introns_b is True:
                orfs_introns_b = True
            else:
                orfs_introns_b = False
            out = [name, gene, model, i, gene_pos, model_pos, length, orfs_b, introns_b, orfs_introns_b]
            out.append('|'.join(iseq))
            out.append('|'.join(['|'.join(orf) for orf in orfs]))
            out.append('|'.join(['|'.join(intron) for intron in introns]))
            print('\t'.join([str(i) for i in out]), file=itable)
    itable.close()
1156483
def _cnvkit_prioritize(sample, genes, allele_file, metrics_file):
    
    mdf = pd.read_table(metrics_file)
    mdf.columns = [x.lower() for x in mdf.columns]
    if len(genes) > 0:
        mdf = mdf[mdf["gene"].str.contains("|".join(genes))]
    mdf = mdf[["chromosome", "start", "end", "gene", "log2", "ci_hi", "ci_lo"]]
    adf = pd.read_table(allele_file)
    if len(genes) > 0:
        adf = adf[adf["gene"].str.contains("|".join(genes))]
    if "cn1" in adf.columns and "cn2" in adf.columns:
        adf = adf[["chromosome", "start", "end", "cn", "cn1", "cn2"]]
    else:
        adf = adf[["chromosome", "start", "end", "cn"]]
    df = pd.merge(mdf, adf, on=["chromosome", "start", "end"])
    df = df[df["cn"] != 2]
    if len(df) > 0:
        def passes(row):
            spread = abs(row["ci_hi"] - row["ci_lo"])
            return spread < 0.25
        df["passes"] = df.apply(passes, axis=1)
    df.insert(0, "sample", [sample] * len(df))
    return df
176247
def load_sequences_to_reference(self, sc=None, force_rerun=False):
        """"""
        log.info('Loading sequences to reference GEM-PRO...')
        from random import shuffle
        g_ids = [g.id for g in self.reference_gempro.functional_genes]
        shuffle(g_ids)

        def _load_sequences_to_reference_gene_sc(g_id, outdir=self.sequences_by_gene_dir,
                                                 g_to_pickle=self.gene_protein_pickles,
                                                 strain_infodict=self.strain_infodict,
                                                 orth_matrix=self.df_orthology_matrix, force_rerun=force_rerun):
            """Load orthologous strain sequences to reference Protein object, save as new pickle"""
            import ssbio.utils
            import ssbio.io
            from Bio import SeqIO
            import os.path as op

            protein_seqs_pickle_path = op.join(outdir, '{}_protein_withseqs.pckl'.format(g_id))

            if ssbio.utils.force_rerun(flag=force_rerun, outfile=protein_seqs_pickle_path):
                protein_pickle_path = g_to_pickle[g_id]
                protein_pickle = ssbio.io.load_pickle(protein_pickle_path)

                for strain, info in strain_infodict.items():
                    strain_sequences = SeqIO.index(info['genome_path'], 'fasta')
                    strain_gene_functional = info['functional_genes'][g_id]
                    if strain_gene_functional:
                        # Pull the gene ID of the strain from the orthology matrix
                        strain_gene_key = orth_matrix.at[g_id, strain]
                        new_id = '{}_{}'.format(g_id, strain)
                        if protein_pickle.sequences.has_id(new_id):
                            continue
                        protein_pickle.load_manual_sequence(seq=strain_sequences[strain_gene_key],
                                                            ident=new_id,
                                                            set_as_representative=False)
                protein_pickle.save_pickle(outfile=protein_seqs_pickle_path)

            return g_id, protein_seqs_pickle_path

        if sc:
            genes_rdd = sc.parallelize(g_ids)
            result = genes_rdd.map(_load_sequences_to_reference_gene_sc).collect()
        else:
            result = []
            for g in tqdm(g_ids):
                result.append(self._load_sequences_to_reference_gene(g, force_rerun))

        log.info('Storing paths to new Protein objects in self.gene_protein_pickles...')
        updated = []
        for g_id, protein_pickle in result:
            self.gene_protein_pickles[g_id] = protein_pickle
            updated.append(g_id)
        not_updated = set(list(self.gene_protein_pickles.keys())).difference(updated)
        log.info('No change to {} genes, removing from gene_protein_pickles'.format(len(not_updated)))
        log.debug(not_updated)
        for rem in not_updated:
            del self.gene_protein_pickles[rem]
8632
def _load_sequences_to_strain(self, strain_id, force_rerun=False):
        """"""
        gp_seqs_path = op.join(self.model_dir, '{}_gp_withseqs.pckl'.format(strain_id))

        if ssbio.utils.force_rerun(flag=force_rerun, outfile=gp_seqs_path):
            gp_noseqs = ssbio.io.load_pickle(self.strain_infodict[strain_id]['gp_noseqs_path'])
            strain_sequences = SeqIO.index(self.strain_infodict[strain_id]['genome_path'], 'fasta')
            for strain_gene in gp_noseqs.functional_genes:
                # Pull the gene ID of the strain from the orthology matrix
                strain_gene_key = self.df_orthology_matrix.at[strain_gene.id, strain_id]
                # Load into the strain GEM-PRO
                new_id = '{}_{}'.format(strain_gene.id, strain_id)
                if strain_gene.protein.sequences.has_id(new_id):
                    continue
                strain_gene.protein.load_manual_sequence(seq=strain_sequences[strain_gene_key], ident=new_id,
                                                         set_as_representative=True)
            gp_noseqs.save_pickle(outfile=gp_seqs_path)

        return strain_id, gp_seqs_path
304717
def tblNomenclator
      @headers = %w{NomenclatorID GenusNameID SubgenusNameID SpeciesNameID SubspeciesNameID LastUpdate ModifiedBy SuitableForGenus SuitableForSpecies InfrasubspeciesNameID InfrasubKind}
      sql = []   
      i = 1

      # Ugh, move build from here 
      @name_collection.nomenclators.keys.each do |i|
        name =  @name_collection.nomenclators[i]
        genus_id = @genus_names[name[0]]
        genus_id ||= 0
        subgenus_id = @genus_names[name[1]]
        subgenus_id ||= 0
        species_id = @species_names[name[2]]
        species_id ||= 0
        subspecies_id = @species_names[name[3]]
        subspecies_id ||= 0
        variety_id = @species_names[name[4]]
        variety_id ||= 0

        cols = {
          NomenclatorID: i,
          GenusNameID: genus_id, 
          SubgenusNameID: subgenus_id, 
          SpeciesNameID: species_id, 
          SubspeciesNameID: subspecies_id,
          InfrasubspeciesNameID: variety_id,
          InfrasubKind: (variety_id == 0 ? 0 : 2), 
          LastUpdate: @time,  
          ModifiedBy: @authorized_user_id, 
          SuitableForGenus: 0,                      # Set in SF w test
          SuitableForSpecies: 0                     # Set in SF w test
        }
        i += 1
        sql << sql_insert_statement('tblNomenclator', cols) 
      end

      sql.join("\n")
    end
646034
def load_strain(self, strain_id, strain_genome_file):
        
        # logging.disable(logging.WARNING)
        strain_gp = GEMPRO(gem_name=strain_id, genome_path=strain_genome_file, write_protein_fasta_files=False)
        # logging.disable(logging.NOTSET)

        self.strains.append(strain_gp)
        return self.strains.get_by_id(strain_id)
1789
def run(self, args, pipeline_command):
        

        # output that must be moved but not renamed
        consistentNaming = ['alignments/normal_dna_fix_pg_sorted.bam',
                           'alignments/normal_dna_fix_pg_sorted.bam.bai',
                           'alignments/rna_genome_sorted.bam',
                           'alignments/rna_genome_sorted.bam.bai',
                           'alignments/rna_transcriptome.bam',
                           'alignments/tumor_dna_fix_pg_sorted.bam',
                           'alignments/tumor_dna_fix_pg_sorted.bam.bai',
                           'mutations/merged/all_merged.vcf',
                           'rankboost/mhcii_rankboost_concise_results.tsv',
                           'rankboost/mhci_rankboost_concise_results.tsv',
                           ]

        # output that must be renamed as well as moved
        # map of the original name to the final name
        renamingNeeded = {'binding_predictions': 'binding_predictions.tar',
                         'expression': 'expression.tar',
                         'haplotyping': 'haplotyping.tar',
                         'peptides': 'peptides.tar',
                         'rankboost': 'rankboost.tar',
                         'reports': 'reports.tar',
                         'mutations/snpeffed/mutations.vcf': 'all_snpeffed.vcf',
                         'mutations/transgened/mutations.vcf': 'all_transgened.vcf',
                         'mutations/merged': 'merged_perchrom.tar',
                         'mutations/muse': 'muse_perchrom.tar',
                         'mutations/mutect': 'mutect_perchrom.tar',
                         'mutations/radia': 'radia_perchrom.tar',
                         'mutations/somaticsniper': 'somaticsniper_perchrom.tar',
                         'mutations/strelka/snv': 'strelka_snv_perchrom.tar',
                         'mutations/strelka/indel': 'strelka_indel_perchrom.tar'}

        def make_output(output_dir, source_dir):
            """
            :param output_dir: dir to write the output to
            :param source_dir: dir containing the directory structure to be parsed
            :return:
            """
            def make_tar(dir, tar):
                with tarfile.open(tar, "w:gz") as tar:
                    tar.add(dir)
            # the output dir is where the real output directories are written
            protect_outputs = os.listdir(source_dir)
            for protectOut in protect_outputs:
                def getName(fileName):
                    return os.path.join(os.path.join(source_dir, protectOut), fileName)
                # move individual files out
                for fileName in consistentNaming:
                    shutil.copyfile(getName(fileName), os.path.join(output_dir, os.path.basename(fileName)))
                for src, dst in renamingNeeded.iteritems():
                    if dst.endswith('.tar'):
                        make_tar(getName(src), os.path.join(output_dir, dst))
                    else:
                        shutil.copyfile(getName(src), os.path.join(output_dir, dst))
            shutil.rmtree(source_dir)


        # prepare workdir
        mount = self._prepare_mount(args)
        self._workdir = os.path.join(mount, 'Toil-' + self._name)
        # insure the pairs are in the same directory, as protect expects
        # This is made more complicated by the fact CWLTool mounts inputs into random, read-only dirs
        # to get around this we copy all inputs into their own directories that we own
        tumor_dna_dir = os.path.expanduser('~/tumorDNA')
        tumor_rna_dir = os.path.expanduser('~/tumorRNA')
        normal_dna_dir = os.path.expanduser('~/normalDNA')
        os.mkdir(tumor_dna_dir)
        os.mkdir(tumor_rna_dir)
        os.mkdir(normal_dna_dir)
        shutil.copy(args.tumor_dna, tumor_dna_dir)
        shutil.copy(args.tumor_rna, tumor_rna_dir)
        shutil.copy(args.normal_dna, normal_dna_dir)
        shutil.copy(args.tumor_dna2, tumor_dna_dir)
        shutil.copy(args.tumor_rna2, tumor_rna_dir)
        shutil.copy(args.normal_dna2, normal_dna_dir)
        args.tumor_dna = os.path.join(tumor_dna_dir, os.path.basename(args.tumor_dna))
        args.tumor_dna2 = os.path.join(tumor_dna_dir, os.path.basename(args.tumor_dna2))
        args.tumor_rna = os.path.join(tumor_rna_dir, os.path.basename(args.tumor_rna))
        args.tumor_rna2 = os.path.join(tumor_rna_dir, os.path.basename(args.tumor_rna2))
        args.normal_dna = os.path.join(normal_dna_dir, os.path.basename(args.normal_dna))
        args.normal_dna2 = os.path.join(normal_dna_dir, os.path.basename(args.normal_dna2))

        # prepare config
        args_dict = vars(args)
        args_dict['output_dir'] = mount
        self._config = textwrap.dedent(self._config.format(**args_dict))
        self._sample_name = args_dict["sample_name"]
        config_path = os.path.join(self._workdir, 'config')
        command = self._make_prefix(os.path.join(self._workdir, 'jobStore'),
                                          config_path,
                                          self._workdir) + pipeline_command
        if self._resume and args.resume:
            command.append('--restart')
        self._create_workdir(args)
        with open(config_path, 'w') as f:
            f.write(self._config)

        try:
            subprocess.check_call(command)
        except subprocess.CalledProcessError as e:
            print(e, file=sys.stderr)
        finally:
            log.info('Pipeline terminated, changing ownership of output files from root to user.')
            stat = os.stat(self._mount)
            subprocess.check_call(['chown', '-R', '{}:{}'.format(stat.st_uid, stat.st_gid),
                                   self._mount])
            make_output(self._mount, os.path.join(self._mount, 'output'))
            if self._no_clean and args.no_clean:
                log.info('Flag "--no-clean" was used, therefore %s was not deleted.', self._workdir)
            else:
                log.info('Cleaning up temporary directory: %s', self._workdir)
                shutil.rmtree(self._workdir)
451196
def print_seqs(seqs, id2name, name):
    
    orfs = open('%s.orfs.faa' % (name.rsplit('.', 1)[0]), 'w')
    introns = open('%s.introns.fa' % (name.rsplit('.', 1)[0]), 'w')
    insertions = open('%s.insertions.fa' % (name.rsplit('.', 1)[0]), 'w')
    for seq in seqs:
        for i, ins in enumerate(seqs[seq][2], 1):
            model_pos = ins[1]
            if ins[6] is True: # orf(s) in ins[4]
                for orf in ins[4]:
                    orf_info = orf[0].split('>')[1].split()
                    id = orf_info[0].split('_', 1)[0]
                    name = id2name[id]
                    annotation = orf_info[1]
                    strand = orf_info[7]
                    if strand == '1':
                        strand = '+'
                    else:
                        strand = '-'
                    start, stop = sorted([int(orf_info[3]), int(orf_info[5])])
                    header = '>%s insertion::seq=%s type=%s strand=%s gene-pos=%s-%s model-pos=%s'\
                            % (name, i, annotation, strand, start, stop, model_pos)
                    print('\n'.join([header, orf[1]]), file=orfs)
            if ins[7] is True: # intron(s) in ins[5]
                for intron in ins[5]:
                    id, type, strand, pos = intron[0].split('>', 1)[1].split()
                    name = id2name[id.split('_')[0]]
                    header = '>%s insertion::seq=%s type=%s strand=%s gene-pos=%s model-pos=%s'\
                            % (name, i, type, strand, pos, model_pos)
                    print('\n'.join([header, intron[1]]), file=introns)
            insertion = ins[3]
            id, info = insertion[0].split('>')[1].split(' ', 1)
            name = id2name[id.split('_')[0]]
            header = '>%s %s' % (name, info)
            print('\n'.join([header, insertion[1]]), file=insertions)
    insertions.close()        
    orfs.close()
    introns.close()
353074
def parse_args():
    """"""
    parser = argparse.ArgumentParser(
        prog="geneparse-indexer",
        description="Genotype file indexer."
    )

    # IMPUTE2 files
    group = parser.add_argument_group("IMPUTE2 index")
    group.add_argument(
        "--impute2", metavar="IMPUTE2", type=str, nargs="+",
        help="Index an IMPUTE2 genotype file format. The file can be plain "
             "text or bgzipped.",
    )

    # BGEN files
    group = parser.add_argument_group("BGEN index")
    group.add_argument(
        "--bgen", metavar="BGEN", type=str, nargs="+",
        help="Index a BGEN genotype file. This requires 'bgenix' to be in the "
             "PATH.",
    )
    group.add_argument(
        "--legacy", action="store_true",
        help="Index the file using the '-with-rowid' option. This flag "
             "enables compatibility with SQLITE prior to version 3.8.2. See "
             "https://bitbucket.org/gavinband/bgen/wiki/bgenix for more "
             "information.",
    )

    return parser.parse_args()
1047820
def _context_names():
    """"""
    import inspect

    from renku.models import provenance
    from renku.models._jsonld import JSONLDMixin

    for name in dir(provenance):
        cls = getattr(provenance, name)
        if inspect.isclass(cls) and issubclass(cls, JSONLDMixin):
            yield name
89687
def main():
    
    import argparse

    parser = argparse.ArgumentParser(description="Ancestral character reconstruction and visualisation "
                                                 "for rooted phylogenetic trees.", prog='pastml')

    tree_group = parser.add_argument_group('tree-related arguments')
    tree_group.add_argument('-t', '--tree', help="input tree in newick format (must be rooted).",
                            type=str, required=True)

    annotation_group = parser.add_argument_group('annotation-file-related arguments')
    annotation_group.add_argument('-d', '--data', required=True, type=str,
                                  help="annotation file in tab/csv format with the first row "
                                       "containing the column names.")
    annotation_group.add_argument('-s', '--data_sep', required=False, type=str, default='\t',
                                  help="column separator for the annotation table. "
                                       "By default is set to tab, i.e. for a tab-delimited file. "
                                       "Set it to ',' if your file is csv.")
    annotation_group.add_argument('-i', '--id_index', required=False, type=int, default=0,
                                  help="index of the annotation table column containing tree tip names, "
                                       "indices start from zero (by default is set to 0).")

    acr_group = parser.add_argument_group('ancestral-character-reconstruction-related arguments')
    acr_group.add_argument('-c', '--columns', nargs='*',
                           help="names of the annotation table columns that contain characters "
                                "to be analysed. "
                                "If not specified, all columns are considered.",
                           type=str)
    acr_group.add_argument('--prediction_method',
                           choices=[MPPA, MAP, JOINT, DOWNPASS, ACCTRAN, DELTRAN, COPY, ALL, ML, MP],
                           type=str, nargs='*', default=MPPA,
                           help='ancestral character reconstruction (ACR) method, '
                                'can be one of the max likelihood (ML) methods: {ml}, '
                                'one of the max parsimony (MP) methods: {mp}; '
                                'or {copy} to keep the annotated character states as-is without inference. '
                                'One can also specify one of the meta-methods {meta} that would perform ACR '
                                'with multiple methods (all of them for {meta_all}, '
                                'all the ML methods for {meta_ml}, or all the MP methods for {meta_mp}) '
                                'and save/visualise the results as multiple characters '
                                'suffixed with the corresponding method.'
                                'When multiple ancestral characters are specified (see -c, --columns), '
                                'the same method can be used for all of them (if only one method is specified), '
                                'or different methods can be used (specified in the same order as -c, --columns). '
                                'If multiple methods are given, but not for all the characters, '
                                'for the rest of them the default method ({default}) is chosen.'
                           .format(ml=', '.join(ML_METHODS), mp=', '.join(MP_METHODS), copy=COPY, default=MPPA,
                                   meta=', '.join(META_ML_METHODS | {MP}), meta_ml=ML, meta_mp=MP, meta_all=ALL))
    acr_group.add_argument('--forced_joint', action='store_true',
                           help='add {joint} state to the {mppa} state selection '
                                'even if it is not selected by Brier score.'.format(joint=JOINT, mppa=MPPA))
    acr_group.add_argument('-m', '--model', default=F81,
                           choices=[JC, F81, EFT, HKY, JTT],
                           type=str, nargs='*',
                           help='evolutionary model for ML methods (ignored by MP methods). '
                                'When multiple ancestral characters are specified (see -c, --columns), '
                                'the same model can be used for all of them (if only one model is specified), '
                                'or different models can be used (specified in the same order as -c, --columns). '
                                'If multiple models are given, but not for all the characters, '
                                'for the rest of them the default model ({}) is chosen.'.format(F81))
    acr_group.add_argument('--parameters', type=str, nargs='*',
                           help='optional way to fix some of the ML-method parameters '
                                'by specifying files that contain them. '
                                'Should be in the same order '
                                'as the ancestral characters (see -c, --columns) '
                                'for which the reconstruction is to be preformed. '
                                'Could be given only for the first few characters. '
                                'Each file should be tab-delimited, with two columns: '
                                'the first one containing parameter names, '
                                'and the second, named "value", containing parameter values. '
                                'Parameters can include character state frequencies '
                                '(parameter name should be the corresponding state, '
                                'and parameter value - the float frequency value, between 0 and 1),'
                                'and tree branch scaling factor (parameter name {}).'.format(SCALING_FACTOR))

    vis_group = parser.add_argument_group('visualisation-related arguments')
    vis_group.add_argument('-n', '--name_column', type=str, default=None,
                           help="name of the character to be used for node names "
                                "in the compressed map visualisation "
                                "(must be one of those specified via -c, --columns). "
                                "If the annotation table contains only one column it will be used by default.")
    vis_group.add_argument('--date_column', required=False, default=None,
                           help="name of the annotation table column that contains tip dates, "
                                "if specified it is used to add a time slider to the visualisation.",
                           type=str)
    vis_group.add_argument('--tip_size_threshold', type=int, default=REASONABLE_NUMBER_OF_TIPS,
                           help="recursively remove the tips of size less than threshold-th largest tip"
                                "from the compressed map (set to 1e10 to keep all tips). "
                                "The larger it is the less tips will be trimmed.")

    out_group = parser.add_argument_group('output-related arguments')
    out_group.add_argument('-o', '--out_data', required=False, type=str,
                           help="path to the output annotation file with the reconstructed ancestral character states.")
    out_group.add_argument('--work_dir', required=False, default=None, type=str,
                           help="path to the folder where pastml parameter, named tree "
                                "and marginal probability (for marginal ML methods ({}) only) files are to be stored. "
                                "Default is <path_to_input_file>/<input_file_name>_pastml. "
                                "If the folder does not exist, it will be created."
                           .format(', '.join(MARGINAL_ML_METHODS)))
    out_group.add_argument('-p', '--html_compressed', required=False, default=None, type=str,
                           help="path to the output compressed map visualisation file (html).")
    out_group.add_argument('-l', '--html', required=False, default=None, type=str,
                           help="path to the output full tree visualisation file (html).")
    out_group.add_argument('-v', '--verbose', action='store_true',
                           help="print information on the progress of the analysis (to console)")

    parser.add_argument('--version', action='version', version='%(prog)s {version}'.format(version=PASTML_VERSION))

    itol_group = parser.add_argument_group('iTOL-related arguments')
    itol_group.add_argument('--upload_to_itol', action='store_true',
                            help="upload the ACR annotated tree to iTOL (https://itol.embl.de/)")
    itol_group.add_argument('--itol_id', required=False, default=None, type=str,
                            help="iTOL user batch upload ID that enables uploading to your iTOL account "
                                 "(see https://itol.embl.de/help.cgi#batch). "
                                 "If not specified, the tree will not be associated to any account.")
    itol_group.add_argument('--itol_project', required=False, default=None, type=str,
                            help="iTOL project the annotated tree should be associated with "
                                 "(must exist, and --itol_id must be specified). "
                                 "If not specified, the tree will not be associated with any project.")
    itol_group.add_argument('--itol_tree_name', required=False, default=None, type=str,
                            help="name for the tree uploaded to iTOL.")

    params = parser.parse_args()

    pastml_pipeline(**vars(params))