diff --git a/.github/workflows/apptainer.yml b/.github/workflows/apptainer.yml index 281d631..d8f3316 100644 --- a/.github/workflows/apptainer.yml +++ b/.github/workflows/apptainer.yml @@ -79,4 +79,4 @@ jobs: - name: Run Apptainer tests run: | sequana_taxonomy --download toydb - sequana_multitax --input-directory test/data/ --apptainer-prefix ~/.apptainer/cache --databases ~/.config/sequana/kraken2_dbs/toydb && cd multitax && bash multitax.sh + sequana_multitax --input-directory test/data/ --exclude-pattern empty --apptainer-prefix ~/.apptainer/cache --databases ~/.config/sequana/kraken2_dbs/toydb && cd multitax && bash multitax.sh diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..a35305a --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,21 @@ +# Changelog + +All notable changes to this project will be documented in this file. + +## [0.15.1] - 2026-05-21 + +### Added +- Improved test coverage with 8 new CLI tests (up from 3 to 11 total) + - 5 new CLI config verification tests: paired reads, multiple databases, store-unclassified, kraken-confidence, blast options + - Tests follow lora pipeline pattern using yaml.safe_load for config verification + - Tests for paired-end data handling and multiple database configuration + +### Fixed +- Fix multitax.rules: ensure unclassified.fastq output always exists + - Wrapper shell script now touches unclassified.fastq file after sequana_taxonomy completes + - Fixes missing output errors when store_unclassified=False + +## [0.15.0] - 2024-08-07 + +### Added +- Initial stable release with Kraken2 and multiple database support diff --git a/pyproject.toml b/pyproject.toml index 8616aea..5ba276b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "poetry.core.masonry.api" [tool.poetry] name = "sequana-multitax" -version = "0.15.0" +version = "0.15.1" description = "A multi-sample and multi-databases taxonomic analysis using Kraken" authors = ["Sequana Team"] license = "BSD-3" diff --git a/sequana_pipelines/multitax/blast.py b/sequana_pipelines/multitax/blast.py index 6ffe4b8..45def6d 100644 --- a/sequana_pipelines/multitax/blast.py +++ b/sequana_pipelines/multitax/blast.py @@ -10,8 +10,9 @@ # Documentation: http://sequana.readthedocs.io # Contributors: https://github.com/sequana/sequana/graphs/contributors ############################################################################## -import pandas as pd from collections import defaultdict + +import pandas as pd from sequana.taxonomy import Taxonomy @@ -44,10 +45,10 @@ def parse_blast(filename): df["taxids"] = df["taxids"].astype(str) # On vérifie si les 1000 reads initiaux ont été blastés - #len_diff = (Ntotal - len(df["qseqid"].unique())) / 10 + # len_diff = (Ntotal - len(df["qseqid"].unique())) / 10 # Percentage of unclassified reads not in blast - #print(f"percentage of input reads not classified by Blast: {len_diff}% ") + # print(f"percentage of input reads not classified by Blast: {len_diff}% ") # keep only the first taxids if there are several in the same row knowing # they belong to the same lineage @@ -62,16 +63,13 @@ def parse_blast(filename): for qseqid in qseqid_set: bitscore_max = df.loc[(df["qseqid"] == str(qseqid))]["bitscore"].max() df.drop( - df[ - ((df["qseqid"] == str(qseqid)) & (df["bitscore"] < int(bitscore_max))) - ].index, + df[((df["qseqid"] == str(qseqid)) & (df["bitscore"] < int(bitscore_max)))].index, inplace=True, ) return df - def read_taxonomy(filename=None): # Reads ID, parent ID, scientific names # Read the sequana taxonomy file @@ -80,22 +78,54 @@ def read_taxonomy(filename=None): tax = Taxonomy() filename = tax.database - with open(filename, "r") as f: - current_key = None - TAXID = defaultdict(list) + # Handle both old text format (taxonomy.dat) and new CSV format (taxonomy.csv.gz) + if filename.endswith(".csv.gz") or filename.endswith(".csv"): + import gzip - for line in f.readlines(): + open_func = gzip.open if filename.endswith(".gz") else open + mode = "rt" if filename.endswith(".gz") else "r" - if " : " in line: - deb = line.split(" : ") - (key, value) = (deb[0].strip(), deb[1].strip()) - - if key == "ID": - current_key = value - else: - TAXID[current_key].append(value) - - return TAXID + TAXID = defaultdict(list) + with open_func(filename, mode) as f: + header = f.readline().strip().split(",") + id_idx = header.index("id") + parent_idx = header.index("parent") + rank_idx = header.index("rank") + name_idx = header.index("scientific_name") + + for line in f: + line = line.rstrip("\n") + # Split only on first 3 commas, rest is the scientific name + parts = line.split(",", 3) + if len(parts) >= 4: + taxid = str(parts[id_idx]) + parent = str(parts[parent_idx]) + rank = str(parts[rank_idx]) + name = str(parts[name_idx]) + TAXID[taxid] = [parent, rank, name] + return TAXID + else: + import gzip + + open_func = gzip.open if filename.endswith(".gz") else open + mode = "rt" if filename.endswith(".gz") else "r" + + with open_func(filename, mode) as f: + current_key = None + TAXID = defaultdict(list) + + for line in f.readlines(): + + if " : " in line: + deb = line.split(" : ") + (key, value) = (deb[0].strip(), deb[1].strip()) + + if key == "ID": + current_key = value + else: + TAXID[current_key].append(value) + + return TAXID # OPTIONNEL : RECUPERATION DES TAXIDS DEPUIS LES ACCESSION NUMBERS @@ -139,6 +169,7 @@ def taxidstolineage(taxid_set): "order", "class", "phylum", + "domain", # Modern taxonomy uses "domain" instead of "superkingdom" for Eukaryota "superkingdom", ) # Peut se modifier en fonction de ce que l'on veut @@ -152,16 +183,10 @@ def taxidstolineage(taxid_set): # ranks attribue à chaque head_rank son sub_group de ranks taxid_sup = taxid for head_rank in head_ranks: - if ( - TAXID[str(taxid_sup)][1] == head_rank - ): # On est déjà au niveau du head_rank qui nous intéresse - value_rank = TAXID[str(taxid_sup)][ - 2 - ] # On prend le nom scientifique au niveau du head rank + if TAXID[str(taxid_sup)][1] == head_rank: # On est déjà au niveau du head_rank qui nous intéresse + value_rank = TAXID[str(taxid_sup)][2] # On prend le nom scientifique au niveau du head rank taxid_rank = taxid_sup - taxid_sup = TAXID[str(taxid_sup)][ - 0 - ] # On remonte de head rank en head rank + taxid_sup = TAXID[str(taxid_sup)][0] # On remonte de head rank en head rank else: # Il faut prendre en compte que certains rangs ne sont pas classés dans la hiérarchie des rangs mais que l'on peut tout de mâme éventuellement remonter et tomber sur le head_rank d'intérêt. Il s'agit des rangs suivants : 'no_rank', 'clade', 'genotype', 'pathogroup', 'serotype', 'serogroup'. taxid_bis = TAXID[str(taxid_sup)][0] @@ -207,9 +232,14 @@ def get_LCA(filename): df["Order"] = df.apply(lambda row: TAXID2LIN[row.taxids][8], axis=1) df["Class"] = df.apply(lambda row: TAXID2LIN[row.taxids][10], axis=1) df["Phylum"] = df.apply(lambda row: TAXID2LIN[row.taxids][12], axis=1) - df["Superkingdom"] = df.apply(lambda row: TAXID2LIN[row.taxids][14], axis=1) + df["Domain"] = df.apply( + lambda row: TAXID2LIN[row.taxids][14] if len(TAXID2LIN[row.taxids]) > 14 else "None", axis=1 + ) + df["Superkingdom"] = df.apply( + lambda row: TAXID2LIN[row.taxids][16] if len(TAXID2LIN[row.taxids]) > 16 else TAXID2LIN[row.taxids][14], axis=1 + ) - # deal with uncultured bacteria + # deal with uncultured bacteria df["SpeciesName"] = df.apply(lambda row: TAXID2LIN[row.taxids][3], axis=1) Rank_Names = [ @@ -228,15 +258,14 @@ def get_LCA(filename): # We get back the LCA for each subset of duplicated entryies (gp de qseqid) # and store them in SEQID2LCA for qseqid in qseqid_set: + taxid_LCA = "None" for rank in Rank_Names: taxid_unique = df.loc[(df["qseqid"] == str(qseqid))][rank].unique() # deal with uncultured bacteria - name_species_unique = df.loc[(df["qseqid"] == str(qseqid))][ - "SpeciesName" - ].unique() + name_species_unique = df.loc[(df["qseqid"] == str(qseqid))]["SpeciesName"].unique() if len(taxid_unique) == 1: @@ -253,9 +282,7 @@ def get_LCA(filename): # special case of uncultered bacteria --> we keep the bacteria # at the level of superkingdom - elif rank == "Superkingdom" and "uncultured" in str( - name_species_unique - ): + elif rank == "Superkingdom" and "uncultured" in str(name_species_unique): taxid_LCA = "2" break @@ -317,7 +344,15 @@ def remove_duplicates(filename): df["Order"] = df.apply(lambda row: TAXID2LIN[row.taxid_LCA][9], axis=1) df["Class"] = df.apply(lambda row: TAXID2LIN[row.taxid_LCA][11], axis=1) df["Phylum"] = df.apply(lambda row: TAXID2LIN[row.taxid_LCA][13], axis=1) - df["Superkingdom"] = df.apply(lambda row: TAXID2LIN[row.taxid_LCA][15], axis=1) + df["Domain"] = df.apply( + lambda row: TAXID2LIN[row.taxid_LCA][15] if len(TAXID2LIN[row.taxid_LCA]) > 15 else "None", axis=1 + ) + df["Superkingdom"] = df.apply( + lambda row: TAXID2LIN[row.taxid_LCA][17] + if len(TAXID2LIN[row.taxid_LCA]) > 17 + else TAXID2LIN[row.taxid_LCA][15], + axis=1, + ) # On exporte les résultats au format .csv df.to_csv("{}_Results.csv".format(filename), index=False) @@ -326,33 +361,61 @@ def remove_duplicates(filename): def krona(filename, output): - """krona visulation""" + """krona visualization""" + import os df = pd.read_csv(filename, sep=",", header="infer") # output of remove_duplicates + # Count total reads from unclassified_subsample.fasta + fasta_dir = os.path.dirname(filename) + fasta_path = os.path.join(os.path.dirname(fasta_dir), "kraken", "unclassified_subsample.fasta") + total_reads = 0 + if os.path.exists(fasta_path): + with open(fasta_path) as f: + total_reads = sum(1 for line in f if line.startswith(">")) + # We keep only the first ranks - df = df.drop(columns=["qseqid", "bitscore", "taxid_LCA", "SpeciesName", "Strain"]) + cols_to_drop = ["qseqid", "bitscore", "taxid_LCA", "SpeciesName", "Strain"] + available_cols = [c for c in cols_to_drop if c in df.columns] + df = df.drop(columns=available_cols) + + # Use Domain if available (modern taxonomy), otherwise use Superkingdom + top_rank = "Domain" if "Domain" in df.columns else "Superkingdom" # Se the correct format for krona (ktImportText) - df = ( - df.groupby( - ["Superkingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species"] - ) - .size() - .reset_index(name="Count") - ) - df = df.reindex( - columns=[ - "Count", - "Superkingdom", - "Phylum", - "Class", - "Order", - "Family", - "Genus", - "Species", - ] - ) + groupby_cols = [top_rank, "Phylum", "Class", "Order", "Family", "Genus", "Species"] + # Filter to only columns that exist + groupby_cols = [c for c in groupby_cols if c in df.columns] + + df = df.groupby(groupby_cols).size().reset_index(name="Count") + + # Add unclassified reads if we know the total + if total_reads > 0: + classified_reads = df["Count"].sum() + unclassified_count = total_reads - classified_reads + if unclassified_count > 0: + unclassified_row = {col: "Unclassified" for col in groupby_cols} + unclassified_row["Count"] = unclassified_count + df = pd.concat([df, pd.DataFrame([unclassified_row])], ignore_index=True) + + # Add percentage column + total = df["Count"].sum() + df["Percentage"] = (df["Count"] / total * 100).round(2) + + # Reindex with the columns that actually exist + reindex_cols = [ + "Count", + "Percentage", + top_rank, + "Phylum", + "Class", + "Order", + "Family", + "Genus", + "Species", + ] + reindex_cols = [c for c in reindex_cols if c in df.columns] + df = df[reindex_cols] df.to_csv("{}".format(filename), index=False) diff --git a/sequana_pipelines/multitax/config.yaml b/sequana_pipelines/multitax/config.yaml index 29125a7..49d6a3c 100644 --- a/sequana_pipelines/multitax/config.yaml +++ b/sequana_pipelines/multitax/config.yaml @@ -10,10 +10,12 @@ input_readtag: _R[12]_ input_pattern: "*fastq.gz" apptainers: - sequana_tools: "https://zenodo.org/record/18257162/files/sequana_tools_26.1.14.img" graphviz: "https://zenodo.org/record/7928262/files/graphviz_7.0.5.img" multiqc: "https://zenodo.org/record/10205070/files/multiqc_1.16.0.img" - sequana: "https://zenodo.org/record/18257162/files/sequana_tools_26.1.14.img" + sequana: "https://zenodo.org/record/20467327/files/sequana_0.22.0.img" + krona: "https://zenodo.org/record/20451801/files/krona_2.8.1.img" + awk: "https://zenodo.org/record/20451789/files/awk_5.3.2.img" + blast: "https://zenodo.org/record/7848524/files/blast_2.12.0.img" ############################################################################## # sequana_taxonomy diff --git a/sequana_pipelines/multitax/main.py b/sequana_pipelines/multitax/main.py index 17751f1..ba802c1 100755 --- a/sequana_pipelines/multitax/main.py +++ b/sequana_pipelines/multitax/main.py @@ -56,15 +56,17 @@ def update_taxonomy(ctx, param, value): # callback for --databases multiple arguments def check_databases(ctx, param, value): - if value: - # click transform the input databases (tuple) into a string - # we need to convert it back to a tuple before checking the databases - values = ast.literal_eval(value) - for db in values: - if not os.path.exists(db): - click.echo(f"{db} does not exists. Check its path name") - sys.exit(1) - return ast.literal_eval(value) + if not value: + # e.g. with --from-project, --databases may not be provided + return value + # click transform the input databases (tuple) into a string + # we need to convert it back to a tuple before checking the databases + values = ast.literal_eval(value) + for db in values: + if not os.path.exists(db): + click.echo(f"{db} does not exists. Check its path name") + sys.exit(1) + return values @click.command(context_settings=help) @@ -93,7 +95,7 @@ def check_databases(ctx, param, value): type=click.STRING, cls=OptionEatAll, callback=check_databases, - required="--update-taxonomy" not in sys.argv, + required="--update-taxonomy" not in sys.argv and "--from-project" not in sys.argv, help="""Path to a valid Kraken database(s). See sequana_taxonomy standaline to download some. You may use several, in which case, an iterative taxonomy is performed as explained in online sequana @@ -118,10 +120,6 @@ def check_databases(ctx, param, value): ) def main(**options): - if options["from_project"]: - click.echo("--from-project Not yet implemented") - sys.exit(1) - # the real stuff is here manager = SequanaManager(options, NAME) manager.setup() @@ -161,7 +159,7 @@ def fill_do_blast_unclassified(): fill_kraken_confidence() if "--store-unclassified" in sys.argv: fill_store_unclassified() - if "--do_blast_unclassified" in sys.argv: + if "--do-blast-unclassified" in sys.argv: fill_do_blast_unclassified() else: fill_databases() diff --git a/sequana_pipelines/multitax/multitax.rules b/sequana_pipelines/multitax/multitax.rules index 6bfa541..3a4c896 100644 --- a/sequana_pipelines/multitax/multitax.rules +++ b/sequana_pipelines/multitax/multitax.rules @@ -36,10 +36,11 @@ if len(config['sequana_taxonomy']["databases"]) > 1: -if config['blast']['do'] and not config['sequana_taxonomy']['store_unclassified']: - from sequana import logger - logger.error("You set blast section so you must set sequana_taxonomy/store_unclassified to true in the configuration file") - sys.exit(1) +# blast needs unclassified reads, so auto-enable store_unclassified when blast.do is set +_store_unclassified = bool( + config["sequana_taxonomy"].get("store_unclassified", False) + or config["blast"].get("do", False) +) @@ -50,7 +51,8 @@ if config['blast']['do']: expand("{sample}/kraken/unclassified_subsample.fasta", sample=sorted(manager.samples)), expand("{sample}/blast/blast.tsv", sample=sorted(manager.samples)), expand("{sample}/blast/krona.txt", sample=sorted(manager.samples)), - expand("{sample}/blast/text.krona.html", sample=sorted(manager.samples)) + expand("{sample}/blast/text.krona.html", sample=sorted(manager.samples)), + expand("{sample}/.blast_summary_updated", sample=sorted(manager.samples)) else: rule all: input: @@ -60,42 +62,54 @@ else: def get_sequana_taxonomy_unclassified(): if manager.paired: - return "{sample}/kraken/unclassified1.fastq" + return "{sample}/kraken/unclassified.fastq" else: return "{sample}/kraken/unclassified.fastq" # ================================================== Taxonomy # -def get_sequana_taxonomy_taxonomy_dep(): - - # in theory, users should call sequana_taxonomy manually, - # but for the singularity case, the taxonomy.dat file will - # be download by the image in .config/sequana. This is fine - # and works. However, on a cluster, this means access to network - # so it should be run locally. Yet, the rule sequana_taxonomy - # required resources so it should be a cluster node that has no - # network access...if we use singularity, we can therefore - # download the file manually and copy it in .config/sequana/taxonomy - # or let this rule work locally +def _using_apptainer(): + # snakemake <8 exposed workflow.use_singularity / workflow.use_apptainer. + # snakemake >=8 dropped those in favour of + # workflow.deployment_settings.deployment_method (a set of DeploymentMethod enums). if getattr(workflow, "use_singularity", False) or getattr(workflow, "use_apptainer", False): + return True + ds = getattr(workflow, "deployment_settings", None) + methods = getattr(ds, "deployment_method", ()) if ds is not None else () + return any("APPTAINER" in str(m).upper() or "SINGULARITY" in str(m).upper() for m in methods) + + +def get_sequana_taxonomy_taxonomy_dep(): + # The sequana_taxonomy rule needs the NCBI taxonomy database + # (~/.config/sequana/taxonomy.csv.gz). When running inside a container, + # snakemake forces --home to the working directory, so sequana looks for the + # database under /.config/sequana and, if missing, tries to download + # it. On a cluster the resource-heavy sequana_taxonomy rule runs on a compute + # node that usually has no network access, so the download fails. + # + # To avoid that, when apptainer/singularity is used we stage the database into + # the working directory up-front: either by copying the user's existing cache + # (done at parse time on the host, which has the real home and network) or via + # a local download_taxonomy rule (kept local so it runs on the submit node). + if _using_apptainer(): + + os.makedirs(".config/sequana", exist_ok=True) home = Path(".").home() - if (home / Path(".config/sequana/taxonomy.dat")).exists(): - os.makedirs(".config/", exist_ok=True) - os.makedirs(".config/sequana", exist_ok=True) - shell("cp ~/.config/sequana/taxonomy.dat .config/sequana/") + if (home / Path(".config/sequana/taxonomy.csv.gz")).exists(): + shell("cp ~/.config/sequana/taxonomy.csv.gz .config/sequana/") else: rule download_taxonomy: input: - output: ".config/sequana/taxonomy.dat" + output: ".config/sequana/taxonomy.csv.gz" container: - config["apptainers"]["sequana_tools"] + config["apptainers"]["sequana"] shell: """ sequana_taxonomy --update-taxonomy """ - return ".config/sequana/taxonomy.dat" + return ".config/sequana/taxonomy.csv.gz" else: return [] @@ -106,11 +120,6 @@ _confidence = config["sequana_taxonomy"].get("confidence", 0) if _confidence: _taxonomy_options += f" --confidence {_confidence}" _taxonomy_options += f" --level {config['sequana_taxonomy'].get('level', 'INFO')}" -if config["sequana_taxonomy"].get("store_unclassified", False): - if manager.paired: - _taxonomy_options += " --unclassified-out unclassified#.fastq" - else: - _taxonomy_options += " --unclassified-out unclassified.fastq" rule sequana_taxonomy: @@ -127,6 +136,7 @@ rule sequana_taxonomy: config['sequana_taxonomy']['threads'] params: databases=config["sequana_taxonomy"]["databases"], + store_unclassified=_store_unclassified, options=_taxonomy_options log: "{sample}/kraken/sequana_taxonomy.log" @@ -149,7 +159,7 @@ if config['blast']['do']: params: nreads=config['blast']['nreads'] container: - config["apptainers"]["sequana_tools"] + config["apptainers"]["awk"] shell: """ awk '{{if(NR%4==1) {{printf(">%s\\n",substr($0,2));}} else if(NR%4==2) print;}}' {input} | head -n {params.nreads} > {output} || echo "wierd" @@ -165,7 +175,7 @@ if config['blast']['do']: fields="qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore staxids" threads: 4 container: - config["apptainers"]["sequana_tools"] + config["apptainers"]["blast"] shell: """ # requires BLASTDB to be defined ! @@ -193,12 +203,29 @@ if config['blast']['do']: output: "{sample}/blast/text.krona.html" container: - config["apptainers"]["sequana_tools"] + config["apptainers"]["krona"] shell: """ ktImportText {input} -o {output} """ + # Regenerate the sample summary.html to include the blast section + # Required because sequana_taxonomy creates summary.html before blast runs + rule regenerate_summary_with_blast: + input: + blast_csv="{sample}/blast/blast_summary.csv", + krona_html="{sample}/blast/text.krona.html", + html="{sample}/summary.html" + output: + touch("{sample}/.blast_summary_updated") + run: + from sequana.modules_report.kraken import KrakenModule + from sequana.utils import config as cfg + from pathlib import Path + sample_dir = Path(input.html).parent + cfg.output_dir = str(sample_dir) + KrakenModule(sample_dir / "kraken", output_filename="summary.html") + # ================================================== summary plot # rule summary_plot: @@ -307,8 +334,10 @@ rule dot2svg: # -# Those rules takes a couple of seconds so no need for a cluster -localrules: multiqc, rulegraph, summary_plot, download_taxonomy +# Those rules takes a couple of seconds so no need for a cluster. +# download_taxonomy is kept local so it runs on the submit node, which has +# network access (compute nodes generally do not). +localrules: multiqc, rulegraph, summary_plot, regenerate_summary_with_blast, download_taxonomy @@ -326,28 +355,26 @@ onsuccess: from sequana.modules_report.kraken import KrakenModule from sequana.utils.datatables_js import DataTable - intro = """

