diff --git a/CHANGELOG.md b/CHANGELOG.md index 8835df0..f5ffe6c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,54 @@ All notable changes to FunVIP are documented here. This project adheres to [Semantic Versioning](https://semver.org/). +## [1.0.2] - 2026-09-30 + +Faster outgroup selection and tree drawing on large databases, earlier input +validation, and a fix for runs at the default verbosity on Windows and macOS. +Analysis results are unchanged. + +### Performance +- **Outgroup selection groups the search table once.** Preparing the outgroup tasks + built a new `groupby` over the whole concatenated search table for every group, + which is O(groups x rows) and ran in the main process before the worker pool + started. On a 121,133-sequence, 5,266-genus ITS database this took about 2 h with no + log output. The table is now grouped once. At 1,500 groups and 300,000 search hits, + preparing the tasks went from 34.8 s to 0.15 s with the same task list. +- **Single-gene groups interpret and draw each tree once.** With one gene, the + concatenated tree is a copy of the gene tree, but it was interpreted and rendered a + second time. When the tree, the trimmed alignment, the partition and the query and + outgroup lists of both datasets are identical, the concatenated tree now reuses the + gene tree's interpretation, SVG and report rows. Otherwise it is interpreted as + before. On the 121,133-sequence ITS database about half of the visualize step was + this duplicate work. + +### Fixed +- **Crash on Windows and macOS at the default verbosity.** Outgroup selection and tree + interpretation handed the sequence records to their worker processes through a + module-level store that only fork-started workers inherit. Windows and macOS start + workers by spawning, so every run below `--verbose 3` stopped at the outgroup step + with `KeyError: 'list_FI'`. Each worker now receives the records when it starts. + The bundled test presets use `--verbose 3`, which is why `--test` runs passed. +- **Free text in a gene column no longer reaches GenMine.** A cell was treated as an + accession whenever an accession appeared anywhere in it, and the whole cell was + sent to GenMine, which failed on free text and accession lists. A cell is now + downloaded only when it is a single token. Other cells containing an accession stop + the run at input validation, reporting the table, line and column. FASTA-formatted + cells (starting with `>`) are read as sequences even when the header holds an + accession. +- **IDs that collide after non-ASCII replacement stop at input validation.** + `LÖ21-04` and `LO21-04`, or IDs that differ only in a Unicode hyphen, became the + same ID and were merged into one record, or failed later as a colliding genus or + datatype. Each such pair is now reported with both original strings, across DB and + query tables. An ID repeated verbatim is still merged with a warning, as before. +- IDs without any letter or digit (`--`, `...`, a lone dash) are rejected as empty + IDs. + +### Documentation +- The `--outgroupoffset` help and `docs/parameters.md` now state that the value also + discards every search hit at or below it before clustering, besides setting the + bitscore gap between ingroup and outgroup. + ## [1.0.1] - 2026-09-10 ### Fixed @@ -82,5 +130,6 @@ so rebuild the conda environment instead of `pip install --upgrade`. See the fitted regression line. - Removed confirmed-dead / broken code paths. +[1.0.2]: https://github.com/Changwanseo/FunVIP/releases/tag/v1.0.2 [1.0.1]: https://github.com/Changwanseo/FunVIP/releases/tag/v1.0.1 [1.0.0]: https://github.com/Changwanseo/FunVIP/releases/tag/v1.0.0 diff --git a/docs/parameters.md b/docs/parameters.md index fdca8e9..36e0876 100644 --- a/docs/parameters.md +++ b/docs/parameters.md @@ -89,7 +89,7 @@ do not edit by hand. Regenerate after changing the CLI._ | `--allow-innertrimming` | flag | Turn off FunVIP adjustment to not to trim inner alignment columns, default: False | | `--criterion` | str | Model-selection criterion for modeltest, one of [AIC, AICc, BIC]. default : BIC | | `--noavx` | flag | Do not use AVX for RAxML, default: False | -| `--outgroupoffset` | int | outgroupoffset value. Highering this value may select more distant outgroup, default : 20 | +| `--outgroupoffset` | int | Bitscore gap between ingroup and outgroup. Outgroup candidates must score more than this below the weakest ingroup hit, so a higher value selects a more distant outgroup. Search hits with a bitscore at or below this value are also discarded before clustering. default : 20 | | `--nosuspicious` | flag | Do not include suspicious samples for sequence-set. Mostly for metabarcoding analysis. May deduce inaccurate result with problematic database. | | `--terminate` | flag | Terminate FunVIP run when critical error detected | diff --git a/funvip/src/cluster.py b/funvip/src/cluster.py index 1d8dcd1..3495fe4 100644 --- a/funvip/src/cluster.py +++ b/funvip/src/cluster.py @@ -205,11 +205,16 @@ def cluster(FI, df_search, opt): ### Append outgroup to given group-gene dataset by search matrix -# V.list_FI shared with outgroup-append pool workers via fork copy-on-write, set -# before the Pool is created, instead of being pickled into every task tuple. +# V.list_FI shared with outgroup-append pool workers once per worker (inherited +# through fork on Linux, sent to each spawned worker on Windows and macOS) instead +# of being pickled into every task tuple. _OUTGROUP_SHARED = {} +def _share_outgroup_data(list_FI): + _OUTGROUP_SHARED["list_FI"] = list_FI + + def append_outgroup(df_search, gene, group, path, opt): logging.info(f"Appending outgroup on group: {group}, Gene: {gene}") list_FI = _OUTGROUP_SHARED["list_FI"] @@ -433,11 +438,13 @@ def outgroup_append_opt_generator(V, path, opt): # if concatenated analysis is true # concatenated + df_group = None for group in V.dict_dataset: if "concatenated" in V.dict_dataset[group]: try: - df = V.cSR - df_group = df.groupby(df["query_group"]) + if df_group is None: + df = V.cSR + df_group = df.groupby(df["query_group"]) df_group_ = df_group.get_group(group) # Generating outgroup opt for multiprocessing for gene in V.dict_dataset[group]: @@ -538,13 +545,15 @@ def pipe_cluster(V, opt, path): def pipe_append_outgroup(V, path, opt): opt_append_outgroup = outgroup_append_opt_generator(V, path, opt) - # Share V.list_FI with pool workers via fork copy-on-write (set before the Pool) - # instead of pickling the whole list into every task. - _OUTGROUP_SHARED["list_FI"] = V.list_FI + # Share V.list_FI with pool workers once per worker instead of pickling the + # whole list into every task. + _share_outgroup_data(V.list_FI) # run multiprocessing start if opt.verbose < 3: - p = mp.Pool(opt.thread) + p = mp.Pool( + opt.thread, initializer=_share_outgroup_data, initargs=(V.list_FI,) + ) result_append_outgroup = p.starmap(append_outgroup, opt_append_outgroup) p.close() p.join() diff --git a/funvip/src/command.py b/funvip/src/command.py index 60dfa60..ea75b43 100644 --- a/funvip/src/command.py +++ b/funvip/src/command.py @@ -315,7 +315,7 @@ def get_args(self) -> argparse.Namespace: ) group_advanced.add_argument( "--outgroupoffset", - help="outgroupoffset value. Highering this value may select more distant outgroup, default : 20", + help="Bitscore gap between ingroup and outgroup. Outgroup candidates must score more than this below the weakest ingroup hit, so a higher value selects a more distant outgroup. Search hits with a bitscore at or below this value are also discarded before clustering. default : 20", type=int, ) group_advanced.add_argument( diff --git a/funvip/src/tree_interpretation_pipe.py b/funvip/src/tree_interpretation_pipe.py index b86dbb5..5246e3b 100644 --- a/funvip/src/tree_interpretation_pipe.py +++ b/funvip/src/tree_interpretation_pipe.py @@ -7,11 +7,13 @@ from funvip.src.reporter import Singlereport from funvip.src.exceptions import TreeError import traceback -from copy import deepcopy +from copy import copy, deepcopy +from Bio import SeqIO import pandas as pd import re import sys import os +import filecmp import shutil import psutil import logging @@ -22,11 +24,17 @@ ### For single dataset # Input : out, group, gene, V, path, opt # Whole-run FI collections (V.dict_hash_FI / V.list_FI) shared with -# interpretation-pool workers via fork copy-on-write, set before the Pool is -# created, instead of being pickled into every (group, gene) task tuple. +# interpretation-pool workers once per worker (inherited through fork on Linux, +# sent to each spawned worker on Windows and macOS) instead of being pickled into +# every (group, gene) task tuple. _INTERP_SHARED = {} +def _share_interpretation_data(funinfo_dict, funinfo_list): + _INTERP_SHARED["funinfo_dict"] = funinfo_dict + _INTERP_SHARED["funinfo_list"] = funinfo_list + + def pipe_module_tree_interpretation( out, group, @@ -41,7 +49,7 @@ def pipe_module_tree_interpretation( ): # time_start = time() - # Read the whole-run FI collections from the fork-inherited shared store rather + # Read the whole-run FI collections from the per-worker shared store rather # than receiving a freshly pickled copy of the entire universe per task. funinfo_dict = _INTERP_SHARED["funinfo_dict"] funinfo_list = _INTERP_SHARED["funinfo_list"] @@ -647,6 +655,76 @@ def _safe_pipe_module_tree_visualization(*args): return None +def _single_gene_twin(V, path, opt, group): + genes = [g for g in V.dict_dataset[group] if g != "concatenated"] + if len(genes) != 1 or "concatenated" not in V.dict_dataset[group]: + return None + gene = genes[0] + + dataset_gene = V.dict_dataset[group][gene] + dataset_concat = V.dict_dataset[group]["concatenated"] + for attr in ("list_qr_FI", "list_og_FI"): + if [FI.hash for FI in getattr(dataset_gene, attr)] != [ + FI.hash for FI in getattr(dataset_concat, attr) + ]: + return None + + def tree_file(g): + return f"{path.out_tree}/hash/hash_{opt.runname}_{group}_{g}.nwk" + + def alignment_file(g): + return f"{path.out_alignment}/hash/{opt.runname}_hash_trimmed_{group}_{g}.fasta" + + for get_file in (tree_file, alignment_file): + a, b = get_file(gene), get_file("concatenated") + if not ( + os.path.isfile(a) + and os.path.isfile(b) + and filecmp.cmp(a, b, shallow=False) + ): + return None + + partition = V.partition.get(group) + if partition is None or partition.get("order") != [gene]: + return None + try: + width = len(next(SeqIO.parse(alignment_file(gene), "fasta")).seq) + except StopIteration: + return None + if partition.get("len", {}).get(gene) != width: + return None + + return gene + + +def _twin_interpretation(tree_info, path, opt): + group, gene = tree_info.group, tree_info.gene + out = f"{path.out_tree}/{opt.runname}_{group}" + for prefix in (f"{path.out_tree}/hash_{opt.runname}_{group}", out): + shutil.copy( + f"{prefix}_{gene}_original.svg", f"{prefix}_concatenated_original.svg" + ) + shutil.move(f"{out}_concatenated.nwk", f"{out}_concatenated_original.nwk") + shutil.copy(f"{out}_{gene}.nwk", f"{out}_concatenated.nwk") + + twin = copy(tree_info) + twin.gene = "concatenated" + twin.tree_name = f"{path.out_tree}/hash/hash_{opt.runname}_{group}_concatenated.nwk" + return twin + + +def _twin_visualization(tree_info, report_list, path, opt): + out = f"{path.out_tree}/{opt.runname}_{tree_info.group}" + shutil.copy(f"{out}_{tree_info.gene}.svg", f"{out}_concatenated.svg") + + twin_report_list = [] + for report in report_list: + twin_report = copy(report) + twin_report.update_gene("concatenated") + twin_report_list.append(twin_report) + return twin_report_list + + def pipe_tree_interpretation(V, path, opt): # Generate tree_interpretation opt to run # tree_interpretation_opt = [] @@ -661,11 +739,9 @@ def pipe_tree_interpretation(V, path, opt): funinfo_list = V.list_FI hash_dict = V.dict_hash_name - # Share the whole-run FI collections with pool workers via fork copy-on-write - # (must be set before the Pool below is created) instead of pickling them into - # every task; workers read them from _INTERP_SHARED. - _INTERP_SHARED["funinfo_dict"] = funinfo_dict - _INTERP_SHARED["funinfo_list"] = funinfo_list + # Share the whole-run FI collections with pool workers once per worker instead + # of pickling them into every task; workers read them from _INTERP_SHARED. + _share_interpretation_data(funinfo_dict, funinfo_list) # Generate options using generator def generate_interpretation_opt(): @@ -720,23 +796,59 @@ def generate_interpretation_opt(): f"Failed interpreting tree {group} {gene} because no outgroup available" ) - tree_interpretation_opt = generate_interpretation_opt() - - tree_info_list = [] + tree_interpretation_opt = list(generate_interpretation_opt()) + + task_keys = [(option[1], option[2]) for option in tree_interpretation_opt] + task_key_set = set(task_keys) + twin_gene = {} + for group in V.dict_dataset: + gene = _single_gene_twin(V, path, opt, group) + if ( + gene is not None + and (group, gene) in task_key_set + and (group, "concatenated") in task_key_set + ): + twin_gene[group] = gene + if twin_gene: + logging.info( + f"Reusing single gene tree interpretation for the concatenated tree of {len(twin_gene)} groups" + ) + run_opt = [ + option + for option in tree_interpretation_opt + if not (option[1] in twin_gene and option[2] == "concatenated") + ] ## Tree interpretation - outgroup, reconstruction(solve_flat), collapsing if opt.verbose < 3: - with mp.Pool(opt.thread) as p: - tree_info_list.extend( - p.starmap(_safe_pipe_module_tree_interpretation, tree_interpretation_opt) - ) + with mp.Pool( + opt.thread, + initializer=_share_interpretation_data, + initargs=(funinfo_dict, funinfo_list), + ) as p: + run_result = p.starmap(_safe_pipe_module_tree_interpretation, run_opt) else: # non-multithreading mode for debugging - tree_info_list = [ - pipe_module_tree_interpretation(*option) - for option in tree_interpretation_opt - ] + run_result = [pipe_module_tree_interpretation(*option) for option in run_opt] + + interpreted = { + (option[1], option[2]): tree_info + for option, tree_info in zip(run_opt, run_result) + } + twin_of = {} + tree_info_list = [] + for group, gene in task_keys: + if group in twin_gene and gene == "concatenated": + source = interpreted[(group, twin_gene[group])] + if source is None: + tree_info_list.append(None) + continue + twin = _twin_interpretation(source, path, opt) + twin_of[id(twin)] = source + tree_info_list.append(twin) + else: + tree_info_list.append(interpreted[(group, gene)]) # Drop genera that failed interpretation (guarded above) so one bad tree # does not sink the whole run; failures are logged as [TREE INTERPRETATION FAILED]. @@ -761,8 +873,10 @@ def generate_interpretation_opt(): V_dict_hash_FI[key] = (FI.original_id, FI.hash, FI.adjusted_group) # Generate options using generator + visualization_list = [ti for ti in tree_info_list if id(ti) not in twin_of] + def generate_visualization_opt(): - for tree_info in tree_info_list: + for tree_info in visualization_list: yield (tree_info, V.tup_genus, V_dict_hash_FI, path, opt) # Generate visualization option to run @@ -770,16 +884,33 @@ def generate_visualization_opt(): ## Tree visualization if opt.verbose < 3: with mp.Pool(opt.thread) as p: - tree_visualization_result = p.starmap( + visualization_result = p.starmap( _safe_pipe_module_tree_visualization, tree_visualization_opt ) else: # non-multithreading mode for debugging - tree_visualization_result = [ + visualization_result = [ pipe_module_tree_visualization(*option) for option in tree_visualization_opt ] + visualized = { + id(tree_info): report_list + for tree_info, report_list in zip(visualization_list, visualization_result) + } + tree_visualization_result = [] + for tree_info in tree_info_list: + if id(tree_info) not in twin_of: + tree_visualization_result.append(visualized[id(tree_info)]) + else: + source = twin_of[id(tree_info)] + if visualized[id(source)] is None: + tree_visualization_result.append(None) + else: + tree_visualization_result.append( + _twin_visualization(source, visualized[id(source)], path, opt) + ) + # Drop genera that failed visualization (guarded) before flattening tree_visualization_result = [r for r in tree_visualization_result if r is not None] diff --git a/funvip/src/validate_input.py b/funvip/src/validate_input.py index 1e2f149..027dc66 100644 --- a/funvip/src/validate_input.py +++ b/funvip/src/validate_input.py @@ -452,7 +452,7 @@ def _resolve_genmine(): return shutil.which("GenMine") or "GenMine" -def input_table(funinfo_dict, path, opt, table_list, datatype): +def input_table(funinfo_dict, path, opt, table_list, datatype, id_origin=None): # Whether to check if GenMine has run GenMine_flag = 0 string_error = 0 @@ -460,6 +460,9 @@ def input_table(funinfo_dict, path, opt, table_list, datatype): initialize_path(path) # this one is ugly df_list = [] + if id_origin is None: + id_origin = {} + warnings = [] errors = [] @@ -612,9 +615,17 @@ def input_table(funinfo_dict, path, opt, table_list, datatype): ): for n, _ in enumerate(df[gene]): if not (pd.isna(df[gene][n])): - if re.search(regex_genbank, df[gene][n]): - # remove unexpected indents with strip - download_set.add(df[gene][n].strip()) + # remove unexpected indents with strip + cell = str(df[gene][n]).strip() + if cell.startswith(">"): + continue + if re.search(regex_genbank, cell): + if re.fullmatch(r"[A-Za-z0-9_.]+", cell): + download_set.add(cell) + else: + errors.append( + f"In table {table}, {gene} of line {n} is neither a single GenBank accession nor a DNA sequence: {cell[:80]!r}" + ) # if NCBI accessions detected in sequence part, download it if len(download_set) > 0: @@ -728,10 +739,19 @@ def update_from_GenMine(string): raise Exception # Manage unicode for ID + raw_id = [str(x) for x in df["id"]] df["id"] = df["id"].apply( lambda x: manage_unicode(str(x), column="ID/Accession") ) + for n, key in enumerate(df["id"]): + if key in id_origin and id_origin[key][0] != raw_id[n]: + errors.append( + f"In table {table}, id {raw_id[n]!r} of line {n} becomes {key!r} after replacing non-ASCII characters, colliding with id {id_origin[key][0]!r} in table {id_origin[key][1]}. Please rename one of them" + ) + elif key not in id_origin: + id_origin[key] = (raw_id[n], table) + # To prevent errors on genus / spcies column # if this function operates with query mode, there might be no genus or species column if "genus" in df.columns: @@ -744,7 +764,7 @@ def update_from_GenMine(string): # Empty id check empty_error = [] for n, acc in enumerate(df["id"]): - if df["id"][n].strip() == "" or df["id"][n].strip() == "-": + if not any(c.isalnum() for c in df["id"][n]): empty_error.append(n) if len(empty_error) > 0: @@ -866,7 +886,7 @@ def update_from_GenMine(string): return funinfo_dict, GenMine_flag, warnings, errors -def db_input(funinfo_dict, opt, path) -> list: +def db_input(funinfo_dict, opt, path, id_origin=None) -> list: # Get DB input logging.info(f"Input DB list: {opt.db}") @@ -881,6 +901,7 @@ def db_input(funinfo_dict, opt, path) -> list: opt=opt, table_list=opt.db, datatype="db", + id_origin=id_origin, ) for warning in sorted(list(warnings)): @@ -919,7 +940,7 @@ def db_input(funinfo_dict, opt, path) -> list: return funinfo_dict, GenMine_flag -def query_input(funinfo_dict, opt, path): +def query_input(funinfo_dict, opt, path, id_origin=None): query_fasta = [ file for file in opt.query @@ -940,6 +961,7 @@ def query_input(funinfo_dict, opt, path): opt=opt, table_list=query_table, datatype="query", + id_origin=id_origin, ) funinfo_dict, fasta_warnings, fasta_errors = input_fasta( path=path, @@ -981,15 +1003,16 @@ def query_input(funinfo_dict, opt, path): # combined db and query input def data_input(V, R, opt, path): funinfo_dict = {} + id_origin = {} # get database input funinfo_dict, GenMine_flag_db = db_input( - funinfo_dict=funinfo_dict, opt=opt, path=path + funinfo_dict=funinfo_dict, opt=opt, path=path, id_origin=id_origin ) # get query input funinfo_dict, GenMine_flag_query = query_input( - funinfo_dict=funinfo_dict, opt=opt, path=path + funinfo_dict=funinfo_dict, opt=opt, path=path, id_origin=id_origin ) if GenMine_flag_db != 0 or GenMine_flag_query != 0: diff --git a/pyproject.toml b/pyproject.toml index b5a0ba2..e607ca1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "FunVIP" -version = "1.0.1" +version = "1.0.2" description = "Fungal Validation & Identification Pipeline" authors = [{name = "Changwan Seo", email = "wan101010@snu.ac.kr"}] urls = { "Homepage" = "https://github.com/Changwanseo/FunVIP" } diff --git a/tests/test_cluster.py b/tests/test_cluster.py index 5a9b33e..d0d4d0a 100644 --- a/tests/test_cluster.py +++ b/tests/test_cluster.py @@ -1,7 +1,9 @@ -"""Tests for funvip.src.cluster._majority_top_group (the group tie-break).""" +"""Tests for funvip.src.cluster: the group tie-break and outgroup task generation.""" +from types import SimpleNamespace + import pandas as pd -from funvip.src.cluster import _majority_top_group +from funvip.src.cluster import _majority_top_group, outgroup_append_opt_generator def _df(rows): @@ -31,3 +33,43 @@ def test_full_tie_returns_highest_bitscore_hit(): assert _majority_top_group(d) in {"A", "B"} # first row after the descending sort is the tiebreak assert _majority_top_group(d) == d.iloc[0]["subject_group"] + + +def test_outgroup_tasks_group_the_search_table_once(monkeypatch): + calls = [] + groupby = pd.DataFrame.groupby + + def counting_groupby(self, *args, **kwargs): + calls.append(1) + return groupby(self, *args, **kwargs) + + monkeypatch.setattr(pd.DataFrame, "groupby", counting_groupby) + cSR = pd.DataFrame( + { + "qseqid": ["HS1HE", "HS2HE", "HS3HE", "HS4HE"], + "bitscore": [300.0, 250.0, 200.0, 150.0], + "query_group": ["A", "B", "A", "C"], + } + ) + V = SimpleNamespace( + cSR=cSR, + dict_dataset={ + "A": {"its": None, "concatenated": None}, + "B": {"its": None, "tef": None, "concatenated": None}, + "C": {"its": None}, + "D": {"its": None, "concatenated": None}, + }, + ) + + tasks = outgroup_append_opt_generator(V, None, None) + + assert len(calls) == 1 + assert [(t[2], t[1]) for t in tasks] == [ + ("A", "its"), + ("A", "concatenated"), + ("B", "its"), + ("B", "tef"), + ("B", "concatenated"), + ] + assert list(tasks[0][0]["qseqid"]) == ["HS1HE", "HS3HE"] + assert list(tasks[2][0]["qseqid"]) == ["HS2HE"] diff --git a/tests/test_tree_interpretation_pipe.py b/tests/test_tree_interpretation_pipe.py new file mode 100644 index 0000000..8065805 --- /dev/null +++ b/tests/test_tree_interpretation_pipe.py @@ -0,0 +1,133 @@ +"""Tests for reusing a single gene tree interpretation for the concatenated tree.""" +import os +from types import SimpleNamespace + +import pytest + +try: + from funvip.src import tree_interpretation_pipe as tip +except Exception as e: + pytest.skip(f"tree_interpretation_pipe unavailable: {e}", allow_module_level=True) + +from funvip.src.reporter import Singlereport + +ALIGNMENT = ">HS1HE\nACGT-ACGT\n>HS2HE\nACGTTACGT\n" +TREE = "(HS1HE:0.1,HS2HE:0.2);\n" + + +def _fi(n): + return SimpleNamespace(hash=f"HS{n}HE") + + +def _dataset(qr, og): + return SimpleNamespace(list_qr_FI=qr, list_og_FI=og, list_db_FI=[]) + + +def _setup(tmp_path, genes=("its",), concat_qr=None, concat_alignment=ALIGNMENT): + path = SimpleNamespace( + out_tree=(tmp_path / "07_Tree").as_posix(), + out_alignment=(tmp_path / "05_Alignment").as_posix(), + ) + os.makedirs(f"{path.out_tree}/hash") + os.makedirs(f"{path.out_alignment}/hash") + opt = SimpleNamespace(runname="R") + + qr, og = [_fi(1)], [_fi(2)] + dict_dataset = {"G": {g: _dataset(qr, og) for g in genes}} + dict_dataset["G"]["concatenated"] = _dataset( + qr if concat_qr is None else concat_qr, og + ) + + alignments = {g: ALIGNMENT for g in genes} + alignments["concatenated"] = concat_alignment + for g, text in alignments.items(): + with open(f"{path.out_tree}/hash/hash_R_G_{g}.nwk", "w") as f: + f.write(TREE) + with open(f"{path.out_alignment}/hash/R_hash_trimmed_G_{g}.fasta", "w") as f: + f.write(text) + + width = len(ALIGNMENT.split("\n")[1]) + V = SimpleNamespace( + dict_dataset=dict_dataset, + partition={"G": {"len": {genes[0]: width}, "order": [genes[0]]}}, + ) + return V, path, opt + + +def test_single_gene_group_with_identical_inputs_is_reused(tmp_path): + V, path, opt = _setup(tmp_path) + assert tip._single_gene_twin(V, path, opt, "G") == "its" + + +def test_multi_gene_group_is_not_reused(tmp_path): + V, path, opt = _setup(tmp_path, genes=("its", "tef")) + assert tip._single_gene_twin(V, path, opt, "G") is None + + +def test_differing_query_list_is_not_reused(tmp_path): + V, path, opt = _setup(tmp_path, concat_qr=[_fi(1), _fi(3)]) + assert tip._single_gene_twin(V, path, opt, "G") is None + + +def test_differing_alignment_is_not_reused(tmp_path): + V, path, opt = _setup(tmp_path, concat_alignment=ALIGNMENT.replace("ACGT-", "ACGTA")) + assert tip._single_gene_twin(V, path, opt, "G") is None + + +def test_partition_not_covering_the_alignment_is_not_reused(tmp_path): + V, path, opt = _setup(tmp_path) + V.partition["G"]["len"]["its"] -= 1 + assert tip._single_gene_twin(V, path, opt, "G") is None + + +def test_twin_interpretation_copies_outputs_and_renames_gene(tmp_path): + out_tree = (tmp_path / "07_Tree").as_posix() + os.makedirs(out_tree) + path = SimpleNamespace(out_tree=out_tree) + opt = SimpleNamespace(runname="R") + files = { + "hash_R_G_its_original.svg": "hash svg", + "R_G_its_original.svg": "svg", + "R_G_its.nwk": "interpreted tree", + "R_G_concatenated.nwk": "tree before interpretation", + } + for name, text in files.items(): + with open(f"{out_tree}/{name}", "w") as f: + f.write(text) + source = SimpleNamespace(group="G", gene="its", tree_name="its tree") + + twin = tip._twin_interpretation(source, path, opt) + + def read(name): + with open(f"{out_tree}/{name}") as f: + return f.read() + + assert read("hash_R_G_concatenated_original.svg") == "hash svg" + assert read("R_G_concatenated_original.svg") == "svg" + assert read("R_G_concatenated_original.nwk") == "tree before interpretation" + assert read("R_G_concatenated.nwk") == "interpreted tree" + assert (twin.gene, source.gene) == ("concatenated", "its") + assert twin.tree_name.endswith("hash/hash_R_G_concatenated.nwk") + + +def test_twin_visualization_copies_svg_and_relabels_reports(tmp_path): + out_tree = (tmp_path / "07_Tree").as_posix() + os.makedirs(out_tree) + with open(f"{out_tree}/R_G_its.svg", "w") as f: + f.write("") + report = Singlereport() + report.hash, report.gene, report.species_assigned = "HS1HE", "its", "Genus species" + source = SimpleNamespace(group="G", gene="its") + + twin_reports = tip._twin_visualization( + source, + [report], + SimpleNamespace(out_tree=out_tree), + SimpleNamespace(runname="R"), + ) + + with open(f"{out_tree}/R_G_concatenated.svg") as f: + assert f.read() == "" + assert [r.gene for r in twin_reports] == ["concatenated"] + assert report.gene == "its" + assert vars(twin_reports[0]) == {**vars(report), "gene": "concatenated"} diff --git a/tests/test_validate_input.py b/tests/test_validate_input.py new file mode 100644 index 0000000..410b40f --- /dev/null +++ b/tests/test_validate_input.py @@ -0,0 +1,177 @@ +"""Tests for input-table validation in funvip.src.validate_input.input_table.""" +import csv +import os +from types import SimpleNamespace + +import pytest + +from funvip.src import validate_input +from funvip.src.validate_input import input_table + +SEQ = "ACGT" * 30 + + +def _path(tmp_path): + for d in ("db", "query", "GenMine"): + (tmp_path / d).mkdir(exist_ok=True) + return SimpleNamespace( + genusdb=os.path.join( + os.path.dirname(validate_input.__file__), "..", "data", "genus_line.txt" + ), + out_db=str(tmp_path / "db"), + out_query=str(tmp_path / "query"), + GenMine=str(tmp_path / "GenMine"), + ) + + +def _opt(email=""): + return SimpleNamespace(level="genus", gene=["its"], email=email, tableformat="csv") + + +def _table(tmp_path, name, rows): + p = tmp_path / name + with open(p, "w", newline="", encoding="utf-8") as f: + w = csv.writer(f, quoting=csv.QUOTE_ALL) + w.writerow(["id", "genus", "species", "its"]) + w.writerows(rows) + return p.as_posix() + + +def _run(tmp_path, tables, datatype="db", email="", shared=None): + kwargs = {} + if shared is not None: + kwargs["id_origin"] = shared[0] + funinfo_dict, genmine_flag, warnings, errors = input_table( + funinfo_dict={} if shared is None else shared[1], + path=_path(tmp_path), + opt=_opt(email), + table_list=tables, + datatype=datatype, + **kwargs, + ) + return funinfo_dict, genmine_flag, warnings, [e for e in errors if e] + + +@pytest.fixture +def no_genmine(monkeypatch): + def fail(*args, **kwargs): + raise AssertionError("GenMine must not be called") + + monkeypatch.setattr(validate_input.subprocess, "call", fail) + + +def test_single_token_accessions_are_sent_to_genmine(tmp_path, monkeypatch): + # Formats the strict patterns do not cover (2 letters + 8 digits, RefSeq with + # 7 digits) must keep reaching GenMine as before. + monkeypatch.setattr(validate_input, "sleep", lambda s: None) + monkeypatch.setattr(validate_input.subprocess, "call", lambda *a, **k: 1) + table = _table( + tmp_path, + "db.csv", + [ + ("A1", "Aspergillus", "terreus", "MN123456"), + ("A2", "Aspergillus", "terreus", "PQ12345678"), + ("A3", "Aspergillus", "terreus", "NR_1234567.1"), + ("A4", "Aspergillus", "terreus", " MN654321.1 "), + ("A5", "Aspergillus", "terreus", SEQ), + ], + ) + with pytest.raises(Exception): + _run(tmp_path, [table], email="test@example.org") + + with open(tmp_path / "GenMine" / "Accessions.txt") as f: + sent = sorted(line.strip() for line in f) + assert sent == ["MN123456", "MN654321.1", "NR_1234567.1", "PQ12345678"] + + +def test_free_text_and_accession_lists_stop_at_validation(tmp_path, no_genmine): + table = _table( + tmp_path, + "db.csv", + [ + ("A1", "Aspergillus", "terreus", SEQ), + ("A2", "Aspergillus", "terreus", "ITS from MN123456 and MN123457"), + ("A3", "Aspergillus", "terreus", "MN123456, MN123457"), + ("A4", "Aspergillus", "terreus", "MN123456;MN123457"), + ], + ) + _, genmine_flag, _, errors = _run(tmp_path, [table]) + + assert genmine_flag == 0 + rejected = [e for e in errors if "neither a single GenBank accession" in e] + assert len(rejected) == 3 + assert any("its of line 1 " in e and "ITS from MN123456" in e for e in rejected) + + +def test_fasta_cell_with_accession_header_is_read_as_sequence(tmp_path, no_genmine): + table = _table( + tmp_path, + "db.csv", + [("A1", "Aspergillus", "terreus", f">MN123456.1 Aspergillus terreus\n{SEQ}")], + ) + funinfo_dict, genmine_flag, _, errors = _run(tmp_path, [table]) + + assert errors == [] + assert genmine_flag == 0 + assert funinfo_dict["A1"].seq["its"] == SEQ + + +@pytest.mark.parametrize( + "first, second", + [("LÖ21-04", "LO21-04"), ("GB‐0065942", "GB-0065942")], +) +def test_ids_colliding_after_unicode_replacement_stop_at_validation( + tmp_path, no_genmine, first, second +): + table = _table( + tmp_path, + "db.csv", + [ + (first, "Aspergillus", "terreus", SEQ), + (second, "Aspergillus", "terreus", "ACGT" * 25), + ], + ) + _, _, _, errors = _run(tmp_path, [table]) + + collisions = [e for e in errors if "colliding with id" in e] + assert len(collisions) == 1 + assert repr(first) in collisions[0] and repr(second) in collisions[0] + + +def test_ids_colliding_across_db_and_query_tables(tmp_path, no_genmine): + shared = ({}, {}) + db = _table(tmp_path, "db.csv", [("LO21-04", "Aspergillus", "terreus", SEQ)]) + query = _table(tmp_path, "query.csv", [("LÖ21-04", "", "", SEQ)]) + + _, _, _, db_errors = _run(tmp_path, [db], shared=shared) + _, _, _, query_errors = _run(tmp_path, [query], datatype="query", shared=shared) + + assert db_errors == [] + assert any("colliding with id 'LO21-04' in table" in e for e in query_errors) + + +def test_identical_id_in_two_tables_is_still_merged(tmp_path, no_genmine): + first = _table(tmp_path, "db1.csv", [("A1", "Aspergillus", "terreus", SEQ)]) + second = _table(tmp_path, "db2.csv", [("A1", "Aspergillus", "terreus", SEQ)]) + + funinfo_dict, _, warnings, errors = _run(tmp_path, [first, second]) + + assert errors == [] + assert list(funinfo_dict) == ["A1"] + assert any("Duplicate id A1" in w for w in warnings) + + +def test_ids_without_letters_or_digits_are_rejected(tmp_path, no_genmine): + table = _table( + tmp_path, + "db.csv", + [ + ("A1", "Aspergillus", "terreus", SEQ), + ("--", "Aspergillus", "terreus", SEQ), + ("–", "Aspergillus", "terreus", SEQ), + ("...", "Aspergillus", "terreus", SEQ), + ], + ) + _, _, _, errors = _run(tmp_path, [table]) + + assert any("Empty id found, line [1, 2, 3]" in e for e in errors)