diff --git a/.github/workflows/cibuildwheel.yml b/.github/workflows/cibuildwheel.yml index 9857c11..83ace72 100644 --- a/.github/workflows/cibuildwheel.yml +++ b/.github/workflows/cibuildwheel.yml @@ -26,7 +26,7 @@ jobs: env: CIBW_ENVIRONMENT: "BRANCH_NAME=$(cat branch_name.tmp)" CIBW_BUILD_VERBOSITY: 3 - CIBW_BUILD: cp36-manylinux_x86_64 cp37-manylinux_x86_64 cp38-manylinux_x86_64 cp39-manylinux_x86_64 cp310-manylinux_x86_64 + CIBW_BUILD: cp36-manylinux_x86_64 cp37-manylinux_x86_64 cp38-manylinux_x86_64 cp39-manylinux_x86_64 cp310-manylinux_x86_64 cp311-manylinux_x86_64 cp36-macosx_x86_64 cp37-macosx_x86_64 cp38-macosx_x86_64 cp38-macosx_universal2 cp38-macosx_arm64 cp39-macosx_x86_64 cp39-macosx_universal2 cp39-macosx_arm64 CIBW_SKIP: pp* *-manylinux_{aarch64,ppc64le,s390x} CIBW_ARCHS_LINUX: x86_64 CIBW_MANYLINUX_X86_64_IMAGE: quay.io/pypa/manylinux2014_x86_64:latest diff --git a/kSpider_version.py b/kSpider_version.py index 3b27a79..3f6e494 100644 --- a/kSpider_version.py +++ b/kSpider_version.py @@ -52,5 +52,5 @@ def get_version(): else: version_tag = dev_version - # print(f"[DEBUG VERSION] = {version_tag}") + print(f"[DEBUG VERSION] = {version_tag}") return version_tag diff --git a/pykSpider/kSpider2/ks_clustering.py b/pykSpider/kSpider2/ks_clustering.py index 6a8f11b..e888238 100644 --- a/pykSpider/kSpider2/ks_clustering.py +++ b/pykSpider/kSpider2/ks_clustering.py @@ -1,37 +1,50 @@ from __future__ import division from collections import defaultdict -import itertools -import sys import os -import sqlite3 import click from kSpider2.click_context import cli -import glob +import rustworkx as rx +from tqdm import tqdm -class kClusters: +class Clusters: + + distance_to_col = { + "min_cont": 3, + "avg_cont": 4, + "max_cont": 5, + "ani": 6, + } - source = [] - target = [] - source2 = [] - target2 = [] seq_to_kmers = dict() names_map = dict() - components = defaultdict(set) - def __init__(self, logger_obj, index_prefix, cut_off_threshold): + def __init__(self, logger_obj, index_prefix, cut_off_threshold, dist_type): + self.index_prefix = index_prefix self.Logger = logger_obj + self.dist_type = dist_type + self.edges_batch_number = 10_000_000 self.names_file = index_prefix + ".namesMap" self.cut_off_threshold = cut_off_threshold self.seqToKmers_file = index_prefix + "_kSpider_seqToKmersNo.tsv" self.pairwise_file = index_prefix + "_kSpider_pairwise.tsv" - self.uncovered_seqs = set() + self.output = index_prefix + \ + f"_kSpider_clusters_{cut_off_threshold}%.tsv" self.shared_kmers_threshold = 200 - self.seq_to_clusterid = dict() - self.max_cluster_id = 0 self.Logger.INFO("Loading TSV pairwise file") self.load_seq_to_kmers(self.seqToKmers_file) self.tsv_get_namesmap() + if dist_type not in self.distance_to_col: + logger_obj.ERROR("unknown distance!") + self.dist_col = self.distance_to_col[dist_type] + + if dist_type == "ani": + if not os.path.exists(self.index_prefix + "_kSpider_pairwise.ani_col.tsv"): + logger_obj.ERROR(f"ANI was selected, but the ani file {self.index_prefix}_kSpider_pairwise.ani_col.tsv was not found!") + + self.graph = rx.PyGraph() + self.nodes_indeces = self.graph.add_nodes_from( + list(self.names_map.keys())) def load_seq_to_kmers(self, tsv): with open(tsv) as KMER_COUNT: @@ -40,13 +53,6 @@ def load_seq_to_kmers(self, tsv): seq_ID, no_of_kmers = tuple(line.strip().split('\t')[1:]) self.seq_to_kmers[int(seq_ID)] = int(no_of_kmers) - def ids_to_names(self, cluster): - new_cluster = [] - for _id in cluster: - new_cluster.append(self.names_map[int(_id)]) - - return new_cluster - def tsv_get_namesmap(self): with open(self.names_file, 'r') as namesMap: next(namesMap) # skip the header @@ -54,148 +60,78 @@ def tsv_get_namesmap(self): row = row.strip().split() self.names_map[int(row[0])] = row[1] - def tsv_build_graph(self): + def construct_graph(self): + batch_counter = 0 + edges_tuples = [] + print("[i] constructing graph") with open(self.pairwise_file, 'r') as pairwise_tsv: next(pairwise_tsv) # skip header - - for row in pairwise_tsv: - row = row.strip().split() - seq1 = int(row[1]) - seq2 = int(row[2]) - shared_kmers = int(row[3]) - containment = 0.0 - - min_seq = float( - min(self.seq_to_kmers[seq1], self.seq_to_kmers[seq2])) - containment = shared_kmers / min_seq - - if containment < self.cut_off_threshold: - continue - - if shared_kmers < self.shared_kmers_threshold: - self.source2.append(seq1) - self.target2.append(seq2) - - elif shared_kmers >= self.shared_kmers_threshold: - self.source.append(seq1) - self.target.append(seq2) - - # # For covering clusters with single sequence - uncovered_seqs_1 = set(self.names_map.keys()) - \ - set(self.source).union(set(self.target)) - for seq in uncovered_seqs_1: - self.uncovered_seqs.add(seq) - - # OR: - # for i in range(1, len(self.names_map) + 1, 1): - # self.source.append(i) - # self.target.append(i) - - def clustering(self): - registers = defaultdict(lambda: None) - - def find(x): - l = registers[x] - if l is not None: - l = find(l) - registers[x] = l - return l - return x - - def union(x, y): - lx, ly = find(x), find(y) - if lx != ly: - registers[lx] = ly - - for i in range(len(self.source)): - union(self.source.pop(), self.target.pop()) - - for x in registers: - self.components[find(x)].add(x) - - temp_components = self.components.copy() - self.components.clear() - - for cluster_id, (k, v) in enumerate(temp_components.items(), 1): - self.components[cluster_id] = set(v) - for seq in v: - self.seq_to_clusterid[seq] = cluster_id - - temp_components.clear() - self.post_clustering() - - def post_clustering(self): - registers2 = defaultdict(lambda: None) - local_components = defaultdict(set) - covered_seqs = set() - - def find(x): - l = registers2[x] - if l is not None: - l = find(l) - registers2[x] = l - return l - return x - - def union(x, y): - lx, ly = find(x), find(y) - if lx != ly: - registers2[lx] = ly - - for i in range(len(self.source2)): - union(self.source2.pop(), self.target2.pop()) - - for x in registers2: - local_components[find(x)].add(x) - - self.components = dict(self.components) - - covered_clusters = set() - - for cluster2_id, (k, v) in enumerate(local_components.items(), 1): - - for seq in v: - covered_seqs.add(seq) - - for seq in v: - if seq in self.seq_to_clusterid: - cluster_id = self.seq_to_clusterid[seq] - to_be_added = set() - - for i in v: - if i not in self.seq_to_clusterid: - to_be_added.add(i) - - self.components[cluster_id] = self.components[cluster_id].union( - to_be_added) - covered_clusters.add(k) - continue - - self.uncovered_seqs = self.uncovered_seqs - covered_seqs - uncovered_clusters = set(local_components.keys()) - covered_clusters - max_id = len(self.components) - for i, unc in enumerate(uncovered_clusters, 1): - max_id += 1 - self.components[max_id] = local_components[unc] - - for seq in self.uncovered_seqs: - max_id += 1 - self.components[max_id] = {seq} - - def export_kCluster(self): - kCluster_file_name = f"kSpider_{self.cut_off_threshold:.2f}%_" - kCluster_file_name += os.path.basename( - self.pairwise_file).split(".")[0] - kCluster_file_name += ".clusters.tsv" - - with open(kCluster_file_name, 'w') as kClusters: - kClusters.write("kClust_id\tseqs_ids\n") - for cluster_id, (k, v) in enumerate(self.components.items(), 1): - kClusters.write( - f"{cluster_id}\t{'|'.join(self.ids_to_names(v))}\n") - - self.Logger.INFO(f"Total Number Of Clusters: {cluster_id}") + if self.dist_type == "ani": + with open(self.index_prefix + "_kSpider_pairwise.ani_col.tsv") as ani_col_file: + for row in pairwise_tsv: + row = row.strip().split('\t') + seq1 = int(row[0]) - 1 + seq2 = int(row[1]) - 1 + distance = float(next(ani_col_file).strip()) * 100 + + # don't make graph edge + if distance < self.cut_off_threshold: + continue + + if batch_counter < self.edges_batch_number: + batch_counter += 1 + edges_tuples.append((seq1, seq2, distance)) + else: + self.graph.add_edges_from(edges_tuples) + batch_counter = 0 + edges_tuples.clear() + + else: + if len(edges_tuples): + self.graph.add_edges_from(edges_tuples) + else: + for row in pairwise_tsv: + row = row.strip().split('\t') + seq1 = int(row[0]) - 1 + seq2 = int(row[1]) - 1 + distance = float(row[self.dist_col]) * 100 + + # don't make graph edge + if distance < self.cut_off_threshold: + continue + + if batch_counter < self.edges_batch_number: + batch_counter += 1 + edges_tuples.append((seq1, seq2, distance)) + else: + self.graph.add_edges_from(edges_tuples) + batch_counter = 0 + edges_tuples.clear() + + else: + if len(edges_tuples): + self.graph.add_edges_from(edges_tuples) + + def cluster_graph(self): + self.connected_components = rx.connected_components(self.graph) + self.Logger.INFO( + f"number of clusters: {len(self.connected_components)}") + single_components = 0 + retworkx_export = self.index_prefix + \ + f"_kSpider_graph_{self.cut_off_threshold}%.json" + # and {self.output} ...") + self.Logger.INFO(f"writing {retworkx_export}") + # rx.node_link_json(self.graph, path = retworkx_export) + with open(self.output, 'w') as CLUSTERS: + for component in self.connected_components: + # uncomment to exclude single genome clusters from exporting + # if len(component) == 1: + # single_components += 1 + # continue + named_component = [self.names_map[node + 1] + for node in component] + CLUSTERS.write(','.join(named_component) + '\n') """ @@ -211,15 +147,14 @@ def export_kCluster(self): @cli.command(name="cluster", help_priority=4) @click.option('-c', '--cutoff', required=False, type=click.FloatRange(0, 1, clamp=False), default=0.0, show_default=True, help="cluster sequences with (containment > cutoff)") @click.option('-i', '--index-prefix', "index_prefix", required=True, type=click.STRING, help="Index file prefix") +@click.option('-d', '--dist-type', "distance_type", required=False, default="max_cont", show_default=True, type=click.STRING, help="select from ['min_containment', 'avg_containment', 'max_containment', 'ani']") @click.pass_context -def main(ctx, index_prefix, cutoff): +def main(ctx, index_prefix, cutoff, distance_type): """Sequence clustering.""" - - kCl = kClusters(logger_obj=ctx.obj, - index_prefix=index_prefix, cut_off_threshold=cutoff) + cutoff = float(cutoff) * 100 + kCl = Clusters(logger_obj=ctx.obj, index_prefix=index_prefix, + cut_off_threshold=cutoff, dist_type=distance_type) ctx.obj.INFO("Building the main graph...") - kCl.tsv_build_graph() + kCl.construct_graph() ctx.obj.INFO("Clustering...") - kCl.clustering() - ctx.obj.INFO("Exporting ...") - kCl.export_kCluster() + kCl.cluster_graph() diff --git a/pykSpider/kSpider2/ks_dataset_indexing.py b/pykSpider/kSpider2/ks_dataset_indexing.py index 348146c..7cf8c6e 100644 --- a/pykSpider/kSpider2/ks_dataset_indexing.py +++ b/pykSpider/kSpider2/ks_dataset_indexing.py @@ -9,9 +9,9 @@ from glob import glob -@cli.command(name="index", help_priority=2) -@click.option('--dir', "sketches_dir", required = True, help="Sketches directory (must contain only the sketches)") -@click.option('-k', '--kmer-size', "kSize", required=False, default = 0, type=click.INT, help="kmer size (only if using --sourmash)") +@cli.command(name="index", help_priority=2) +@click.option('--dir', "sketches_dir", required=True, help="Sketches directory (must contain only the sketches)") +@click.option('-k', '--kmer-size', "kSize", required=False, default=0, type=click.INT, help="kmer size (only if using --sourmash)") @click.option('--sourmash', "sourmash", is_flag=True, show_default=True, default=False, help="use sourmash sigs instead of kProcessor") @click.pass_context def main(ctx, sketches_dir, sourmash, kSize): @@ -20,20 +20,23 @@ def main(ctx, sketches_dir, sourmash, kSize): """ if not os.path.exists(sketches_dir): ctx.obj.ERROR(f"{sketches_dir} does not exist!") - + if sourmash: if not kSize: ctx.obj.ERROR(f"must select kSize when using --sourmash") - ctx.obj.INFO(f"Indexing sourmash sigs in {sketches_dir} with kSize={kSize}.") + ctx.obj.INFO( + f"Indexing sourmash sigs in {sketches_dir} with kSize={kSize}.") kSpider_internal.sourmash_sigs_indexing(sketches_dir, kSize) + ctx.obj.SUCCESS("DONE!") - else: + else: all_extra = list(glob(f"{sketches_dir}/*extra")) - all_sketches_phmap = glob(f"{sketches_dir}/*phmap") - all_sketches_mqf = glob(f"{sketches_dir}/*mqf") - + all_sketches_phmap = glob(f"{sketches_dir}/*phmap") + all_sketches_mqf = glob(f"{sketches_dir}/*mqf") + if len(all_extra) != (len(all_sketches_phmap) + len(all_sketches_mqf)): ctx.obj.ERROR(f"Inconsistent sketches files.") - + ctx.obj.INFO(f"Indexing sketches in {sketches_dir}.") kSpider_internal.index_datasets(sketches_dir) + ctx.obj.SUCCESS("DONE!") diff --git a/pykSpider/kSpider2/ks_export.py b/pykSpider/kSpider2/ks_export.py index f753eec..3a97917 100644 --- a/pykSpider/kSpider2/ks_export.py +++ b/pykSpider/kSpider2/ks_export.py @@ -40,23 +40,45 @@ def get_newick(node, parent_dist, leaf_names, newick='') -> str: @cli.command(name="export", help_priority=5) @click.option('-i', '--index-prefix', required=True, type=click.STRING, help="Index file prefix") -@click.option('--dist-mat', "distance_matrix", is_flag=True, help="Convert pairwise matrix to NxN distance matrix", default=False) +# @click.option('--dist-mat', "distance_matrix", is_flag=True, help="Convert pairwise matrix to NxN distance matrix", default=False) @click.option('--newick', "newick", is_flag=True, help="Convert pairwise (containment) matrix to newick format", default=False) -@click.option('--containment', "containment", is_flag=True, help="Use max %containment instead of number of kmers", default=True) +@click.option('-d', '--dist-type', "distance_type", required=False, default="max_cont", show_default=True, type=click.STRING, help="select from ['min_containment', 'avg_containment', 'max_containment', 'ani']") +@click.option('-o', "overwritten_output", default="na", required=False, type=click.STRING, help="custom output file name prefix") @click.pass_context -def main(ctx, index_prefix, containment, newick, distance_matrix): +def main(ctx, index_prefix, newick, distance_type, overwritten_output): """ Export kSpider pairwise to multiple formats. """ - + index_basename = os.path.basename(index_prefix) kSpider_pairwise_tsv = f"{index_prefix}_kSpider_pairwise.tsv" namesMap_file = f"{index_prefix}.namesMap" seqToKmers_tsv = f"{index_prefix}_kSpider_seqToKmersNo.tsv" + + LOGGER = ctx.obj + + distance_to_col = { + "min_cont": 3, + "avg_cont": 4, + "max_cont": 5, + "ani": 99 + } + + if distance_type not in distance_to_col: + LOGGER.ERROR("unknown distance!") + + dist_col = distance_to_col[distance_type] + if dist_col == "ani": + with open(kSpider_pairwise_tsv, 'r') as pairwise_tsv: + if "ani" not in next(pairwise_tsv).lower(): + LOGGER.ERROR("ANI was selected but was not found in the pairwise file.\nPlease, run kSpider pairwise --extend_with_ani -i script") + + + # Check for existing pairwise file for _file in [kSpider_pairwise_tsv, namesMap_file, seqToKmers_tsv]: if not os.path.exists(_file): - ctx.obj.ERROR("File {_file} is not found.") + LOGGER.ERROR(f"File {_file} is not found.") """ # Load kmer count per record @@ -83,59 +105,70 @@ def main(ctx, index_prefix, containment, newick, distance_matrix): """Parse kSpider's pairwise """ + distances = dict() labeled_out = f"kSpider_{index_basename}_pairwise.tsv" distmatrix_out = f"kSpider_{index_basename}_distmat.tsv" newick_out = f"kSpider_{index_basename}.newick" - - with open(kSpider_pairwise_tsv) as PAIRWISE, open(labeled_out, 'w') as NEW: - # Skip header - third_column = "shared_kmers" if not containment else "max_containment" - ctx.obj.INFO(f"Writing pairwise matrix to {labeled_out}") - NEW.write(f"grp1\tgrp2\t{third_column}\n") - next(PAIRWISE) - for line in PAIRWISE: - line = (line.strip().split('\t')[1:]) - origin_grp1 = line[0] - origin_grp2 = line[1] - shared_kmers = int(line[2]) - grp1 = namesMap_dict[origin_grp1] - grp2 = namesMap_dict[origin_grp2] - min_seq = float( - min(seq_to_kmers[origin_grp1], seq_to_kmers[origin_grp2])) - max_containment = (shared_kmers / min_seq) - value = shared_kmers - if containment: - value = max_containment - if distance_matrix or newick: - distances[(grp1, grp2)] = max_containment - - NEW.write(f"{grp1}\t{grp2}\t{value}\n") - - if distance_matrix or newick: - unique_ids = sorted(set([x for y in distances.keys() for x in y])) - df = pd.DataFrame(index=unique_ids, columns=unique_ids) - for k, v in distances.items(): - df.loc[k[0], k[1]] = 1-v - df.loc[k[1], k[0]] = 1-v - - df = df.fillna(0) - if distance_matrix: - ctx.obj.INFO(f"Writing distance matrix to {distmatrix_out}") - df.to_csv(distmatrix_out, sep='\t') - if newick: - if not distance_matrix: - df.to_csv(distmatrix_out, sep='\t') - loaded_df = pd.read_csv(distmatrix_out, sep='\t') - # os.remove(distmatrix_out) - ctx.obj.INFO(f"Writing newick to {newick_out}.") - names = list(loaded_df.columns[1:]) - dist = loaded_df[loaded_df.columns[1:]].to_numpy() - Z = linkage(dist, 'complete') - tree = to_tree(Z, False) - - newick = get_newick(tree, tree.dist, names) - with open(newick_out, 'w') as NW: - NW.write(newick) - - ctx.obj.SUCCESS("Done.") + + if overwritten_output != "na": + labeled_out = f"{overwritten_output}_pairwise.tsv" + distmatrix_out = f"{overwritten_output}_distmat.tsv" + newick_out = f"{overwritten_output}.newick" + + if distance_type == "ani": + with open(kSpider_pairwise_tsv) as PAIRWISE, open(labeled_out, 'w') as NEW, open(index_prefix + "_kSpider_pairwise.ani_col.tsv") as ANI: + ctx.obj.INFO(f"Writing pairwise matrix to {labeled_out}") + NEW.write(f"grp1\tgrp2\t{distance_type}\n") + # Skip header + next(PAIRWISE) + next(ANI) + for line in PAIRWISE: + line = (line.strip().split('\t')) + origin_grp1 = line[0] + origin_grp2 = line[1] + grp1 = namesMap_dict[origin_grp1] + grp2 = namesMap_dict[origin_grp2] + dist_metric = float(next(ANI).strip()) + distances[(grp1, grp2)] = dist_metric + NEW.write(f"{grp1}\t{grp2}\t{dist_metric}\n") + + else: + with open(kSpider_pairwise_tsv) as PAIRWISE, open(labeled_out, 'w') as NEW: + ctx.obj.INFO(f"Writing pairwise matrix to {labeled_out}") + NEW.write(f"grp1\tgrp2\t{distance_type}\n") + # Skip header + next(PAIRWISE) + for line in PAIRWISE: + line = (line.strip().split('\t')) + origin_grp1 = line[0] + origin_grp2 = line[1] + grp1 = namesMap_dict[origin_grp1] + grp2 = namesMap_dict[origin_grp2] + dist_metric = float(line[dist_col]) + distances[(grp1, grp2)] = dist_metric + NEW.write(f"{grp1}\t{grp2}\t{dist_metric}\n") + + unique_ids = sorted(set([x for y in distances.keys() for x in y])) + df = pd.DataFrame(index=unique_ids, columns=unique_ids) + for k, v in distances.items(): + df.loc[k[0], k[1]] = 1-v + df.loc[k[1], k[0]] = 1-v + + df = df.fillna(0) + LOGGER.INFO(f"Writing distance matrix to {distmatrix_out}") + df.to_csv(distmatrix_out, sep='\t') + + if newick: + loaded_df = pd.read_csv(distmatrix_out, sep='\t') + LOGGER.INFO(f"Writing newick to {newick_out}.") + names = list(loaded_df.columns[1:]) + dist = loaded_df[loaded_df.columns[1:]].to_numpy() + Z = linkage(dist, 'single') + tree = to_tree(Z, False) + + newick = get_newick(tree, tree.dist, names) + with open(newick_out, 'w') as NW: + NW.write(newick) + + LOGGER.SUCCESS("Done.") diff --git a/pykSpider/kSpider2/ks_pairwise.py b/pykSpider/kSpider2/ks_pairwise.py index 97d6bfb..4755e6e 100644 --- a/pykSpider/kSpider2/ks_pairwise.py +++ b/pykSpider/kSpider2/ks_pairwise.py @@ -7,18 +7,72 @@ from kSpider2.click_context import cli + @cli.command(name="pairwise", help_priority=3) @click.option('-i', '--index-prefix', required=True, type=click.STRING, help="Index file prefix") +@click.option('--estimate-ani', "ani", is_flag=True, show_default=True, default=False, help="estimate ANI and write result in a new file with single column") @click.option('-t', '--threads', "user_threads", default=1, required=False, type=int, help="number of cores") +@click.option('-s', '--scale', "sourmash_scale", required=False, default=0, type=int, help="scale used in creating sourmash sigs (only when using --estimate-ani)") @click.pass_context -def main(ctx, index_prefix, user_threads): +def main(ctx, index_prefix, user_threads, ani, sourmash_scale): """ Generate containment pairwise matrix. """ + if not ani: + ctx.obj.INFO( + f"Constructing the containment pairwise matrix using {user_threads} cores.") + kSpider_internal.pairwise(index_prefix, user_threads) + ctx.obj.SUCCESS("Done.") + else: + if user_threads > 1: + ctx.obj.WARNING("sorry, current ANI estimation does not allow multithreading") + + if not sourmash_scale: + ctx.obj.ERROR("estimating ANI requires to provide --scale value") + + # Detect kmer size + kSize = None + with open(f"{index_prefix}.extra") as extra: + kSize = int(next(extra)) + + from sourmash.distance_utils import containment_to_distance + pairwise_file = index_prefix + "_kSpider_pairwise.tsv" + ani_col = index_prefix + "_kSpider_pairwise.ani_col.tsv" + seqToKmers_tsv = f"{index_prefix}_kSpider_seqToKmersNo.tsv" + + + """ + # Load kmer count per record + """ + id_to_kmer_count = dict() + with open(seqToKmers_tsv) as KMER_COUNT: + next(KMER_COUNT) + for line in KMER_COUNT: + seq_ID, no_of_kmers = tuple(line.strip().split('\t')[1:]) + id_to_kmer_count[int(seq_ID)] = int(no_of_kmers) + + + + with open(pairwise_file) as PAIRWISE, open(ani_col, 'w') as ANI_COL: + next(PAIRWISE) + ANI_COL.write("avg_ani\n") + for or_line in PAIRWISE: + line = or_line.strip().split('\t') + shared_kmers = int(line[2]) + if shared_kmers < 5: + continue + id_1 = int(line[0]) + id_2 = int(line[1]) + min_containment = float(line[3]) + max_containment = float(line[5]) + # print("Debug:") + # print( + # (min_containment, kSize, sourmash_scale, id_to_kmer_count[]*sourmash_scale) + # ) + ani_1_in_2 = containment_to_distance(min_containment, kSize, sourmash_scale, n_unique_kmers= id_to_kmer_count[id_2]*sourmash_scale).ani + ani_2_in_1 = containment_to_distance(max_containment, kSize, sourmash_scale, n_unique_kmers= id_to_kmer_count[id_1]*sourmash_scale).ani + avg_ani = (ani_1_in_2 + ani_2_in_1) / 2.0 + new_line = f"{avg_ani}\n" + ANI_COL.write(new_line) - ctx.obj.INFO( - f"Constructing the containment pairwise matrix using {user_threads} cores.") - - kSpider_internal.pairwise(index_prefix, user_threads) - - ctx.obj.SUCCESS("Done.") + \ No newline at end of file diff --git a/setup.py b/setup.py index b869874..22ac1da 100644 --- a/setup.py +++ b/setup.py @@ -176,6 +176,7 @@ class CustomBuild(build): "Programming Language :: Python :: 3.8", "Programming Language :: Python :: 3.9", "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", ] from setuptools.command.build_py import build_py @@ -204,6 +205,8 @@ def run(self): 'pandas', 'scipy', 'numpy', + 'rustworkx', + 'tqdm' ], include_package_data=True, entry_points=''' diff --git a/src/pairwise.cpp b/src/pairwise.cpp index 47ea7e8..f53dc10 100644 --- a/src/pairwise.cpp +++ b/src/pairwise.cpp @@ -188,13 +188,31 @@ namespace kSpider { std::ofstream myfile; myfile.open(index_prefix + "_kSpider_pairwise.tsv"); - myfile << "bin_1" << '\t' << "bin_2" << '\t' << "shared_kmers" << '\t' << "max_containment" << '\n'; + myfile + << "source_1" + << "\tsource_2" + << "\tshared_kmers" + << "\tmin_containment" + << "\tavg_containment" + << "\tmax_containment" + << '\n'; uint64_t line_count = 0; for (const auto& edge : edges) { - float max_containment = (float)edge.second / min(groupID_to_kmerCount[edge.first.first], groupID_to_kmerCount[edge.first.second]); - myfile << edge.first.first << '\t' << edge.first.second << '\t' << edge.second << '\t' << max_containment << '\n'; + float cont_1_in_2 = (float)edge.second / groupID_to_kmerCount[edge.first.second]; + float cont_2_in_1 = (float)edge.second / groupID_to_kmerCount[edge.first.first]; + float min_containment = min(cont_1_in_2, cont_2_in_1); + float avg_containment = (cont_1_in_2 + cont_2_in_1) / 2.0; + float max_containment = max(cont_1_in_2, cont_2_in_1); + + myfile + << edge.first.first + << '\t' << edge.first.second + << '\t' << edge.second + << '\t' << min_containment + << '\t' << avg_containment + << '\t' << max_containment + << '\n'; } myfile.close(); - } } \ No newline at end of file diff --git a/src/sourmash_indexing.cpp b/src/sourmash_indexing.cpp index 89e7531..b79df5b 100644 --- a/src/sourmash_indexing.cpp +++ b/src/sourmash_indexing.cpp @@ -50,8 +50,16 @@ namespace kSpider { void sourmash_sigs_indexing(string sigs_dir, int selective_kSize) { kDataFrame* frame; + while (sigs_dir.size() > 0 && sigs_dir[sigs_dir.size() - 1] == '/') sigs_dir.erase(sigs_dir.size() - 1, 1); + std::string dir_prefix = sigs_dir.substr(sigs_dir.find_last_of("/\\") + 1); + while (dir_prefix.size() > 0 && dir_prefix[dir_prefix.size() - 1] == '/') { + dir_prefix.erase(dir_prefix.size() - 1, 1); + } + + cout << "dir_prefix: " << dir_prefix << endl; + flat_hash_map namesMap; string names_fileName = sigs_dir; @@ -252,7 +260,7 @@ namespace kSpider { } } - cout << " saved_kmers(~" << frame->size() << ")." << endl << endl; + cout << " saved_kmers(~" << frame->size() << ")." << endl; cout << " colors(~" << legend->size() << ")." << endl << endl; break;