diff --git a/.github/workflows/build_and_test.yml b/.github/workflows/build_and_test.yml index 7dcd868..a888cec 100644 --- a/.github/workflows/build_and_test.yml +++ b/.github/workflows/build_and_test.yml @@ -38,21 +38,6 @@ jobs: run: | ln -s ~/kraken2_db_8 kraken2_db - - name: make sylph db # Needed so tests can run quickly - run: | - cd ~ - if [ ! -d "gatekeeper_sylph_test_data" ]; then - mkdir gatekeeper_sylph_test_data - cp /mnt/portal-sp3dev-knowledge/sylph/gtdb-r220-c1000-dbv1.syldb gatekeeper_sylph_test_data - cp /mnt/portal-sp3dev-knowledge/sylph/virus_imgvr_c200_v0.3.0.syldb gatekeeper_sylph_test_data - cp /mnt/portal-sp3dev-knowledge/sylph/fungi-refseq-2024-07-25-c200-v0.3.syldb gatekeeper_sylph_test_data - cp /mnt/portal-sp3dev-knowledge/sylph/*.tsv* gatekeeper_sylph_test_data - fi - - - name: symlink sylph db - run: | - ln -s ~/gatekeeper_sylph_test_data sylph_dbs - # Run pytest in container to use conda environment - name: pytest in docker run: | diff --git a/.gitignore b/.gitignore index 4143746..7483ef2 100644 --- a/.gitignore +++ b/.gitignore @@ -22,8 +22,6 @@ kraken2_db* results report*.html -sylph_dbs - # Python *__pycache__ diff --git a/CHANGELOG.md b/CHANGELOG.md index db8c17e..0ac5c0d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,8 +1,6 @@ ## New -- Code for generating sylph databases has been altered to incorporate ani_groups (based on taxid). This is only set for groups, not singletons -- when making sylph databases the taxonomy tsv now have a header and gets the ani_group column -- code for making virus sylph db now filters for human viruses +- Remove Sylph ## 2.4.0 diff --git a/README.md b/README.md index d5cb55d..ef78ade 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,6 @@ This pipeline should be used to perform check quality and trimming in sequence d * Nextflow (tested on version 23.04) * dependencies in [env.yml](env.yml) * Kraken2 [Standard dataset](https://benlangmead.github.io/aws-indexes/k2) (tested with Standard-8 from 2023/06/05). -* For sylph: pre-built indexes and associated metadata. ## Running Pipeline ## Parameters: @@ -14,11 +13,9 @@ Parameters: - **input_dir** = A directory containing input fastq.gz files - **seq_platform** = Sequencing platform `illumina` or `ont` - **kraken2_db_path** = Path to the kraken2 dataset -- **workflow** = `sylph`, `flu`, or `myco`. `speciated` is deprecated. +- **workflow** = `flu`, or `myco`. `speciated` is deprecated. - **mem_map** (Optional) = Use Kraken 2 memory mapping (requires preloading kraken2 db `cp -r /dev/shm`) `true` or `false` - **publish_dir** (Optional) If set will output files to directory provided. -- **sylph_dbs** = Path to sylph database(s). See [here](https://github.com/bluenote-1577/sylph/wiki/Pre%E2%80%90built-databases) -- **sylph_db_taxonomy_files** = Path taxonomy file(s) for sylph db. Should be tsv (optionally gzipped) See parameters with: ```bash @@ -56,56 +53,6 @@ cp -r kraken_index /dev/shm/ kraken2 --db /dev/shm/kraken_index --memory-mapping ... ``` -## Sylph - -Sylph is a k-mer based classification system that processes all the reads in one go, rather than classifying per read. -It can accept multiple databases, so to incorporate fungi or viruses an extra database can be provided rather than needing to rebuild a single db. -To use multiple databases it is easiest to set it in the `nextflow.config`: -```bash -params { - sylph_dbs = ["knowledge_bucket/sylph/refseq_bacteria_2025-08-04/refseq_bacteria_db_c200.syldb", - "knowledge_bucket/sylph/refseq_fungi_2025-08-04/refseq_fungi_db_c200.syldb", - "knowledge_bucket/sylph/refseq_viral_2025-08-04/refseq_viral_db_c100.syldb"] - sylph_db_taxonomy_files = ["knowledge_bucket/sylph/refseq_bacteria_2025-08-04/taxonomy.tsv.gz", - "knowledge_bucket/sylph/refseq_fungi_2025-08-04/taxonomy.tsv.gz", - "knowledge_bucket/sylph/refseq_viral_2025-08-04/taxonomy.tsv.gz"] -} -``` - -The taxonomy file needs to fit the format required for [sylph-tax](https://sylph-docs.github.io/sylph-tax-custom-taxonomies/). -As such it must be a two column TSV where the first column if the name of the fasta file, and the second column is the taxonomy string. -From testing: -- Can add a header without causing issues -- Can add additional columns - -The taxonomy format for sylph-tax is the same as provided in the downloads from [GTDB](https://data.gtdb.aau.ecogenomic.org/releases/release226/226.0/) - -### Parameters -There are two parameters used in the Sylph process: -- `sylph_query_ani_threshold`: sets the minimum ani threshold from the query step. Sylph authors recommend 90. -- `sylph_profile_ani_threshold`: sets the minimum ani threshold from the profile step. Sylph authors recommend 95. - -### Databases for use -Preference is to go with NCBI over GTDB as it is more familiar to clinical users (Though GTDB likely works better). -Likewise with viruses the focus should be on human pathogenic viruses rather than phages etc. - -In the knowledge bucket there are a variety of sylph databates. -To use one in this pipeline requires the database itself plus a tsv with the taxonomy metadata. - -See [sylph_references](sylph_references/README.md) for details of creating databases from RefSeq (NCBI). - -Note: For database to be used in competitive mapping downstream it also requires a folder of the reference genomes to make a manifest from. -As such the pre-built Sylph dbs are less useful except for GTDB which matches certain releases. - -### Data for tests -Sylph tests expect a folder `sylph_dbs` containing the sylph [pre-built databases](https://github.com/bluenote-1577/sylph/wiki/Pre%E2%80%90built-databases) for testing. -These are included in the knowledge bucket in the sylph directory. Specifically uses: -- `gtdb-r220-c1000-dbv1.syldb` and `gtdb_r220_metadata.tsv.gz` -- `virus_imgvr_c200_v0.3.0.syldb` and `virus_imgvr_4.1_metadata.tsv.gz` - -The c1000 means that the k-mers are sampled every 1000bp making it smaller than c200, and so useful for testing. -The metadata/taxonomy file can be optionally gzipped. - ## Development/Install The python side of this repo forms a python project. diff --git a/env.yml b/env.yml index e5c9c77..bc08a2a 100644 --- a/env.yml +++ b/env.yml @@ -11,5 +11,3 @@ dependencies: - just=1.40.0 - bracken=2.9 - pytest=8.3.2 -- sylph=0.8.0 -- sylph-tax=1.1.2 diff --git a/main.nf b/main.nf index 8ff1eb2..4db0431 100644 --- a/main.nf +++ b/main.nf @@ -4,7 +4,6 @@ include { fastp } from './process/quality_check.nf' include { kraken2 ; summarise_taxa ; bracken_summarise_taxa } from './process/kraken2.nf' include { filter_myco ; filter_flu } from './process/filter_reads.nf' include { has_enough_reads } from './process/has_enough_reads.nf' -include { sylph } from './process/sylph.nf' // Parameters params.input_paired_suffix = "*_{1,2}.fastq.gz" @@ -15,8 +14,6 @@ params.gatekeeper_illumina_min_num_reads = 10000 params.gatekeeper_ont_min_num_reads = 1000 params.gatekeeper_confidence_threshold = 0.00 params.gatekeeper_taxa_summary_threshold = 0.5 -params.sylph_query_ani_threshold = 90 -params.sylph_profile_ani_threshold = 95 workflow { if (params.help) { @@ -27,12 +24,10 @@ workflow { Parameters: ------------------------------------------------------------------------ - --workflow workflow to run: speciate, myco, sylph, speciate_detailed + --workflow workflow to run: speciate, myco, speciate_detailed --input_dir Path to directory containing input fastq files --seq_platform sequencing platform: illumina or ont --kraken2_db_path Path to the kraken dataset directory - --sylph_db Path to the sylph database - --sylph_db_taxonomy Path to the sylph database taxonomy metadata Optional: --mem_map If 'true' use memory mapping for kraken2. See README for details @@ -43,12 +38,12 @@ workflow { } if (params.seq_platform == 'ont') { - input_files = Channel.fromPath("${params.input_dir}/${params.input_single_suffix}", checkIfExists: true) + input_files = channel.fromPath("${params.input_dir}/${params.input_single_suffix}", checkIfExists: true) .ifEmpty { error("cannot find any reads matching ${params.input_single_suffix} in ${params.input_dir}") } .map { it -> tuple(it.simpleName, it) } } else if (params.seq_platform == 'illumina') { - input_files = Channel.fromFilePairs( + input_files = channel.fromFilePairs( "${params.input_dir}/${params.input_paired_suffix}", flat: false, checkIfExists: true, @@ -57,12 +52,12 @@ workflow { .ifEmpty { error("cannot find any reads matching ${params.input_paired_suffix} in ${params.input_dir}") } } - kraken2_db_path = Channel.fromPath(params.kraken2_db_path).first() + kraken2_db_path = channel.fromPath(params.kraken2_db_path).first() check_seq_platform(params.seq_platform, input_files) if (params.workflow == 'speciate') { - log.info("Speciate workflow is deprecated, consider sylph instead") + log.info("Speciate workflow is deprecated, consider sylph instead within competitive mapping") gatekeeper_speciate(input_files, kraken2_db_path, params.seq_platform) } else if (params.workflow == 'myco') { @@ -71,9 +66,6 @@ workflow { else if (params.workflow == 'flu') { gatekeeper_flu(input_files, kraken2_db_path, params.seq_platform) } - else if (params.workflow == 'sylph') { - gatekeeper_sylph(input_files, params.sylph_dbs, params.sylph_db_taxonomy_files, params.seq_platform) - } else { error("Invalid workflow: ${params.workflow}") } @@ -169,28 +161,6 @@ workflow gatekeeper_flu { gatekeeper_report = filter_flu.out.gatekeeper_report_json } -workflow gatekeeper_sylph { - take: - fastq_files - sylph_dbs - sylph_db_taxonomy_files - seq_platform - - main: - sylph( - fastq_files, - sylph_dbs, - sylph_db_taxonomy_files, - seq_platform, - params.sylph_query_ani_threshold, - params.sylph_profile_ani_threshold, - ) - - emit: - sylph_query = sylph.out.sylph_query - sylph_report = sylph.out.sylph_report - taxonomy_report = sylph.out.taxonomy_report -} def check_seq_platform(seq_platform, fastq_files) { def seq_platforms = ['ont', 'illumina'] diff --git a/minimal_files.txt b/minimal_files.txt index eb7e7af..b50571c 100644 --- a/minimal_files.txt +++ b/minimal_files.txt @@ -7,5 +7,3 @@ process/templates/filter_myco.sh process/templates/filter_flu.sh process/templates/run_fastp.sh process/templates/run_kraken.sh -process/sylph.nf -process/templates/run_sylph.sh diff --git a/process/sylph.nf b/process/sylph.nf deleted file mode 100644 index 76de1b4..0000000 --- a/process/sylph.nf +++ /dev/null @@ -1,36 +0,0 @@ -process sylph { - publishDir "${params.publish_dir}", enabled: params.publish_dir != "", mode: "copy", saveAs: { filename -> sample_name + "_" + filename } - maxRetries 5 - - container { - params.test_container_gatekeeper == "" ? params.container_prefix + '/gpas/gatekeeper_pipeline:37bfc21' : params.test_container_gatekeeper - } - cpus 4 - memory { 14.GB * task.attempt } - - pod label: "name", value: "gatekeeper_pipeline:sylph" - pod label: "sample_id", value: "${params.sample_id}" - pod label: "run_id", value: "${params.run_id}" - - input: - tuple val(sample_name), path(fqs) - // Use fixed name to avoid name conflicts. - // Can use multiple databases - path "database??.syldb" - path "taxonomies?/*" - val seq_platform - val query_ani_threshold - val profile_ani_threshold - - output: - tuple val(sample_name), path(sylph_query), emit: sylph_query - tuple val(sample_name), path(sylph_report), emit: sylph_report - tuple val(sample_name), path(taxonomy_report), emit: taxonomy_report - - script: - sylph_query = "sylph_query.tsv" - sylph_report = "sylph_report.tsv" - taxonomy_report = "taxonomy_report.tsv" - - template('run_sylph.sh') -} diff --git a/process/templates/run_sylph.sh b/process/templates/run_sylph.sh deleted file mode 100644 index 965a50f..0000000 --- a/process/templates/run_sylph.sh +++ /dev/null @@ -1,34 +0,0 @@ - -if [ $seq_platform == 'ont' ] -then - sylph sketch -t ${task.cpus} -c 100 -r $fqs -elif [ $seq_platform == 'illumina' ] -then - sylph sketch -t ${task.cpus} -c 100 \ - -1 ${fqs[0]} -2 ${fqs[1]} -fi -# either way a file ending .sylsp is created - -sylph query -t ${task.cpus} -m ${query_ani_threshold} \ - -o ${sylph_query} database*.syldb *.sylsp - -sylph profile -u -t ${task.cpus} -m ${profile_ani_threshold} \ - -o ${sylph_report} database*.syldb *.sylsp - - -if [[ \$(wc -l <${sylph_report}) -le 1 ]] -then - echo "Sylph failed to generate a profile report" - echo "Generating a dummy taxonomy report" - if [ $seq_platform == 'ont' ] - then - echo -e "#SampleID\t${fqs}" > ${taxonomy_report} - elif [ $seq_platform == 'illumina' ] - then - echo -e "#SampleID\t${fqs[0]}" > ${taxonomy_report} - fi - echo -e "clade_name\trelative_abundance\tsequence_abundance\tANI (if strain-level)\tCoverage (if strain-level)" >> ${taxonomy_report} -else - sylph-tax taxprof ${sylph_report} -t taxonomies*/* -o taxo_ - cat taxo_* > ${taxonomy_report} -fi diff --git a/pyproject.toml b/pyproject.toml index 62135af..66c4a28 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -55,5 +55,4 @@ version_files = [ "process/kraken2.nf", "process/quality_check.nf", "process/filter_reads.nf", - "process/sylph.nf", ] diff --git a/sylph_references/.gitignore b/sylph_references/.gitignore deleted file mode 100644 index 6744e18..0000000 --- a/sylph_references/.gitignore +++ /dev/null @@ -1,14 +0,0 @@ -assembly_summary_refseq.txt -ncbi_taxonomy/ -virushostdb.daily.tsv - -*download/ -*download.sh -*metadata.csv -*rsync_file_list.txt -*taxonomy.tsv -*taxonomy.tsv.gz -*ani_edge_list.tsv -*.syldb -*genome_paths.tsv -genomes/ diff --git a/sylph_references/README.md b/sylph_references/README.md deleted file mode 100644 index 4428e76..0000000 --- a/sylph_references/README.md +++ /dev/null @@ -1,120 +0,0 @@ -# Making Reference databases for Sylph - -## Currently available databases -These are available in the sp3dev knowledge bucket - -Sylph pre built (options)[https://github.com/bluenote-1577/sylph/wiki/Pre%E2%80%90built-databases]: -- GTDB-220 release from 24/04/24 using c200 or c1000 `gtdb-r220-c200-dbv1.syldb` -- RefSeq fungi from 25/07/24 using c200 `fungi-refseq-2024-07-25-c200-v0.3.syldb` -- IMG/VR4.1 high-confidence vOTU reps using c200 `virus_imgvr_c200_v0.3.0.syldb` - -RefSeq databases built on 4/8/2025 using method below. - -TODO: Look into IMG/VR4.1 for making custom viral database - ---- - -## Dependencies: -```bash -conda env create -f env.yml -conda activate ncbi_refseq -``` - -## Components of a database -- a folder called `genomes` containing all the ref fastas -- A file called `genome_paths.tsv` with paths to the genomes based on their filenames -- a `.syldb` which is used by sylph and built from the genomes -- `taxonomy.tsv` is a table going from ref file name to taxonomy -- `metadata.csv` has the metadata for chosen genomes. The parts used by comp mapping is `accession`, `taxonomy`, `ani_group`. So the name used for the reference will be the species in the taxonomy column. - -## Walk through making databases -This walks through the steps to produce sylph databases from NCBI. -Note that it takes some time to download all the references! - -To make reference we need: -- The list of reference genomes [assembly_summary_refseq.tsv](https://ftp.ncbi.nlm.nih.gov/genomes/ASSEMBLY_REPORTS/). -- The taxonomy [taxdump.tar.gz](https://ftp.ncbi.nih.gov/pub/taxonomy/). This should be extracted to a folder. -Both of these are updated regularly, so note when they were collected and aim to be at the same time. -- Viral-host database file [virushostdb.daily.tsv](https://www.genome.jp/ftp/db/virushostdb/). Only human viruses will be used. - -```bash -wget https://ftp.ncbi.nlm.nih.gov/genomes/ASSEMBLY_REPORTS/assembly_summary_refseq.txt -# equivalently -wget https://ftp.ncbi.nlm.nih.gov/genomes/refseq/assembly_summary_refseq.txt - -# For explanations -wget https://ftp.ncbi.nlm.nih.gov/genomes/README_assembly_summary.txt - -mkdir ncbi_taxonomy -wget https://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz -tar -xvf taxdump.tar.gz -C ncbi_taxonomy -rm taxdump.tar.gz - -wget https://www.genome.jp/ftp/db/virushostdb/virushostdb.daily.tsv -``` - -### Selecting reference genomes - -The reference genomes can then be selected: -```bash -python3 select_references.py assembly_summary_refseq.txt ncbi_taxonomy virushostdb.daily.tsv -``` - -For each group (bacteria/archaea, fungi, viral) it will produce a folder containing: -- csv with metadata on all the references -- a file to download all genomes -- a taxonomy file used by sylph to produce taxonomy summary. - -Note: some references do not seem to have a taxonomy, this seems to be for very recently published refs. -You may need to manually set these. Look online for that reference and fill in the taxonomy value accordingly. Or try again tomorrow. - -### Downloading reference genomes - -You can then download the samples (likely with screen) using: -```bash -cd refseq_bacteria -mkdir -p genomes -rsync -avP --files-from=rsync_file_list.txt rsync://ftp.ncbi.nlm.nih.gov/genomes ./genomes -``` - -This runs a single rsync command for all the files. You should just be able to run it again if connection drops (etc). - -At this point you should have a directory with all the reference genomes - -### Making resources for species identification pipeline -To make the sylph database. Can lower c to make more sensitive (suggested for viruses is c=100) -```bash -find genomes -type f > sylph_file_list.txt -sylph sketch -l sylph_file_list.txt -o refseq_bacteria_db -t 30 -c 200 -rm sylph_file_list.txt -``` -to produce refseq_bacteria_db.syldb - -The `taxonomy.tsv` acts as the metadata for sylph. Can be gzipped to be smaller: -```bash -gzip -c taxonomy.tsv > taxonomy.tsv.gz -``` - -Will also require a `genome_paths.tsv` which is used by dynamic competitive mapping: -```bash -find genomes -type f | while read file_path; do - file_name=$(basename $file_path) - dir=$(dirname $file_path) - echo -e "$file_name\t$dir/" >> genome_paths.tsv -done -``` - - -### Making ANI groups - -Assign genomes to ANI groups. This adds "ani_group" and "in_group" to the metadata csv. -```bash -cd refseq_bacteria -python3 ../assign_ani_groups.py -g genomes -m metadata.csv -t taxonomy.tsv --cpus 16 --ani_threshold 97 -``` -This also produces `ani_edge_list.tsv` which has the actual distances. - -The ani_groups are based on the lowest taxid within each group. -If taxid's are not present than just uses row index. -An ani group is only given if there is more than one genome in the group. -Both metadata.csv and taxonomy.tsv get the ani_group column added. diff --git a/sylph_references/assign_ani_groups.py b/sylph_references/assign_ani_groups.py deleted file mode 100644 index 39fbb9b..0000000 --- a/sylph_references/assign_ani_groups.py +++ /dev/null @@ -1,129 +0,0 @@ -import argparse -import os -import subprocess -from pathlib import Path - -import networkx as nx -import pandas as pd - - -def set_ani_distances( - genome_dir: str, - metadata: pd.DataFrame, - ani_threshold: float, - cpus: int, -) -> pd.DataFrame: - """Get ANI distances between genomes using skani""" - with open("tmp_file_list.txt", "w", encoding="utf-8") as f: - for file in Path(genome_dir).rglob("*.fna.gz"): - f.write(str(file) + "\n") - - command = ( - f"skani triangle -E -t {cpus} --medium -l tmp_file_list.txt > ani_edge_list.tsv" - ) - subprocess.run(command, shell=True, check=True) - ani_df = pd.read_csv("ani_edge_list.tsv", sep="\t") - os.remove("tmp_file_list.txt") - - # filter ANI dataframe based on threshold - ani_df = ani_df[ani_df["ANI"] >= ani_threshold] - - # get assembly accessions from file paths - ani_df["Ref"] = ( - ani_df["Ref_file"].str.replace("_genomic.fna.gz", "").str.split("/").str[-1] - ) - ani_df["Query"] = ( - ani_df["Query_file"].str.replace("_genomic.fna.gz", "").str.split("/").str[-1] - ) - - # Ensure that both Ref and Query are in metadata - ani_df = ani_df[ - ani_df["Ref"].isin(metadata["accession"]) - & ani_df["Query"].isin(metadata["accession"]) - ] - - # build graph from filtered data - graph = nx.from_pandas_edgelist(ani_df, "Ref", "Query") - - # Set initial ani_group as taxid (or row index if no taxid) - if "taxid" in metadata.columns: - metadata["ani_group"] = metadata["taxid"] - else: - metadata["ani_group"] = metadata.index - metadata["in_group"] = False - - # Assign ani_group based on connected components - for component in nx.connected_components(graph): - if len(component) <= 1: - continue - ani_group = min( - metadata.loc[metadata["accession"].isin(component), "ani_group"] - ) - metadata.loc[metadata["accession"].isin(component), "ani_group"] = ani_group - metadata.loc[metadata["accession"].isin(component), "in_group"] = True - - # If not in a group set ani_group to nan - metadata["ani_group"] = metadata["ani_group"].astype("Int64") - metadata.loc[~metadata["in_group"], "ani_group"] = pd.NA - - front_cols = ["accession", "ani_group", "in_group"] - metadata = metadata[ - front_cols + [col for col in metadata.columns if col not in front_cols] - ] - - # sort by ani_group and accession - metadata.sort_values(by=["ani_group", "accession"], inplace=True) - - return metadata - - -def main(): - parser = argparse.ArgumentParser( - description="Run skani triangle to create ani groups." - ) - parser.add_argument( - "-g", - "--genome_dir", - required=True, - help="Directory containing gzipped fasta files", - ) - parser.add_argument( - "-m", "--metadata", required=True, help="CSV file with metadata" - ) - parser.add_argument( - "-t", "--taxonomy", required=True, help="TSV file with taxonomy" - ) - parser.add_argument( - "-c", "--cpus", type=int, default=8, help="Number of CPUs to use for skani" - ) - parser.add_argument( - "-a", - "--ani_threshold", - type=float, - default=95, - help="ANI threshold for filtering (default: 95)", - ) - args = parser.parse_args() - - metadata = pd.read_csv(args.metadata) - - metadata = set_ani_distances( - genome_dir=args.genome_dir, - metadata=metadata, - ani_threshold=args.ani_threshold, - cpus=args.cpus, - ) - - metadata.to_csv(args.metadata, index=False) - - taxonomy = pd.read_csv(args.taxonomy, sep="\t", usecols=["accession", "taxonomy"]) - taxonomy = taxonomy.merge( - metadata[["accession", "ani_group"]], on="accession", how="left" - ) - taxonomy.to_csv(args.taxonomy, sep="\t", index=False) - - print("ANI groups assigned and metadata updated.") - - -if __name__ == "__main__": - main() diff --git a/sylph_references/check_for_missing.py b/sylph_references/check_for_missing.py deleted file mode 100644 index f43955b..0000000 --- a/sylph_references/check_for_missing.py +++ /dev/null @@ -1,36 +0,0 @@ -import argparse -import os - -import pandas as pd - - -def main(): - parser = argparse.ArgumentParser(description="Check for missing genome files.") - parser.add_argument("input_dir", help="Directory containing .fna.gz files") - parser.add_argument("input_csv", help="CSV file with accession column") - args = parser.parse_args() - - df = pd.read_csv(args.input_csv) - if "assembly_accession" not in df.columns: - print("Error: 'assembly_accession' column not found in CSV.") - return - - missing = [] - for accession in df["assembly_accession"]: - filename = f"{accession}_genomic.fna.gz" - filepath = os.path.join(args.input_dir, filename) - if not os.path.isfile(filepath): - missing.append(accession) - - if not missing: - print("All genome files are present.") - - missing_df = df[df["assembly_accession"].isin(missing)] - missing_df.to_csv("missing_genomes.csv", index=False) - print( - f"Missing genome files saved to 'missing_genomes.csv'. Found {len(missing)} missing files." - ) - - -if __name__ == "__main__": - main() diff --git a/sylph_references/env.yml b/sylph_references/env.yml deleted file mode 100644 index dc5fb75..0000000 --- a/sylph_references/env.yml +++ /dev/null @@ -1,14 +0,0 @@ -name: ncbi_refseq -channels: - - conda-forge - - bioconda - - defaults -dependencies: - - "python>=3.11,<4.0" - - pandas>=2.2.0 - - sylph=0.8.0 - - skani - - pip - - pip: - - git+https://github.com/GlobalPathogenAnalysisService/kraken2_analysis.git - - networkx[default] diff --git a/sylph_references/make_dbs.sh b/sylph_references/make_dbs.sh deleted file mode 100644 index c739235..0000000 --- a/sylph_references/make_dbs.sh +++ /dev/null @@ -1,30 +0,0 @@ -# Script for doing the steps from downloading refs onwards - -threads=100 - -for group in fungi bacteria viral -do - cd refseq_${group} - mkdir -p genomes - rsync -avP --files-from=rsync_file_list.txt rsync://ftp.ncbi.nlm.nih.gov/genomes ./genomes - - c=200 - if [ "$group" == "viral" ]; then - c=100 - fi - find genomes -type f > sylph_file_list.txt - sylph sketch -l sylph_file_list.txt -o refseq_${group}_db -t ${threads} -c ${c} - rm sylph_file_list.txt - - gzip -c taxonomy.tsv > taxonomy.tsv.gz - - find genomes -type f | while read file_path; do - file_name=$(basename $file_path) - dir=$(dirname $file_path) - echo -e "$file_name\t$dir/" >> genome_paths.tsv - done - - python3 ../assign_ani_groups.py -g genomes/ -m metadata.csv --cpus ${threads} -a 97 - - cd .. -done diff --git a/sylph_references/select_references.py b/sylph_references/select_references.py deleted file mode 100644 index 7dcea86..0000000 --- a/sylph_references/select_references.py +++ /dev/null @@ -1,209 +0,0 @@ -import argparse -import os - -import pandas as pd -from k2_analysis.kraken_tree import Tree, ncbi_taxonomy_to_tree - -PROVIDED_COLS = [ - "assembly_accession", - "bioproject", - "biosample", - "wgs_master", - "refseq_category", - "taxid", - "species_taxid", - "organism_name", - "infraspecific_name", - "isolate", - "version_status", - "assembly_level", - "release_type", - "genome_rep", - "seq_rel_date", - "asm_name", - "asm_submitter", - "gbrs_paired_asm", - "paired_asm_comp", - "ftp_path", - "excluded_from_refseq", - "relation_to_type_material", - "asm_not_live_date", - "assembly_type", - "group", - "genome_size", - "genome_size_ungapped", - "gc_percent", - "replicon_count", - "scaffold_count", - "contig_count", - "annotation_provider", - "annotation_name", - "annotation_date", - "total_gene_count", - "protein_coding_gene_count", - "non_coding_gene_count", - "pubmed_id", -] - -SELECTED_COLS = [ - "assembly_accession", - "bioproject", - "biosample", - "refseq_category", - "taxid", - "species_taxid", - "organism_name", - "isolate", - "version_status", - "assembly_level", - "genome_rep", - "relation_to_type_material", - "group", - "genome_size", - "contig_count", - "annotation_provider", - "ftp_path", -] - -LEVEL_CODES = { - "superkingdom": "d", - "phylum": "p", - "class": "c", - "order": "o", - "family": "f", - "genus": "g", - "species": "s", -} - - -def get_taxonomy(node: Tree | None) -> list[str]: - if node is None: - return ["unknown_taxa"] - - if node.level_code in ["R"] or node.name == "cellular organisms": - return [] - - level_code = LEVEL_CODES.get(node.level_code, "x") - - taxa_string = level_code + "__" + str(node.name) - return get_taxonomy(node.parent) + [taxa_string] - - -def main(): - # parse arguments - parser = argparse.ArgumentParser( - description="Select reference genomes from refseq metadata" - ) - parser.add_argument("metadata", type=str, help="Path to the metadata file") - parser.add_argument( - "ncbi_taxonomy_dir", - type=str, - help="Path to directory with names.dmp and nodes.dmp", - ) - parser.add_argument( - "virushostdb", - type=str, - help="Path to the viral-host database file (virushostdb.daily.tsv)", - ) - args = parser.parse_args() - - df = pd.read_csv( - args.metadata, sep="\t", skiprows=2, names=PROVIDED_COLS, usecols=SELECTED_COLS - ) - - print("groups: ", df["group"].unique()) - - df = df[df["group"].isin(["bacteria", "archaea", "fungi", "viral"])] - df = df[df["version_status"] == "latest"] - refs = df[ - (df["refseq_category"] == "reference genome") | (df["group"] == "viral") - ].copy() - - missing_ftp = refs[refs["ftp_path"].isna()] - if not missing_ftp.empty: - print( - f"Warning: {len(missing_ftp)} samples have no ftp_path and will be excluded from the download list" - ) - missing_ftp.to_csv("missing_ftp.csv", index=False) - print("Missing ftp paths saved to missing_ftp.csv") - - refs = refs[~refs["ftp_path"].isna()] - - def get_rsync_path(ftp_path): - accession = ftp_path.split("/")[-1] - return f"{ftp_path.replace('https://ftp.ncbi.nlm.nih.gov/genomes/', '')}/{accession}_genomic.fna.gz" - - refs["rsync_path"] = refs.apply( - lambda row: get_rsync_path(row["ftp_path"]), - axis=1, - ) - - refs["accession"] = refs["ftp_path"].str.split("/").str[-1] - refs = refs[["accession"] + [c for c in refs.columns if c not in ["accession"]]] - - taxonomy = ncbi_taxonomy_to_tree( - f"{args.ncbi_taxonomy_dir}/nodes.dmp", f"{args.ncbi_taxonomy_dir}/names.dmp" - ) - taxa_to_node = {taxon.taxid: taxon for taxon in taxonomy.to_list()} - print("Taxonomy loaded") - refs["taxonomy"] = refs["taxid"].apply( - lambda x: ";".join(get_taxonomy(taxa_to_node.get(x, None))) - ) - - unknown_taxa = refs[refs["taxonomy"] == "unknown_taxa"] - if not unknown_taxa.empty: - print( - f"Warning: {len(unknown_taxa)} samples have no taxonomy please set manually." - ) - unknown_taxa.to_csv("unknown_taxa.csv", index=False) - print("Unknown taxa saved to unknown_taxa.csv") - - def save_outputs(df_to_save: pd.DataFrame, prefix: str): - os.makedirs(f"refseq_{prefix}", exist_ok=True) - df_to_save.drop(columns=["rsync_path", "ftp_path"], inplace=False).to_csv( - f"refseq_{prefix}/metadata.csv", index=False - ) - - df_to_save[["rsync_path"]].to_csv( - f"refseq_{prefix}/rsync_file_list.txt", index=False, header=False - ) - - df_to_save[["accession", "taxonomy"]].to_csv( - f"refseq_{prefix}/taxonomy.tsv", sep="\t", index=False, header=True - ) - - bacteria = refs[refs["group"].isin(["bacteria", "archaea"])] - save_outputs(bacteria, "bacteria") - - fungi = refs[refs["group"].isin(["fungi"])] - save_outputs(fungi, "fungi") - - ## Viruses - - viralhostdb = pd.read_csv( - args.virushostdb, - sep="\t", - usecols=["virus tax id", "virus name", "refseq id", "DISEASE", "host tax id"], - ) - viralhostdb = viralhostdb[ - viralhostdb["host tax id"] == 9606 - ] # Filter for human hosts - viralhostdb.rename( - columns={ - "virus tax id": "taxid", - "virus name": "virus_name", - "refseq id": "refseq_ids", - "DISEASE": "disease", - }, - inplace=True, - ) - viralhostdb.drop(columns=["host tax id"], inplace=True) - - # viruses are do not have refseq category - viral = refs[refs["group"].isin(["viral"])] - viral = viral.merge(viralhostdb, on="taxid", how="inner") - save_outputs(viral, "viral") - - -if __name__ == "__main__": - main() diff --git a/tests/nextflow/samples/clinic/sylph_report.tsv b/tests/nextflow/samples/clinic/sylph_report.tsv deleted file mode 100644 index 09d632d..0000000 --- a/tests/nextflow/samples/clinic/sylph_report.tsv +++ /dev/null @@ -1,4 +0,0 @@ -Sample_file Genome_file Taxonomic_abundance Sequence_abundance Adjusted_ANI True_cov ANI_5-95_percentile Eff_lambda Lambda_5-95_percentile Median_cov Mean_cov_geq1 Containment_ind Naive_ANI kmers_reassigned Contig_name -test_clinic_sample_1.fastq.gz gtdb_genomes_reps_r220/database/GCF/002/086/515/GCF_002086515.1_genomic.fna.gz 97.6912 97.1919 99.99 47.385 NA-NA HIGH NA-NA 36 36.637 4279/4287 99.99 0 NZ_MVHH01000001.1 Mycolicibacter arupensis strain DSM 44942 Mycobacterium_arupense_DSM44942_contig1, whole genome shotgun sequence -test_clinic_sample_1.fastq.gz gtdb_genomes_reps_r220/database/GCF/025/215/495/GCF_025215495.1_genomic.fna.gz 1.6754 1.4966 99.68 0.813 99.39-99.92 0.617 0.56-0.69 1 1.346 1600/3773 97.27 47 NZ_CP104302.1 Mycolicibacterium brumae strain ATCC 51384 chromosome, complete genome -test_clinic_sample_1.fastq.gz gtdb_genomes_reps_r220/database/GCF/000/195/955/GCF_000195955.2_genomic.fna.gz 0.6334 0.6258 99.55 0.307 98.94-100.06 0.233 0.20-0.28 1 1.129 785/4253 94.70 9 NC_000962.3 Mycobacterium tuberculosis H37Rv, complete genome diff --git a/tests/nextflow/samples/clinic/sylph_taxonomy_report.tsv b/tests/nextflow/samples/clinic/sylph_taxonomy_report.tsv deleted file mode 100644 index 64e4b23..0000000 --- a/tests/nextflow/samples/clinic/sylph_taxonomy_report.tsv +++ /dev/null @@ -1,14 +0,0 @@ -#SampleID test_clinic_sample_1.fastq.gz Taxonomies_used:['taxonomies1/gtdb_r220_metadata.tsv.gz'] -clade_name relative_abundance sequence_abundance ANI (if strain-level) Coverage (if strain-level) -d__Bacteria 99.99999999999999 99.3143 NA NA -d__Bacteria|p__Actinomycetota 99.99999999999999 99.3143 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes 99.99999999999999 99.3143 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales 99.99999999999999 99.3143 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae 99.99999999999999 99.3143 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae|g__Mycobacterium 99.99999999999999 99.3143 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae|g__Mycobacterium|s__Mycobacterium arupense 97.6912 97.1919 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae|g__Mycobacterium|s__Mycobacterium brumae 1.6754 1.4966 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae|g__Mycobacterium|s__Mycobacterium tuberculosis 0.6334 0.6258 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae|g__Mycobacterium|s__Mycobacterium arupense|t__GCF_002086515.1 97.6912 97.1919 99.99 47.385 -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae|g__Mycobacterium|s__Mycobacterium brumae|t__GCF_025215495.1 1.6754 1.4966 99.68 0.813 -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae|g__Mycobacterium|s__Mycobacterium tuberculosis|t__GCF_000195955.2 0.6334 0.6258 99.55 0.307 diff --git a/tests/nextflow/samples/clinic/with_virus_sylph_report.tsv b/tests/nextflow/samples/clinic/with_virus_sylph_report.tsv deleted file mode 100644 index 7b6939f..0000000 --- a/tests/nextflow/samples/clinic/with_virus_sylph_report.tsv +++ /dev/null @@ -1,7 +0,0 @@ -Sample_file Genome_file Taxonomic_abundance Sequence_abundance Adjusted_ANI True_cov ANI_5-95_percentile Eff_lambda Lambda_5-95_percentile Median_cov Mean_cov_geq1 Containment_ind Naive_ANI kmers_reassigned Contig_name -test_clinic_sample_1.fastq.gz imgvr_reps.fna 26.0540 0.4168 100.00 50.018 NA-NA HIGH NA-NA 38 37.136 81/81 100.00 0 IMGVR_UViG_3300048835_000078|3300048835|Ga0498919_0008194 -test_clinic_sample_1.fastq.gz imgvr_reps.fna 24.6828 1.0613 98.55 47.385 NA-NA HIGH NA-NA 36 36.156 141/222 98.55 81 IMGVR_UViG_2636415716_000002|2636415716|2636450064|114-49182 -test_clinic_sample_1.fastq.gz gtdb_genomes_reps_r220/database/GCF/002/086/515/GCF_002086515.1_genomic.fna.gz 24.6828 96.0894 99.95 47.385 NA-NA HIGH NA-NA 36 36.668 4223/4287 99.95 56 NZ_MVHH01000001.1 Mycolicibacter arupensis strain DSM 44942 Mycobacterium_arupense_DSM44942_contig1, whole genome shotgun sequence -test_clinic_sample_1.fastq.gz imgvr_reps.fna 23.9971 0.3342 100.00 46.069 NA-NA HIGH NA-NA 35 35.463 67/67 100.00 0 IMGVR_UViG_2636415716_000003|2636415716|2636450084|42096-57988 -test_clinic_sample_1.fastq.gz gtdb_genomes_reps_r220/database/GCF/025/215/495/GCF_025215495.1_genomic.fna.gz 0.4233 1.4796 99.68 0.813 99.39-99.92 0.617 0.56-0.69 1 1.346 1600/3773 97.27 47 NZ_CP104302.1 Mycolicibacterium brumae strain ATCC 51384 chromosome, complete genome -test_clinic_sample_1.fastq.gz gtdb_genomes_reps_r220/database/GCF/000/195/955/GCF_000195955.2_genomic.fna.gz 0.1600 0.6187 99.55 0.307 98.94-100.06 0.233 0.20-0.28 1 1.129 785/4253 94.70 9 NC_000962.3 Mycobacterium tuberculosis H37Rv, complete genome diff --git a/tests/nextflow/samples/clinic/with_virus_sylph_taxonomy_report.tsv b/tests/nextflow/samples/clinic/with_virus_sylph_taxonomy_report.tsv deleted file mode 100644 index 55a22b0..0000000 --- a/tests/nextflow/samples/clinic/with_virus_sylph_taxonomy_report.tsv +++ /dev/null @@ -1,25 +0,0 @@ -#SampleID test_clinic_sample_1.fastq.gz Taxonomies_used:['taxonomies1/gtdb_r220_metadata.tsv.gz', 'taxonomies2/virus_imgvr_4.1_metadata.tsv.gz'] -clade_name relative_abundance sequence_abundance ANI (if strain-level) Coverage (if strain-level) -r__Duplodnaviria 74.7339 1.8123 NA NA -d__Bacteria 25.2661 98.1877 NA NA -r__Duplodnaviria|k__Heunggongvirae 74.7339 1.8123 NA NA -d__Bacteria|p__Actinomycetota 25.2661 98.1877 NA NA -r__Duplodnaviria|k__Heunggongvirae|p__Uroviricota 74.7339 1.8123 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes 25.2661 98.1877 NA NA -r__Duplodnaviria|k__Heunggongvirae|p__Uroviricota|c__Caudoviricetes 74.7339 1.8123 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales 25.2661 98.1877 NA NA -r__Duplodnaviria|k__Heunggongvirae|p__Uroviricota|c__Caudoviricetes| 74.7339 1.8123 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae 25.2661 98.1877 NA NA -r__Duplodnaviria|k__Heunggongvirae|p__Uroviricota|c__Caudoviricetes|| 74.7339 1.8123 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae|g__Mycobacterium 25.2661 98.1877 NA NA -r__Duplodnaviria|k__Heunggongvirae|p__Uroviricota|c__Caudoviricetes||| 74.7339 1.8123 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae|g__Mycobacterium|s__Mycobacterium arupense 24.6828 96.0894 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae|g__Mycobacterium|s__Mycobacterium brumae 0.4233 1.4796 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae|g__Mycobacterium|s__Mycobacterium tuberculosis 0.16 0.6187 NA NA -r__Duplodnaviria|k__Heunggongvirae|p__Uroviricota|c__Caudoviricetes|||| 74.7339 1.8123 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae|g__Mycobacterium|s__Mycobacterium arupense|t__GCF_002086515.1 24.6828 96.0894 99.95 47.385 -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae|g__Mycobacterium|s__Mycobacterium brumae|t__GCF_025215495.1 0.4233 1.4796 99.68 0.813 -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae|g__Mycobacterium|s__Mycobacterium tuberculosis|t__GCF_000195955.2 0.16 0.6187 99.55 0.307 -r__Duplodnaviria|k__Heunggongvirae|p__Uroviricota|c__Caudoviricetes|||||t__IMGVR_UViG_3300048835_000078 26.054 0.4168 100.0 50.018 -r__Duplodnaviria|k__Heunggongvirae|p__Uroviricota|c__Caudoviricetes|||||t__IMGVR_UViG_2636415716_000002 24.6828 1.0613 98.55 47.385 -r__Duplodnaviria|k__Heunggongvirae|p__Uroviricota|c__Caudoviricetes|||||t__IMGVR_UViG_2636415716_000003 23.9971 0.3342 100.0 46.069 diff --git a/tests/nextflow/samples/ont_full/sylph_report.tsv b/tests/nextflow/samples/ont_full/sylph_report.tsv deleted file mode 100644 index be08132..0000000 --- a/tests/nextflow/samples/ont_full/sylph_report.tsv +++ /dev/null @@ -1,2 +0,0 @@ -Sample_file Genome_file Taxonomic_abundance Sequence_abundance Adjusted_ANI True_cov ANI_5-95_percentile Eff_lambda Lambda_5-95_percentile Median_cov Mean_cov_geq1 Containment_ind Naive_ANI kmers_reassigned Contig_name -reads.fastq.gz gtdb_genomes_reps_r220/database/GCF/000/195/955/GCF_000195955.2_genomic.fna.gz 100.0000 96.3226 99.96 15.070 NA-NA HIGH NA-NA 9 9.079 4196/4253 99.96 0 NC_000962.3 Mycobacterium tuberculosis H37Rv, complete genome diff --git a/tests/nextflow/samples/ont_full/sylph_taxonomy_report.tsv b/tests/nextflow/samples/ont_full/sylph_taxonomy_report.tsv deleted file mode 100644 index 9da3f8c..0000000 --- a/tests/nextflow/samples/ont_full/sylph_taxonomy_report.tsv +++ /dev/null @@ -1,10 +0,0 @@ -#SampleID reads.fastq.gz Taxonomies_used:['taxonomies1/gtdb_r220_metadata.tsv.gz'] -clade_name relative_abundance sequence_abundance ANI (if strain-level) Coverage (if strain-level) -d__Bacteria 100.0 96.3226 NA NA -d__Bacteria|p__Actinomycetota 100.0 96.3226 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes 100.0 96.3226 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales 100.0 96.3226 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae 100.0 96.3226 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae|g__Mycobacterium 100.0 96.3226 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae|g__Mycobacterium|s__Mycobacterium tuberculosis 100.0 96.3226 NA NA -d__Bacteria|p__Actinomycetota|c__Actinomycetes|o__Mycobacteriales|f__Mycobacteriaceae|g__Mycobacterium|s__Mycobacterium tuberculosis|t__GCF_000195955.2 100.0 96.3226 99.96 15.07 diff --git a/tests/nextflow/sylph.nf.test b/tests/nextflow/sylph.nf.test deleted file mode 100644 index 8e1122b..0000000 --- a/tests/nextflow/sylph.nf.test +++ /dev/null @@ -1,175 +0,0 @@ -nextflow_process { - name "Test Process Sylph" - script "process/sylph.nf" - process "sylph" - - test("Illumina - sylph should run on a sample") { - tag 'illumina' - when { - process{ - """ - input[0] = ['sample_test', [file("$projectDir/tests/nextflow/samples/clinic/test_clinic_sample_1.fastq.gz"), file("$projectDir/tests/nextflow/samples/clinic/test_clinic_sample_2.fastq.gz")]] - input[1] = [file("$projectDir/sylph_dbs/gtdb-r220-c1000-dbv1.syldb")] - input[2] = [file("$projectDir/sylph_dbs/gtdb_r220_metadata.tsv.gz")] - input[3] = 'illumina' - input[4] = 90 - input[5] = 95 - """ - } - } - - then { - assert process.success - - with(process.out.sylph_report.get(0)) { - assert get(1).contains("sylph_report.tsv") - assert (path(get(1)).md5 == path("$projectDir/tests/nextflow/samples/clinic/sylph_report.tsv").md5) - } - with(process.out.taxonomy_report.get(0)) { - assert get(1).contains("taxonomy_report.tsv") - assert (path(get(1)).md5 == path("$projectDir/tests/nextflow/samples/clinic/sylph_taxonomy_report.tsv").md5) - } - } - } - - test("Illumina - sylph should run with edited taxonomy metadata") { - tag 'illumina' - when { - process{ - """ - input[0] = ['sample_test', [file("$projectDir/tests/nextflow/samples/clinic/test_clinic_sample_1.fastq.gz"), file("$projectDir/tests/nextflow/samples/clinic/test_clinic_sample_2.fastq.gz")]] - input[1] = [file("$projectDir/sylph_dbs/gtdb-r220-c1000-dbv1.syldb")] - input[2] = [file("$projectDir/tests/nextflow/other_files/myco_taxonomy_edited.gz")] - input[3] = 'illumina' - input[4] = 90 - input[5] = 95 - """ - } - } - - then { - assert process.success - - with(process.out.sylph_report.get(0)) { - assert get(1).contains("sylph_report.tsv") - assert (path(get(1)).md5 == path("$projectDir/tests/nextflow/samples/clinic/sylph_report.tsv").md5) - } - with(process.out.taxonomy_report.get(0)) { - assert get(1).contains("taxonomy_report.tsv") - assert path(get(1)).readLines().drop(1) == path("$projectDir/tests/nextflow/samples/clinic/sylph_taxonomy_report.tsv").readLines().drop(1) - } - } - } - - test("Illumina - sylph should run on an empty sample") { - tag 'illumina' - when { - process{ - """ - input[0] = ['sample_test', [file("$projectDir/tests/nextflow/samples/empty_samples/reads_1.fastq.gz"), file("$projectDir/tests/nextflow/samples/empty_samples/reads_2.fastq.gz")]] - input[1] = [file("$projectDir/sylph_dbs/gtdb-r220-c1000-dbv1.syldb")] - input[2] = [file("$projectDir/sylph_dbs/gtdb_r220_metadata.tsv.gz")] - input[3] = 'illumina' - input[4] = 90 - input[5] = 95 - """ - } - } - - then { - assert process.success - - with(process.out.sylph_report.get(0)) { - assert get(1).contains("sylph_report.tsv") - } - with(process.out.taxonomy_report.get(0)) { - assert get(1).contains("taxonomy_report.tsv") - } - } - } - - test("Illumina - sylph should run on multiple dbs") { - tag 'illumina' - when { - process{ - """ - input[0] = ['sample_test', [file("$projectDir/tests/nextflow/samples/clinic/test_clinic_sample_1.fastq.gz"), file("$projectDir/tests/nextflow/samples/clinic/test_clinic_sample_2.fastq.gz")]] - input[1] = [file("$projectDir/sylph_dbs/gtdb-r220-c1000-dbv1.syldb"), file("$projectDir/sylph_dbs/virus_imgvr_c200_v0.3.0.syldb")] - input[2] = [file("$projectDir/sylph_dbs/gtdb_r220_metadata.tsv.gz"), file("$projectDir/sylph_dbs/virus_imgvr_4.1_metadata.tsv.gz")] - input[3] = 'illumina' - input[4] = 90 - input[5] = 95 - """ - } - } - - then { - assert process.success - - with(process.out.sylph_report.get(0)) { - assert get(1).contains("sylph_report.tsv") - assert (path(get(1)).md5 == path("$projectDir/tests/nextflow/samples/clinic/with_virus_sylph_report.tsv").md5) - } - with(process.out.taxonomy_report.get(0)) { - assert get(1).contains("taxonomy_report.tsv") - assert (path(get(1)).md5 == path("$projectDir/tests/nextflow/samples/clinic/with_virus_sylph_taxonomy_report.tsv").md5) - } - } - } - - test("ONT - should run on a sample") { - tag 'ont' - when { - process{ - """ - input[0] = ['sample_test', file("$projectDir/tests/nextflow/samples/ont_full/reads.fastq.gz")] - input[1] = [file("$projectDir/sylph_dbs/gtdb-r220-c1000-dbv1.syldb")] - input[2] = [file("$projectDir/sylph_dbs/gtdb_r220_metadata.tsv.gz")] - input[3] = 'ont' - input[4] = 90 - input[5] = 95 - """ - } - } - - then { - assert process.success - - with(process.out.sylph_report.get(0)) { - assert get(1).contains("sylph_report.tsv") - assert (path(get(1)).md5 == path("$projectDir/tests/nextflow/samples/ont_full/sylph_report.tsv").md5) - } - with(process.out.taxonomy_report.get(0)) { - assert get(1).contains("taxonomy_report.tsv") - assert (path(get(1)).md5 == path("$projectDir/tests/nextflow/samples/ont_full/sylph_taxonomy_report.tsv").md5) - } - } - } - - test("ONT - sylph should run on an empty sample") { - tag 'ont' - when { - process{ - """ - input[0] = ['sample_test', file("$projectDir/tests/nextflow/samples/empty_samples/reads_1.fastq.gz")] - input[1] = [file("$projectDir/sylph_dbs/gtdb-r220-c1000-dbv1.syldb")] - input[2] = [file("$projectDir/sylph_dbs/gtdb_r220_metadata.tsv.gz")] - input[3] = 'ont' - input[4] = 90 - input[5] = 95 - """ - } - } - - then { - assert process.success - - with(process.out.sylph_report.get(0)) { - assert get(1).contains("sylph_report.tsv") - } - with(process.out.taxonomy_report.get(0)) { - assert get(1).contains("taxonomy_report.tsv") - } - } - } -}