General Information

{}""".format(manager.get_html_summary()) - intro += """
This pipeline summarizes the taxonomic analysis - made on {} samples. Here below you can find a summary - of the analysis. The first plot (left) shows the proportion of - reads classified by the databases in each found kingdom (Eukaryota, - Bacteria, Viruses, Archea). If several databases were used, - you should also see a plot (right panel) with proportion of reads - classified in each database (rather than by kingdom). - The per-sample analysis are made with kraken. A multi-sample - report is also available as + intro = """

General Information

+
This pipeline summarizes the taxonomic analysis + made on {} samples. Below you can find a summary of the analysis. + The first plot (left) shows the proportion of reads classified by + the databases in each found kingdom (Eukaryota, Bacteria, Viruses, + Archaea). If several databases were used, you should also see a plot + (right panel) with the proportion of reads classified in each database + (rather than by kingdom). Per-sample analyses are performed with Kraken. + A multi-sample report is also available as a sequana multiqc report. - Individual report are also browsable (Go to the - Individual taxonomic - reports section here below).
""".format(len(manager.samples)) + Individual reports are also browsable (see the + Individual taxonomic reports section below).
""".format(len(manager.samples)) intro += """
""" intro += """

Analysis overview

-

Proportion of reads in virus, bacteria, human and unclassified categories is - shown here below. The underlying data can be downloaded from - sequana_kraken_summary.json file - or here below as a CSV file.

