Add standalone command-line version

standalone/README.md

@@ -0,0 +1,56 @@
+# SINCOPA — Standalone Command-Line Version
+
+SINCOPA detects selective sweeps in bacteria based on **S**imilarity of **IN**con**CO**ngruent **PA**tterns. It identifies regions of a DNA multiple sequence alignment where allelic patterns are highly incongruent with the species phylogeny, consistent with a transfer-mediated selective sweep.
+
+SINCOPA is also available as a [web server](https://sincopa.tau.ac.il/).
+
+## Requirements
+
+- Python ≥ 3.6
+- conda
+- git ≥ 2.25
+
+## Installation
+
+```bash
+git clone --filter=blob:none --sparse https://github.com/orenavram/SINCOPA.git
+cd SINCOPA
+git sparse-checkout set standalone
+cd standalone
+conda env create -f environment.yml
+conda activate sincopa
+```
+
+## Usage
+
+```bash
+python sincopa.py <alignment.fasta> <tree.newick> <output_dir> [--window_size 50]
+```
+
+| Argument | Description |
+|---|---|
+| `alignment.fasta` | DNA multiple sequence alignment in FASTA format |
+| `tree.newick` | Species phylogeny in Newick format |
+| `output_dir` | Directory where results will be written |
+| `--window_size` | Sliding window size (default: 50) |
+
+**Input requirements:**
+- Alignment must be ≥ 300 bp after trimming
+- The species tree must contain all taxa present in the alignment
+- The tree should be reconstructed from independent data, not from the input alignment
+
+## Output
+
+| File | Description |
+|---|---|
+| `homoplasy.txt` | Parsimony-based homoplasy score per alignment column |
+| `sweeps_scores.txt` | S* score per sliding window |
+| `sweeps_summary.txt` | Summary statistics including peak S* score |
+| `sweeps_scores.png` | S* score plot across the alignment |
+| `done.txt` | Created upon successful completion |
+
+## Quick test
+
+```bash
+python sincopa.py test/example.fasta test/example.newick test/out --window_size 50
+```

standalone/adjust_tree_to_msa.py

@@ -0,0 +1,86 @@
+import os
+import subprocess
+import re
+import logging
+
+from auxiliaries import get_tree_labels
+
+logging.basicConfig(level=logging.INFO)
+logger = logging.getLogger('main')
+
+
+def remove_bootstrap_values(in_tree_path, out_tree_path):
+    with open(in_tree_path) as f:
+        tree_as_str = f.read()
+    tree_as_str = re.sub(r'\)\d+:', '):', tree_as_str)
+    with open(out_tree_path, 'w') as f:
+        f.write(tree_as_str)
+
+
+def prune_tree(msa_path, tree_to_prune_path, tmp_dir, output_tree_path,
+               remove_taxa_script_path=None):
+    # Use local path if not specified
+    if remove_taxa_script_path is None:
+        remove_taxa_script_path = os.path.join(os.path.dirname(__file__), "bin/removeTaxa")
+
+    logger.debug(f'Pruning tree for {msa_path}')
+    msa_name = os.path.split(msa_path)[1]
+
+    # get list of taxa in the full tree
+    tree_taxa = get_tree_labels(tree_to_prune_path)
+
+    # get list of taxa in msa
+    with open(msa_path) as f:
+        msa = f.read()
+
+    msa_taxa = re.findall(r'>(\S+)\r?\n', msa)
+    list_of_taxa_to_remove = [taxon for taxon in tree_taxa if taxon not in msa_taxa]
+
+    # list of names to prune
+    list_of_taxa_to_remove_path = f'{tmp_dir}/{msa_name}.txt'
+    with open(list_of_taxa_to_remove_path, 'w') as f:
+        f.write('\n'.join(list_of_taxa_to_remove))
+
+    cmd_for_pruning = f'{remove_taxa_script_path} {tree_to_prune_path} {list_of_taxa_to_remove_path} {output_tree_path}'
+    logger.info(f'Fetching pruning command:\n{cmd_for_pruning}')
+    subprocess.run(cmd_for_pruning, shell=True)
+
+
+def fix_tree(msa_path, tree_to_adjust, tmp_dir, output_tree_path):
+
+    tree_without_bootstrap_path = tree_to_adjust+'.no_bootstrap'
+
+    remove_bootstrap_values(tree_to_adjust, tree_without_bootstrap_path)
+
+    prune_tree(msa_path, tree_without_bootstrap_path, tmp_dir, output_tree_path)
+
+    logger.info(f'Tree was adjusted successfully and was saved to {output_tree_path}')
+
+    return output_tree_path
+
+
+if __name__ == '__main__':
+        from sys import argv
+        print(f'Starting {argv[0]}. Executed command is:\n{" ".join(argv)}')
+
+        import argparse
+        parser = argparse.ArgumentParser()
+        parser.add_argument('msa_path',
+                            help='A path to an MSA file to which the tree should be adjusted',
+                            type=lambda path: path if os.path.exists(path) else parser.error(f'{path} does not exist!'))
+        parser.add_argument('tree_to_adjust',
+                            help='A path to a species tree that contains (at least) all the species in the input MSA',
+                            type=lambda path: path if os.path.exists(path) else parser.error(f'{path} does not exist!'))
+        parser.add_argument('tmp_dir',
+                            help='A path to a folder in which a txt file (with the same name as the msa_file) will be'
+                                 'created. All the species that do not appear in the msa (and thus will be removed) '
+                                 'will be written to the file that was created',
+                            type=lambda path: path if os.path.exists(path) else parser.error(f'{path} does not exist!'))
+        parser.add_argument('output_tree_path',
+                            help='A path to a file in which the pruned tree will be written',
+                            type=lambda path: path if os.path.exists(os.path.split(path)[0]) else parser.error(
+                                f'output folder {os.path.split(path)[0]} does not exist!'))
+        args = parser.parse_args()
+
+        fix_tree(args.msa_path, args.tree_to_adjust,
+                 args.tmp_dir, args.output_tree_path)

standalone/auxiliaries.py

@@ -0,0 +1,86 @@
+from Bio import Phylo
+import logging
+
+logging.basicConfig(level=logging.INFO)
+logger = logging.getLogger('main')
+
+
+def get_tree_labels(tree_to_prune_path):
+    i = 0
+    logger.info(f'Extracting leaf labels from {tree_to_prune_path}')
+    while True:
+        i += 1
+        logger.info(f'Iteration #{i}')
+        if i > 100:
+            logger.info(f'Failed 100 times to extract leaf labels from {tree_to_prune_path}!')
+            raise AssertionError(f'Failed 100 times to extract leaf labels from {tree_to_prune_path}!')
+        try:
+            tree = Phylo.read(tree_to_prune_path, 'newick')
+            result = [leaf.name for leaf in tree.get_terminals()]
+            break
+        except:
+            pass
+
+    return result
+
+
+def fail(error_msg, error_file_path):
+    with open(error_file_path, 'w') as error_f:
+        error_f.write(error_msg + '\n')
+    raise Exception(error_msg)
+
+
+def update_html(html_path, src, dst):
+    # The initial file exists (generate by the cgi) so we can read and parse it.
+    with open(html_path) as f:
+        html_content = f.read()
+    html_content = html_content.replace(src, dst)
+    with open(html_path, 'w') as f:
+        f.write(html_content)
+
+
+def append_to_html(html_path, new_content):
+    with open(html_path) as f:
+        html_content = f.read()
+    html_content += new_content
+    with open(html_path, 'w') as f:
+        f.write(html_content)
+
+
+def load_header2sequences_dict(fasta_path, get_length=False, upper_sequence=False):
+    header_to_sequence_dict = {}
+    seq_length = 0
+
+    with open(fasta_path) as f:
+        header = f.readline().lstrip('>').rstrip()
+        sequence = ''
+        for line in f:
+            line = line.rstrip()
+            if line.startswith('>'):
+                seq_length = len(sequence)
+                if upper_sequence:
+                    header_to_sequence_dict[header] = sequence.upper()
+                else:
+                    # leave untouched
+                    header_to_sequence_dict[header] = sequence
+                header = line.lstrip('>')
+                sequence = ''
+            else:
+                sequence += line
+
+        # don't forget last record!!
+        if sequence != '':
+
+            if upper_sequence:
+                header_to_sequence_dict[header] = sequence.upper()
+            else:
+                # leave untouched
+                header_to_sequence_dict[header] = sequence
+
+    if get_length:
+        return header_to_sequence_dict, seq_length
+    else:
+        return header_to_sequence_dict
+
+
+

standalone/compute_homoplasy.py

@@ -0,0 +1,80 @@
+import os
+import logging
+
+logging.basicConfig(level=logging.INFO)
+logger = logging.getLogger('main')
+
+
+def fix_column(column, strain2sequence, legal_chars='ACGT-'):
+    col = [strain2sequence[strain][column] for strain in strain2sequence
+           if strain2sequence[strain][column] in legal_chars]
+
+    major_allele = max(set(col), key=col.count)
+    for strain in strain2sequence:
+        if strain2sequence[strain][column] not in legal_chars:
+            logger.info(f'Replacing column #{column} of {strain} from {strain2sequence[strain][column]} to {major_allele}')
+            strain2sequence[strain] = strain2sequence[strain][:column] + major_allele + strain2sequence[strain][column+1:]
+
+
+def compute_homoplasy(msa_path, tree_path, control_path, output_path,
+                      MPReconstruct_script_path=None):
+    # Use local path if not specified
+    if MPReconstruct_script_path is None:
+        MPReconstruct_script_path = os.path.join(os.path.dirname(__file__), "bin/MPreconstruct")
+
+    # make sure that an appropriate gcc module is loaded. E.g., gcc/gcc-6.2.0
+    import subprocess
+
+    # prepare a control file for the c++ code:
+    with open(control_path, 'w') as f:
+        f.write(f'''### parameters for MPreconstruct
+
+## in and out files
+_treefile {tree_path}
+_seqfile {msa_path}
+_logfile {os.path.splitext(control_path)[0]}.log
+_outfile {os.path.splitext(control_path)[0]}.out  
+
+## types are: nuc, amino, threeState, integer
+_alphabetType nuc
+
+## alphabet: a, c, g, t, and -
+_alphabetSize 5
+
+## types are: file,fitch,diff,diffSquare
+_costMatrixType fitch''')
+
+    # what to execute in the shell:
+    cmd = f'{MPReconstruct_script_path} {control_path} {output_path}'
+
+    logger.info(f'Computing homoplasy. Executed command is:\n{cmd}')
+    subprocess.run(cmd, shell=True)
+
+
+if __name__ == '__main__':
+        from sys import argv
+        print(f'Starting {argv[0]}. Executed command is:\n{" ".join(argv)}')
+
+        import argparse
+        parser = argparse.ArgumentParser()
+        parser.add_argument('msa_path',
+                            help='A path to an MSA file to compute homoplasy',
+                            type=lambda path: path if os.path.exists(path) else parser.error(f'{path} does not exist!'))
+        parser.add_argument('tree_path',
+                            help='A path to a(n adjusted) phylogenetic tree',
+                            type=lambda path: path if os.path.exists(path) else parser.error(f'{path} does not exist!'))
+        parser.add_argument('control_path',
+                            help='A path to control file in which the MPReconstruct parameters will be written to',
+                            type=lambda path: path if os.path.exists(os.path.split(path)[0]) else parser.error(
+                                f'output folder {os.path.split(path)[0]} does not exist!'))
+        parser.add_argument('output_path',
+                            help='A path in which the computed homoplasy will be written to',
+                            type=lambda path: path if os.path.exists(os.path.split(path)[0]) else parser.error(
+                                f'output folder {os.path.split(path)[0]} does not exist!'))
+        parser.add_argument('-v', '--verbose', help='Increase output verbosity', action='store_true')
+
+        args = parser.parse_args()
+
+        compute_homoplasy(args.msa_path, args.tree_path, args.control_path, args.output_path)
+
+