diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 0dee238..3be8ff8 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -2,16 +2,15 @@ name: Tests on: push: - branches: [ main ] + branches: [main] pull_request: - branches: [ main ] - + branches: [main] jobs: Formatting: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v6 - name: Formatting uses: github/super-linter@v5 env: @@ -23,33 +22,33 @@ jobs: Linting: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v3 - - name: Lint workflow - uses: snakemake/snakemake-github-action@v1 - with: - directory: . - snakefile: workflow/Snakefile - args: "--lint" + - uses: actions/checkout@v6 + - name: Lint workflow + uses: snakemake/snakemake-github-action@v2 + with: + directory: . + snakefile: workflow/Snakefile + args: "--lint" # Testing: - # runs-on: ubuntu-latest - # needs: - # - Linting - # - Formatting +# runs-on: ubuntu-latest +# needs: +# - Linting +# - Formatting # steps: - # - uses: actions/checkout@v3 +# - uses: actions/checkout@v3 # - name: Test workflow - # uses: snakemake/snakemake-github-action@v1 - # with: - # directory: .tests - # snakefile: workflow/Snakefile - # args: "--use-conda --show-failed-logs --cores 3 --conda-cleanup-pkgs cache --all-temp" +# uses: snakemake/snakemake-github-action@v1 +# with: +# directory: .tests +# snakefile: workflow/Snakefile +# args: "--use-conda --show-failed-logs --cores 3 --conda-cleanup-pkgs cache --all-temp" # - name: Test report - # uses: snakemake/snakemake-github-action@v1 - # with: - # directory: .test - # snakefile: workflow/Snakefile - # args: "--report report.zip" +# uses: snakemake/snakemake-github-action@v1 +# with: +# directory: .test +# snakefile: workflow/Snakefile +# args: "--report report.zip" diff --git a/README.md b/README.md index 21adc6a..f0c7b89 100644 --- a/README.md +++ b/README.md @@ -1,121 +1,78 @@ # ResMAG - name pending -[![Snakemake](https://img.shields.io/badge/snakemake-≥6.3.0-brightgreen.svg)](https://snakemake.github.io) +[![Snakemake](https://img.shields.io/badge/snakemake-≥9.0.0-brightgreen.svg)](https://snakemake.github.io) [![GitHub actions status](https://github.com///workflows/Tests/badge.svg?branch=main)](https://github.com///actions?query=branch%3Amain+workflow%3ATests) -ResMAG is a state-of-the-art and user-friendly Snakemake workflow designed for the analysis of metagenomic data. It integrates multiple bioinformatics tools and algorithms to facilitate key steps in metagenome analysis, including bin refinement, metagenome-assembled genome (MAG) reconstruction, taxonomic classification of MAGs, and identification of antibiotic resistance genes.
+ResMAG is a state-of-the-art and user-friendly Snakemake workflow designed for the analysis of metagenomic data. It integrates multiple bioinformatics tools and algorithms to facilitate key steps in metagenome analysis, including quality control, assembling, bin refinement, metagenome-assembled genome (MAG) reconstruction and taxonomic classification. ResMAG has a special focus on highly diverse samples, such as wastewater, and the identification of antibiotic resistance genes.
-### Key Features +--- -**Binning Techniques**: Employ a collection of five state-of-the-art binning tools to partition metagenomic contigs into individual bins, allowing for comprehensive and accurate analysis.
+## Key Features -**MAG Reconstruction**: Utilize cutting-edge algorithms to reconstruct high-quality metagenome-assembled genomes (MAGs) from sequencing data.
+**Binning Techniques**: Employ a collection of two state-of-the-art binning tools to partition metagenomic contigs into individual bins, allowing for comprehensive and accurate analysis.
-**Taxonomic Classification**: Apply advanced taxonomic classification methods to assign taxonomic labels to MAGs and identify the microbial community composition within the metagenomic samples.
+**MAG Reconstruction**: Utilize cutting-edge algorithms to reconstruct high-quality metagenome-assembled genomes (MAGs), especially from highly diverse microbial communities, like wastewater.
+ +**Taxonomic Classification**: Apply advanced taxonomic classification methods to assign taxonomic labels to reads, contigs and MAGs and identify the microbial community composition within the metagenomic samples.
**Antibiotic Resistance Gene Identification**: Perform in-depth analysis to detect and characterize antibiotic resistance genes within the metagenomic data, providing valuable insights into antimicrobial resistance profiles.
**Performance Refinement**: Continuously optimize the pipeline by incorporating the latest advancements in metagenomics research, ensuring the highest accuracy and efficiency in metagenomic data analysis.
-### Overview -```mermaid -%%{init: { - 'theme':'base', - 'themeVariables': { - 'secondaryColor': '#fff', - 'tertiaryColor': '#fff', - 'tertiaryBorderColor' : '#fff'} - }}%% - -flowchart TB; - - subgraph " " - direction TB - - %% Nodes - A[/short reads/] - B["QC
fastp and fastQC"] - C["Filtering human reads
minimap2"] - L["optional
Filtering host reads
minimap2"] - N["taxonomic classification
Kaiju"] - D["Assembly
MegaHIT"] - E["Binning
MetaBAT and MetaCoAG"] - F["Bin refinement
DAS Tool"] - G["Bin quality control
CheckM2"] - M[/MAGs/] - H["Resistance analysis
CARD"] - I["taxonomic classification
GTDBTk"] - J[/MultiQC report/] - K[/Assembly summary/] - - %% input & output node design - classDef in_output fill:#fff,stroke:#cde498,stroke-width:4px - class A,J,K,M in_output - %% rule node design - classDef rule fill:#cde498,stroke:#000 - class B,C,D,E,F,G,H,I,L,N rule - - %% Node links - A --> B - B --> C - B --- J - C --> L - L --> D - C --> N - L --> N - C --> D - D --> E - D --- K - E --> F - F --> G - G --"filtering"--> M - M --- H - M --- I - - end +--- -``` +## Overview -## Usage +--- -### Preparations +## Preparations -#### Clone workflow -To prepare the workflow clone it to your desired working folder via git or your preferred IDE -``` -git clone https://github.com/IKIM-Essen/metagenomics_res.git -``` +Create a snakemake environment using [mamba](https://mamba.readthedocs.io/en/latest/) via: -#### Download GTDB -The GTDB needs to be downloaded and decompressed, it requires about 110 Gb. -1. Change to the cloned workflow directory -2. Create a new folder `resources/gtdb/` and change to this directory -3. Download the latest version of GTDB - **or** - if you have already downloaded a version of GTDB move the `gtdbtk_data.tar.gz` file to `resources/gtdb/` -4. Unarchive the downloaded file -5. After successful step 4: the archive can be removed + ```mamba create -c conda-forge -c bioconda -n snakemake snakemake snakemake-storage-plugin-fs``` -``` -wget https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/auxillary_files/gtdbtk_package/full_package/gtdbtk_data.tar.gz -tar xvzf gtdbtk_data.tar.gz -``` +For installation details, see the [instructions in the Snakemake documentation](https://snakemake.readthedocs.io/en/stable/getting_started/installation.html). -#### Install Snakemake -Create a snakemake environment using [mamba](https://mamba.readthedocs.io/en/latest/) via: +### Download GTDB +The GTDB needs to be downloaded and decompressed, it requires about 140 Gb. - ```mamba create -c conda-forge -c bioconda -n snakemake snakemake=7.32.3``` +1. Change to the directory where the GTDB should be stored +2. Download the latest or your desired version of GTDB + Please make sure this version is compatible with GTDB-tk version `2.6.1` + ``` + wget https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/auxillary_files/gtdbtk_package/full_package/gtdbtk_data.tar.gz + ``` +3. Decompress the downloaded archive + ``` + tar xzf gtdbtk_data.tar.gz + ``` +4. After successful step 3: the archive can be removed +5. Please specify the path to your decompressed GTDB in the config file (see [Configuring workflow](#configuring-the-workflow)) -For installation details, see the [instructions in the Snakemake documentation](https://snakemake.readthedocs.io/en/stable/getting_started/installation.html). +### Create UniCARD database +Please follow the steps outlined in the [UniCARD repository](https://github.com/IKIM-Essen/uniCARD) to create a UniCARD database for resistance gene annotation. +And specify the path to your database in the config file (see [Configuring workflow](#configuring-the-workflow)) + +--- -#### Configuring workflow +## Usage + +Please obtain a copy of this workflow by [cloning](https://docs.github.com/en/repositories/creating-and-managing-repositories/cloning-a-repository) this repository. +``` +git clone https://github.com/IKIM-Essen/metagenomics_res.git +``` + +### Configuring the workflow 1. Edit the `config/config.yaml` file: - Specify a project name (`project-name`) - Specify filtering options for human reads (`human-filtering`) - Specify host filtering options, if you have a non-human host (`host-filtering`) - - Specify options for GTDB database (see [Download GTDB](#Download-GTDB)) -2. Provide sample information in the `config/pep/samples.csv` file while keeping the header and the format: + - Specify options for different databases: + - GTDB database needs to be downloaded before (see [Download GTDB](#download-gtdb)) + - UniCARD database needs to be create before (see [Create UniCARD database](#create-unicard-database)) + - other databases (kaiju, CheckM2, CARD, genomad) can be given as a local path or downloaded when running the pipeline +2. Provide sample information in the `config/pep/samples.csv` file while keeping the header and the format as: ``` sample_name,fq1,fq2 @@ -131,47 +88,66 @@ Test your configuration by performing a dry-run via ```snakemake --use-conda -n``` Executing the workflow: -```snakemake --use-conda --cores $N --rerun-incomplete``` +```snakemake --use-conda --cores $N -k``` using `$N` cores. It is recommended to use all available cores. -The usage of this workflow is described in the [Snakemake Workflow Catalog](https://snakemake.github.io/snakemake-workflow-catalog/?usage=%2F). +--- ## Output +--- + ## Tools A list of tools used in the pipeline: | Tool | Link | | --------- | ------------------------------------------------- | +| CARD | https://doi.org/10.1093/nar&gkac920 | | CheckM2 | https://doi.org/10.1038/s41592-023-01940-w | | CoverM | https://github.com/wwood/CoverM | | DAS Tool | https://doi.org/10.1038/s41564-018-0171-1 | +| DIAMOND | https://doi.org/10.1038/s41592-021-01101-x | | fastp | https://doi.org/10.1093/bioinformatics/bty560 | | FastQC | www.bioinformatics.babraham.ac.uk/projects/fastqc | +| geNomad | https://doi.org/10.1038/s41587-023-01953-y | +| GTDB-Tk | https://doi.org/10.1093/bioinformatics/btz848 | | Kaiju | https://doi.org/10.1038/ncomms11257 | | MEGAHIT | https://doi.org/10.1093/bioinformatics/btv033 | | MetaBAT | http://dx.doi.org/10.7717/peerj.1165 | | MetaCoAG | https://doi.org/10.1101/2021.09.10.459728 | | minimap2 | https://doi.org/10.1093/bioinformatics/bty191 | | MultiQC | www.doi.org/10.1093/bioinformatics/btw354 | -| pandas | pandas.pydata.org | +| pprodigal | https://github.com/sjaenick/pprodigal | +| Rust-Bio | https://doi.org/10.1093/bioinformatics/btv573 | | samtools | https://doi.org/10.1093/gigascience/giab008 | | Snakemake | www.doi.org/10.12688/f1000research.29032.1 | +| UniCARD | https://github.com/IKIM-Essen/uniCARD | + +--- + +## Contributions -## Contact Information +Pull requests and feature suggestions are very welcome! Feel free to fork and submit improvements. +For any other questions, or feedback, please contact the project maintainer Josefa Welling (@josefawelling) at josefa.welling@uk-essen.de. -For any questions, or feedback, please contact the project maintainer at josefa.welling@uk-essen.de. We appreciate your input and support in using and improving ResMAG. +We appreciate your input and support in using and improving ResMAG. + +--- ## Acknowledgements We would like to express our gratitude towards Katharina Block, Adrian Doerr, Miriam Balzer, Alexander Thomas, Johannes Köster, Ann-Kathrin Doerr and the IKIM who have contributed to the development and testing of ResMAG. Their valuable insights and feedback have been helpful throughout the creation of the workflow. +--- + ## Citation -A paper is on its way. If you use ResMAG in your work before the paper, then please consider citing this GitHub. +A paper is on its way. If you use ResMAG in your work, don't forget to give credits to the authors by citing the URL of this repository. + +--- ## License -ResMAG is released under the [BSD-2 Clause](https://www.open-xchange.com/hubfs/2_Clause_BSD_License.pdf?hsLang=en). Please review the license file for more details. \ No newline at end of file +ResMAG is licensed under the [BSD-2 Clause](https://www.open-xchange.com/hubfs/2_Clause_BSD_License.pdf?hsLang=en). \ No newline at end of file diff --git a/config/config.yaml b/config/config.yaml index a864f70..61b3a1a 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -3,14 +3,6 @@ pepfile: config/pep/config.yaml ## Please change to a name describing your project ## All results can be found under results/project-name/ project-name: "test" -# -run-date: "2025-06-26" - -data-handling: - # path to store data within the workflow - data: data/ - # path where databases and reference genomes are stored - resources: resources/ quality-criteria: # minimal length of acceptable reads @@ -23,11 +15,11 @@ adapter-seqs: "-a CTGTCTCTTATACACATCT -g AGATGTGTATAAGAGACAG" # reads are always filtered for human contamination human-filtering: # if you have a human reference genome on your device, - # you can set use-local to True and specify the path to the human reference genome - use-local: False + # you can set use-shared to True and specify the path to the human reference genome + use-shared: True # path to a locally stored human reference genome - local-path: /groups/ds/databases_refGenomes/refGenomes/latest/human/GCA_000001405.29_GRCh38.p14_genomic.fna.gz - # if use-local = False the reference genome is downloaded via the following url + local-path: /groups/ds/databases_refGenomes/refGenomes/latest/human/GRCh38_latest_genomic.fna.gz + # if use-shared = False the reference genome is downloaded via the following url download-path: https://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/annotation/GRCh38_latest/refseq_identifiers/GRCh38_latest_genomic.fna.gz # if filtering for non-human host is required, @@ -42,35 +34,46 @@ host-filtering: # minimum contig length used during assembly step min-contig-length: 300 -das-tool: - # cores to use for das_tool - threads: 64 - binner-list: ["metabat", "metacoag"] - MAG-criteria: min-completeness: 50.00 max-contamination: 30.00 ### Handling of different databases +## The GTDB has about 140 Gb +# please specify the folder where the decompressed database is stored +gtdb: /groups/ds/databases_refGenomes/databases/latest/gtdb/release226/ + kaiju: + use-shared: True + # path of kaiju DB .fmi file + # if use-shared False: DB is downloaded and saved to this path + fmi-file: /groups/ds/databases_refGenomes/databases/latest/kaiju/refseq_nr/unpacked_2024-08-13/kaiju_db_refseq_nr.fmi download: https://kaiju-idx.s3.eu-central-1.amazonaws.com/2024/kaiju_db_refseq_nr_2024-08-13.tgz - fmi-file: kaiju_db/kaiju_db_refseq_nr.fmi -checkm2: CheckM2_database/uniref100.KO.1.dmnd +checkm2: + use-shared: True + #if use-shared False a folder named CheckM2_database will be created under the given path + db-folder: /groups/ds/databases_refGenomes/databases/latest/checkM2/ -## The GTDB has about 100 Gb -## if use-local is set to False this will be downloaded during the pipeline run -gtdb: - use-local: True - # if use-local is set to True, please specify the folder where the decompressed database is stored - # this path is expected to lay under the data-handling resources folder - db-folder: gtdb/release226/ +# folder of genomad DB +genomad: + use-shared: True + # a folder named genomad_db is expected under the given path + # if use-shared False this folder will be created + db-folder: /groups/ds/databases_refGenomes/databases/latest/ card: - version: v4.0.0 - dbfile: card.json + use-shared: True + db-file: /groups/ds/databases_refGenomes/databases/latest/card/card.json + version: v4.0.1 url: https://card.mcmaster.ca/latest/data +uniCARD: + #path to the fasta.gz file of UniCARD database + db-file: /groups/ds/databases_refGenomes/databases/latest/uniCARD/uniCARD.fasta.gz + # path to CARDs hierarchy file + hierarchy-json: /groups/ds/databases_refGenomes/databases/latest/uniCARD/CARD_hierarchy_v4.0.0.json + ## string term used for formatting output tables tablular-config: '/>github<\/a>/a \\t\t\t\n\t\t\t
  • \n\t\t\t\t' diff --git a/config/multiqc_config.yaml b/config/multiqc_config.yaml index 24f4723..005aa10 100644 --- a/config/multiqc_config.yaml +++ b/config/multiqc_config.yaml @@ -31,7 +31,8 @@ fn_clean_exts: - "_L002_R2_001" - "_S" - ".1" - + - ".2" + # customising general Statistics table_columns_visible: Reads Quality Control: @@ -46,4 +47,4 @@ table_columns_visible: after_filtering_q30_bases: False after_filtering_gc_content: False pct_surviving: True - pct_adapter: False \ No newline at end of file + pct_adapter: False diff --git a/workflow/Snakefile b/workflow/Snakefile index f8e0c36..3d0b47e 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -5,7 +5,7 @@ from snakemake.utils import min_version -min_version("6.3.0") +min_version("9.0.0") configfile: "config/config.yaml" @@ -20,7 +20,6 @@ report: "report/workflow.rst" include: "rules/common.smk" include: "rules/qc.smk" include: "rules/host_filtering.smk" -include: "rules/trimming.smk" include: "rules/assembly.smk" include: "rules/metabat.smk" include: "rules/metacoag.smk" @@ -28,6 +27,7 @@ include: "rules/das_tool.smk" include: "rules/bin_qc.smk" include: "rules/classify.smk" include: "rules/analysis.smk" +include: "rules/resistance.smk" include: "rules/report.smk" @@ -35,7 +35,7 @@ rule all: input: # read classification output expand( - "results/{project}/output/classification/reads/{sample}/{sample}_{level}_kaiju_summary.tsv", + "results/{project}/output/classification/reads/{sample}/{sample}_{level}_abundance.tsv", project=get_project(), sample=get_samples(), level=["species", "genus", "family", "order", "class", "phylum"], @@ -63,48 +63,24 @@ rule all: project=get_project(), sample=get_samples(), ), - # resistance analysis - expand( - "results/{project}/output/ARGs/reads/{sample}/{sample}_read_ARGs.csv", - project=get_project(), - sample=get_samples(), - ), - expand( - "results/{project}/output/ARGs/assembly/{sample}/{sample}_assembly_ARGs.csv", - project=get_project(), - sample=get_samples(), - ), # overall report expand( "results/{project}/output/report/report_{project}.zip", project=get_project(), ), + # resistance analysis expand( - "results/{project}/megahit/{sample}_cleanup.done", - project=get_project(), - sample=get_samples(), - ), - expand( - "results/{project}/binning/metabat/{sample}_cleanup.done", - project=get_project(), - sample=get_samples(), - ), - expand( - "results/{project}/binning/metacoag/{sample}_cleanup.done", + "results/{project}/output/resistance/CARD/reads/{sample}/{sample}.gene_mapping_data.txt", project=get_project(), sample=get_samples(), ), expand( - "results/{project}/binning/das_tool/{sample}_cleanup.done", + "results/{project}/output/resistance/uniCARD/{project}_drug_class_abundance_all_possible.csv", project=get_project(), - sample=get_samples(), ), - - -""" + #contig classification expand( - "results/{project}/output/ARGs/mags/{sample}/all_mags.done", + "results/{project}/output/classification/assembly/{sample}/kaiju.out.gz", project=get_project(), sample=get_samples(), ), -""" diff --git a/workflow/envs/card.yaml b/workflow/envs/card.yaml index 88f740f..1c69e2a 100644 --- a/workflow/envs/card.yaml +++ b/workflow/envs/card.yaml @@ -3,5 +3,4 @@ channels: - bioconda - nodefaults dependencies: - - rgi = 6.0.3 - - kma = 1.4.9 \ No newline at end of file + - rgi = 6.0.5 diff --git a/workflow/envs/checkm2.yaml b/workflow/envs/checkm2.yaml index a9d75be..40d8f09 100644 --- a/workflow/envs/checkm2.yaml +++ b/workflow/envs/checkm2.yaml @@ -3,4 +3,4 @@ channels: - bioconda - nodefaults dependencies: - - checkm2 = 1.0.1 \ No newline at end of file + - checkm2 = 1.1.0 diff --git a/workflow/envs/coverm.yaml b/workflow/envs/coverm.yaml index 4ef59df..3fc7948 100644 --- a/workflow/envs/coverm.yaml +++ b/workflow/envs/coverm.yaml @@ -3,4 +3,4 @@ channels: - bioconda - nodefaults dependencies: - - coverm = 0.6.1 \ No newline at end of file + - coverm = 0.7.0 diff --git a/workflow/envs/das_tool.yaml b/workflow/envs/das_tool.yaml index bcb2764..3c32a4e 100644 --- a/workflow/envs/das_tool.yaml +++ b/workflow/envs/das_tool.yaml @@ -2,4 +2,4 @@ channels: - conda-forge - bioconda dependencies: - - das_tool=1.1.6 + - das_tool=1.1.7 diff --git a/workflow/envs/diamond.yaml b/workflow/envs/diamond.yaml new file mode 100644 index 0000000..42d39a8 --- /dev/null +++ b/workflow/envs/diamond.yaml @@ -0,0 +1,6 @@ +channels: + - conda-forge + - bioconda + - nodefaults +dependencies: + - diamond = 2.1.12 diff --git a/workflow/envs/fastqc.yaml b/workflow/envs/fastqc.yaml new file mode 100644 index 0000000..2e917d0 --- /dev/null +++ b/workflow/envs/fastqc.yaml @@ -0,0 +1,6 @@ +channels: + - conda-forge + - bioconda + - nodefaults +dependencies: + - fastqc = 0.12.1 diff --git a/workflow/envs/genomad.yaml b/workflow/envs/genomad.yaml index d6589cf..3424746 100644 --- a/workflow/envs/genomad.yaml +++ b/workflow/envs/genomad.yaml @@ -1,6 +1,6 @@ -channels: +channels: - conda-forge - bioconda - nodefaults dependencies: - - genomad = 1.8.0 + - genomad = 1.11.2 diff --git a/workflow/envs/gtdbtk.yaml b/workflow/envs/gtdbtk.yaml index 3144b50..46b47f4 100644 --- a/workflow/envs/gtdbtk.yaml +++ b/workflow/envs/gtdbtk.yaml @@ -3,5 +3,6 @@ channels: - bioconda - nodefaults dependencies: - - gtdbtk = 2.4.1 - - numpy = 1.23.1 + - gtdbtk = 2.6.1 + # fix because gtdbtk doesn't work with python 3.14 + - python = 3.13 diff --git a/workflow/envs/metabat.yaml b/workflow/envs/metabat.yaml index 47d94d6..a31cac3 100644 --- a/workflow/envs/metabat.yaml +++ b/workflow/envs/metabat.yaml @@ -3,4 +3,4 @@ channels: - bioconda - nodefaults dependencies: - - metabat2 = 2.15 \ No newline at end of file + - metabat2 = 2.18 diff --git a/workflow/envs/metacoag.yaml b/workflow/envs/metacoag.yaml index 6a3e0d4..e8c0241 100644 --- a/workflow/envs/metacoag.yaml +++ b/workflow/envs/metacoag.yaml @@ -3,5 +3,5 @@ channels: - bioconda - nodefaults dependencies: - - metacoag = 1.2.1 + - metacoag = 1.2.2 - biopython = 1.81 diff --git a/workflow/envs/minimap2.yaml b/workflow/envs/minimap2.yaml index 30b228a..f625c9c 100644 --- a/workflow/envs/minimap2.yaml +++ b/workflow/envs/minimap2.yaml @@ -3,8 +3,8 @@ channels: - bioconda - nodefaults dependencies: - - minimap2 = 2.24 - - samtools = 1.16.1 + - minimap2 = 2.30 + - samtools = 1.23 - pip - pip: - - Command \ No newline at end of file + - Command diff --git a/workflow/envs/pprodigal.yaml b/workflow/envs/pprodigal.yaml new file mode 100644 index 0000000..dad1fe8 --- /dev/null +++ b/workflow/envs/pprodigal.yaml @@ -0,0 +1,6 @@ +channels: + - conda-forge + - bioconda + - nodefaults +dependencies: + - pprodigal = 1.0.1 diff --git a/workflow/envs/python.yaml b/workflow/envs/python.yaml index bf746a3..1682687 100644 --- a/workflow/envs/python.yaml +++ b/workflow/envs/python.yaml @@ -3,9 +3,9 @@ channels: - bioconda - anaconda dependencies: - - altair = 5.2.0 - - biopython = 1.82 + - altair = 6.0.0 + - biopython = 1.85 - distinctipy = 1.3.4 - - numpy = 1.26.3 - - matplotlib = 3.8.2 - - pandas = 2.1.4 + - numpy = 2.4.1 + - matplotlib = 3.10.8 + - pandas = 2.3.3 diff --git a/workflow/envs/rbt.yaml b/workflow/envs/rbt.yaml index c11823a..ebc035f 100644 --- a/workflow/envs/rbt.yaml +++ b/workflow/envs/rbt.yaml @@ -3,4 +3,4 @@ channels: - bioconda - nodefaults dependencies: - - rust-bio-tools = 0.42.0 \ No newline at end of file + - rust-bio-tools = 0.42.2 diff --git a/workflow/envs/snakemake.yaml b/workflow/envs/snakemake.yaml index 612d3a1..9551cc7 100644 --- a/workflow/envs/snakemake.yaml +++ b/workflow/envs/snakemake.yaml @@ -2,4 +2,5 @@ channels: - conda-forge - bioconda dependencies: - - snakemake=7.32.4 \ No newline at end of file + - snakemake=9.16.3 + - snakemake-storage-plugin-fs=1.1.3 diff --git a/workflow/envs/unix.yaml b/workflow/envs/unix.yaml index f148464..fc89a8a 100644 --- a/workflow/envs/unix.yaml +++ b/workflow/envs/unix.yaml @@ -6,4 +6,4 @@ dependencies: - git - sed - curl - - pigz \ No newline at end of file + - gzip diff --git a/workflow/report/host_plot.rst b/workflow/report/host_plot.rst deleted file mode 100644 index 1388dd2..0000000 --- a/workflow/report/host_plot.rst +++ /dev/null @@ -1,2 +0,0 @@ -A barplot for each sample, showing the percentage of human reads and if given also the percentage of reads of another host. -The reads are classified by mapping against a reference genome using `minimap2 `_ and `samtools `_. \ No newline at end of file diff --git a/workflow/rules/analysis.smk b/workflow/rules/analysis.smk index 1f59532..6a5262e 100644 --- a/workflow/rules/analysis.smk +++ b/workflow/rules/analysis.smk @@ -1,222 +1,96 @@ -# Plasmid analysis -rule load_genomad_DB: +# gene identification +rule gene_identification: + input: + contigs=get_assembly, output: - file=get_genomad_DB_file(), - params: - folder=lambda wildcards, output: Path(output.file).parent.parent, + faa="results/{project}/output/proteins/{sample}/{sample}_proteins.faa", + fna="results/{project}/output/proteins/{sample}/{sample}_nucleotides.fna", + gff="results/{project}/output/proteins/{sample}/{sample}_annotations.gff", log: - "logs/load_genomad_DB.log", + "logs/{project}/proteins/{sample}.log", conda: - "../envs/genomad.yaml" + "../envs/pprodigal.yaml" + threads: 32 shell: - "genomad download-database {params.folder}/ > {log} 2>&1" + "pprodigal -i {input.contigs} -o {output.gff} -a {output.faa} " + "-d {output.fna} -p meta --tasks {threads} > {log} 2>&1" + + +rule gzip_proteins: + input: + faa=rules.gene_identification.output.faa, + fna=rules.gene_identification.output.fna, + gff=rules.gene_identification.output.gff, + output: + faa="results/{project}/output/proteins/{sample}/{sample}_proteins.faa.gz", + fna="results/{project}/output/proteins/{sample}/{sample}_nucleotides.fna.gz", + gff="results/{project}/output/proteins/{sample}/{sample}_annotations.gff.gz", + log: + "logs/{project}/proteins/{sample}_gzip.log", + priority: 1 + conda: + "../envs/unix.yaml" + threads: 20 + shell: + "gzip {input.faa} {input.fna} {input.gff} > {log} 2>&1" + + +# Plasmid analysis + +if not config["genomad"]["use-shared"]: + + rule load_genomad_DB: + output: + folder=get_genomad_DB_folder(), + file=get_genomad_DB_file(), + log: + "logs/load_genomad_DB.log", + conda: + "../envs/genomad.yaml" + params: + res_folder=lambda wildcards, output: Path(output.folder).parent, + shell: + "genomad download-database {params.res_folder}/ > {log} 2>&1" rule genomad_run: input: - db=rules.load_genomad_DB.output.file, + db=get_genomad_DB_folder(), asmbl=rules.gzip_assembly.output, output: - outdir=temp(directory("results/{project}/genomad/{sample}/")), plasmid_tsv=temp( "results/{project}/genomad/{sample}/{sample}_summary/{sample}_plasmid_summary.tsv" ), virus_tsv=temp( "results/{project}/genomad/{sample}/{sample}_summary/{sample}_virus_summary.tsv" ), - params: - db_folder=lambda wildcards, input: Path(input.db).parent, log: "logs/{project}/plasmids/{sample}_run.log", - threads: 64 conda: "../envs/genomad.yaml" + threads: 64 + params: + outdir=lambda wildcards, output: Path(output.plasmid_tsv).parent.parent, shell: "genomad end-to-end --cleanup -t {threads} " - "{input.asmbl} {output.outdir}/ " - "{params.db_folder}/ > {log} 2>&1" + "{input.asmbl} {params.outdir}/ " + "{input.db}/ > {log} 2>&1" rule move_genomad_output: input: - folder=rules.genomad_run.output.outdir, plasmid_tsv=rules.genomad_run.output.plasmid_tsv, virus_tsv=rules.genomad_run.output.virus_tsv, output: plasmid_tsv="results/{project}/output/plasmids/{sample}/{sample}_plasmid_summary.tsv", virus_tsv="results/{project}/output/plasmids/{sample}/{sample}_virus_summary.tsv", - params: - outdir=lambda wildcards, output: Path(output.plasmid_tsv).parent, - sum_folder=lambda wildcards, input: Path(input.plasmid_tsv).parent, log: "logs/{project}/plasmids/{sample}_move_output.log", - threads: 64 conda: "../envs/unix.yaml" - shell: - "scp {params.sum_folder}/* {params.outdir}/ > {log} 2>&1" - - -# Resistance analysis -rule download_CARD_data: - output: - json=get_card_db_file(), - params: - download=config["card"]["url"], - folder=lambda wildcards, output: Path(output.json).parent, - filename=lambda wildcards, output: Path(output.json).name, - log: - "logs/CARD_data_download.log", - threads: 30 - conda: - "../envs/unix.yaml" - shell: - "(cd {params.folder} && " - "wget {params.download} && " - "tar -xvf data ./{params.filename}) > {log} 2>&1" - - -rule CARD_load_DB_for_reads: - input: - get_card_db_file(), - output: - temp(touch("results/CARD_load_DB_for_reads.done")), - log: - "logs/CARD_load_DB_for_reads.log", - conda: - "../envs/card.yaml" - shell: - "rgi clean --local && " - "rgi load --card_json {input} --local > {log} 2>&1" - - -rule CARD_annotation: - input: - json=get_card_db_file(), - load=rules.CARD_load_DB_for_reads.output, - output: - temp(touch("results/CARD_annotation.done")), - params: - folder=lambda wildcards, input: Path(input.json).parent, - file=lambda wildcards, input: Path(input.json).name, - ann=get_card_annotation_file(), - log: - "logs/CARD_annotation.log", - conda: - "../envs/card.yaml" - shell: - "(cd {params.folder}/ && " - "rgi card_annotation -i {params.file}) && " - "rgi load -i {input.json} --card_annotation {params.ann} --local > {log} 2>&1" - - -rule CARD_read_run: - input: - fa=get_filtered_gz_fastqs, - db=rules.CARD_annotation.output, - output: - txt="results/{project}/output/ARGs/reads/{sample}/{sample}.gene_mapping_data.txt", - bam=temp( - "results/{project}/output/ARGs/reads/{sample}/{sample}.sorted.length_100.bam" - ), - params: - folder=lambda wildcards, output: Path(output.txt).parent, - log: - "logs/{project}/ARGs/reads/{sample}.log", - threads: 64 - conda: - "../envs/card.yaml" - shell: - "rgi bwt -1 {input.fa[0]} -2 {input.fa[1]} " - "-o {params.folder}/{wildcards.sample} --local " - "--clean -n {threads} > {log} 2>&1" - - -rule CARD_read_sample_summary: - input: - txt=rules.CARD_read_run.output.txt, - output: - csv="results/{project}/output/ARGs/reads/{sample}/{sample}_read_ARGs.csv", - params: - case="reads", - log: - "logs/{project}/ARGs/reads/{sample}.log", - threads: 2 - conda: - "../envs/python.yaml" - script: - "../scripts/arg_summary_sample.py" - - -# updates CARD database to use for contigs instead of reads -# read based classification must be finished before -rule CARD_load_DB: - input: - db=get_card_db_file(), - read_args=expand( - "results/{project}/output/ARGs/reads/{sample}/{sample}.gene_mapping_data.txt", - sample=get_samples(), - project=get_project(), - ), - output: - touch("results/CARD_load_DB.done"), - log: - "logs/CARD_load_DB.log", - conda: - "../envs/card.yaml" - shell: - "rgi clean --local && " - "rgi load --card_json {input.db} --local > {log} 2>&1" - - -rule CARD_assembly_run: - input: - fa=rules.gzip_assembly.output, - db=rules.CARD_load_DB.output, - output: - txt="results/{project}/output/ARGs/assembly/{sample}/{sample}.txt", - json="results/{project}/output/ARGs/assembly/{sample}/{sample}.json", + threads: 32 params: - path_wo_ext=lambda wildcards, output: Path(output.txt).with_suffix(""), - log: - "logs/{project}/ARGs/assembly/{sample}.log", - threads: 64 - conda: - "../envs/card.yaml" + outdir=lambda wildcards, output: Path(output.plasmid_tsv).parent, + sum_folder=lambda wildcards, input: Path(input.plasmid_tsv).parent, shell: - "rgi main -i {input.fa} -o {params.path_wo_ext} " - "-t contig -a DIAMOND --low_quality --local " - "-n {threads} --clean > {log} 2>&1" - - -use rule CARD_read_sample_summary as CARD_assembly_sample_summary with: - input: - txt=rules.CARD_assembly_run.output.txt, - json=rules.CARD_assembly_run.output.json, - output: - csv="results/{project}/output/ARGs/assembly/{sample}/{sample}_assembly_ARGs.csv", - params: - case="assembly", - log: - "logs/{project}/ARGs/assembly/{sample}.log", - - -use rule CARD_assembly_run as CARD_mag_run with: - input: - fa="results/{project}/output/fastas/{sample}/mags/{binID}.fa.gz", - db=rules.CARD_load_DB.output, - output: - txt="results/{project}/output/ARGs/mags/{sample}/{binID}.txt", - params: - path_wo_ext=lambda wildcards, output: Path(output.txt).with_suffix(""), - log: - "logs/{project}/ARGs/mags/{sample}/{binID}.log", - - -rule wrap_mag_ARGs: - input: - get_mag_ARGs, - output: - touch("results/{project}/output/ARGs/mags/{sample}/all_mags.done"), - log: - "logs/{project}/ARGs/mags/{sample}/all_mags.log", + "scp {params.sum_folder}/*.tsv {params.sum_folder}/*.json {params.outdir}/ > {log} 2>&1" diff --git a/workflow/rules/assembly.smk b/workflow/rules/assembly.smk index 3c9b10a..6596587 100644 --- a/workflow/rules/assembly.smk +++ b/workflow/rules/assembly.smk @@ -7,11 +7,13 @@ rule megahit: output: contigs=temp("results/{project}/megahit/{sample}/final.contigs.fa"), outdir=temp(directory("results/{project}/megahit/{sample}/")), - log="results/{project}/report_prerequisites/assembly/{sample}_megahit.log", - done=touch("results/{project}/megahit/{sample}.done"), + log="results/{project}/output/report/prerequisites/assembly/{sample}_megahit.log", params: threshold=get_contig_length_threshold(), threads: 64 + priority: 2 + resources: + heavy=2, log: "logs/{project}/assembly/{sample}_megahit.log", conda: @@ -29,7 +31,7 @@ rule map_to_assembly: fastqs=get_filtered_gz_fastqs, output: bam=temp( - "results/{project}/report_prerequisites/assembly/{sample}_reads_mapped.bam" + "results/{project}/output/report/prerequisites/assembly/{sample}_reads_mapped.bam" ), threads: 64 log: @@ -47,7 +49,7 @@ rule index_assembly_alignment: rules.map_to_assembly.output.bam, output: bai=temp( - "results/{project}/report_prerequisites/assembly/{sample}_reads_mapped.bam.bai" + "results/{project}/output/report/prerequisites/assembly/{sample}_reads_mapped.bam.bai" ), threads: 20 log: @@ -63,7 +65,7 @@ rule reads_mapped_assembly: bam=rules.map_to_assembly.output.bam, bai=rules.index_assembly_alignment.output.bai, output: - bai="results/{project}/report_prerequisites/assembly/{sample}_reads_mapped.txt", + bai="results/{project}/output/report/prerequisites/assembly/{sample}_reads_mapped.txt", threads: 20 log: "logs/{project}/assembly/{sample}_mapping_reads.log", @@ -79,24 +81,33 @@ rule gzip_assembly: contigs=get_assembly, output: "results/{project}/output/fastas/{sample}/{sample}.fa.gz", - threads: 64 + threads: 20 + priority: 1 log: "logs/{project}/assembly/{sample}_gzip.log", conda: "../envs/unix.yaml" shell: - "pigz -c {input.contigs} > {output} 2> {log}" + "gzip -c {input.contigs} > {output} 2> {log}" rule assembly_summary: input: qc_csv=rules.qc_summary.output.csv, asbl=expand( - "results/{{project}}/report_prerequisites/assembly/{sample}_megahit.log", + "results/{{project}}/output/report/prerequisites/assembly/{sample}_megahit.log", sample=get_samples(), ), mapped=expand( - "results/{{project}}/report_prerequisites/assembly/{sample}_reads_mapped.txt", + "results/{{project}}/output/report/prerequisites/assembly/{sample}_reads_mapped.txt", + sample=get_samples(), + ), + csv_bins=expand( + "results/{{project}}/output/report/{sample}/{sample}_bin_summary.csv", + sample=get_samples(), + ), + csv_mags=expand( + "results/{{project}}/output/report/{sample}/{sample}_mags_summary.csv", sample=get_samples(), ), output: @@ -114,13 +125,15 @@ use rule qc_summary_report as assembly_report with: input: rules.assembly_summary.output.vis_csv, output: - report( - directory("results/{project}/output/report/all/assembly/"), - htmlindex="index.html", - category="3. Assembly results", - labels={ - "sample": "all samples", - }, + temp( + report( + directory("results/{project}/output/report/all/assembly/"), + htmlindex="index.html", + category="3. Assembly results", + labels={ + "sample": "all samples", + }, + ) ), params: pin_until="sample", @@ -130,18 +143,3 @@ use rule qc_summary_report as assembly_report with: pattern=config["tablular-config"], log: "logs/{project}/report/assembly_rbt_csv.log", - - -# remove megahit intermediate results when all dependent results are produced -rule cleanup_megahit_output: - input: - # folder to remove - asmbl_folder=rules.megahit.output.outdir, - #dependent results - gz_asmbl=rules.gzip_assembly.output, - asmbl_summary=rules.assembly_summary.output.csv, - binning_done="results/{project}/binning/das_tool/{sample}_run.done", - output: - touch("results/{project}/megahit/{sample}_cleanup.done"), - log: - "logs/{project}/assembly/{sample}_cleanup.log", diff --git a/workflow/rules/bin_qc.smk b/workflow/rules/bin_qc.smk index c2a47be..519b926 100644 --- a/workflow/rules/bin_qc.smk +++ b/workflow/rules/bin_qc.smk @@ -1,15 +1,17 @@ ## bin QC -rule checkm2_DB_download: - output: - dbfile=get_checkm2_db(), #"{}/{}".format(config["data-handling"]["resources"], config["checkm2"]), - params: - direct=lambda wildcards, output: Path(output.dbfile).parent.parent, - log: - "logs/checkm2_DB_download.log", - conda: - "../envs/checkm2.yaml" - shell: - "checkm2 database --download --path {params.direct} > {log} 2>&1" +if not config["checkm2"]["use-shared"]: + + rule checkm2_DB_download: + output: + dbfile=get_checkm2_db(), + log: + "logs/checkm2_DB_download.log", + conda: + "../envs/checkm2.yaml" + params: + direct=lambda wildcards, output: Path(output.dbfile).parent.parent, + shell: + "checkm2 database --download --path {params.direct} > {log} 2>&1" rule checkm2_run: @@ -17,19 +19,19 @@ rule checkm2_run: bins="results/{project}/output/fastas/{sample}/bins/", dbfile=get_checkm2_db(), output: - outdir=temp(directory("results/{project}/qc/checkm2/{sample}/")), stats="results/{project}/output/report/{sample}/checkm2_quality_report.tsv", - params: - outname="quality_report.tsv", log: "logs/{project}/checkm2/{sample}.log", - threads: 24 conda: "../envs/checkm2.yaml" + threads: 24 + params: + outdir="results/{project}/qc/checkm2/{sample}/", shell: "(checkm2 predict -x fa.gz --threads {threads} --force " - "--input {input.bins}/ --output-directory {output.outdir}/ && " - "cp {output.outdir}/{params.outname} {output.stats}) > {log} 2>&1" + "--input {input.bins}/ --output-directory {params.outdir}/ " + "--database_path {input.dbfile} && " + "cp {params.outdir}/quality_report.tsv {output.stats}) > {log} 2>&1" rule bin_summary_sample: @@ -37,21 +39,18 @@ rule bin_summary_sample: tool="results/{project}/output/report/{sample}/{sample}_DASTool_summary.tsv", checkm="results/{project}/output/report/{sample}/checkm2_quality_report.tsv", gtdb="results/{project}/output/classification/bins/{sample}/{sample}.summary.tsv", - #args="results/{project}/output/ARGs/bins/{sample}/all_bins.done", output: csv_bins="results/{project}/output/report/{sample}/{sample}_bin_summary.csv", - csv_checkm="results/{project}/output/report/{sample}/{sample}_checkm2_summary.csv", - csv_dastool="results/{project}/output/report/{sample}/{sample}_DASTool_summary.csv", csv_tax="results/{project}/output/report/{sample}/{sample}_bin_taxonomy.csv", csv_mags="results/{project}/output/report/{sample}/{sample}_mags_summary.csv", - params: - max_cont=config["MAG-criteria"]["max-contamination"], - min_comp=config["MAG-criteria"]["min-completeness"], log: "logs/{project}/bin_summary/{sample}.log", - threads: 4 conda: "../envs/python.yaml" + threads: 4 + params: + max_cont=config["MAG-criteria"]["max-contamination"], + min_comp=config["MAG-criteria"]["min-completeness"], script: "../scripts/bin_summary_sample.py" @@ -60,150 +59,66 @@ use rule qc_summary_report as bin_sample_report with: input: "results/{project}/output/report/{sample}/{sample}_bin_summary.csv", output: - report( - directory("results/{project}/output/report/{sample}/bin/"), - htmlindex="index.html", - category="4. Binning results", - subcategory="4.1 Summary", - labels={"sample": "{sample}"}, + temp( + report( + directory("results/{project}/output/report/{sample}/bin/"), + htmlindex="index.html", + category="4. Binning results", + subcategory="4.1 Summary", + labels={"sample": "{sample}"}, + ) ), - params: - pin_until="bin", - styles="resources/report/tables/", - name="{sample}_bin_summary", - header="Bin summary for sample {sample}", - pattern=config["tablular-config"], log: "logs/{project}/report/{sample}/bin_rbt_csv.log", - - -use rule qc_summary_report as dastool_report with: - input: - "results/{project}/output/report/{sample}/{sample}_DASTool_summary.csv", - output: - report( - directory("results/{project}/output/report/{sample}/dastool/"), - htmlindex="index.html", - category="4. Binning results", - subcategory="4.2 Quality control", - labels={ - "sample": "{sample}", - "tool": "DAS Tool", - }, - ), params: pin_until="bin", styles="resources/report/tables/", - name="{sample}_DASTool_summary", - header="DAS Tool summary for sample {sample}", - pattern=config["tablular-config"], - log: - "logs/{project}/report/{sample}/dastool_rbt_csv.log", - - -use rule qc_summary_report as checkm2_report with: - input: - "results/{project}/output/report/{sample}/{sample}_checkm2_summary.csv", - output: - report( - directory("results/{project}/output/report/{sample}/checkm2/"), - htmlindex="index.html", - category="4. Binning results", - subcategory="4.2 Quality control", - labels={ - "sample": "{sample}", - "tool": "CheckM 2", - }, - ), - params: - pin_until="bin", - styles="resources/report/tables/", - name="{sample}_CheckM2_summary", - header="CheckM2 summary for sample {sample}", + name="{sample}_bin_summary", + header="Bin summary for sample {sample}", pattern=config["tablular-config"], - log: - "logs/{project}/report/{sample}/checkm2_rbt_csv.log", use rule qc_summary_report as taxonomy_report with: input: "results/{project}/output/report/{sample}/{sample}_bin_taxonomy.csv", output: - report( - directory("results/{project}/output/report/{sample}/taxonomy/"), - htmlindex="index.html", - category="4. Binning results", - subcategory="4.3 Taxonomy classification", - labels={"sample": "{sample}"}, + temp( + report( + directory("results/{project}/output/report/{sample}/taxonomy/"), + htmlindex="index.html", + category="4. Binning results", + subcategory="4.3 Taxonomy classification", + labels={"sample": "{sample}"}, + ) ), + log: + "logs/{project}/report/{sample}/taxonomy_rbt_csv.log", params: pin_until="bin", styles="resources/report/tables/", name="{sample}_taxonomy_summary", header="Taxonomy summary for sample {sample}", pattern=config["tablular-config"], - log: - "logs/{project}/report/{sample}/taxonomy_rbt_csv.log", use rule qc_summary_report as mag_report with: input: "results/{project}/output/report/{sample}/{sample}_mags_summary.csv", output: - report( - directory("results/{project}/output/report/{sample}/mags/"), - htmlindex="index.html", - category="5. Taxonomic classification", - subcategory="5.1 MAGs classification", - labels={"sample": "{sample}"}, + temp( + report( + directory("results/{project}/output/report/{sample}/mags/"), + htmlindex="index.html", + category="5. Taxonomic classification", + subcategory="5.1 MAGs classification", + labels={"sample": "{sample}"}, + ) ), + log: + "logs/{project}/report/{sample}/mag_rbt_csv.log", params: pin_until="MAG", styles="resources/report/tables/", name="{sample}_MAG_summary", header="MAG summary for sample {sample}", pattern=config["tablular-config"], - log: - "logs/{project}/report/{sample}/mag_rbt_csv.log", - - -rule bin_summary_all: - input: - csv_mags=expand( - "results/{{project}}/output/report/{sample}/{sample}_mags_summary.csv", - sample=get_samples(), - ), - csv_bins=expand( - "results/{{project}}/output/report/{sample}/{sample}_bin_summary.csv", - sample=get_samples(), - ), - output: - "results/{project}/output/report/all/binning_summary.csv", - log: - "logs/{project}/bin_summary/all.log", - threads: 4 - conda: - "../envs/python.yaml" - script: - "../scripts/bin_summary_all.py" - - -use rule qc_summary_report as bin_all_report with: - input: - "results/{project}/output/report/all/binning_summary.csv", - output: - report( - directory("results/{project}/output/report/all/binning/"), - htmlindex="index.html", - category="4. Binning results", - subcategory="4.1 Summary", - labels={"sample": "all"}, - ), - params: - pin_until="sample", - styles="resources/report/tables/", - name="bin_summary", - header="Bin summary for all samples", - pattern=config["tablular-config"], - log: - "logs/{project}/report/all_bin_rbt_csv.log", diff --git a/workflow/rules/classify.smk b/workflow/rules/classify.smk index 8c29097..c785336 100644 --- a/workflow/rules/classify.smk +++ b/workflow/rules/classify.smk @@ -1,20 +1,24 @@ -rule download_kaiju: - output: - db_files=get_kaiju_files(), - params: - download=config["kaiju"]["download"], - db_folder=lambda wildcards, output: Path(output.db_files[0]).parent, - log: - "logs/kaiju_DB_download.log", - conda: - "../envs/unix.yaml" - shell: - "(mkdir -p {params.db_folder} && " - "wget -c {params.download} -O - | " - "tar -zxv -C {params.db_folder}) > {log} 2>&1" +# if there is no local database version to use, it is downloaded +if not config["kaiju"]["use-shared"]: + + rule download_kaiju: + output: + db_files=get_kaiju_files(), + params: + download=config["kaiju"]["download"], + db_folder=lambda wildcards, output: Path(output.db_files[0]).parent, + log: + "logs/kaiju_DB_download.log", + conda: + "../envs/unix.yaml" + shell: + "(mkdir -p {params.db_folder} && " + "wget -c {params.download} -O - | " + "tar -zxv -C {params.db_folder}) > {log} 2>&1" -rule run_kaiju: +# classification of reads +rule run_kaiju_reads: input: db_files=get_kaiju_files(), fastqs=get_filtered_gz_fastqs, @@ -24,7 +28,7 @@ rule run_kaiju: evalue=0.00001, threads: 64 log: - "logs/{project}/kaiju/{sample}_run.log", + "logs/{project}/kaiju/{sample}_reads_run.log", conda: "../envs/kaiju.yaml" shell: @@ -33,15 +37,15 @@ rule run_kaiju: "-j {input.fastqs[1]} -o {output.report} > {log} 2>&1" -rule kaiju_table: +rule kaiju_table_reads: input: - report=rules.run_kaiju.output.report, + report=rules.run_kaiju_reads.output.report, db_files=get_kaiju_files(), output: - summary="results/{project}/output/classification/reads/{sample}/{sample}_{level}_kaiju_summary.tsv", + summary="results/{project}/output/classification/reads/{sample}/{sample}_{level}_abundance.tsv", threads: 2 log: - "logs/{project}/kaiju/{sample}_{level}_table.log", + "logs/{project}/kaiju/{sample}_reads_{level}_table.log", conda: "../envs/kaiju.yaml" shell: @@ -51,23 +55,23 @@ rule kaiju_table: # -p to print full taxonomy -rule kaiju2krona: +rule kaiju2krona_reads: input: db_files=get_kaiju_files(), - report=rules.run_kaiju.output.report, + report=rules.run_kaiju_reads.output.report, output: krona=temp( "results/{project}/output/classification/reads/{sample}/kaiju.out.krona" ), html=report( - "results/{project}/output/report/{sample}/{sample}_kaiju.out.html", + "results/{project}/output/report/{sample}/{sample}_reads_kaiju.out.html", htmlindex="index.html", category="5. Taxonomic classification", subcategory="5.2 Read classification", labels={"sample": "{sample}"}, ), log: - "logs/{project}/kaiju/{sample}_2krona.log", + "logs/{project}/kaiju/{sample}_reads_2krona.log", conda: "../envs/kaiju.yaml" shell: @@ -76,86 +80,144 @@ rule kaiju2krona: "ktImportText -o {output.html} {output.krona}) > {log} 2>&1" -if config["gtdb"]["use-local"]: +# classification of contigs +rule run_kaiju_contigs: + input: + db_files=get_kaiju_files(), + assembly=get_gz_assembly, + output: + report=temp( + "results/{project}/output/classification/assembly/{sample}/kaiju.out" + ), + threads: 60 + log: + "logs/{project}/contig_classification/{sample}_kaiju_run.log", + conda: + "../envs/kaiju.yaml" + shell: + "kaiju -z {threads} -t {input.db_files[0]} -f {input.db_files[1]} " + "-i {input.assembly} -o {output.report} > {log} 2>&1" - rule prepare_gtdb: - output: - done=touch("results/GTDB_prep.done"), - params: - db_folder=get_gtdb_folder(), - threads: 1 - log: - "logs/GTDB_prep.log", - conda: - "../envs/gtdbtk.yaml" - shell: - "(conda env config vars set " - "GTDBTK_DATA_PATH='{params.db_folder}') > {log} 2>&1" - rule gtdbtk_classify_wf: - input: - bins=rules.gzip_bins.output.bins, - db_prep=rules.prepare_gtdb.output.done, - output: - json="results/{project}/output/classification/bins/{sample}/gtdbtk.json", - outdir=temp( - directory( - "results/{project}/output/classification/bins/{sample}/gtdbtk/" - ) - ), - params: - db_folder=get_gtdb_folder(), - clf_outdir=lambda wildcards, output: Path(output.json).parent, - json=lambda wildcards, output: Path(output.json).name, - threads: 64 - log: - "logs/{project}/gtdbtk/{sample}_classify.log", - conda: - "../envs/gtdbtk.yaml" - shell: - "(gtdbtk classify_wf --prefix {wildcards.sample} -x fa.gz " - "--mash_db {params.db_folder}mash_db/ --cpus 40 " - "--genome_dir {input.bins}/ --out_dir {output.outdir}/ && " - "cp {output.outdir}/{params.json} {params.clf_outdir}/) > {log} 2>&1" +rule gzip_kaiju_contigs: + input: + report=rules.run_kaiju_contigs.output.report, + output: + report="results/{project}/output/classification/assembly/{sample}/kaiju.out.gz", + threads: 20 + priority: 1 + log: + "logs/{project}/contig_classification/{sample}_kaiju_gzip.log", + conda: + "../envs/unix.yaml" + shell: + "gzip -k {input.report} > {log} 2>&1" -else: - rule download_GTDB: - output: - done=touch("results/gtdbtk_download_DB.done"), - log: - "logs/gtdbtk_download_DB.log", - threads: 64 - conda: - "../envs/gtdbtk.yaml" - shell: - "download-db.sh > {log} 2>&1" +use rule kaiju_table_reads as kaiju_table_contigs with: + input: + report=rules.run_kaiju_contigs.output.report, + db_files=get_kaiju_files(), + output: + summary="results/{project}/output/classification/assembly/{sample}/{sample}_{level}_abundance.tsv", + threads: 2 + log: + "logs/{project}/contig_classification/{sample}_{level}_table.log", - rule gtdbtk_classify_wf: - input: - bins=rules.gzip_bins.output.bins, - load=rules.download_GTDB.output.done, - output: - json="results/{project}/output/classification/bins/{sample}/gtdbtk.json", - outdir=temp( - directory( - "results/{project}/output/classification/bins/{sample}/gtdbtk/" - ) - ), - params: - clf_outdir=lambda wildcards, output: Path(output.json).parent, - json=lambda wildcards, output: Path(output.json).name, - threads: 64 - log: - "logs/{project}/gtdbtk/{sample}_classify.log", - conda: - "../envs/gtdbtk.yaml" - shell: - "(gtdbtk classify_wf --prefix {wildcards.sample} -x fa.gz " - "--cpus {threads} --genome_dir {input.bins}/ --out_dir {output.outdir}/ && " - "cp {output.outdir}/{params.json} {params.clf_outdir}/) > {log} 2>&1" + +use rule kaiju2krona_reads as kaiju2krona_contigs with: + input: + report=rules.run_kaiju_contigs.output.report, + db_files=get_kaiju_files(), + output: + krona=temp( + "results/{project}/output/classification/assembly/{sample}/kaiju.out.krona" + ), + html=report( + "results/{project}/output/report/{sample}/{sample}_contigs_kaiju.out.html", + htmlindex="index.html", + category="5. Taxonomic classification", + subcategory="5.3 Contig classification", + labels={"sample": "{sample}"}, + ), + log: + "logs/{project}/contig_classification/{sample}_2krona.log", + + +rule extract_contig_headers: + input: + assembly=get_gz_assembly, + output: + headers="results/{project}/output/classification/assembly/{sample}/headers.txt", + threads: 10 + log: + "logs/{project}/contig_classification/{sample}_header.log", + conda: + "../envs/unix.yaml" + shell: + "(zgrep '>' {input.assembly} > {output.headers}) > {log} 2>&1" + + +rule contig_classification: + input: + out=rules.run_kaiju_contigs.output.report, + genus="results/{project}/output/classification/assembly/{sample}/{sample}_genus_abundance.tsv", + species="results/{project}/output/classification/assembly/{sample}/{sample}_species_abundance.tsv", + headers=rules.extract_contig_headers.output.headers, + output: + csv="results/{project}/output/classification/assembly/{sample}/{sample}_classified_contigs.csv", + log: + "logs/{project}/contig_classification/{sample}_classified_contigs.log", + threads: 5 + conda: + "../envs/python.yaml" + script: + "../scripts/contig_classification.py" + + +rule prepare_gtdb: + output: + done=temp(touch("results/GTDB_prep.done")), + params: + db_folder=get_gtdb_folder(), + threads: 1 + log: + "logs/GTDB_prep.log", + conda: + "../envs/gtdbtk.yaml" + shell: + "(conda env config vars set " + "GTDBTK_DATA_PATH='{params.db_folder}') > {log} 2>&1" + + +rule gtdbtk_classify_wf: + input: + bins=rules.gzip_bins.output.bins, + db_prep=rules.prepare_gtdb.output.done, + output: + json="results/{project}/output/classification/bins/{sample}/gtdbtk.json", + outdir=temp( + directory("results/{project}/output/classification/bins/{sample}/gtdbtk/") + ), + params: + clf_outdir=lambda wildcards, output: Path(output.json).parent, + json=lambda wildcards, output: Path(output.json).name, + threads: 40 + # to ensure only one of these commands are run at the same time + resources: + heavy=1, + log: + "logs/{project}/gtdbtk/{sample}_classify.log", + conda: + "../envs/gtdbtk.yaml" + shell: + "(gtdbtk classify_wf --prefix {wildcards.sample} -x fa.gz " + "--cpus {threads} --pplacer_cpus {threads} " + "--genome_dir {input.bins}/ --out_dir {output.outdir}/ && " + "cp {output.outdir}/{params.json} {params.clf_outdir}/) > {log} 2>&1" +# combines bacterial and archaeal rule gtdb_summary: input: json=rules.gtdbtk_classify_wf.output.json, diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index abed525..7169886 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -4,14 +4,6 @@ import os configfile: "config/config.yaml" -def get_data_path(): - return config["data-handling"]["data"] - - -def get_resource_path(): - return config["data-handling"]["resources"] - - def get_project(): return config["project-name"] @@ -22,17 +14,10 @@ def get_samples(): def get_fastqs(wildcards): file_r1 = pep.sample_table.loc[wildcards.sample]["fq1"] - folder = str(Path(file_r1).parent) - filename_r1 = Path(file_r1).name - filename_r2 = Path(pep.sample_table.loc[wildcards.sample]["fq2"]).name - return [folder, filename_r1, filename_r2] - - -def get_local_fastqs(wildcards): - path = get_data_path() + file_r2 = pep.sample_table.loc[wildcards.sample]["fq2"] return ( - "{data}{{project}}/{{sample}}_R1.fastq.gz".format(data=path), - "{data}{{project}}/{{sample}}_R2.fastq.gz".format(data=path), + file_r1, + file_r2, ) @@ -60,7 +45,7 @@ def get_prefiltered_fastqs(wildcards): def get_host_map_statistics(wildcards): if config["host-filtering"]["do-host-filtering"]: logs = expand( - "results/{{project}}/report_prerequisites/qc/filter_host_{sample}.log", + "results/{{project}}/output/report/prerequisites/qc/filter_host_{sample}.log", sample=get_samples(), ) return logs @@ -68,38 +53,17 @@ def get_host_map_statistics(wildcards): return [] -def get_human_ref(): - if config["human-filtering"]["use-local"]: - path = config["human-filtering"]["local-path"] - else: - path = config["human-filtering"]["download-path"] - filename = path.split("/")[-1] - local_ref = "".join([get_resource_path(), "ref_genome/", filename]) - return local_ref - - -def get_human_local_folder(): - path = config["human-filtering"]["local-path"] - folder = Path(path).parent - return folder - - -def get_checkm2_db(): - file = "{}{}".format(get_resource_path(), config["checkm2"]) - return file - - def get_filtered_fastqs(wildcards): return [ - "results/{project}/filtered/fastqs/{sample}_R1.fastq", - "results/{project}/filtered/fastqs/{sample}_R2.fastq", + "results/{project}/output/filtered_reads/{sample}_R1.fastq", + "results/{project}/output/filtered_reads/{sample}_R2.fastq", ] def get_filtered_gz_fastqs(wildcards): return [ - "results/{project}/filtered/fastqs/{sample}_R1.fastq.gz", - "results/{project}/filtered/fastqs/{sample}_R2.fastq.gz", + "results/{project}/output/filtered_reads/{sample}_R1.fastq.gz", + "results/{project}/output/filtered_reads/{sample}_R2.fastq.gz", ] @@ -111,8 +75,45 @@ def get_gz_assembly(wildcards): return "results/{project}/output/fastas/{sample}/{sample}.fa.gz" +def get_proteins(wildcards): + return ("results/{project}/output/proteins/{sample}/{sample}_proteins.faa.gz",) + + +## binning parameters +def get_contig_length_threshold(): + return config["min-contig-length"] + + +def bins_for_sample(wildcards): + binfolder = "results/{project}/binning/das_tool/{sample}/{sample}_DASTool_bins/" + bin_list = [b for b in os.listdir(binfolder) if b.endswith(".fa")] + if len(bin_list) > 0: + return True + else: + return False + + +def get_human_ref_download(): + return config["human-filtering"]["download-path"] + + +def get_human_ref(): + if config["human-filtering"]["use-shared"]: + path = config["human-filtering"]["local-path"] + return path + else: + path = get_human_ref_download() + filename = path.split("/")[-1] + local_ref = "".join(["resources/ref_genome/", filename]) + return local_ref + + +def get_host_ref(): + return config["host-filtering"]["ref-genome"] + + def get_kaiju_files(): - file = "".join([get_resource_path(), config["kaiju"]["fmi-file"]]) + file = config["kaiju"]["fmi-file"] path = str(Path(file).parent) fmi = Path(file).name names = ["nodes.dmp", fmi, "names.dmp"] @@ -120,49 +121,59 @@ def get_kaiju_files(): return files -## binning parameters -def get_contig_length_threshold(): - return config["min-contig-length"] +def get_checkm2_db(): + path = "".join( + [config["checkm2"]["db-folder"], "CheckM2_database/uniref100.KO.1.dmnd"] + ) + return path -def get_binners(): - return config["das-tool"]["binner-list"] +def get_gtdb_folder(): + return config["gtdb"] -def get_all_contig2bin_files(wildcards): - binners = get_binners() - file_list = [ - "".join( - [ - "results/{project}/output/contig2bins/{sample}/", - binner, - "_contig2bin.tsv", - ] - ) - for binner in binners - ] - return file_list +def get_genomad_DB_folder(): + path = "".join([config["genomad"]["db-folder"], "genomad_db/"]) + return path -## reads in binner control file and returns list with paths to contig2bin files -## and a list with name of the binners that produced results -def get_paths_binner(wildcards): - file = "results/{}/binning/das_tool/binner_control_{}.csv".format( - wildcards.project, wildcards.sample - ) - lines = open(file).readlines() - paths = str(lines[0].rstrip("\n")) - binner = str(lines[1].rstrip("\n")) - return paths, binner +def get_genomad_DB_file(): + path = "".join([get_genomad_DB_folder(), "names.dmp"]) + return path -def bins_for_sample(wildcards): - if len(get_paths_binner[0]) > 0: - return True - else: - return False +def get_card_db_file(): + return config["card"]["db-file"] + + +def get_card_annotation_file(): + version = config["card"]["version"] + filename = "".join(["card_database_", version, ".fasta"]) + folder = str(Path(get_card_db_file()).parent) + + path = "/".join([folder, filename]) + return path + + +def get_uniCARD_db(): + return config["uniCARD"]["db-file"] + + +def get_CARD_hierarchy(): + return config["uniCARD"]["hierarchy-json"] + + +def get_uniCARD_db_wo_ext(): + db_wo_ext = str(Path(get_uniCARD_db()).with_suffix("").with_suffix("")) + return db_wo_ext + + +def get_unicard_dmnd(): + return "".join([get_uniCARD_db_wo_ext(), ".dmnd"]) +""" + def get_mag_fa(wildcards): folder = "results/{}/output/fastas/{}/mags/".format( wildcards.project, wildcards.sample @@ -175,14 +186,6 @@ def get_mag_fa(wildcards): return files -def get_binIDs_for_sample(wildcards): - folder = "results/{}/output/fastas/{}/mags/".format( - wildcards.project, wildcards.sample - ) - binIDs = [binID for binID in os.listdir(folder) if binID.endswith("fa.gz")] - return binIDs - - def get_mag_ARGs(wildcards): bin_fastas = (get_mag_fa(wildcards),) bin_fastas = list(bin_fastas)[0] @@ -194,24 +197,4 @@ def get_mag_ARGs(wildcards): os.path.join(folder, binID.replace(".fa.gz", ".txt")) for binID in binIDs ] return arg_files - - -def get_gtdb_folder(): - path = "".join([get_resource_path(), config["gtdb"]["db-folder"]]) - return path - - -def get_genomad_DB_file(): - path = "".join([get_resource_path(), "genomad_db/names.dmp"]) - return path - - -def get_card_db_file(): - path = "".join([get_resource_path(), "CARD_db/", config["card"]["dbfile"]]) - return path - - -def get_card_annotation_file(): - version = config["card"]["version"] - path = "".join([get_resource_path(), "CARD_db/card_database_", version, ".fasta"]) - return path +""" diff --git a/workflow/rules/das_tool.smk b/workflow/rules/das_tool.smk index aa1f5ae..81238bd 100644 --- a/workflow/rules/das_tool.smk +++ b/workflow/rules/das_tool.smk @@ -5,7 +5,7 @@ rule postprocess_metabat: input: outdir=rules.metabat.output.outdir, output: - "results/{project}/output/contig2bins/{sample}/metabat_contig2bin.tsv", + tsv="results/{project}/output/contig2bins/{sample}/metabat_contig2bin.tsv", params: binner="metabat", prefix="bin", @@ -21,9 +21,8 @@ rule postprocess_metabat: rule postprocess_metacoag: input: c2bin="results/{project}/binning/metacoag/{sample}/contig_to_bin.tsv", - folder="results/{project}/binning/metacoag/{sample}/", output: - "results/{project}/output/contig2bins/{sample}/metacoag_contig2bin.tsv", + tsv="results/{project}/output/contig2bins/{sample}/metacoag_contig2bin.tsv", threads: 12 log: "logs/{project}/contig2bins/{sample}/postprocess_metacoag.log", @@ -31,30 +30,14 @@ rule postprocess_metacoag: "../envs/unix.yaml" shell: "(awk '{{print $1 \"\t\" $NF}}' {input.c2bin} | " - "sed 's/len=[0-9]*,/metacoag_/g' > {output}) > {log} 2>&1" - - -## tests if binner created bins and writes file for DAS Tool -rule binner_control: - input: - get_all_contig2bin_files, - output: - temp("results/{project}/binning/das_tool/binner_control_{sample}.csv"), - params: - get_binners(), - log: - "logs/{project}/das_tool/{sample}/binner_control.log", - conda: - "../envs/unix.yaml" - script: - "../scripts/binner_control.py" + "sed 's/len=[0-9]*,/metacoag_/g' > {output.tsv}) > {log} 2>&1" rule dastool_run: input: - binc=rules.binner_control.output, + metacoag=rules.postprocess_metacoag.output.tsv, + metabat=rules.postprocess_metabat.output.tsv, contigs=get_assembly, - asml_folder=rules.megahit.output.outdir, output: summary=temp( "results/{project}/binning/das_tool/{sample}/{sample}_DASTool_summary.tsv" @@ -67,22 +50,22 @@ rule dastool_run: "results/{project}/binning/das_tool/{sample}/{sample}_DASTool_bins/" ) ), - outdir=temp(directory("results/{project}/binning/das_tool/{sample}/")), - done=touch("results/{project}/binning/das_tool/{sample}_run.done"), params: - path_bin_list=get_paths_binner, + outdir=lambda wildcards, output: Path(output.bins).parent, threshold=0.001, threads: 64 + resources: + heavy=2, log: "logs/{project}/das_tool/{sample}/das_tool_run.log", conda: "../envs/das_tool.yaml" shell: "DAS_Tool --debug --write_bins " - "-i {params.path_bin_list[0]} " - "-l {params.path_bin_list[1]} " + "-i {input.metabat},{input.metacoag} " + "-l metabat,metacoag " "-c {input.contigs} " - "-o {output.outdir}/{wildcards.sample} " + "-o {params.outdir}/{wildcards.sample} " "--score_threshold {params.threshold} " "--threads={threads} " "> {log} 2>&1 " @@ -113,13 +96,14 @@ if bins_for_sample: output: bins=directory("results/{project}/output/fastas/{sample}/bins/"), done=touch("results/{project}/binning/das_tool/{sample}_bins.done"), - threads: 64 + threads: 10 + priority: 1 log: "logs/{project}/bins/{sample}/gz_bins.log", conda: "../envs/unix.yaml" shell: - "(pigz -k {input.bins}/*.fa && " + "(gzip -k {input.bins}/*.fa && " "mkdir -p {output.bins}/ && " "mv {input.bins}/*.fa.gz {output.bins}/ ) > {log} 2>&1" @@ -136,14 +120,3 @@ if bins_for_sample: "../envs/python.yaml" script: "../scripts/move_MAGs.py" - - rule cleanup_dastool_output: - input: - folder=rules.dastool_run.output.outdir, - bins=rules.gzip_bins.output.done, - move=rules.move_dastool_output.output.done, - output: - done=touch("results/{project}/binning/das_tool/{sample}_cleanup.done"), - threads: 2 - log: - "logs/{project}/das_tool/{sample}/cleanup.log", diff --git a/workflow/rules/host_filtering.smk b/workflow/rules/host_filtering.smk index 238d47f..e896199 100644 --- a/workflow/rules/host_filtering.smk +++ b/workflow/rules/host_filtering.smk @@ -1,40 +1,20 @@ from pathlib import Path - -if config["human-filtering"]["use-local"]: - - rule copy_local_human_ref: - output: - fasta=get_human_ref(), - params: - local=get_human_local_folder(), - folder=lambda wildcards, output: Path(output.fasta).parent, - file=lambda wildcards, output: Path(output.fasta).name, - log: - "logs/human_ref_local_copy.log", - group: - "refGenome_depended" - conda: - "../envs/unix.yaml" - shell: - "(mkdir -p {params.folder} && " - "tar cpfz - -C {params.local} {params.file} | " - "(cd {params.folder} ; tar xpfz -)) > {log} 2>&1" - -else: +# Download human reference genome if not local file is given +if not config["human-filtering"]["use-shared"]: rule download_human_ref: output: fasta=get_human_ref(), - params: - download=config["human-filtering"]["download-path"], - folder=lambda wildcards, output: Path(output.fasta).parent, log: "logs/human_ref_download.log", group: "refGenome_depended" conda: "../envs/unix.yaml" + params: + download=get_human_ref_download(), + folder=lambda wildcards, output: Path(output.fasta).parent, shell: "(mkdir -p {params.folder} && " "cd {params.folder} && " @@ -47,11 +27,11 @@ rule map_to_human: ref=get_human_ref(), output: bam=temp("results/{project}/human_filtering/alignments/{sample}.bam"), - threads: 64 log: "logs/{project}/human_filtering/map_to_human_{sample}.log", conda: "../envs/minimap2.yaml" + threads: 64 shell: "(minimap2 -a -xsr -t {threads} {input.ref} {input.fastqs} | " "samtools view -bh | " @@ -63,11 +43,11 @@ rule index_human_alignment: rules.map_to_human.output.bam, output: bai=temp("results/{project}/human_filtering/alignments/{sample}.bam.bai"), - threads: 20 log: "logs/{project}/human_filtering/index_human_alignment_{sample}.log", conda: "../envs/minimap2.yaml" + threads: 20 shell: "samtools index {input} > {log} 2>&1" @@ -79,15 +59,15 @@ rule filter_human: output: filtered=temp( expand( - "results/{{project}}/filtered/fastqs/{{sample}}_{read}.fastq", + "results/{{project}}/output/filtered_reads/{{sample}}_{read}.fastq", read=["R1", "R2"], ) ), - threads: 64 log: - "results/{project}/report_prerequisites/qc/filter_human_{sample}.log", + "results/{project}/output/report/prerequisites/qc/{sample}_filter_human.log", conda: "../envs/minimap2.yaml" + threads: 64 shell: "(samtools fastq --threads {threads} -F 3584 -f 77 " "-o {output.filtered[0]} {input.bam} && " @@ -97,16 +77,17 @@ rule filter_human: rule gzip_filtered_reads: input: - "results/{project}/filtered/fastqs/{sample}_{read}.fastq", + "results/{project}/output/filtered_reads/{sample}_{read}.fastq", output: - "results/{project}/filtered/fastqs/{sample}_{read}.fastq.gz", + "results/{project}/output/filtered_reads/{sample}_{read}.fastq.gz", log: "logs/{project}/human_filtering/gzip_{sample}_{read}.log", - threads: 64 + priority: 1 conda: "../envs/unix.yaml" + threads: 20 shell: - "pigz -k {input} > {log} 2>&1" + "gzip -k {input} > {log} 2>&1" if config["host-filtering"]["do-host-filtering"]: @@ -116,23 +97,21 @@ if config["host-filtering"]["do-host-filtering"]: use rule map_to_human as map_to_host with: input: fastqs=get_trimmed_fastqs, - ref=config["host-filtering"]["ref-genome"], + ref=get_host_ref(), output: bam=temp("results/{project}/host_filtering/alignments/{sample}.bam"), - params: - ref=config["host-filtering"]["ref-genome"], - threads: 20 log: "logs/{project}/host_filtering/map_to_host_{sample}.log", + threads: 20 use rule index_human_alignment as index_host_alignment with: input: rules.map_to_host.output.bam, output: bai=temp("results/{project}/host_filtering/alignments/{sample}.bam.bai"), - threads: 3 log: "logs/{project}/host_filtering/index_host_alignment_{sample}.log", + threads: 3 rule filter_host: input: @@ -145,90 +124,13 @@ if config["host-filtering"]["do-host-filtering"]: read=["R1", "R2"], ) ), - threads: 64 log: - "results/{project}/report_prerequisites/qc/filter_host_{sample}.log", + "results/{project}/output/report/prerequisites/qc/{sample}_filter_host.log", conda: "../envs/minimap2.yaml" + threads: 64 shell: "(samtools fastq -F 3584 -f 77 {input.bam} | " - "pigz -c > {output.filtered[0]} && " + "gzip -c > {output.filtered[0]} && " "samtools fastq -F 3584 -f 141 {input.bam} | " - "pigz -c > {output.filtered[1]}) > {log} 2>&1" - - -## TODO -## change to using kaiju output or just to present human contamination -## all reports are done on human (+ optional other different host) filtered -"""rule diversity_summary: - input: - reports=expand( - "results/{{project}}/output/classification/reads/{sample}/{sample}_kraken2_report.tsv", - sample=get_samples(), - ), - jsons=expand( - "results/{{project}}/report_prerequisites/qc/{sample}.fastp.json", - sample=get_samples(), - ), - human_logs=expand( - "results/{{project}}/report_prerequisites/qc/filter_human_{sample}.log", - sample=get_samples(), - ), - host_logs=get_host_map_statistics, - output: - csv="results/{project}/output/report/all/diversity_summary.csv", - log: - "logs/{project}/kraken2/summary.log", - params: - other_host=config["host-filtering"]["do-host-filtering"], - hostname=config["host-filtering"]["host-name"], - threads: 2 - conda: - "../envs/python.yaml" - script: - "../scripts/diversity_summary.py" - - -use rule qc_summary_report as diversity_summary_report with: - input: - rules.diversity_summary.output.csv, - output: - report( - directory("results/{project}/output/report/all/diversity_summary/"), - htmlindex="index.html", - caption="../report/kraken.rst", - category="2. Species diversity", - labels={ - "sample": "all samples", - }, - ), - params: - pin_until="sample", - styles="resources/report/tables/", - name="diversity_summary", - header="Diversity summary based on mapping to host genome(s) and Kraken2", - pattern=config["tablular-config"], - log: - "logs/{project}/report/kraken2_rbt_csv.log", - - -rule create_host_plot: - input: - csv=rules.diversity_summary.output.csv, - output: - html=report( - "results/{project}/output/report/all/host_contamination.html", - caption="../report/host_plot.rst", - category="2. Species diversity", - labels={"sample": "all samples"}, - ), - params: - other_host=config["host-filtering"]["do-host-filtering"], - hostname=config["host-filtering"]["host-name"], - log: - "logs/{project}/report/host_plot.log", - conda: - "../envs/python.yaml" - script: - "../scripts/plot_host.py" -""" + "gzip -c > {output.filtered[1]}) > {log} 2>&1" diff --git a/workflow/rules/metabat.smk b/workflow/rules/metabat.smk index 7d8328c..44e467b 100644 --- a/workflow/rules/metabat.smk +++ b/workflow/rules/metabat.smk @@ -36,14 +36,3 @@ rule metabat: "-m 1500 -s 100000 " "--minCorr 95 --minContigByCorr 300 --minSamples 1 " "-t {threads} -v -l -o {output}/{params.prefix} > {log} 2>&1" - - -rule cleanup_metabat_output: - input: - folder=rules.metabat.output.outdir, - dastool="results/{project}/binning/das_tool/{sample}/{sample}_DASTool_summary.tsv", - output: - done=touch("results/{project}/binning/metabat/{sample}_cleanup.done"), - threads: 2 - log: - "logs/{project}/metabat/{sample}/cleanup.log", diff --git a/workflow/rules/metacoag.smk b/workflow/rules/metacoag.smk index 49319a4..b5fc240 100644 --- a/workflow/rules/metacoag.smk +++ b/workflow/rules/metacoag.smk @@ -63,12 +63,10 @@ rule metacoag_run: contigs=get_assembly, gfa="results/{project}/binning_prep/{sample}/assembly_tree.gfa", abd="results/{project}/binning_prep/{sample}/abundance_metacoag.tsv", + #assembly folder needs to be there + folder=rules.megahit.output.outdir, output: out_tsv=temp("results/{project}/binning/metacoag/{sample}/contig_to_bin.tsv"), - folder=temp(directory("results/{project}/binning/metacoag/{sample}/")), - intermediate=temp( - "results/{project}/binning/metacoag/{sample}final.contigs.fa.normalized_contig_tetramers.pickle" - ), params: outdir=lambda wildcards, output: Path(output.out_tsv).parent, threads: 64 @@ -82,14 +80,3 @@ rule metacoag_run: "--output {params.outdir} --min_length 300 " "--bin_mg_threshold 0.2 --min_bin_size 100000 " "--nthreads {threads} > {log} 2>&1" - - -rule cleanup_metacoag_output: - input: - folder=rules.metacoag_run.output.folder, - dastool="results/{project}/binning/das_tool/{sample}/{sample}_DASTool_summary.tsv", - output: - done=touch("results/{project}/binning/metacoag/{sample}_cleanup.done"), - threads: 2 - log: - "logs/{project}/metacoag/{sample}/cleanup.log", diff --git a/workflow/rules/qc.smk b/workflow/rules/qc.smk index b00345e..e377b4b 100644 --- a/workflow/rules/qc.smk +++ b/workflow/rules/qc.smk @@ -1,68 +1,104 @@ ## read QC +# fastp in paired-end mode for Illumina paired-end data +# version in this wrapper: fastp=1.0.1 +rule fastp: + input: + sample=get_fastqs, + output: + trimmed=temp( + [ + "results/{project}/trimmed/fastp/{sample}.1.fastq.gz", + "results/{project}/trimmed/fastp/{sample}.2.fastq.gz", + ] + ), + html=temp("results/{project}/trimmed/fastp/{sample}.html"), + json="results/{project}/output/report/prerequisites/qc/{sample}.fastp.json", + log: + "logs/{project}/fastp/{sample}.log", + threads: 16 + params: + adapters=get_adapters, + extra="--qualified_quality_phred {phred} --length_required {minlen}".format( + phred=(config["quality-criteria"]["min-PHRED"]), + minlen=(config["quality-criteria"]["min-length-reads"]), + ), + wrapper: + "v7.1.0/bio/fastp" + + rule fastqc: input: - get_trimmed_fastqs, + rules.fastp.output.trimmed, output: - html=temp("results/{project}/qc/fastqc/{sample}_trimmed.html"), - zip=temp("results/{project}/qc/fastqc/{sample}_trimmed_fastqc.zip"), + zip1=temp( + "results/{project}/trimmed/fastp/{sample}.1_fastqc.zip", + ), + zip2=temp( + "results/{project}/trimmed/fastp/{sample}.2_fastqc.zip", + ), log: "logs/{project}/fastqc/{sample}.log", - wrapper: - "v1.23.5/bio/fastqc" + conda: + "../envs/fastqc.yaml" + threads: 4 + resources: + mem_mb=1024, + shell: + "fastqc --memory {resources.mem_mb} --threads {threads} " + "--format fastq --quiet {input} > {log} 2>&1" +# version in this wrapper: multiqc=1.33 rule multiqc: input: expand( [ - "results/{{project}}/qc/fastqc/{sample}_trimmed_fastqc.zip", - "results/{{project}}/report_prerequisites/qc/{sample}.fastp.json", + "results/{{project}}/trimmed/fastp/{sample}.{read}_fastqc.zip", + "results/{{project}}/output/report/prerequisites/qc/{sample}.fastp.json", ], sample=get_samples(), + read=["1", "2"], ), + config="config/multiqc_config.yaml", output: report( - "results/{project}/output/report/all/multiqc.html", - htmlindex="multiqc.html", + "results/{project}/output/report/all/multiqc_{project}.html", + htmlindex="multiqc_{project}.html", category="1. Quality control", labels={"sample": "all samples"}, ), - "results/{project}/output/report/all/multiqc_data.zip", - params: - extra=( - "--zip-data-dir " - "--config config/multiqc_config.yaml " - "--title 'Results for data from {project} project'" - ), - use_input_files_only=True, + "results/{project}/output/report/all/multiqc_{project}.zip", log: "logs/{project}/multiqc.log", + params: + extra=("--title 'Results for data from {project} project'"), + use_input_files_only=True, wrapper: - "v3.3.1/bio/multiqc" + "v8.1.1/bio/multiqc" rule qc_summary: input: # move these to report_prerequistes jsons=expand( - "results/{{project}}/report_prerequisites/qc/{sample}.fastp.json", + "results/{{project}}/output/report/prerequisites/qc/{sample}.fastp.json", sample=get_samples(), ), human_logs=expand( - "results/{{project}}/report_prerequisites/qc/filter_human_{sample}.log", + "results/{{project}}/output/report/prerequisites/qc/{sample}_filter_human.log", sample=get_samples(), ), host_logs=get_host_map_statistics, output: csv="results/{project}/output/report/all/quality_summary.csv", vis_csv=temp("results/{project}/output/report/all/quality_summary_visual.csv"), - params: - other_host=config["host-filtering"]["do-host-filtering"], - hostname=config["host-filtering"]["host-name"], log: "logs/{project}/report/qc_summary.log", conda: "../envs/python.yaml" + params: + other_host=config["host-filtering"]["do-host-filtering"], + hostname=config["host-filtering"]["host-name"], script: "../scripts/quality_summary.py" @@ -71,24 +107,26 @@ rule qc_summary_report: input: rules.qc_summary.output.vis_csv, output: - report( - directory("results/{project}/output/report/all/quality_summary/"), - htmlindex="index.html", - category="1. Quality control", - labels={ - "sample": "all samples", - }, + temp( + report( + directory("results/{project}/output/report/all/quality_summary/"), + htmlindex="index.html", + category="1. Quality control", + labels={ + "sample": "all samples", + }, + ) ), + log: + "logs/{project}/report/qc_summary_rbt_csv.log", + conda: + "../envs/rbt.yaml" params: pin_until="sample", styles="resources/report/tables/", name="quality_summary", header="Quality summary based on fastp report and mapping to host genome(s)", pattern=config["tablular-config"], - log: - "logs/{project}/report/qc_summary_rbt_csv.log", - conda: - "../envs/rbt.yaml" shell: "rbt csv-report {input} --pin-until {params.pin_until} {output} && " "(sed -i '{params.pattern} {params.header}' " diff --git a/workflow/rules/report.smk b/workflow/rules/report.smk index 149cbc1..421644f 100644 --- a/workflow/rules/report.smk +++ b/workflow/rules/report.smk @@ -7,25 +7,15 @@ rule snakemake_report: input: # 1. Quality control report_input, - "results/{project}/output/report/all/multiqc.html", + "results/{project}/output/report/all/multiqc_{project}.html", # 2. Species diversity "results/{project}/output/report/all/quality_summary/", # 3. Assembly results "results/{project}/output/report/all/assembly/", - # 4. Binning results - "results/{project}/output/report/all/binning/", expand( "results/{{project}}/output/report/{sample}/bin/", sample=get_samples(), ), - expand( - "results/{{project}}/output/report/{sample}/checkm2/", - sample=get_samples(), - ), - expand( - "results/{{project}}/output/report/{sample}/dastool/", - sample=get_samples(), - ), expand( "results/{{project}}/output/report/{sample}/taxonomy/", sample=get_samples(), @@ -36,7 +26,11 @@ rule snakemake_report: sample=get_samples(), ), expand( - "results/{{project}}/output/report/{sample}/{sample}_kaiju.out.html", + "results/{{project}}/output/report/{sample}/{sample}_reads_kaiju.out.html", + sample=get_samples(), + ), + expand( + "results/{{project}}/output/report/{sample}/{sample}_contigs_kaiju.out.html", sample=get_samples(), ), output: @@ -50,6 +44,5 @@ rule snakemake_report: "../envs/snakemake.yaml" shell: "snakemake --nolock --report {output} --report-stylesheet {params.style} " - "> {log} 2>&1" #"{params.for_testing} " diff --git a/workflow/rules/resistance.smk b/workflow/rules/resistance.smk new file mode 100644 index 0000000..08dfcef --- /dev/null +++ b/workflow/rules/resistance.smk @@ -0,0 +1,272 @@ +# Resistance analysis +if not config["card"]["use-shared"]: + + rule download_CARD_data: + output: + json=get_card_db_file(), + log: + "logs/CARD_data_download.log", + conda: + "../envs/unix.yaml" + threads: 30 + params: + download=config["card"]["url"], + folder=lambda wildcards, output: Path(output.json).parent, + filename=lambda wildcards, output: Path(output.json).name, + shell: + "(cd {params.folder} && " + "wget {params.download} && " + "tar -xvf data ./{params.filename}) > {log} 2>&1" + + +rule CARD_load_DB_for_reads: + input: + get_card_db_file(), + output: + temp(touch("results/CARD_load_DB_for_reads.done")), + log: + "logs/CARD_load_DB_for_reads.log", + conda: + "../envs/card.yaml" + shell: + "rgi clean --local && " "rgi load --card_json {input} --local > {log} 2>&1" + + +rule CARD_annotation: + input: + json=get_card_db_file(), + load=rules.CARD_load_DB_for_reads.output, + output: + done=temp(touch("results/CARD_annotation.done")), + ann=get_card_annotation_file(), + log: + "logs/CARD_annotation.log", + conda: + "../envs/card.yaml" + params: + folder=lambda wildcards, input: Path(input.json).parent, + file=lambda wildcards, input: Path(input.json).name, + shell: + "(cd {params.folder}/ && " + "rgi card_annotation -i {params.file}) && " + "rgi load -i {input.json} --card_annotation {output.ann} --local > {log} 2>&1" + + +rule CARD_read_run: + input: + fa=get_filtered_gz_fastqs, + db=rules.CARD_annotation.output, + output: + txt="results/{project}/output/resistance/CARD/reads/{sample}/{sample}.gene_mapping_data.txt", + log: + "logs/{project}/ARGs/reads/{sample}.log", + conda: + "../envs/card.yaml" + threads: 64 + params: + folder=lambda wildcards, output: Path(output.txt).parent, + shell: + "rgi bwt -1 {input.fa[0]} -2 {input.fa[1]} " + "-o {params.folder}/{wildcards.sample} --local " + "--clean -n {threads} > {log} 2>&1" + + +rule uniCARD_makeDB: + input: + db=get_uniCARD_db(), + output: + dmdb=get_unicard_dmnd(), + log: + "logs/unicard_makeDB.log", + conda: + "../envs/diamond.yaml" + threads: 60 + params: + filename=lambda wildcards, input: Path(input.db).name, + filename_wo_ext=Path(get_uniCARD_db_wo_ext()).name, + folder=lambda wildcards, input: Path(input.db).parent, + shell: + "(cd {params.folder} && diamond makedb --in {params.filename} -d {params.filename_wo_ext}) > {log} 2>&1" + + +rule uniCARD_run: + input: + faa=get_proteins, + dmdb=get_unicard_dmnd(), + output: + tsv=temp( + "results/{project}/output/resistance/uniCARD/assembly/{sample}_out.tsv" + ), + log: + "logs/{project}/uniCARD/{sample}.log", + # ensure priority over potential contig classification + priority: 1 + conda: + "../envs/diamond.yaml" + threads: 60 + # to ensure only one of these commands are run at the same time + resources: + heavy=1, + params: + db_wo_ext=lambda w, input: os.path.splitext(input.dmdb)[0], + max_targets=15, + shell: + "diamond blastp -q {input.faa} -d {params.db_wo_ext} " + "-p {threads} --max-target-seqs {params.max_targets} " + "-o {output.tsv} > {log} 2>&1" + + +rule gz_uniCARD_tsv: + input: + unicard_in=rules.uniCARD_run.output.tsv, + output: + tsv=temp( + "results/{project}/output/resistance/uniCARD/assembly/{sample}_out.tsv.gz" + ), + log: + "logs/{project}/uniCARD/{sample}_gz.log", + conda: + "../envs/unix.yaml" + threads: 16 + shell: + "(gzip -k {input.unicard_in}) > {log} 2>&1" + + +rule filter_uniCARD: + input: + tsv=rules.gz_uniCARD_tsv.output.tsv, + json=get_CARD_hierarchy(), + output: + csv="results/{project}/output/resistance/uniCARD/assembly/{sample}.csv", + log: + "logs/{project}/uniCARD/{sample}_filter.log", + conda: + "../envs/python.yaml" + threads: 20 + params: + # filter window for filtering UniCARD results, default: 1% + filter_window=0.01, + shell: + "python workflow/scripts/uniCARD_filtering.py --infile {input.tsv} " + "--card_hierarchy {input.json} --outfile {output.csv} " + "--filter_window {params.filter_window} > {log} 2>&1" + + +rule resistance_abundance_per_sample: + input: + csv=rules.filter_uniCARD.output.csv, + json=get_CARD_hierarchy(), + output: + abd_class="results/{project}/output/resistance/uniCARD/per_sample_resistance_abundance/{sample}_drug_class_abundance.csv", + abd_ARGs="results/{project}/output/resistance/uniCARD/per_sample_resistance_abundance/{sample}_ARGs_abundance.csv", + log: + "logs/{project}/uniCARD/{sample}_abundance.log", + conda: + "../envs/python.yaml" + threads: 5 + script: + "../scripts/uniCARD_abundance_sample.py" + + +rule resistance_abundance_all: + input: + contig_files=expand( + "results/{{project}}/output/classification/assembly/{sample}/{sample}_classified_contigs.csv", + sample=get_samples(), + ), + abd_class_files=expand( + "results/{{project}}/output/resistance/uniCARD/per_sample_resistance_abundance/{sample}_drug_class_abundance.csv", + sample=get_samples(), + ), + abd_ARGs_files=expand( + "results/{{project}}/output/resistance/uniCARD/per_sample_resistance_abundance/{sample}_ARGs_abundance.csv", + sample=get_samples(), + ), + output: + abd_class_possible="results/{project}/output/resistance/uniCARD/{project}_drug_class_abundance_all_possible.csv", + abd_class_present="results/{project}/output/resistance/uniCARD/{project}_drug_class_abundance_all_present.csv", + abd_ARGs_possible="results/{project}/output/resistance/uniCARD/{project}_ARGs_abundance_all_possible.csv", + abd_ARGs_present="results/{project}/output/resistance/uniCARD/{project}_ARGs_abundance_all_present.csv", + log: + "logs/{project}/uniCARD/all_abundances.log", + conda: + "../envs/python.yaml" + threads: 5 + script: + "../scripts/uniCARD_abundance.py" + + +""" +# updates CARD database to use for contigs instead of reads +# read based classification must be finished before +rule CARD_load_DB: + input: + db=get_card_db_file(), + read_args=expand( + "results/{project}/output/resistance/CARD/reads/{sample}/{sample}.gene_mapping_data.txt", + sample=get_samples(), + project=get_project(), + ), + output: + touch("results/CARD_load_DB.done"), + log: + "logs/CARD_load_DB.log", + conda: + "../envs/card.yaml" + shell: + "rgi clean --local && " + "rgi load --card_json {input.db} --local > {log} 2>&1" + + +rule CARD_assembly_run: + input: + faa=rules.gzip_proteins.output.faa, + db=rules.CARD_load_DB.output, + output: + txt="results/{project}/output/resistance/CARD/assembly/{sample}/{sample}.txt", + json="results/{project}/output/resistance/CARD/assembly/{sample}/{sample}.json", + params: + path_wo_ext=lambda wildcards, output: Path(output.txt).with_suffix(""), + log: + "logs/{project}/ARGs/assembly/{sample}.log", + threads: 64 + conda: + "../envs/card.yaml" + shell: + "rgi main -i {input.faa} -o {params.path_wo_ext} " + "-t protein -a DIAMOND --local " + "-n {threads} --clean > {log} 2>&1" + + +use rule CARD_read_sample_summary as CARD_assembly_sample_summary with: + input: + txt=rules.CARD_assembly_run.output.txt, + json=rules.CARD_assembly_run.output.json, + output: + csv="results/{project}/output/resistance/CARD/assembly/{sample}/{sample}_assembly_ARGs.csv", + params: + case="assembly", + log: + "logs/{project}/ARGs/assembly/{sample}.log", + + +use rule CARD_assembly_run as CARD_mag_run with: + input: + fa="results/{project}/output/fastas/{sample}/mags/{binID}.fa.gz", + db=rules.CARD_load_DB.output, + output: + txt="results/{project}/output/resistance/CARD/mags/{sample}/{binID}.txt", + params: + path_wo_ext=lambda wildcards, output: Path(output.txt).with_suffix(""), + log: + "logs/{project}/ARGs/mags/{sample}/{binID}.log", + + +rule wrap_mag_ARGs: + input: + get_mag_ARGs, + output: + touch("results/{project}/output/ARGs/mags/{sample}/all_mags.done"), + log: + "logs/{project}/ARGs/mags/{sample}/all_mags.log", +""" diff --git a/workflow/rules/trimming.smk b/workflow/rules/trimming.smk deleted file mode 100644 index 2cf0384..0000000 --- a/workflow/rules/trimming.smk +++ /dev/null @@ -1,48 +0,0 @@ -RAW_DATA_PATH = get_data_path() - -""" -# copy files to local -rule copy_fastq: - output: - raw1=f"{RAW_DATA_PATH}{{project}}/{{sample}}_R1.fastq.gz", - raw2=f"{RAW_DATA_PATH}{{project}}/{{sample}}_R2.fastq.gz", - params: - outdir=lambda wildcards, output: Path(output.raw1).parent, - # returns a list of folder and the filenames for R1 and R2 reads - fastqs=get_fastqs, - threads: 20 - log: - "logs/{project}/copy_data/{sample}.log", - conda: - "../envs/unix.yaml" - shell: - "(mkdir -p {params.outdir} && " - "tar cpfz - -C {params.fastqs[0]}/ {params.fastqs[1]} {params.fastqs[2]} | " - "(cd {params.outdir} ; tar xpfz -)) > {log} 2>&1" -""" - - -# fastp in paired-end mode for Illumina paired-end data -rule fastp: - input: - sample=get_local_fastqs, - output: - trimmed=temp( - [ - "results/{project}/trimmed/fastp/{sample}.1.fastq.gz", - "results/{project}/trimmed/fastp/{sample}.2.fastq.gz", - ] - ), - html=temp("results/{project}/trimmed/fastp/{sample}.html"), - json="results/{project}/report_prerequisites/qc/{sample}.fastp.json", - params: - adapters=get_adapters, - extra="--qualified_quality_phred {phred} --length_required {minlen}".format( - phred=(config["quality-criteria"]["min-PHRED"]), - minlen=(config["quality-criteria"]["min-length-reads"]), - ), - log: - "logs/{project}/fastp/{sample}.log", - threads: 16 - wrapper: - "v2.6.0/bio/fastp" diff --git a/workflow/scripts/assembly_summary.py b/workflow/scripts/assembly_summary.py index ee5fbda..535e0a3 100644 --- a/workflow/scripts/assembly_summary.py +++ b/workflow/scripts/assembly_summary.py @@ -1,22 +1,27 @@ import pandas as pd import sys +import os sys.stderr = open(snakemake.log[0], "w") asbl_logs = snakemake.input.asbl txts = snakemake.input.mapped qc_csv = snakemake.input.qc_csv +csv_mags = snakemake.input.csv_mags +csv_bins = snakemake.input.csv_bins -df=pd.read_csv(qc_csv,index_col="sample") -samples=df.index.to_list() +df = pd.read_csv(qc_csv, index_col="sample") +samples = df.index.to_list() summary_dict = {} for sample in samples: summary_sample_dict = {} - summary_sample_dict["#reads_after_filtering"] = df.loc[sample,"#reads_after_filtering"] + summary_sample_dict["#reads_after_filtering"] = df.loc[ + sample, "#reads_after_filtering" + ] for asbl_log in asbl_logs: - if asbl_log.rfind(sample) >= 0: + if asbl_log.rfind(f"{sample}_megahit") >= 0: with open(asbl_log, "r") as a_log: for line in a_log: if line.find("total") >= 0: @@ -26,30 +31,46 @@ sep = sep.split() if "contigs" in sep: colname = "#contigs" - value=int(sep[0]) + value = int(sep[0]) elif "total" in sep: colname = "assembly_size(bp)" - value=int(sep[1]) + value = int(sep[1]) elif "max" in sep: colname = "longest_contig(bp)" - value=int(sep[1]) + value = int(sep[1]) elif "avg" in sep: colname = "average_contig_size(bp)" - value=int(sep[1]) + value = int(sep[1]) elif "N50" in sep: colname = "N50(bp)" - value=int(sep[1]) - + value = int(sep[1]) + summary_sample_dict[colname] = value for txt in txts: - if txt.rfind(sample) >= 0: + if txt.rfind(f"{sample}_reads_mapped") >= 0: with open(txt) as t: - mapped=int(t.readline()) + mapped = int(t.readline()) summary_sample_dict["#assembled_reads"] = mapped - summary_sample_dict["%assembled_reads"] = round(((mapped/summary_sample_dict["#reads_after_filtering"]) * 100),4) + summary_sample_dict["%assembled_reads"] = round( + ((mapped / summary_sample_dict["#reads_after_filtering"]) * 100), 2 + ) break + for bin_file in csv_bins: + if bin_file.rfind(sample) >= 0: + + bin_df = pd.read_csv(bin_file) + summary_sample_dict["#bins"] = len(bin_df) + + mag_file = [ + file + for file in csv_mags + if os.path.basename(os.path.dirname(file)) == sample + ][0] + mag_df = pd.read_csv(mag_file) + summary_sample_dict["#MAGs"] = len(mag_df) + summary_dict[sample] = summary_sample_dict summary_df = pd.DataFrame.from_dict(summary_dict, orient="index") @@ -61,19 +82,22 @@ summary_df.to_csv(snakemake.output.csv) + ## edit decimal separator for visual of table def format_int_with_commas(x): - #Formats an integer with commas as thousand separators + # Formats an integer with commas as thousand separators return f"{x:,}" + def format_float_with_percent(x): - #Formats a float with percent sign + # Formats a float with percent sign return f"{x:.2f}%" -cols=[col for col in summary_df.columns.to_list() if col.find("%")>=0] -summary_df[cols]=summary_df[cols].map(format_float_with_percent) -cols=[col for col in summary_df.columns.to_list() if col.find("%")<0] -summary_df[cols]=summary_df[cols].map(format_int_with_commas) +cols = [col for col in summary_df.columns.to_list() if col.find("%") >= 0] +summary_df[cols] = summary_df[cols].map(format_float_with_percent) + +cols = [col for col in summary_df.columns.to_list() if col.find("%") < 0] +summary_df[cols] = summary_df[cols].map(format_int_with_commas) summary_df.to_csv(snakemake.output.vis_csv) diff --git a/workflow/scripts/bin_summary_all.py b/workflow/scripts/bin_summary_all.py deleted file mode 100644 index 57ff66d..0000000 --- a/workflow/scripts/bin_summary_all.py +++ /dev/null @@ -1,32 +0,0 @@ -import pandas as pd -import sys -import os - -sys.stderr = open(snakemake.log[0], "w") -csv_mags = snakemake.input.csv_mags -csv_bins = snakemake.input.csv_bins - -summary_dict = {} -for bin_file in csv_bins: - sample = os.path.basename(os.path.dirname(bin_file)) - - summary_sample_dict = {} - - mag_file=[file for file in csv_mags if os.path.basename(os.path.dirname(file)) == sample][0] - mag_df=pd.read_csv(mag_file) - summary_sample_dict["# MAGs"] = len(mag_df) - - bin_df=pd.read_csv(bin_file) - summary_sample_dict["# bins"] = len(bin_df) - - summary_sample_dict["least contigs"] = bin_df["contigs"].min() - - summary_sample_dict["highest N50"] = f'{(bin_df["contig_N50"].max()):,}' - - summary_dict[sample]=summary_sample_dict - -summary_df = pd.DataFrame.from_dict(summary_dict, orient="index") -summary_df.index.name = "sample" -summary_df.sort_index(inplace=True) - -summary_df.to_csv(snakemake.output[0]) diff --git a/workflow/scripts/bin_summary_sample.py b/workflow/scripts/bin_summary_sample.py index a880292..d7acfc9 100644 --- a/workflow/scripts/bin_summary_sample.py +++ b/workflow/scripts/bin_summary_sample.py @@ -4,33 +4,15 @@ sys.stderr = open(snakemake.log[0], "w") ## input files -in_dastool = ( - snakemake.input.tool -) # "ResMAG/results/autobrewer/das_tool/ABS_24_07/ABS_24_07_DASTool_summary.tsv" -in_checkm = ( - snakemake.input.checkm -) # "ResMAG/results/autobrewer/qc/checkm2/ABS_24_07/quality_report.tsv" -in_gtdb = ( - snakemake.input.gtdb -) # "ResMAG/results/autobrewer/classification/ABS_24_07/ABS_24_07.bac120.summary.tsv" +in_dastool = snakemake.input.tool +in_checkm = snakemake.input.checkm +in_gtdb = snakemake.input.gtdb ## output files ### csv -csv_path_mags = ( - snakemake.output.csv_mags -) # "ResMAG/results/autobrewer/report/ABS_24_07/mags_summary.csv" # -csv_path_bins = ( - snakemake.output.csv_bins -) # "ResMAG/results/autobrewer/report/ABS_24_07/bin_summary.csv" -csv_path_tax = ( - snakemake.output.csv_tax -) # "ResMAG/results/autobrewer/report/ABS_24_07/bin_taxonomy.csv" -csv_path_checkm = ( - snakemake.output.csv_checkm -) # "ResMAG/results/autobrewer/report/ABS_24_07/checkm_summary.csv" -csv_path_dastool = ( - snakemake.output.csv_dastool -) # "ResMAG/results/autobrewer/report/ABS_24_07/DASTool_summary.csv" +csv_path_mags = snakemake.output.csv_mags +csv_path_bins = snakemake.output.csv_bins +csv_path_tax = snakemake.output.csv_tax ## params max_cont = snakemake.params.max_cont @@ -46,7 +28,7 @@ def save_csv_table(csv_path, summary_df): tool_df = pd.read_table(in_dastool) tool_df.drop(["bin_set"], axis=1, inplace=True) -# tool_df.rename({"bin":"bin"},axis=1,inplace=True) + tool_df.set_index("bin", inplace=True) col_order = [ "bin_score", @@ -64,9 +46,6 @@ def save_csv_table(csv_path, summary_df): del col_order[:2] tool_red_df = tool_df.drop(col_order, axis=1) -save_csv_table(csv_path_dastool, tool_red_df) - - checkm_df = pd.read_table(in_checkm) rm_list = ["Translation_Table_Used", "Additional_Notes"] checkm_df.drop(rm_list, axis=1, inplace=True) @@ -93,8 +72,6 @@ def save_csv_table(csv_path, summary_df): del col_order[:5] checkm_red_df = checkm_df.drop(col_order, axis=1) -save_csv_table(csv_path_checkm, checkm_red_df) - gtdb_df = pd.read_table(in_gtdb) gtdb_df.rename({"user_genome": "bin"}, axis=1, inplace=True) @@ -144,23 +121,3 @@ def save_csv_table(csv_path, summary_df): ] mags_df.index.name = "MAG" save_csv_table(csv_path_mags, mags_df) - - -""" -if we want to sort by mean of deviation from min completeness & max contamination - -mean_dict={} -for binid in bins_df.index: - cont=bins_df.at[binid,"contamination"] - cont_per=round(((max_contamination - cont)/max_contamination), 4) - comp=bins_df.at[binid,"completeness"] - comp_per=round((comp/100),4) - mean_dict[binid]=round(((cont_per + comp_per)/2),4) - -to_sort=pd.DataFrame.from_dict(mean_dict,orient="index",columns=['mean']) - -# index list sorted by mean used to reindex original df -new_ind=to_sort.sort_values(["mean"],ascending = [False]).index.to_list() -sorted_all=bins_df.reindex(new_ind) - -""" diff --git a/workflow/scripts/binner_control.py b/workflow/scripts/binner_control.py deleted file mode 100644 index a0e38cb..0000000 --- a/workflow/scripts/binner_control.py +++ /dev/null @@ -1,32 +0,0 @@ -import os -import sys - -sys.stderr = open(snakemake.log[0], "w") - -def check_files(files, binner_ls): - path_ls = [] - results_binner = binner_ls.copy() - - for binner in binner_ls: - for file_path in files: - if file_path.find(binner) >= 0: - # check if file exists and is not empty - if os.path.isfile(file_path) and os.path.getsize(file_path) == 0: - results_binner.remove(binner) - break - else: - path_ls.append(file_path) - break - - with open(output_csv, 'w') as f_out: - f_out.write("{paths}\n{binners}\n".format(paths=(",".join(path_ls)),binners=(",".join(results_binner)))) - - - -# Get the input and output file paths from command-line arguments -files = snakemake.input -binner_ls = snakemake.params[0] -output_csv = snakemake.output[0] - -# Call the function with the provided file paths -check_files(files, binner_ls) \ No newline at end of file diff --git a/workflow/scripts/contig_classification.py b/workflow/scripts/contig_classification.py new file mode 100644 index 0000000..6f3611e --- /dev/null +++ b/workflow/scripts/contig_classification.py @@ -0,0 +1,75 @@ +import pandas as pd +import re +import sys + +sys.stderr = open(snakemake.log[0], "w") + +kaiju_out = snakemake.input.out +kaiju_species = snakemake.input.species +kaiju_genus = snakemake.input.genus +headers = snakemake.input.headers + +classification = snakemake.output.csv + + +## Prepare classification and length information for all contigs classified to at least genus level ## + +df_kaiju_raw = pd.read_table(kaiju_out, names=["classified", "contig", "taxid"]) +kaiju_df = df_kaiju_raw[df_kaiju_raw["classified"] != "U"].copy() + +species_df = pd.read_csv(kaiju_species, sep="\t") +# remove rows for unclassified and virus +species_df = species_df.iloc[:-3].copy() + +genus_df = pd.read_csv(kaiju_genus, sep="\t") +# remove rows for unclassified and virus +genus_df = genus_df.iloc[:-3].copy() + +# make taxon name to taxon id dictionary for all occurring species and genera +species_map = dict(zip(species_df["taxon_id"], species_df["taxon_name"])) +genus_map = dict(zip(genus_df["taxon_id"], genus_df["taxon_name"])) + +# kaiju_df is the original DataFrame with only classified contigs +kaiju_df["species_name"] = kaiju_df["taxid"].map(species_map) +kaiju_df["genus_name"] = kaiju_df["taxid"].map(genus_map) + +# removing all not classified to species or genus level +kaiju_df_red = kaiju_df[ + ~(kaiju_df["species_name"].isna() & kaiju_df["genus_name"].isna()) +].copy() +# adding info about the level of classification +kaiju_df_red["level"] = ( + kaiju_df_red["species_name"].notna().map({True: "species", False: "genus"}) +) +kaiju_df_red["classification"] = kaiju_df_red["species_name"].combine_first( + kaiju_df_red["genus_name"] +) +kaiju_df_red = kaiju_df_red.drop( + columns=["species_name", "genus_name", "classified", "taxid"] +) + +# get contig length +records = [] + +with open(headers) as f: + for line in f: + # Example header: >k141_570066 flag=0 multi=1.0000 len=450 + contig_match = re.search(r"^>(\S+)", line) + length_match = re.search(r"len=(\d+)", line) + + if contig_match and length_match: + contig = contig_match.group(1) + length = int(length_match.group(1)) + records.append((contig, length)) + +# Convert to DataFrame +length_df = pd.DataFrame(records, columns=["contig", "contig_len"]) + +# merge contig classification with length information +merged_df = kaiju_df_red.merge(length_df, on="contig", how="left") +cols = merged_df.columns.tolist() +cols.insert(1, cols.pop(cols.index("contig_len"))) +merged_df = merged_df[cols] + +# Write out contig classification file +merged_df.to_csv(classification, index=False) diff --git a/workflow/scripts/diversity_summary.py b/workflow/scripts/diversity_summary.py deleted file mode 100644 index cea0b21..0000000 --- a/workflow/scripts/diversity_summary.py +++ /dev/null @@ -1,106 +0,0 @@ -import pandas as pd -import json -import sys - -sys.stderr = open(snakemake.log[0], "w") - -jsons = snakemake.input.jsons -human_logs = snakemake.input.human_logs -if snakemake.params.other_host: - hostname = snakemake.params.hostname - host_logs = snakemake.input.host_logs -reports = snakemake.input.reports -tax_ids = snakemake.params.taxid_dict -outfile = snakemake.output.csv - -# header for kraken report file -header_ls = ["prct", "total_reads", "lvl_reads", "lvl", "tax_id", "name"] -results_dict = {} - -for report in reports: - sample = report[report.rfind("/") + 1 : report.rfind("_kraken2_report")] - sample_results_dict = {} - - for json_file in jsons: - if json_file.find(sample) >= 0: - with open(json_file, "r") as read_file: - jdata = json.load(read_file) - # fastp json shows reads instead of read pairs as in the kraken reports - total_pre_host_filt = int( - (jdata["summary"]["after_filtering"]["total_reads"]) / 2 - ) - break - sample_results_dict["#total_reads"] = total_pre_host_filt - - if snakemake.params.other_host: - for host_log in host_logs: - if host_log.find(sample) >= 0: - with open(host_log, "r") as file: - for line in file.readlines(): - if line.find("processed") >= 0: - non_host_reads = int(line.split()[2]) - no_reads = total_pre_host_filt - non_host_reads - prct_reads = ( - int(no_reads) / int(total_pre_host_filt) - ) * 100 - break - - sample_results_dict[f"%{hostname}"] = "%.3f" % prct_reads - sample_results_dict[f"#reads_{hostname}"] = no_reads - break - - for human_log in human_logs: - if human_log.find(sample) >= 0: - with open(human_log, "r") as file: - for line in file.readlines(): - if line.find("processed") >= 0: - non_human_reads = int(line.split()[2]) - - if snakemake.params.other_host: - no_reads = ( - total_pre_host_filt - non_host_reads - non_human_reads - ) - - else: - no_reads = total_pre_host_filt - non_human_reads - - prct_reads = (int(no_reads) / int(total_pre_host_filt)) * 100 - break - - - sample_results_dict[f"%human"] = "%.3f" % prct_reads - sample_results_dict[f"#reads_human"] = no_reads - break - - report_df = pd.read_table(report, names=header_ls) - - for ref in tax_ids: - line = report_df.loc[report_df["tax_id"] == tax_ids[ref]] - - no_reads = line.iloc[0]["total_reads"] - prct_reads = (int(no_reads) / int(total_pre_host_filt)) * 100 - - sample_results_dict[f"%{ref}"] = "%.3f" % prct_reads - sample_results_dict[f"#reads_{ref}"] = no_reads - - # add unclassified reads - line = report_df.loc[report_df["tax_id"] == 0] - if not line.empty: - no_reads = line.iloc[0]["total_reads"] - prct_reads = (int(no_reads) / int(total_pre_host_filt)) * 100 - - else: - no_reads = 0 - prct_reads = 0 - - sample_results_dict["%unclassified"] = "%.3f" % prct_reads - sample_results_dict["#reads_unclassified"] = no_reads - - results_dict[sample] = sample_results_dict - - -results_df = pd.DataFrame.from_dict(results_dict, orient="index") -results_df.index.name = "sample" - -results_df = results_df.sort_index() -results_df.to_csv(outfile) diff --git a/workflow/scripts/metabat_postprocess.py b/workflow/scripts/metabat_postprocess.py index a38d269..eac3999 100644 --- a/workflow/scripts/metabat_postprocess.py +++ b/workflow/scripts/metabat_postprocess.py @@ -42,14 +42,19 @@ def append_string_to_file(filepath, append_string, mod_ext): for filename in os.listdir(directory) if filename.endswith(mod_ext) ] - - print("starting to concatenate files") - for modfile in modfiles: - # Concatenate all modified files into a single file - cmd = f"cat {modfile} >> {output_tsv}" + if len(modfiles) == 0: + print("no metabat bins") + # create an empty file + cmd = f"touch {output_tsv}" + os.system(cmd) + else: + print("starting to concatenate files") + for modfile in modfiles: + # Concatenate all modified files into a single file + cmd = f"cat {modfile} >> {output_tsv}" + os.system(cmd) + + print("all files concatenated") + # remove all modified files after concatenating + cmd = f"rm {directory}/*{mod_ext}" os.system(cmd) - - print("all files concatenated") - # remove all modified files after concatenating - cmd = f"rm {directory}/*{mod_ext}" - os.system(cmd) diff --git a/workflow/scripts/plot_host.py b/workflow/scripts/plot_host.py deleted file mode 100644 index 87602cc..0000000 --- a/workflow/scripts/plot_host.py +++ /dev/null @@ -1,142 +0,0 @@ -import pandas as pd -import altair as alt -import sys - -## write to log file -sys.stderr = open(snakemake.log[0], "w") - -## input file -csv_in = snakemake.input.csv - -#input parameter -other_host = snakemake.params.other_host - -## output file -host_percentage_html = snakemake.output.html - -## variables -color_red = "#e03e3e" -color_green = "#6aa84f" - -## prepare dataframe -df=pd.read_csv(csv_in) - -human_cont_df=pd.DataFrame() -human_cont_df["sample"]=df["sample"] -human_cont_df["human"]=df["%human"].divide(100) - -if other_host: - hostname=snakemake.params.hostname - human_cont_df["host"]=df[f"%{hostname}"].divide(100) - - -slider = alt.binding_range( - min=0, max=100, step=0.5, name="maximum acceptable:" -) -# selector = alt.param(name='SelectorName', value=50, bind=slider) -selector = alt.selection_point( - name="SelectorName", fields=["max_contamination"], bind=slider, value=50 -) - -# plot for human percentage -if other_host: - human_base = ( - alt.Chart(human_cont_df) - .encode( - alt.X("human:Q") - .axis(format="%", labelFontSize=12, titleFontSize=12) - .title("Percentage of human reads") - .scale(domain=[0, 1]), - alt.Y("sample:N") - .axis(labelFontSize=12, titleFontSize=12), - ) - .add_params(selector) - .properties(width="container") - .interactive() - ) - -else: - title="Share of human reads among all reads (after QC)" - title_object = alt.TitleParams(title, anchor='middle',fontSize=14) - - human_base = ( - alt.Chart(human_cont_df, title=title_object) - .encode( - alt.X("human:Q") - .axis(format="%", labelFontSize=12, titleFontSize=12) - .title("Percentage of human reads") - .scale(domain=[0, 1]), - alt.Y("sample:N") - .axis(labelFontSize=12, titleFontSize=12), - ) - .add_params(selector) - .properties(width="container") - .interactive() - ) - - -human_bars = human_base.mark_bar().encode( - color=alt.condition( - (alt.datum.human * 100) >= selector.max_contamination, - alt.value(color_red), - alt.value(color_green), - ) -) - -human_text = human_base.mark_text( - align="center", - baseline="middle", - dx=30, - fontSize=12, -).encode( - text=alt.Text("human:Q", format=".2%"), -) - -human_full = human_bars + human_text - -#if second host -if other_host: - # plot for other host percentage - host_base = ( - alt.Chart(human_cont_df) - .encode( - alt.X("host:Q") - .axis(format="%", labelFontSize=12, titleFontSize=12) - .title(f"Percentage of {hostname} reads") - .scale(domain=[0, 1]), - alt.Y("sample:N") - .axis(labelFontSize=12, titleFontSize=12), - ) - .add_params(selector) - .properties(width="container") - .interactive() - ) - - host_bars = host_base.mark_bar().encode( - color=alt.condition( - (alt.datum.host * 100) >= selector.max_contamination, - alt.value(color_red), - alt.value(color_green), - ) - ) - - host_text = host_base.mark_text( - align="center", - baseline="middle", - dx=30, - fontSize=12, - ).encode( - text=alt.Text("host:Q", format=".2%"), - ) - host_full = host_bars + host_text - - title=f"Share of human and {hostname} reads among all reads (after QC)" - title_object = alt.TitleParams(title, anchor='middle',fontSize=14) - - all_chart=alt.vconcat(host_full, human_full, title=title_object) - all_chart.save(host_percentage_html) - -else: - human_full.save(host_percentage_html) - - diff --git a/workflow/scripts/quality_summary.py b/workflow/scripts/quality_summary.py index cd9e0f2..6f282aa 100644 --- a/workflow/scripts/quality_summary.py +++ b/workflow/scripts/quality_summary.py @@ -11,8 +11,8 @@ hostname=snakemake.params.hostname host_logs = snakemake.input.host_logs -outfile = snakemake.output.csv #"/local/work/josefa/ResMAG/results/lanuv/output/report/all/seq_summary.csv" -outfile_vis = snakemake.output.vis_csv #"/local/work/josefa/ResMAG/results/lanuv/output/report/all/seq_summary_vis.csv" +outfile = snakemake.output.csv +outfile_vis = snakemake.output.vis_csv results_dict = {} @@ -29,21 +29,21 @@ total_pre_host_filt = int(jdata["summary"]["after_filtering"]["total_reads"]) sample_results_dict["#reads_afterQC"] = total_pre_host_filt - sample_results_dict["%reads_filtered_by_quality"] = round(((1 - (total_pre_host_filt/reads_before)) * 100),4) + sample_results_dict["%reads_filtered_by_quality"] = round(((1 - (total_pre_host_filt/reads_before)) * 100),2) sample_results_dict["#bp_afterQC"] = int(jdata["summary"]["after_filtering"]["total_bases"]) - sample_results_dict["%Q30_reads_afterQC"] = float(jdata["summary"]["after_filtering"]["q30_rate"]) * 100 + sample_results_dict["%Q30_reads_afterQC"] = round((float(jdata["summary"]["after_filtering"]["q30_rate"]) * 100),2) if snakemake.params.other_host: for host_log in host_logs: - if host_log.find(sample) >= 0: + if host_log.find(f'{sample}_filter_host') >= 0: with open(host_log, "r") as file: for line in file.readlines(): if line.find("processed") >= 0: non_host_reads = int(line.split()[2]) * 2 no_reads = total_pre_host_filt - non_host_reads - prct_reads = round(((int(no_reads) / int(total_pre_host_filt)) * 100),4) + prct_reads = round(((int(no_reads) / int(total_pre_host_filt)) * 100),2) break sample_results_dict[f"#{hostname}_reads"] = no_reads @@ -52,7 +52,7 @@ for human_log in human_logs: - if human_log.find(sample) >= 0: + if human_log.find(f'{sample}_filter_human') >= 0: with open(human_log, "r") as file: for line in file.readlines(): if line.find("processed") >= 0: @@ -61,7 +61,7 @@ no_reads = non_host_reads - non_human_reads else: no_reads = total_pre_host_filt - non_human_reads - prct_reads = round(((int(no_reads) / int(total_pre_host_filt)) * 100),4) + prct_reads = round(((int(no_reads) / int(total_pre_host_filt)) * 100),2) break sample_results_dict[f"#human_reads"] = no_reads @@ -69,7 +69,7 @@ break sample_results_dict["#reads_after_filtering"] = non_human_reads - sample_results_dict["%reads_filtered_total"] = round(((1 - (non_human_reads/reads_before)) * 100),4) + sample_results_dict["%reads_filtered_total"] = round(((1 - (non_human_reads/reads_before)) * 100),2) results_dict[sample] = sample_results_dict diff --git a/workflow/scripts/uniCARD_abundance.py b/workflow/scripts/uniCARD_abundance.py new file mode 100644 index 0000000..3235bf4 --- /dev/null +++ b/workflow/scripts/uniCARD_abundance.py @@ -0,0 +1,52 @@ +import os +import pandas as pd +import sys + +sys.stderr = open(snakemake.log[0], "w") + + +contig_files = snakemake.input.contig_files +abd_class_files = snakemake.input.abd_class_files +abd_ARGs_files = snakemake.input.abd_ARGs_files + + +## Combining contig and ARG information for all samples +classes_all_samples = [] +names_all_samples = [] +for contig_file in contig_files: + sample_id = os.path.basename(os.path.dirname(contig_file)) + + classification_df = pd.read_csv(contig_file) + abd_class_file = [file for file in abd_class_files if file.find(sample_id) >= 0][0] + abd_classes_df = pd.read_csv(abd_class_file) + abd_ARGs_file = [file for file in abd_ARGs_files if file.find(sample_id) >= 0][0] + abd_names_df = pd.read_csv(abd_ARGs_file) + + merged_classes_df = classification_df.merge( + abd_classes_df, on="contig", how="inner" + ) + merged_classes_df.insert(0, "sampleID", sample_id) + + classes_all_samples.append(merged_classes_df) + + merged_names_df = classification_df.merge(abd_names_df, on="contig", how="inner") + merged_names_df.insert(0, "sampleID", sample_id) + + names_all_samples.append(merged_names_df) + + +final_classes_df = pd.concat(classes_all_samples, ignore_index=True) +final_classes_df.to_csv(snakemake.output.abd_class_possible, index=False) + +# in antibiotic classes per sample file all antibiotic classes in CARD are included +# All that do not appear in at least one sample are removed: +final_classes_df_red = final_classes_df.loc[:, (final_classes_df != 0).any(axis=0)] +final_classes_df_red.to_csv(snakemake.output.abd_class_present, index=False) + +final_names_df = pd.concat(names_all_samples, ignore_index=True) +final_names_df.to_csv(snakemake.output.abd_ARGs_possible, index=False) + +# in ARG per sample file all ARGs in CARD are included +# All that do not appear in at least one sample are removed: +final_names_df_red = final_names_df.loc[:, (final_names_df != 0).any(axis=0)] +final_names_df_red.to_csv(snakemake.output.abd_ARGs_present, index=False) diff --git a/workflow/scripts/uniCARD_abundance_sample.py b/workflow/scripts/uniCARD_abundance_sample.py new file mode 100644 index 0000000..1f5553a --- /dev/null +++ b/workflow/scripts/uniCARD_abundance_sample.py @@ -0,0 +1,95 @@ +import pandas as pd +import json +import sys + +sys.stderr = open(snakemake.log[0], "w") + +card_drugs = snakemake.input.json +args_file = snakemake.input.csv +# abundance outfiles +abd_class = snakemake.output.abd_class +abd_ARGs = snakemake.output.abd_ARGs + +## Write drug classes and ARO names abundance: +with open(card_drugs) as f: + data = json.load(f) + +# Convert to DataFrame +df_card = pd.DataFrame.from_dict(data, orient="index") + +# Split semicolon-separated values into lists +df_card["classes"] = df_card["classes"].str.split("; ") +df_card["antibiotics"] = df_card["antibiotics"].str.split("; ") + +# Flatten and get unique values +all_classes = pd.Series( + [item for sublist in df_card["classes"].dropna() for item in sublist] +).unique() +all_classes.sort() +all_antibiotics = pd.Series( + [item for sublist in df_card["antibiotics"].dropna() for item in sublist if item] +).unique() +all_aros = df_card["ARO_name"].unique() +all_aros.sort() + + +##Create classes and ARG names abundance files +df_args = pd.read_csv(args_file) +df_args_for_classes = df_args.copy() +# Step 1: Preprocess the classes column into lists +df_args_for_classes["class_list"] = ( + df_args_for_classes["classes"] + .fillna("") + .apply(lambda x: [c.strip() for c in x.split(";")]) +) + +# Step 2: Create one-hot columns per class +for c in all_classes: + df_args_for_classes[c] = df_args_for_classes["class_list"].apply( + lambda lst: int(c in lst) + ) + +# Step 3: Group by contig +agg_dict = {"ARO_ID": "count"} +agg_dict.update( + {c: "sum" for c in all_classes} # was max before +) # "sum" counts how many times the class appears + + +class_abundance_df = df_args_for_classes.groupby("contig").agg(agg_dict).reset_index() +# Write out abundance file +class_abundance_df = class_abundance_df.rename(columns={"ARO_ID": "#ARGs"}) +class_abundance_df.to_csv(abd_class, index=False) + +df_args["present"] = 1 + +# Step 2: Pivot to get binary presence/absence for all ARO_name per contig +aro_matrix = df_args.pivot_table( + index="contig", + columns="ARO_name", + values="present", + aggfunc="sum", # <-- changed from "max" to "sum" + fill_value=0, +) + +# Step 3: Count ARGs (genes) per contig +gene_counts = df_args.groupby("contig").size().rename("#ARGs") + +# Step 4: Join gene counts and presence matrix +names_abundance_df = pd.concat([gene_counts, aro_matrix], axis=1).reset_index() + +# Find missing AROs +missing_aros = [aro for aro in all_aros if aro not in names_abundance_df.columns] + +# Create a DataFrame with 0s for missing AROs (same index as result) +missing_df = pd.DataFrame(0, index=names_abundance_df.index, columns=missing_aros) + +# Concatenate in one go — avoids fragmentation +names_abundance_df = pd.concat([names_abundance_df, missing_df], axis=1) + +# Optional: reorder columns — num_genes first, then sorted AROs +aro_columns = sorted([a for a in all_aros if a in names_abundance_df.columns]) +names_abundance_df = names_abundance_df[["contig", "#ARGs"] + aro_columns] + +# Write out abundance file +names_abundance_df.to_csv(abd_ARGs, index=False) diff --git a/workflow/scripts/uniCARD_filtering.py b/workflow/scripts/uniCARD_filtering.py new file mode 100644 index 0000000..0e62acd --- /dev/null +++ b/workflow/scripts/uniCARD_filtering.py @@ -0,0 +1,171 @@ +import json +import pandas as pd +import argparse + +# identifier for CARD derived sequences +WT_PREFIX = "CARD_WT:" +CARD_PREFIX = "CARD:" + +# columns to reduce DIAMOND output +DIAMOND_COLS = [ + "query", + "target", + "%identity", + "length", + "#mismatches", + "#gaps", + "query_start", + "query_end", + "target_start", + "target_end", + "evalue", + "bitscore", +] + + +def write_filtered_ARGs(card_classes, unicard_in, outfile, filter_window): + + try: + # read in CARD drugs and classes json file + with open(card_classes) as f: + data = json.load(f) + + card_df = pd.DataFrame.from_dict(data, orient="index") + card_df.index.name = "ARO_ID" + card_df = card_df.reset_index() + + if unicard_in.endswith(".gz"): + df_unicard = pd.read_table( + unicard_in, names=DIAMOND_COLS, compression="gzip" + ) + else: + df_unicard = pd.read_table(unicard_in, names=DIAMOND_COLS) + + ## top 10 ## + # Step 1: Keep only 10 best hits per query + df_top10 = df_unicard[df_unicard.groupby("query").cumcount() < 10] + + # Step 2: Keep all rows for a query **only if at least one Target starts with CARD_PREFIX, + # so wild-type CARD entries are not kept (starting with WT_PREFIX) + df_top10 = df_top10[ + df_top10["query"].isin( + df_top10[df_top10["target"].str.startswith(CARD_PREFIX)]["query"] + ) + ] + + ## filter besthit CARD entries + + # filter for the first row per query + df_firsts = df_top10[df_top10.groupby("query").cumcount() < 1] + # 1. check if best hit is CARD wild-type -> remove from df + df_firsts_reduced = df_firsts[ + ~df_firsts["query"].isin( + df_firsts[df_firsts["target"].str.startswith(WT_PREFIX)]["query"] + ) + ] + # 2. from remaining save those, which are CARD hits + df_bests = df_firsts_reduced[ + df_firsts_reduced["query"].isin( + df_firsts_reduced[ + df_firsts_reduced["target"].str.startswith(CARD_PREFIX) + ]["query"] + ) + ] + + # remove first hit CARD queries from first row per query df + df_firsts_remaining = df_firsts_reduced[ + ~df_firsts_reduced["query"].isin(df_bests["query"]) + ] + + df_filtered = df_top10[df_top10["query"].isin(df_firsts_remaining["query"])] + + filtered_card_rows = [] + # go through all remaining best hits and compare to CARD hits + for idx, row in df_firsts_remaining.iterrows(): + # get all rows with same query to compare + compare_df = df_filtered[df_filtered["query"] == row["query"]] + compare_df = compare_df[ + compare_df["target"].str.startswith(CARD_PREFIX) + & ~compare_df["target"].str.contains("_snp_") + ] + # check for each row in this comparison df if the CARD entry has same values than first match + for idx_card, row_card in compare_df.iterrows(): + if ( + row["%identity"] == row_card["%identity"] + and row["bitscore"] == row_card["bitscore"] + ): + filtered_card_rows.append(row_card) + break + # additional filtering in default 1% window below best hit + elif (row["%identity"] * (1 - filter_window)) <= row_card[ + "%identity" + ] and (row["bitscore"] * (1 - filter_window)) <= row_card["bitscore"]: + filtered_card_rows.append(row_card) + break + + filtered_card_df = pd.DataFrame(filtered_card_rows) + + df_bests = pd.concat([df_bests, filtered_card_df]) + + # Adding contig column + df_bests.insert( + 1, "contig", df_top10["query"].map(lambda x: x.rsplit("_", 1)[0]) + ) + # extract ARO_ID, for variants it's eg 'CARD:3003394_snp_R123S' + df_bests.insert( + 3, + "ARO_ID", + df_top10["target"].map( + lambda x: x.split(CARD_PREFIX, 1)[-1].split("_", 1)[0] + ), + ) + + # Adding CARD drug info: + merged_df = pd.merge(df_bests, card_df, on="ARO_ID", how="left") + + # Write besthits to output file + merged_df.to_csv(outfile, index=False) + + except Exception as e: + print(f"Error filtering UniCARD results: {e}") + raise + + +def main(): + parser = argparse.ArgumentParser(description="Filter UniCARD results") + + parser.add_argument( + "--infile", + required=True, + help="DIAMOND output tsv file with UniCARD results to filter, may be gz", + ) + parser.add_argument( + "--card_hierarchy", required=True, help="Path to CARDs hierarchy json file" + ) + parser.add_argument( + "--outfile", + required=True, + help="Output csv file to store final filtered ARGs", + ) + parser.add_argument( + "--filter_window", + type=float, + default=0.01, + help="Percent window for filtering (default 0.01 (1%%))", + ) + + args = parser.parse_args() + + infile = args.infile + outfile = args.outfile + card_hierarchy = args.card_hierarchy + filter_window = args.filter_window + + print("Filtering UniCARD diamond results...") + write_filtered_ARGs(card_hierarchy, infile, outfile, filter_window) + + print(f"Filtering done! Results saved to: {outfile}") + + +if __name__ == "__main__": + main()