+

The proportion of reads in virus, bacteria, human, and unclassified categories + is shown below. The underlying data can be downloaded from the + sequana_kraken_summary.json file + or as a CSV file below.

""" img = SequanaReport.png_to_embedded_png("dummy_self", "outputs/proportion_kraken.png", @@ -373,6 +400,9 @@ in sequana_kraken_dbs_summary. df.reset_index(inplace=True) df['link'] = links df = df.rename({"index": "sample"}, axis=1) + # Sanitize column names to avoid DataTables errors with special characters (colons, dots, etc) + # Dots break DataTables because it interprets them as nested-object accessors. + df.columns = [str(col).replace(":", "_").replace("/", "_").replace(".", "_") for col in df.columns] datatable = DataTable(df, 'kraken', index=False) datatable.datatable.datatable_options = {'paging': 'false', 'buttons': ['copy', 'csv'], @@ -394,6 +424,9 @@ in sequana_kraken_dbs_summary. df.reset_index(inplace=True) df['link'] = links df = df.rename({"index": "sample"}, axis=1) + # Sanitize column names to avoid DataTables errors with special characters (colons, dots, etc) + # Dots break DataTables because it interprets them as nested-object accessors. + df.columns = [str(col).replace(":", "_").replace("/", "_").replace(".", "_") for col in df.columns] datatable = DataTable(df, 'kraken_kingdom', index=False) datatable.datatable.datatable_options = {'paging': 'false', 'buttons': ['copy', 'csv'], diff --git a/test/data/paired/data_R1_.fastq.gz b/test/data/paired/data_R1_.fastq.gz new file mode 100644 index 0000000..22fdbfc Binary files /dev/null and b/test/data/paired/data_R1_.fastq.gz differ diff --git a/test/data/paired/data_R2_.fastq.gz b/test/data/paired/data_R2_.fastq.gz new file mode 100644 index 0000000..22fdbfc Binary files /dev/null and b/test/data/paired/data_R2_.fastq.gz differ diff --git a/test/test_main.py b/test/test_main.py index 8b7046b..0b0dad4 100644 --- a/test/test_main.py +++ b/test/test_main.py @@ -1,8 +1,12 @@ import os +import shutil import subprocess import sys import tempfile +from pathlib import Path +import pytest +import yaml from click.testing import CliRunner from sequana_pipelines.multitax.main import main @@ -11,6 +15,8 @@ sharedir = f"{test_dir}/data" krakendb = f"{test_dir}/data/krakendb" +paireddir = f"{test_dir}/data/paired" +simpledir = f"{test_dir}/data/simple" def test_standalone_subprocess(): @@ -50,3 +56,122 @@ def test_standalone_script(): def test_version(): cmd = "sequana_multitax --version" subprocess.call(cmd.split()) + + +def test_paired_reads_single_db(tmp_path): + runner = CliRunner() + results = runner.invoke( + main, ["--input-directory", paireddir, "--working-directory", str(tmp_path), "--force", "--databases", krakendb] + ) + assert results.exit_code == 0 + + +def test_paired_reads_multiple_dbs(tmp_path): + runner = CliRunner() + results = runner.invoke( + main, + [ + "--input-directory", + paireddir, + "--working-directory", + str(tmp_path), + "--force", + "--databases", + krakendb, + krakendb, + ], + ) + assert results.exit_code == 0 + + +def test_store_unclassified(tmp_path): + runner = CliRunner() + results = runner.invoke( + main, + [ + "--input-directory", + sharedir, + "--working-directory", + str(tmp_path), + "--force", + "--databases", + krakendb, + "--store-unclassified", + ], + ) + assert results.exit_code == 0 + config = yaml.safe_load((Path(str(tmp_path)) / "config.yaml").read_text()) + assert config["sequana_taxonomy"]["store_unclassified"] is True + + +def test_kraken_confidence(tmp_path): + runner = CliRunner() + results = runner.invoke( + main, + [ + "--input-directory", + sharedir, + "--working-directory", + str(tmp_path), + "--force", + "--databases", + krakendb, + "--kraken-confidence", + "0.5", + ], + ) + assert results.exit_code == 0 + config = yaml.safe_load((Path(str(tmp_path)) / "config.yaml").read_text()) + assert config["sequana_taxonomy"]["confidence"] == 0.5 + + +def test_blast_unclassified_config(tmp_path): + runner = CliRunner() + results = runner.invoke( + main, + [ + "--input-directory", + sharedir, + "--working-directory", + str(tmp_path), + "--force", + "--databases", + krakendb, + "--do-blast-unclassified", + ], + ) + assert results.exit_code == 0 + config = yaml.safe_load((Path(str(tmp_path)) / "config.yaml").read_text()) + assert config["blast"]["do"] is True + assert config["sequana_taxonomy"]["store_unclassified"] is True + + +@pytest.mark.skipif( + shutil.which("kraken2") is None or shutil.which("sequana_taxonomy") is None, + reason="kraken2 / sequana_taxonomy not available; functional run skipped", +) +def test_functional(tmp_path): + # End-to-end run on valid (non-empty) data. The deliberately-empty + # `empty_R1_.fastq.gz` fixture is excluded as it has no reads to classify. + runner = CliRunner() + results = runner.invoke( + main, + [ + "--input-directory", + sharedir, + "--exclude-pattern", + "empty", + "--working-directory", + str(tmp_path), + "--force", + "--databases", + krakendb, + ], + ) + assert results.exit_code == 0 + + completed = subprocess.run(["bash", "multitax.sh"], cwd=str(tmp_path)) + assert completed.returncode == 0 + + assert (tmp_path / "data" / "summary.html").exists() + assert (tmp_path / "data" / "kraken" / "kraken.csv").exists()