diff --git a/standalone/README.md b/standalone/README.md new file mode 100644 index 0000000..00605c7 --- /dev/null +++ b/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 [--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 +``` diff --git a/standalone/adjust_tree_to_msa.py b/standalone/adjust_tree_to_msa.py new file mode 100644 index 0000000..84c0fdf --- /dev/null +++ b/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) diff --git a/standalone/auxiliaries.py b/standalone/auxiliaries.py new file mode 100644 index 0000000..badf86d --- /dev/null +++ b/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 + + + diff --git a/standalone/bin/MPreconstruct b/standalone/bin/MPreconstruct new file mode 100755 index 0000000..929da3a Binary files /dev/null and b/standalone/bin/MPreconstruct differ diff --git a/standalone/bin/removeTaxa b/standalone/bin/removeTaxa new file mode 100755 index 0000000..078c479 Binary files /dev/null and b/standalone/bin/removeTaxa differ diff --git a/standalone/compute_homoplasy.py b/standalone/compute_homoplasy.py new file mode 100644 index 0000000..7244871 --- /dev/null +++ b/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) + + diff --git a/standalone/compute_sweeps_score.py b/standalone/compute_sweeps_score.py new file mode 100644 index 0000000..83fa8ab --- /dev/null +++ b/standalone/compute_sweeps_score.py @@ -0,0 +1,354 @@ +import logging +import sys + +from auxiliaries import load_header2sequences_dict + +logging.basicConfig(level=logging.INFO) +logger = logging.getLogger('main') + +import re +import os +import matplotlib +matplotlib.use('Agg') + +import matplotlib.pyplot as plt +from sklearn import metrics +from numpy import mean, median, argmax +from collections import Counter +from itertools import combinations + + +def get_sequences_from_fasta(msa_path): + # return a list of the sequences without inner "\n"s + sequences = re.split('@', re.sub('\r?\n', '', re.sub('>.+\r?\n', '@', open(msa_path).read())))[1:] + return sequences + + +def get_pairwise_distance(seq1, seq2): + return float(sum([c1 != c2 for c1, c2 in zip(seq1, seq2)])) + + +def get_average_pairwise_distance_and_pi(header2sequence, msa_length): + number_of_species = len(header2sequence) + print(number_of_species) + num_of_pairs = number_of_species * (number_of_species-1) / 2 # number_of_species choose 2 + total_relative_pairwise_distance = 0 + pi = 0 + + sequence2count = {} + for header in header2sequence: + sequence2count[header] = sequence2count.get(header2sequence[header], 0) + 1 + + sequences2frequency = {sequence: sequence2count[sequence]/number_of_species for sequence in sequence2count} + + sequence_pairs = set(combinations(sequence2count, 2)) + for sequence_pair in sequence_pairs: + apd = get_pairwise_distance(*sequence_pair)/msa_length + total_relative_pairwise_distance += apd + pi += 2 * sequences2frequency[sequence_pair[0]] * sequences2frequency[sequence_pair[1]] * apd + + logger.info(f'pi is {pi}') + logger.info(f'total relative apd is {total_relative_pairwise_distance / num_of_pairs}') + return total_relative_pairwise_distance / num_of_pairs, pi + + +def get_alignment_columns_from_dict_msa(header2sequence, msa_length): + msa_columns = [] + + # iterate over msa columns + for j in range(msa_length): + # extract coulmn j + column = ''.join([header2sequence[header][j] for header in header2sequence]) + msa_columns.append(column) + + return msa_columns, len(column) + + +def load_homoplasy(homoplasy_path): + # return a list of ints indicating whether the corresponding columns in an alignment are homoplasy or not + with open(homoplasy_path) as f: + homoplasy = [int(x) for x in f.read().split()] + + return homoplasy + + +# calculate the Hamming similarity of two given columns +def calculate_hamming_similarity(col1, col2): + + chars1 = sorted(set(col1), key=col1.count, reverse=True) + chars2 = sorted(set(col2), key=col2.count, reverse=True) + + for i in range(len(chars1)): + col1 = col1.replace(chars1[i], str(i)) + + for i in range(len(chars2)): + col2 = col2.replace(chars2[i], str(i)) + + return sum(c1 == c2 for c1, c2 in zip(col1, col2)) / len(col1) + + +# calculate the Eli similarity of two given columns +def calculate_greedy_association(col1, col2): + res = 0 + char_to_num = {'A': 0, 'C': 1, 'G': 2, 'T': 3, '-': 4} + pairs_occurences = [[0] * 5, [0] * 5, [0] * 5, [0] * 5, [0] * 5] + for i in range(len(col1)): + pairs_occurences[char_to_num[col1[i]]][char_to_num[col2[i]]] += 1 + + counts = Counter(col1) + sortedCounts = [x[0] for x in sorted(counts.items(), key=lambda x: x[1], reverse=True)] + for char in sortedCounts: + maxOccurenceIndex = argmax(pairs_occurences[char_to_num[char]]) + res += pairs_occurences[char_to_num[char]][maxOccurenceIndex] + + for i in range(len(char_to_num.keys())): + pairs_occurences[i][maxOccurenceIndex] = -1 + + return res / len(col1) + + +# def calculate_greedy_association2(col1, col2, chars='ACGT-'): +# res = 0 +# char2char2charcount = {char: dict.fromkeys(chars, 0) for char in chars} +# +# print(char2char2charcount) +# for ch1, ch2 in zip(col1, col2): +# char2char2charcount[ch1][ch2] = char2char2charcount[ch1].get(ch2, 0) + 1 +# +# counts = Counter(col1) +# sorted_chars_by_frequency = [char for char in sorted(counts, key=counts.get, reverse=True)] +# for char in sorted_chars_by_frequency: +# max_occurence_mate_pair = max(char2char2charcount[char], key=char2char2charcount[char].get) +# res += char2char2charcount[char][max_occurence_mate_pair] +# +# char2char2charcount[char] = {char: dict.fromkeys(chars, -1) for char in chars} +# +# return res / len(col1) + + +def calculate_symmetric_greedy_association(col1, col2): + # calculate a symmetric greedy association of two given columns + # assert calculate_greedy_association(col1, col2) == calculate_greedy_association2(col1, col2) + # assert calculate_greedy_association(col2, col1) == calculate_greedy_association2(col2, col1) + return (calculate_greedy_association(col1, col2) + calculate_greedy_association(col2, col1)) / 2 + + +def get_ith_window(iterable, i, window_size): + return iterable[i: i + window_size] + + +def calculate_total_edge_contribution(edge_column, similarity_function, + window_columns_without_edge_column, window_homoplasy_without_edge_column): + + # calculate the contribution of an edge column to the total window score + edge_contributions = 0 + + for column, homoplasy in zip(window_columns_without_edge_column, window_homoplasy_without_edge_column): + if homoplasy: + # otherwise, ignore this column + edge_contributions += similarity_function(edge_column, column) + + return edge_contributions + + +def get_first_window_score(window_columns, window_homoplasy, similarity_function): + # calculate a window score. + # inefficient for second window (etc...) since the every window overlaps with the previous one + window_score = 0 + + homoplasious_window_columns = [column for i, column in enumerate(window_columns) if window_homoplasy[i]] + for i in range(len(homoplasious_window_columns) - 1): + for j in range(i + 1, len(homoplasious_window_columns)): + window_score += similarity_function(homoplasious_window_columns[i], homoplasious_window_columns[j]) + + return window_score + + +def get_next_window_score(prev_window_score, similarity_function, + prev_window_columns, prev_window_homoplasy, + window_columns, window_homoplasy): + + leftmost_edge_contribution = rightmost_edge_contribution = 0 + + if prev_window_homoplasy[0]: + # previous leftmost edge has homoplasy; otherwise, no contribution... + edge_column = prev_window_columns[0] + leftmost_edge_contribution = calculate_total_edge_contribution(edge_column, similarity_function, + prev_window_columns[1:], prev_window_homoplasy[1:]) + + if window_homoplasy[-1]: + # current rightmost edge has homoplasy; otherwise, no contribution... + edge_column = window_columns[-1] + rightmost_edge_contribution = calculate_total_edge_contribution(edge_column, similarity_function, + window_columns[:-1], window_homoplasy[:-1]) + + return prev_window_score - leftmost_edge_contribution + rightmost_edge_contribution + + +def get_contig_scores(msa_columns, msa_length, homoplasy, window_size, similarity_function, normalization_factor): + + # setting similarity function + if similarity_function == 'symmetric_greedy': + similarity_function = calculate_symmetric_greedy_association + elif similarity_function == 'greedy': + similarity_function = calculate_greedy_association + elif similarity_function == 'mutual': + similarity_function = metrics.mutual_info_score + else: + similarity_function = calculate_hamming_similarity + + # calculates scores for every window in msa + number_of_windows = msa_length - window_size + 1 + logger.info(f'Total number of windows is {number_of_windows}...') + + # scores vector initialization + window_scores = [0] * (msa_length - window_size + 1) + + # first window initialization: + logger.info(f'Calculating score for window #0 with {similarity_function.__name__} as similarity function') + window_columns = get_ith_window(msa_columns, 0, window_size) + window_homoplasy = get_ith_window(homoplasy, 0, window_size) + window_scores[0] = get_first_window_score(window_columns, window_homoplasy, similarity_function) + + # sliding window over the next windows + for i in range(1, number_of_windows): + if i % 100 == 0: + logger.info(f'Calculating score for window #{i}...') + prev_window_columns = window_columns + prev_window_homoplasy = window_homoplasy + window_columns = get_ith_window(msa_columns, i, window_size) + window_homoplasy = get_ith_window(homoplasy, i, window_size) + + window_scores[i] = get_next_window_score(window_scores[i - 1], similarity_function, + prev_window_columns, prev_window_homoplasy, + window_columns, window_homoplasy) + + # normalizing each score by the answer of possible pairs in a window + # normalization must be done AFTER the all the scores were obtained since in order to compute each score + # efficiently, current (unnormalized) score uses previous(unnormalized) score! + normalized_window_scores = [score * normalization_factor for score in window_scores] + return normalized_window_scores + + +def write_summary(msa_name, scores, window_size, header2sequence, msa_length, number_of_sequences, + meta_output_path, mode='w'): + + max_score = max(scores) + index_of_max = scores.index(max_score) + 1 # start from 1 rather than 0 + mean_score = mean(scores) + median_score = median(scores) + max_mean_division = -1 if mean_score == 0 else max_score / mean_score + max_median_division = -1 if median_score == 0 else max_score / median_score + relative_location_of_peak = index_of_max / msa_length + centrality = min(1 - relative_location_of_peak, relative_location_of_peak) * 2 # scaling from 0 to 1 + apd, pi = get_average_pairwise_distance_and_pi(header2sequence, msa_length) + + above05 = above25 = above50 = above75 = above95 = 0 + if max_score > 0.95: + above05 = above25 = above50 = above75 = above95 = 1 + elif max_score > 0.75: + above05 = above25 = above50 = above75 = 1 + elif max_score > 0.5: + above05 = above25 = above50 = 1 + elif max_score > 0.25: + above05 = above25 = 1 + elif max_score > 0.05: + above05 = 1 + + meta_data = [msa_name,max_score,number_of_sequences,centrality,msa_length,window_size,index_of_max, + mean_score,median_score,max_mean_division,max_median_division,relative_location_of_peak, + apd,pi,above95,above75,above50,above25,above05] + + with open(meta_output_path, mode) as f: + f.write(','.join([meta_data[0]] + [f'{abs(score):.4f}' for score in meta_data[1:]]) + '\n') + + +def plot_scores(scores, msa_length, window_size, output_path, similarity_function, + title='', x_label='Window #', y_label='SINCOPA Score', epsilon=0.05): + + # create plot + plt.plot(range(msa_length - window_size + 1), scores) + + # set axes limits + # plt.xlim(-1, msa_length - window_size + 1) + plt.ylim(0 - epsilon, 1 + epsilon) + + # set labels + plt.title(title) + plt.xlabel(x_label) + plt.ylabel(y_label) + + if similarity_function == 'mutual': + # Semi-log for mutual info. Irrelelvant for symmetric_greedy. + plt.yscale('log') + + # save + plt.savefig(output_path) + plt.close() + + +def compute_sweeps_score(msa_path, homoplasy_path, scores_output_path, meta_output_path, plot_path, window_size, + similarity_function='symmetric_greedy', stats_writing_mode='w'): + + logger.info('Starting to compute_sweeps_score...') + + header2sequence, msa_length = load_header2sequences_dict(msa_path, get_length=True, upper_sequence=True) + + msa_columns, number_of_sequences = get_alignment_columns_from_dict_msa(header2sequence, msa_length) + msa_name = os.path.split(msa_path)[-1] + homoplasy = load_homoplasy(homoplasy_path) + + assert msa_length == len(homoplasy), 'MSA length is inconsistent with homoplasy length ' \ + '(they should be both of the same length).' + + scores = get_contig_scores(msa_columns, msa_length, homoplasy, window_size, similarity_function, + normalization_factor=2/(window_size * (window_size-1))) # 1 / window_size choose 2 + + # write scores to file + with open(scores_output_path, 'w') as f: + f.write('\n'.join([f'{abs(score):.4f}' for score in scores]) + '\n') + + write_summary(msa_name, scores, window_size, header2sequence, msa_length, number_of_sequences, + meta_output_path, stats_writing_mode) + + plot_scores(scores, msa_length, window_size, plot_path, similarity_function, + title=f'{msa_name}\n(across {number_of_sequences} sequences)') + + +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('homoplasy_path', help='A path to the homoplasy computation') + # homoplasy path does not exist in the case where a species tree is provided + parser.add_argument('scores_output_path', + help='A path in which the window scores 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('meta_output_path', + help='A path in which the meta output 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('plot_path', + help='A path in which the scores disdribution will be plotted 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('window_size', type=int, + help='The size of a window to which a score will be computed') + parser.add_argument('--similarity_function', default='symmetric_greedy', + help='The type of similarity to be computed between msa columns', + choices=['symmetric_greedy', 'greedy', 'mutual', 'hamming']) + + parser.add_argument('-v', '--verbose', help='Increase output verbosity', action='store_true') + + args = parser.parse_args() + + compute_sweeps_score(args.msa_path, args.homoplasy_path, args.scores_output_path, args.meta_output_path, + args.plot_path, args.window_size, args.similarity_function) + + diff --git a/standalone/environment.yml b/standalone/environment.yml new file mode 100644 index 0000000..84e2af0 --- /dev/null +++ b/standalone/environment.yml @@ -0,0 +1,10 @@ +name: sincopa +channels: + - conda-forge + - defaults +dependencies: + - python>=3.6 + - biopython + - numpy + - matplotlib + - scikit-learn diff --git a/standalone/fix_msa.py b/standalone/fix_msa.py new file mode 100644 index 0000000..6929ad3 --- /dev/null +++ b/standalone/fix_msa.py @@ -0,0 +1,92 @@ +import os +import logging + +from auxiliaries import load_header2sequences_dict + +logging.basicConfig(level=logging.INFO) +logger = logging.getLogger('main') + + +def fix_ambiguous_columns(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'Ambiguous character {strain2sequence[strain][column]} was detected. Replacing column #{column} of {strain} from {strain2sequence[strain][column]} to {major_allele}') + strain2sequence[strain] = strain2sequence[strain][:column] + major_allele + strain2sequence[strain][column+1:] + + +def get_first_column_without_gap(strain2sequence, indexes): + for col_index in indexes: + for strain in strain2sequence: + if strain2sequence[strain][col_index] == '-': + break + else: + # no gap was found! + return col_index + + +def trim_msa(strain2sequence, msa_length): + first_left_col_without_gaps = get_first_column_without_gap(strain2sequence, range(msa_length)) # left side + first_right_col_without_gaps = get_first_column_without_gap(strain2sequence, range(msa_length-1, -1, -1)) # right side + + for strain in strain2sequence: + strain2sequence[strain] = strain2sequence[strain][first_left_col_without_gaps: 1+first_right_col_without_gaps] + + return 1+first_right_col_without_gaps - first_left_col_without_gaps + + +def fix_msa(msa_path, output_path, minimal_length=300): + strain2sequence, msa_length = load_header2sequences_dict(msa_path, get_length=True, upper_sequence=True) + + for strain in strain2sequence: + if len(strain2sequence[strain]) != msa_length: + raise ValueError(f'Illegal MSA. Not all sequences are of the same length. E.g., {strain} sequence length is {len(strain2sequence[strain])} where others are of length {msa_length}.') + + logger.info(f'Trimming {os.path.split(msa_path)[-1]}') + logger.info(f'MSA length before trimming is {msa_length}bps') + trimmed_msa_length = trim_msa(strain2sequence, msa_length) + logger.info(f'MSA length after trimming is {trimmed_msa_length}bps (a total of {msa_length - trimmed_msa_length}bps were trimmed)') + + if trimmed_msa_length < minimal_length: + logger.error(f'Discarding too short MSA (needs to be at least {minimal_length}bps wide).') + return + + for col in range(trimmed_msa_length): + fix_ambiguous_columns(col, strain2sequence) + + fixed_alignment = '' + for strain in strain2sequence: + fixed_alignment += f'>{strain}\n{strain2sequence[strain]}\n' + + with open(output_path, 'w') as f: + f.write(fixed_alignment) + + logger.info(f'MSA was fixed and stored successfully at {output_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('input_msa_path', + help='A path to an MSA file that should be fixed. Fixation includes: ' + '(1) ambiguous chars removal and (2) msa edges trimming.', + type=lambda path: path if os.path.exists(path) else parser.error(f'{path} does not exist!')) + parser.add_argument('output_msa_path', + help='A path to the fixed MSA', + 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('--drop-shorter', default=300, help='A shorter MSA will be discarded.', + type=lambda x: int(x) if int(x) > 0 else parser.error(f'Minimal number of upstream basepairs should be positive!')) + parser.add_argument('-v', '--verbose', help='Increase output verbosity', action='store_true') + + args = parser.parse_args() + + fix_msa(args.input_msa_path, args.output_msa_path, args.drop_shorter) + + diff --git a/standalone/sincopa.py b/standalone/sincopa.py new file mode 100644 index 0000000..78ee77f --- /dev/null +++ b/standalone/sincopa.py @@ -0,0 +1,172 @@ +import os +import logging +import Bio.SeqUtils +from time import sleep +from compute_homoplasy import compute_homoplasy +from compute_sweeps_score import compute_sweeps_score +from fix_msa import fix_msa +from adjust_tree_to_msa import fix_tree +from auxiliaries import * + +logging.basicConfig(level=logging.INFO) +logger = logging.getLogger('main') + + +def verify_fasta_format(fasta_path): + logger.info('Validating FASTA format') + Bio.SeqUtils.IUPACData.ambiguous_dna_letters += 'U-' + legal_chars = set(Bio.SeqUtils.IUPACData.ambiguous_dna_letters.lower() + Bio.SeqUtils.IUPACData.ambiguous_dna_letters) + with open(fasta_path) as f: + line_number = 0 + try: + line = f.readline() + line_number += 1 + if not line.startswith('>'): + return f'Illegal FASTA format. First line starts with "{line[0]}" instead of ">".' + previous_line_was_header = True + putative_end_of_file = False + curated_content = f'>{line[1:]}'.replace("|", "_") + for line in f: + line_number += 1 + line = line.strip() + if not line: + if not putative_end_of_file: + putative_end_of_file = line_number + continue + if putative_end_of_file: + return f'Illegal FASTA format. Line {putative_end_of_file} in MSA is empty.' + if line.startswith('>'): + if previous_line_was_header: + return f'Illegal FASTA format. MSA contains an empty record. Both lines {line_number-1} and {line_number} start with ">".' + else: + previous_line_was_header = True + curated_content += f'>{line[1:]}\n'.replace("|", "_") + continue + else: + previous_line_was_header = False + for c in line: + if c not in legal_chars: + return f'Illegal FASTA format. Line {line_number} contains illegal DNA character "{c}".' + curated_content += f'{line}\n' + except UnicodeDecodeError as e: + logger.info(e.args) + line_number += 1 + return f'Illegal FASTA format. Line {line_number} contains non-ASCII character(s).' + with open(fasta_path, 'w') as f: + f.write(curated_content) + + +def verify_newick_format(tree_path): + logger.info('Validating NEWICK format') + pass + + +def verify_msa_is_consistent_with_tree(msa_path, tree_path): + logger.info('Validating MSA is consistent with the species tree') + tree_strains = get_tree_labels(tree_path) + logger.info(f'Phylogenetic tree contains the following strains:\n{tree_strains}') + with open(msa_path) as f: + logger.info(f'Checking MSA...') + for line in f: + if line.startswith('>'): + strain = line.lstrip('>').rstrip('\n') + if strain not in tree_strains: + msg = f'{strain} appears in the MSA but not in the tree. ' \ + f'Please make sure the tree contains all species in the MSA.' + logger.error(msg) + return msg + else: + logger.info(f'{strain} appears in tree!') + + +def validate_input(msa_path, tree_path, error_path): + logger.info('Validating input...') + error_msg = verify_fasta_format(msa_path) + if error_msg: + fail(error_msg, error_path) + error_msg = verify_newick_format(tree_path) + if error_msg: + fail(error_msg, error_path) + error_msg = verify_msa_is_consistent_with_tree(msa_path, tree_path) + if error_msg: + fail(error_msg, error_path) + + +def fix_input(msa_path, tree_path, output_dir, tmp_dir): + tree_name = os.path.split(tree_path)[-1] + adjusted_tree = f'{output_dir}/{os.path.splitext(tree_name)[0]}_fixed{os.path.splitext(tree_name)[-1]}' + fix_tree(msa_path, tree_path, tmp_dir, adjusted_tree) + msa_name = os.path.split(msa_path)[-1] + fixed_msa_path = f'{output_dir}/{os.path.splitext(msa_name)[0]}_fixed{os.path.splitext(msa_name)[-1]}' + fix_msa(msa_path, fixed_msa_path) + return fixed_msa_path, adjusted_tree + + +def sincopa(msa_path, tree_path, window_size, output_dir, tmp_dir): + homplasy_path = f'{output_dir}/homoplasy.txt' + control_file_path = os.path.join(tmp_dir, 'control.txt') + sweeps_scores_path = f'{output_dir}/sweeps_scores.txt' + sweeps_plot_path = f'{output_dir}/sweeps_scores.png' + sweeps_summary_path = f'{output_dir}/sweeps_summary.txt' + done_path = f'{output_dir}/done.txt' + + header = 'msa_name,max_score,number_of_sequences,centrality,msa_length,window_size,index_of_max,' \ + 'mean_score,median_score,max_mean_division,max_median_division,relative_location_of_peak,' \ + 'apd,pi,above95,above75,above50,above25,above05' + with open(sweeps_summary_path, 'w') as f: + f.write(f'{header}\n') + + compute_homoplasy(msa_path, tree_path, control_file_path, homplasy_path) + + compute_sweeps_score(msa_path, homplasy_path, sweeps_scores_path, sweeps_summary_path, + sweeps_plot_path, window_size, stats_writing_mode='a') + + with open(done_path, 'w'): + pass + + +def main(msa_path, tree_path, window_size, output_dir_path, html_path=None): + error_path = f'{output_dir_path}/error.txt' + try: + os.makedirs(output_dir_path, exist_ok=True) + tmp_dir = f'{os.path.split(msa_path)[0]}/tmp' + os.makedirs(tmp_dir, exist_ok=True) + validate_input(msa_path, tree_path, error_path) + msa_path, tree_path = fix_input(msa_path, tree_path, output_dir_path, tmp_dir) + sincopa(msa_path, tree_path, window_size, output_dir_path, tmp_dir) + logger.info('SUCCEEDED = True') + except Exception as e: + logger.info(f'SUCCEEDED = False') + logger.error(str(e)) + import traceback + traceback.print_exc() + + +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('input_msa_path', + help='A path to a DNA MSA file to look for sweeps.', + type=lambda path: path if os.path.exists(path) else parser.error(f'{path} does not exist!')) + parser.add_argument('input_tree_path', + help='A path to a background species tree that contains (at least) all the species in the ' + 'input MSA. The tree should be reconstructed by external data and not by the MSA provided.', + type=lambda path: path if os.path.exists(path) else parser.error(f'{path} does not exist!')) + parser.add_argument('output_dir_path', + help='A path to a folder in which the sweeps analysis will be written.', + type=lambda path: path.rstrip('/')) + parser.add_argument('--window_size', type=int, default=50, + help='The size of a window to which a score will be computed') + parser.add_argument('-v', '--verbose', help='Increase output verbosity', action='store_true') + + args = parser.parse_args() + + if args.verbose: + logging.basicConfig(level=logging.DEBUG) + else: + logging.basicConfig(level=logging.INFO) + + main(args.input_msa_path, args.input_tree_path, args.window_size, args.output_dir_path) diff --git a/standalone/test/example.fasta b/standalone/test/example.fasta new file mode 100644 index 0000000..ec5740d --- /dev/null +++ b/standalone/test/example.fasta @@ -0,0 +1,8 @@ +>oren +AGGCGTCCCAGCTTGTTGCGGTGAGAGCTTTTCACTTAGACGCGAATTACCACATGTAATACGAACACCTCAGCGAGGTCTACTTGCCTCGGGCCCCTGGTTTTGGAGGCCGTTATCTTACCCGCGAGGCGAGCATATCACATTACCATCGAAACCGTTGTCCGGGATTATCTACCCATGTCGGCAGTATTAAGCCATTAACTGGTGTAAATCTTCATAGTTCGGGTAGTGTTGGCTCAGGACCAACCTTTAGGCCAATACCGCCCACGGAGGAGTTTCACGCCGCGGCGGGCATATCGA +>tal +TGATTATCATTTAACTGTGAGGCTTTAGCCAAAGCATACCAGCCTGACGTCGTACAAGTGCCACGCTGTTTAACGAGGTCTACTTGCCTCGGGCCCCTGGTTTTGGAGGCCGTTATCTTACCCGCGAGGCGAGCATATCACATTACGTAAATGCTAGCCGCGACCCCAGGAAATTAGAAATGGTGTCTCTTGGAGACGGTGCATACGTTCACACCGCCGCCAAAGAGCTTTATGCCTTGTGTATAGCTGAATCCACGGTTTATAAGCGCAGGAAGCAACGAGATCCGCTTTTAGACTACC +>eli +ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGAGGTCTACTTGCCTCGGGCCCCTGGTTTTGGAGGCCGTTATCTTACCCGCGAGGCGAGCATATCACATTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGT +>michal +AGCCTTCGATACTAATTGAGCATCTTGAATGCTAAGTCAAACTAGTCCGATGAGTGCCACTTACCCGTACAGCCGAGGTCTACTTGCCTCGGGCCCCTGGTTTTGGAGGCCGTTATCTTACCCGCGAGGCGAGCATATCACATTACCATAACCGCGGGCAAGACGCAGCCGAACCTGAATCGCCCAAAGCCATGTCGAGTGTGGATTTCATACCGTCGGGACGTTCTTAGAATTGTGCAATTAGTATATCCTATCCGGTTCCATCCGAAATTGGAGTGAAGGGTTAAGGATATATGAGCG diff --git a/standalone/test/example.newick b/standalone/test/example.newick new file mode 100755 index 0000000..5aef6e3 --- /dev/null +++ b/standalone/test/example.newick @@ -0,0 +1 @@ +(tal:6.0,eli:5.0,(oren:3.0,michal:1.0));