diff --git a/.github/ISSUE_TEMPLATE/bug_report.yml b/.github/ISSUE_TEMPLATE/bug_report.yml index 5e1ed5d..b58cfe9 100644 --- a/.github/ISSUE_TEMPLATE/bug_report.yml +++ b/.github/ISSUE_TEMPLATE/bug_report.yml @@ -2,53 +2,53 @@ name: Bug report description: Report something that is broken or incorrect labels: bug body: -- type: textarea - id: description - attributes: - label: Description of the bug - description: A clear and concise description of what the bug is. - validations: - required: true -- type: textarea - id: command_used - attributes: - label: Command used and terminal output - description: Steps to reproduce the behaviour. Please paste the command you used - to launch the pipeline and the output from your terminal. - render: console - placeholder: '$ nextflow run ... - - - Some output where something broke - - ' -- type: textarea - id: files - attributes: - label: Relevant files - description: 'Please drag and drop the relevant files here. Create a `.zip` archive - if the extension is not allowed. - - Your verbose log file `.nextflow.log` is often useful _(this is a hidden file - in the directory where you launched the pipeline)_ as well as custom Nextflow - configuration files. - - ' -- type: textarea - id: system - attributes: - label: System information - description: '* Nextflow version _(eg. 23.04.0)_ - - * Hardware _(eg. HPC, Desktop, Cloud)_ - - * Executor _(eg. slurm, local, awsbatch)_ - - * Container engine: _(e.g. Docker, Singularity, Conda, Podman, Shifter, Charliecloud, - or Apptainer)_ - - * OS _(eg. CentOS Linux, macOS, Linux Mint)_ - - * Version of jlab/refbasedassemblereval _(eg. 1.1, 1.5, 1.8.2)_ - - ' + - type: textarea + id: description + attributes: + label: Description of the bug + description: A clear and concise description of what the bug is. + validations: + required: true + - type: textarea + id: command_used + attributes: + label: Command used and terminal output + description: Steps to reproduce the behaviour. Please paste the command you used + to launch the pipeline and the output from your terminal. + render: console + placeholder: "$ nextflow run ... + + + Some output where something broke + + " + - type: textarea + id: files + attributes: + label: Relevant files + description: "Please drag and drop the relevant files here. Create a `.zip` archive + if the extension is not allowed. + + Your verbose log file `.nextflow.log` is often useful _(this is a hidden file + in the directory where you launched the pipeline)_ as well as custom Nextflow + configuration files. + + " + - type: textarea + id: system + attributes: + label: System information + description: "* Nextflow version _(eg. 23.04.0)_ + + * Hardware _(eg. HPC, Desktop, Cloud)_ + + * Executor _(eg. slurm, local, awsbatch)_ + + * Container engine: _(e.g. Docker, Singularity, Conda, Podman, Shifter, Charliecloud, + or Apptainer)_ + + * OS _(eg. CentOS Linux, macOS, Linux Mint)_ + + * Version of jlab/refbasedassemblereval _(eg. 1.1, 1.5, 1.8.2)_ + + " diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 4b214e9..9890e78 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -14,30 +14,3 @@ env: concurrency: group: "${{ github.workflow }}-${{ github.event.pull_request.number || github.ref }}" cancel-in-progress: true - -jobs: - test: - name: Run pipeline with test data - # Only run on push if this is the nf-core dev branch (merged PRs) - if: "${{ github.event_name != 'push' || (github.event_name == 'push' && github.repository == 'jlab/refbasedassemblereval') }}" - runs-on: ubuntu-latest - strategy: - matrix: - NXF_VER: - - "23.04.0" - - "latest-everything" - steps: - - name: Check out pipeline code - uses: actions/checkout@v3 - - - name: Install Nextflow - uses: nf-core/setup-nextflow@v1 - with: - version: "${{ matrix.NXF_VER }}" - - - name: Run pipeline with test data - # TODO nf-core: You can customise CI pipeline run tests as required - # For example: adding multiple test runs with different parameters - # Remember that you can parallelise this by using strategy.matrix - run: | - nextflow run ${GITHUB_WORKSPACE} -profile test,docker --outdir ./results diff --git a/.github/workflows/linting.yml b/.github/workflows/linting.yml index b8bdd21..0b6eb95 100644 --- a/.github/workflows/linting.yml +++ b/.github/workflows/linting.yml @@ -66,43 +66,3 @@ jobs: Thanks again for your contribution! repo-token: ${{ secrets.GITHUB_TOKEN }} allow-repeats: false - - nf-core: - runs-on: ubuntu-latest - steps: - - name: Check out pipeline code - uses: actions/checkout@v3 - - - name: Install Nextflow - uses: nf-core/setup-nextflow@v1 - - - uses: actions/setup-python@v4 - with: - python-version: "3.11" - architecture: "x64" - - - name: Install dependencies - run: | - python -m pip install --upgrade pip - pip install nf-core - - - name: Run nf-core lint - env: - GITHUB_COMMENTS_URL: ${{ github.event.pull_request.comments_url }} - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} - GITHUB_PR_COMMIT: ${{ github.event.pull_request.head.sha }} - run: nf-core -l lint_log.txt lint --dir ${GITHUB_WORKSPACE} --markdown lint_results.md - - - name: Save PR number - if: ${{ always() }} - run: echo ${{ github.event.pull_request.number }} > PR_number.txt - - - name: Upload linting log file artifact - if: ${{ always() }} - uses: actions/upload-artifact@v3 - with: - name: linting-logs - path: | - lint_log.txt - lint_results.md - PR_number.txt diff --git a/.gitignore b/.gitignore index 90061f3..ad2f7a2 100644 --- a/.gitignore +++ b/.gitignore @@ -26,4 +26,5 @@ plotting/*.csv* plotting/data_cache plotting/envs/ plotting/envs_newtest/ -plotting/archive \ No newline at end of file +plotting/archive +plotting/outdir/ \ No newline at end of file diff --git a/.nf-core.yml b/.nf-core.yml index 5eeaef5..f39ce35 100644 --- a/.nf-core.yml +++ b/.nf-core.yml @@ -1,32 +1,32 @@ lint: files_exist: - - CODE_OF_CONDUCT.md - - assets/nf-core-refbasedassemblereval_logo_light.png - - docs/images/nf-core-refbasedassemblereval_logo_light.png - - docs/images/nf-core-refbasedassemblereval_logo_dark.png - - .github/ISSUE_TEMPLATE/config.yml - - .github/workflows/awstest.yml - - .github/workflows/awsfulltest.yml - - conf/igenomes.config - - conf/igenomes.config + - CODE_OF_CONDUCT.md + - assets/nf-core-refbasedassemblereval_logo_light.png + - docs/images/nf-core-refbasedassemblereval_logo_light.png + - docs/images/nf-core-refbasedassemblereval_logo_dark.png + - .github/ISSUE_TEMPLATE/config.yml + - .github/workflows/awstest.yml + - .github/workflows/awsfulltest.yml + - conf/igenomes.config + - conf/igenomes.config files_unchanged: - - CODE_OF_CONDUCT.md - - assets/nf-core-refbasedassemblereval_logo_light.png - - docs/images/nf-core-refbasedassemblereval_logo_light.png - - docs/images/nf-core-refbasedassemblereval_logo_dark.png - - .github/ISSUE_TEMPLATE/bug_report.yml + - CODE_OF_CONDUCT.md + - assets/nf-core-refbasedassemblereval_logo_light.png + - docs/images/nf-core-refbasedassemblereval_logo_light.png + - docs/images/nf-core-refbasedassemblereval_logo_dark.png + - .github/ISSUE_TEMPLATE/bug_report.yml multiqc_config: - - report_comment + - report_comment nextflow_config: - - manifest.name - - manifest.homePage - - process.cpus - - process.memory - - process.time - - custom_config + - manifest.name + - manifest.homePage + - process.cpus + - process.memory + - process.time + - custom_config repository_type: pipeline template: prefix: jlab skip: - - igenomes - - nf_core_configs + - igenomes + - nf_core_configs diff --git a/CITATIONS.md b/CITATIONS.md index 22d7d4c..8df7905 100644 --- a/CITATIONS.md +++ b/CITATIONS.md @@ -13,51 +13,49 @@ - [FastQC](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) > Andrews, S. (2010). FastQC: A Quality Control Tool for High Throughput Sequence Data [Online]. - > + - [MultiQC](https://pubmed.ncbi.nlm.nih.gov/27312411/) > Ewels P, Magnusson M, Lundin S, Käller M. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics. 2016 Oct 1;32(19):3047-8. doi: 10.1093/bioinformatics/btw354. Epub 2016 Jun 16. PubMed PMID: 27312411; PubMed Central PMCID: PMC5039924. - > + - [Bowtie2](https://link.springer.com/article/10.1186/s13040-014-0034-0) - > LANGDON, William B. Performance of genetic programming optimised Bowtie2 on genome comparison and analytic testing (GCAT) benchmarks.*BioData mining* , 2015, 8. Jg., Nr. 1, S. 1. - > + > LANGDON, William B. Performance of genetic programming optimised Bowtie2 on genome comparison and analytic testing (GCAT) benchmarks._BioData mining_ , 2015, 8. Jg., Nr. 1, S. 1. + - [BEDTools](https://currentprotocols.onlinelibrary.wiley.com/doi/full/10.1002/0471250953.bi1112s47) - > QUINLAN, Aaron R. BEDTools: the Swiss‐army tool for genome feature analysis. *Current protocols in bioinformatics* , 2014, 47. Jg., Nr. 1, S. 11.12. 1-11.12. 34. - > + > QUINLAN, Aaron R. BEDTools: the Swiss‐army tool for genome feature analysis. _Current protocols in bioinformatics_ , 2014, 47. Jg., Nr. 1, S. 11.12. 1-11.12. 34. + - [Salmon](https://www.nature.com/articles/nmeth.4197) - > PATRO, Rob, et al. Salmon provides fast and bias-aware quantification of transcript expression.*Nature methods* , 2017, 14. Jg., Nr. 4, S. 417-419. - > + > PATRO, Rob, et al. Salmon provides fast and bias-aware quantification of transcript expression._Nature methods_ , 2017, 14. Jg., Nr. 4, S. 417-419. + - [SeqKit](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0163962) - > SHEN, Wei, et al. SeqKit: a cross-platform and ultrafast toolkit for FASTA/Q file manipulation. *PloS one* , 2016, 11. Jg., Nr. 10, S. e0163962. - > + > SHEN, Wei, et al. SeqKit: a cross-platform and ultrafast toolkit for FASTA/Q file manipulation. _PloS one_ , 2016, 11. Jg., Nr. 10, S. e0163962. + - [Subread](https://academic.oup.com/nar/article/41/10/e108/1075719?login=true) - > LIAO, Yang; SMYTH, Gordon K.; SHI, Wei. The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote. *Nucleic acids research* , 2013, 41. Jg., Nr. 10, S. e108-e108. - > + > LIAO, Yang; SMYTH, Gordon K.; SHI, Wei. The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote. _Nucleic acids research_ , 2013, 41. Jg., Nr. 10, S. e108-e108. ## Software packaging/containerisation tools - [Anaconda](https://anaconda.com) > Anaconda Software Distribution. Computer software. Vers. 2-2.4.0. Anaconda, Nov. 2016. Web. - > + - [Bioconda](https://pubmed.ncbi.nlm.nih.gov/29967506/) > Grüning B, Dale R, Sjödin A, Chapman BA, Rowe J, Tomkins-Tinch CH, Valieris R, Köster J; Bioconda Team. Bioconda: sustainable and comprehensive software distribution for the life sciences. Nat Methods. 2018 Jul;15(7):475-476. doi: 10.1038/s41592-018-0046-7. PubMed PMID: 29967506. - > + - [BioContainers](https://pubmed.ncbi.nlm.nih.gov/28379341/) > da Veiga Leprevost F, Grüning B, Aflitos SA, Röst HL, Uszkoreit J, Barsnes H, Vaudel M, Moreno P, Gatto L, Weber J, Bai M, Jimenez RC, Sachsenberg T, Pfeuffer J, Alvarez RV, Griss J, Nesvizhskii AI, Perez-Riverol Y. BioContainers: an open-source and community-driven framework for software standardization. Bioinformatics. 2017 Aug 15;33(16):2580-2582. doi: 10.1093/bioinformatics/btx192. PubMed PMID: 28379341; PubMed Central PMCID: PMC5870671. - > + - [Docker](https://dl.acm.org/doi/10.5555/2600239.2600241) > Merkel, D. (2014). Docker: lightweight linux containers for consistent development and deployment. Linux Journal, 2014(239), 2. doi: 10.5555/2600239.2600241. - > + - [Singularity](https://pubmed.ncbi.nlm.nih.gov/28494014/) > Kurtzer GM, Sochat V, Bauer MW. Singularity: Scientific containers for mobility of compute. PLoS One. 2017 May 11;12(5):e0177459. doi: 10.1371/journal.pone.0177459. eCollection 2017. PubMed PMID: 28494014; PubMed Central PMCID: PMC5426675. - > diff --git a/README.md b/README.md index d3fdc9c..9bc28eb 100644 --- a/README.md +++ b/README.md @@ -17,13 +17,13 @@ Recommended usage: -* Create *in silico* datasets using [Marbel](https://anaconda.org/bioconda/marbel), name the output dirs: NAME_microbiome -* Assemble the datasets -* Run the pipeline +- Create _in silico_ datasets using [Marbel](https://anaconda.org/bioconda/marbel), name the output dirs: NAME_microbiome +- Assemble the datasets +- Run the pipeline Fill the required parameters per run and adjust the example config files in `example/tool.config`, `example/dataset.config` and `example/resources.config`. -For each dataset create a dataset.config with the data set name. For outdir parameter chose the same parent dir, if you want to visualise the datasets together, e.g., PATH/group/dataset1, PATH/group/dataset2. +For each dataset create a dataset.config with the data set name. For outdir parameter chose the same parent dir, if you want to visualise the datasets together, e.g., PATH/group/dataset1, PATH/group/dataset2. Run each dataset with: @@ -35,7 +35,7 @@ nextflow run . \ -c example/tool.config ``` -Afterwards there is multiple methods which can visualise the runs. A script which creates all orb files is provided: +Afterwards there is multiple methods which can visualise the runs. A script which creates all orb files is provided: `plotting/plot_orb_figures.py For the dependencies you can use: @@ -52,7 +52,7 @@ python plotting/plot_orb_figures.py 0: assembler_dge.loc[group.index, "de"] = True -cds_contig_dict = assembler_dge.drop_duplicates(subset="gene_name").set_index("gene_name")[["ContigName", "de"]].apply(list, axis=1).to_dict() +cds_contig_dict = ( + assembler_dge.drop_duplicates(subset="gene_name") + .set_index("gene_name")[["ContigName", "de"]] + .apply(list, axis=1) + .to_dict() +) contig_map, de_map = zip(*[cds_contig_dict.get(gene, ["", False]) for gene in gene_summary["ContigName"]]) @@ -47,6 +52,13 @@ cm = confusion_matrix(true_labels, pred_labels, labels=[0, 1]) -linear_cm = {f"{row_prefix}DE_TN": cm[0, 0], f"{row_prefix}DE_FP": cm[0, 1], f"{row_prefix}DE_FN": cm[1, 0], f"{row_prefix}DE_TP": cm[1, 1]} +linear_cm = { + f"{row_prefix}DE_TN": cm[0, 0], + f"{row_prefix}DE_FP": cm[0, 1], + f"{row_prefix}DE_FN": cm[1, 0], + f"{row_prefix}DE_TP": cm[1, 1], +} -pd.DataFrame.from_dict(linear_cm, orient="index", columns=[prefix]).to_csv(f"{row_prefix}{prefix}_linear_cm.tsv", sep="\t") +pd.DataFrame.from_dict(linear_cm, orient="index", columns=[prefix]).to_csv( + f"{row_prefix}{prefix}_linear_cm.tsv", sep="\t" +) diff --git a/bin/calculate_og_confusion.py b/bin/calculate_og_confusion.py index 9941765..fa92910 100755 --- a/bin/calculate_og_confusion.py +++ b/bin/calculate_og_confusion.py @@ -13,8 +13,8 @@ log_name = sys.argv[5] row_prefix = sys.argv[6] -gene_summary = pd.read_csv(gene_summary_file, sep='\t') -assembler_dge = pd.read_csv(assembler_dge_file, sep='\t') +gene_summary = pd.read_csv(gene_summary_file, sep="\t") +assembler_dge = pd.read_csv(assembler_dge_file, sep="\t") gene_summary["cds_de"] = (gene_summary[pval_name] < 0.05) & (gene_summary[log_name].abs() > 1) @@ -24,7 +24,7 @@ cds_contig_dict = dict(zip(assembler_dge["ContigName"], assembler_dge["de"])) -gene_summary["assembler_de"] = gene_summary["ContigName"].apply(lambda x: cds_contig_dict.get(x, False)) +gene_summary["assembler_de"] = gene_summary["ContigName"].apply(lambda x: cds_contig_dict.get(x, False)) true_labels = (gene_summary["cds_de"]).astype(int) pred_labels = (gene_summary["assembler_de"]).astype(int) @@ -32,6 +32,13 @@ cm = confusion_matrix(true_labels, pred_labels, labels=[0, 1]) -linear_cm = {f"{row_prefix}DE_TN": cm[0, 0], f"{row_prefix}DE_FP": cm[0, 1], f"{row_prefix}DE_FN": cm[1, 0], f"{row_prefix}DE_TP": cm[1, 1]} +linear_cm = { + f"{row_prefix}DE_TN": cm[0, 0], + f"{row_prefix}DE_FP": cm[0, 1], + f"{row_prefix}DE_FN": cm[1, 0], + f"{row_prefix}DE_TP": cm[1, 1], +} -pd.DataFrame.from_dict(linear_cm, orient="index", columns=[prefix]).to_csv(f"{row_prefix}{prefix}_linear_cm.tsv", sep="\t") \ No newline at end of file +pd.DataFrame.from_dict(linear_cm, orient="index", columns=[prefix]).to_csv( + f"{row_prefix}{prefix}_linear_cm.tsv", sep="\t" +) diff --git a/bin/calculate_ratios.py b/bin/calculate_ratios.py deleted file mode 100755 index a055397..0000000 --- a/bin/calculate_ratios.py +++ /dev/null @@ -1,22 +0,0 @@ -#!/usr/bin/env python - -import pandas as pd -import sys - -score_over_view = sys.argv[1] -blocks = sys.argv[2] -chimeric_blocks = sys.argv[3] - -scores_df = pd.read_csv(score_over_view, sep='\t', index_col=0) -blocks_df = pd.read_csv(blocks, sep='\t', index_col=0) -chimeric_blocks_df = pd.read_csv(chimeric_blocks, sep='\t', index_col=0) - -n_blocks = blocks_df.shape[0] -n_chimeric_blocks = chimeric_blocks_df.shape[0] - -scores_df.loc["mapped_blocks_ratio"] = scores_df.loc["mapped_contigs"] / n_blocks -scores_df.loc["mapped_chim_blocks_ratio"] = scores_df.loc["chimeric_mapped_contigs"] / n_chimeric_blocks -scores_df.loc["mapped_contigs_ratio"] = (scores_df.loc["mapped_contigs"] + scores_df.loc["chimeric_mapped_contigs"])/ scores_df.loc["total_contigs"] -scores_df.loc["mapped_f1_score"] = (2 * (scores_df.loc["mapped_blocks_ratio"] + scores_df.loc["mapped_chim_blocks_ratio"]) * scores_df.loc["mapped_contigs_ratio"]) / (scores_df.loc["mapped_blocks_ratio"] + scores_df.loc["mapped_chim_blocks_ratio"] + scores_df.loc["mapped_contigs_ratio"]) - -scores_df.to_csv(sys.stdout, sep='\t') \ No newline at end of file diff --git a/bin/calculate_reference_info.py b/bin/calculate_reference_info.py index f09ad8b..e61b616 100755 --- a/bin/calculate_reference_info.py +++ b/bin/calculate_reference_info.py @@ -28,11 +28,17 @@ total_cds_length = sum(cds_lengths) pd.DataFrame({"Block Lengths": blocks["length"]}).to_csv(f"{prefix}_block_lengths.tsv", sep="\t", index=False) -pd.DataFrame({"Chimeric Blocks Lengths":chimeric_blocks["length"]}).to_csv(f"{prefix}_chimeric_block_lengths.tsv", sep="\t", index=False) +pd.DataFrame({"Chimeric Blocks Lengths": chimeric_blocks["length"]}).to_csv( + f"{prefix}_chimeric_block_lengths.tsv", sep="\t", index=False +) pd.DataFrame({"CDS Lengths": cds_lengths}).to_csv(f"{prefix}_cds_lengths.tsv", sep="\t", index=False) cds_summary = pd.concat([pd.Series(cds_lengths).describe(), pd.Series({"total_length": total_cds_length})]) blocks_summary = pd.concat([blocks["length"].describe(), pd.Series({"total_length": total_block_length})]) -chim_blocks_summary = pd.concat([chimeric_blocks["length"].describe(), pd.Series({"total_length": total_chimeric_block_length})]) +chim_blocks_summary = pd.concat( + [chimeric_blocks["length"].describe(), pd.Series({"total_length": total_chimeric_block_length})] +) -pd.DataFrame({'CDS': cds_summary, 'Blocks': blocks_summary, 'Chimeric Blocks': chim_blocks_summary}).to_csv(f"{prefix}_reference_stats.tsv", sep="\t") +pd.DataFrame({"CDS": cds_summary, "Blocks": blocks_summary, "Chimeric Blocks": chim_blocks_summary}).to_csv( + f"{prefix}_reference_stats.tsv", sep="\t" +) diff --git a/bin/calour_calculation.py b/bin/calour_calculation.py index c52a6e6..a455c4e 100755 --- a/bin/calour_calculation.py +++ b/bin/calour_calculation.py @@ -13,7 +13,9 @@ exp = ca.read(data_file=matrix, data_file_type="tsv", normalize=None) exp = exp.filter_sum_abundance(min_abundance) -exp.sample_metadata['Group'] = exp.sample_metadata['_sample_id'].apply(lambda x: [part for part in x.split('_') if 'group' in part][0]) +exp.sample_metadata["Group"] = exp.sample_metadata["_sample_id"].apply( + lambda x: [part for part in x.split("_") if "group" in part][0] +) res = exp.diff_abundance("Group", "group1", "group2") -res.feature_metadata.to_csv(f"{prefix}_calour_full_table.tsv", sep='\t', index=True) +res.feature_metadata.to_csv(f"{prefix}_calour_full_table.tsv", sep="\t", index=True) diff --git a/bin/categorize_contigs.py b/bin/categorize_contigs.py index a06680e..34d02d4 100755 --- a/bin/categorize_contigs.py +++ b/bin/categorize_contigs.py @@ -23,19 +23,27 @@ import polars as pl -all_contig_ids_fl = sys.argv[1] # "/vol/jlab/tlin/no_backup/nextflow_workdir/00/40ed8ff06c8b8b7cc7997b34a4eb99/rnaspades_contigs_ids.txt" +all_contig_ids_fl = sys.argv[ + 1 +] # "/vol/jlab/tlin/no_backup/nextflow_workdir/00/40ed8ff06c8b8b7cc7997b34a4eb99/rnaspades_contigs_ids.txt" -minimap2_categories_fl = sys.argv[2] # "/vol/jlab/tlin/no_backup/nextflow_workdir/42/62113b3dfc94ec10187ea81edf48a0/rnaspades_values.txt" +minimap2_categories_fl = sys.argv[ + 2 +] # "/vol/jlab/tlin/no_backup/nextflow_workdir/42/62113b3dfc94ec10187ea81edf48a0/rnaspades_values.txt" -chimeric_mapped_ids_fl = sys.argv[3] # "/vol/jlab/tlin/no_backup/nextflow_workdir/b7/51144ac49d14bec67d7b2a21adbda3/rnaspades_chimeric_values.txt" +length_filtered_ids_fl = sys.argv[ + 3 +] # "/vol/jlab/tlin/no_backup/nextflow_workdir/00/40ed8ff06c8b8b7cc7997b34a4eb99/rnaspades_length_filtered_ids.txt" -length_filtered_ids_fl = sys.argv[4] # "/vol/jlab/tlin/no_backup/nextflow_workdir/00/40ed8ff06c8b8b7cc7997b34a4eb99/rnaspades_length_filtered_ids.txt" +assembler_mapping_fl = sys.argv[ + 4 +] # "/vol/jlab/tlin/no_backup/nextflow_workdir/3c/f92325783c8fceb5ac01245489cc77/rnaspades_1000000.tsv" -assembler_mapping_fl = sys.argv[5] # "/vol/jlab/tlin/no_backup/nextflow_workdir/3c/f92325783c8fceb5ac01245489cc77/rnaspades_1000000.tsv" +gene_summary_path = sys.argv[ + 5 +] # "/vol/jlab/tlin/marbel_benchmarking_integration/benchmarking_sets_all_sparse_fixed_libsize/moss_microbiome/summary/gene_summary.csv" -gene_summary_path = sys.argv[6] # "/vol/jlab/tlin/marbel_benchmarking_integration/benchmarking_sets_all_sparse_fixed_libsize/moss_microbiome/summary/gene_summary.csv" - -prefix = sys.argv[7] # "rnaspades" +prefix = sys.argv[6] # "rnaspades" all_contigs_ids = pl.read_csv(all_contig_ids_fl, has_header=False, new_columns=["contigs"]) @@ -46,86 +54,81 @@ mapped_ids = [] try: - chimeric_mapped_ids = pl.read_csv(chimeric_mapped_ids_fl, has_header=False, new_columns=["chimeric_mapped"])["chimeric_mapped"].to_list() -except (pl.exceptions.NoDataError, OSError): - chimeric_mapped_ids = [] - -try: - length_filtered_ids = pl.read_csv(length_filtered_ids_fl, has_header=False, new_columns=["l_filtered"])["l_filtered"].to_list() + length_filtered_ids = pl.read_csv(length_filtered_ids_fl, has_header=False, new_columns=["l_filtered"])[ + "l_filtered" + ].to_list() except (pl.exceptions.NoDataError, OSError): length_filtered_ids = [] -assembler_mapping = pl.read_csv(assembler_mapping_fl, separator="\t", has_header=False, new_columns=["contigs", "start", "end", "read_id"]) +assembler_mapping = pl.read_csv( + assembler_mapping_fl, separator="\t", has_header=False, new_columns=["contigs", "start", "end", "read_id"] +) gene_summary = pl.read_csv(gene_summary_path).select(["gene_name", "orthogroup"]) gene_og_dict = dict(zip(gene_summary["gene_name"], gene_summary["orthogroup"])) # Get original gene name and orthogroup from read name -assembler_mapping = assembler_mapping.with_columns([ - pl.col("read_id").str.extract(r"^([^_]+_[^_]+)", 1).alias("origin_gene") -]) +assembler_mapping = assembler_mapping.with_columns( + [pl.col("read_id").str.extract(r"^([^_]+_[^_]+)", 1).alias("origin_gene")] +) # Map gene_og_dict -assembler_mapping = assembler_mapping.with_columns([ - pl.col("origin_gene").map_elements(lambda x: gene_og_dict.get(x, None),return_dtype=pl.String).alias("origin_orthogroup") -]) +assembler_mapping = assembler_mapping.with_columns( + [ + pl.col("origin_gene") + .map_elements(lambda x: gene_og_dict.get(x, None), return_dtype=pl.String) + .alias("origin_orthogroup") + ] +) -contig_aggregation = assembler_mapping.group_by("contigs").agg([ - pl.col("origin_gene").n_unique().alias("origin_gene_nunique"), - pl.col("origin_orthogroup").n_unique().alias("origin_orthogroup_nunique") -]) +contig_aggregation = assembler_mapping.group_by("contigs").agg( + [ + pl.col("origin_gene").n_unique().alias("origin_gene_nunique"), + pl.col("origin_orthogroup").n_unique().alias("origin_orthogroup_nunique"), + ] +) -single_mapped_contigs = contig_aggregation.filter( - pl.col("origin_gene_nunique") == 1 -)["contigs"].to_list() +single_mapped_contigs = contig_aggregation.filter(pl.col("origin_gene_nunique") == 1)["contigs"].to_list() # Step 1: Filter for contigs with multiple origin genes multi_mapped = contig_aggregation.filter(pl.col("origin_gene_nunique") != 1) # multi mapped with multiple ogs -multi_mapped_contigs_multi_og = multi_mapped.filter( - pl.col("origin_orthogroup_nunique") != 1 -)["contigs"].to_list() +multi_mapped_contigs_multi_og = multi_mapped.filter(pl.col("origin_orthogroup_nunique") != 1)["contigs"].to_list() # multimapped but single og -multi_mapped_contigs_single_og = multi_mapped.filter( - pl.col("origin_orthogroup_nunique") == 1 -)["contigs"].to_list() +multi_mapped_contigs_single_og = multi_mapped.filter(pl.col("origin_orthogroup_nunique") == 1)["contigs"].to_list() + +unassigned_contigs = all_contigs_ids.filter(~pl.col("contigs").is_in(mapped_ids)) assembler_contigs = assembler_mapping["contigs"].to_list() -unmapped_contigs = all_contigs_ids.filter( - ~pl.col("contigs").is_in(assembler_contigs) -)["contigs"].to_list() - -contigs_with_cat = all_contigs_ids.with_columns([ - pl.when(pl.col("contigs").is_in(mapped_ids)) - .then(pl.lit("mapped_contigs")) - .when(pl.col("contigs").is_in(chimeric_mapped_ids)) - .then(pl.lit("chimeric_mapped_contigs")) - .when(pl.col("contigs").is_in(length_filtered_ids)) - .then(pl.lit("length_filtered_contigs")) - .when(pl.col("contigs").is_in(unmapped_contigs)) - .then(pl.lit("unmapped_contigs")) - .when(pl.col("contigs").is_in(single_mapped_contigs)) - .then(pl.lit("single_mapped_contigs")) - .when(pl.col("contigs").is_in(multi_mapped_contigs_multi_og)) - .then(pl.lit("multi_mapped_contigs_multi_og")) - .when(pl.col("contigs").is_in(multi_mapped_contigs_single_og)) - .then(pl.lit("multi_mapped_contigs_single_og")) - .otherwise(pl.lit("no_category")) - .alias("category") -]) - -contigs_with_cat_no_mapped = contigs_with_cat.filter( - pl.col("category") != "mapped_contigs" -) + +unmapped_contigs = unassigned_contigs.filter(~pl.col("contigs").is_in(assembler_contigs))["contigs"].to_list() + +contigs_with_cat_no_mapped = unassigned_contigs.with_columns( + [ + pl.when(pl.col("contigs").is_in(length_filtered_ids)) + .then(pl.lit("length_filtered_contigs")) + .when(pl.col("contigs").is_in(unmapped_contigs)) + .then(pl.lit("unmapped_contigs")) + .when(pl.col("contigs").is_in(single_mapped_contigs)) + .then(pl.lit("single_mapped_contigs")) + .when(pl.col("contigs").is_in(multi_mapped_contigs_multi_og)) + .then(pl.lit("multi_mapped_contigs_multi_og")) + .when(pl.col("contigs").is_in(multi_mapped_contigs_single_og)) + .then(pl.lit("multi_mapped_contigs_single_og")) + .otherwise(pl.lit("no_category")) + .alias("category") + ] +) print("before concat") result = pl.concat( - [contigs_with_cat_no_mapped, minimap2_categories.select(["contig", "category"]).rename({"contig":"contigs"})], how="vertical" + [contigs_with_cat_no_mapped, minimap2_categories.select(["contig", "category"]).rename({"contig": "contigs"})], + how="vertical", ) -result.write_csv(f"{prefix}_contigs_categorised.tsv", separator="\t") \ No newline at end of file +result.write_csv(f"{prefix}_contigs_categorised.tsv", separator="\t") diff --git a/bin/check_samplesheet.py b/bin/check_samplesheet.py index 66cf3a8..51bf102 100755 --- a/bin/check_samplesheet.py +++ b/bin/check_samplesheet.py @@ -9,6 +9,7 @@ logger = logging.getLogger() + def read_head(handle, num_lines=10): """Read the specified number of lines from the current position in the file.""" lines = [] @@ -55,11 +56,12 @@ def _validate_sample(self, row): row[self._sample_col] = row[self._sample_col].replace(" ", "_") def _validate_assemblers(self, row): - assemblers = row[self._assemblers_col].split(';') + assemblers = row[self._assemblers_col].split(";") if not assemblers: raise AssertionError("At least one assembler is required.") row[self._assemblers_col] = assemblers + def check_samplesheet(file_in, file_out): required_columns = {"sample", "assemblers"} with file_in.open(newline="") as in_handle: @@ -86,9 +88,10 @@ def check_samplesheet(file_in, file_out): writer = csv.DictWriter(out_handle, header, delimiter=",") writer.writeheader() for row in checker.modified: - row['assemblers'] = ';'.join(row['assemblers']) + row["assemblers"] = ";".join(row["assemblers"]) writer.writerow(row) + def parse_args(argv=None): """Define and immediately parse command line arguments.""" parser = argparse.ArgumentParser( @@ -127,6 +130,7 @@ def main(argv=None): args.file_out.parent.mkdir(parents=True, exist_ok=True) check_samplesheet(args.file_in, args.file_out) + # Rest of the code remains the same. if __name__ == "__main__": diff --git a/bin/check_samplesheet_reads.py b/bin/check_samplesheet_reads.py index 10301b5..27b604c 100755 --- a/bin/check_samplesheet_reads.py +++ b/bin/check_samplesheet_reads.py @@ -2,7 +2,6 @@ """Provide a command line tool to validate and transform tabular samplesheets.""" - import argparse import csv import logging @@ -29,12 +28,7 @@ class RowChecker: ) def __init__( - self, - sample_col="sample", - first_col="fastq_1", - second_col="fastq_2", - single_col="single_end", - **kwargs + self, sample_col="sample", first_col="fastq_1", second_col="fastq_2", single_col="single_end", **kwargs ): """ Initialize the row checker with the expected column names. diff --git a/bin/count_genes_and_ogs_per_contig.py b/bin/count_genes_and_ogs_per_contig.py index a2ba465..f1c73be 100755 --- a/bin/count_genes_and_ogs_per_contig.py +++ b/bin/count_genes_and_ogs_per_contig.py @@ -12,17 +12,17 @@ gene_dict = gene_summary.set_index("gene_name")["orthogroup"].to_dict() -mapping_result["origin_gene"] = mapping_result[3].apply(lambda x: re.sub(r"(.*?_.*?)_.*", r"\1", x)) +mapping_result["origin_gene"] = mapping_result[3].apply(lambda x: re.sub(r"(.*?_.*?)_.*", r"\1", x)) mapping_result["orthogroup"] = mapping_result["origin_gene"].apply(lambda x: gene_dict[x]) mapping_result["contig_origin_gene"] = mapping_result[0].astype(str) + "_" + mapping_result["origin_gene"].astype(str) filtered = mapping_result.drop_duplicates(subset="contig_origin_gene") -#TODO: maybe also add the actual names of the origin cds and the orthogroups +# TODO: maybe also add the actual names of the origin cds and the orthogroups result = mapping_result.groupby(0).agg( origin_gene_count=("origin_gene", "count"), origin_gene_nunique=("origin_gene", "nunique"), - orthogroup_nunique=("orthogroup", "nunique") + orthogroup_nunique=("orthogroup", "nunique"), ) result["origin_genes"] = result.index.map(lambda x: ",".join(filtered[filtered[0] == x]["origin_gene"].to_list())) result["origin_og"] = result.index.map(lambda x: ",".join(set(filtered[filtered[0] == x]["orthogroup"].to_list()))) diff --git a/bin/gather_results.py b/bin/gather_results.py index d2481fa..66be7af 100755 --- a/bin/gather_results.py +++ b/bin/gather_results.py @@ -8,12 +8,11 @@ all_contigs_ids = sys.argv[1] minimap2_categories = sys.argv[2] -chimeric_mapped_ids = sys.argv[3] -length_filtered_ids = sys.argv[4] -assembler_mapping = sys.argv[5] -gene_summary = sys.argv[6] -contig_fasta = sys.argv[7] -prefix = sys.argv[8] +length_filtered_ids = sys.argv[3] +assembler_mapping = sys.argv[4] +gene_summary = sys.argv[5] +contig_fasta = sys.argv[6] +prefix = sys.argv[7] all_contigs_ids = pd.read_csv(all_contigs_ids, header=None).iloc[:, 0] @@ -23,10 +22,6 @@ mapped_ids = pd.read_csv(minimap2_categories, sep="\t", index_col=0)["contig"] except pd.errors.EmptyDataError: mapped_ids = pd.Series([]) -try: - chimeric_mapped_ids = pd.read_csv(chimeric_mapped_ids, header=None).iloc[:, 0] -except pd.errors.EmptyDataError: - chimeric_mapped_ids = pd.Series([]) try: length_filtered_ids = pd.read_csv(length_filtered_ids, header=None).iloc[:, 0] except pd.errors.EmptyDataError: @@ -35,36 +30,34 @@ gene_summary = pd.read_csv(gene_summary, usecols=["gene_name", "orthogroup"]).set_index("gene_name") gene_og_dict = gene_summary.to_dict()["orthogroup"] -# remove ids that are already mapped -chimeric_mapped_ids = chimeric_mapped_ids[~chimeric_mapped_ids.isin(mapped_ids)] +# first gather results -chimeric_mapped_ids_double_mapped = chimeric_mapped_ids[chimeric_mapped_ids.isin(mapped_ids)] +# add full index? -if len(chimeric_mapped_ids_double_mapped) > 0: - print("Some contigs are double mapped, are removed from summary") - print(chimeric_mapped_ids_double_mapped) +# remove ids that are already mapped assembler_mapping = pd.read_csv(assembler_mapping, sep="\t", header=None) -all_contigs_ids = all_contigs_ids[~all_contigs_ids.isin(pd.concat([mapped_ids, chimeric_mapped_ids, length_filtered_ids]))] +all_contigs_ids = all_contigs_ids[~all_contigs_ids.isin(pd.concat([mapped_ids, length_filtered_ids]))] unmapped_contigs = all_contigs_ids[~all_contigs_ids.isin(assembler_mapping[0].unique())] assembler_mapping = assembler_mapping[assembler_mapping[0].isin(all_contigs_ids)] -assembler_mapping["origin_gene"] = assembler_mapping[3].apply(lambda x: re.sub(r"(.*?_.*?)_.*", r"\1", x)) +assembler_mapping["origin_gene"] = assembler_mapping[3].apply(lambda x: re.sub(r"(.*?_.*?)_.*", r"\1", x)) assembler_mapping["origin_orthogroup"] = assembler_mapping["origin_gene"].apply(lambda x: gene_og_dict.get(x)) contig_aggregation = assembler_mapping.groupby(0).agg( - origin_gene_nunique=("origin_gene", "nunique"), - origin_orthogroup_nunique=("origin_orthogroup", "nunique") + origin_gene_nunique=("origin_gene", "nunique"), origin_orthogroup_nunique=("origin_orthogroup", "nunique") ) single_mapped_contigs = contig_aggregation[contig_aggregation["origin_gene_nunique"] == 1].index multi_mapped_contigs = contig_aggregation[contig_aggregation["origin_gene_nunique"] != 1].index -number_ogs_per_multimapped = contig_aggregation[contig_aggregation["origin_gene_nunique"] != 1]["origin_orthogroup_nunique"].value_counts() +number_ogs_per_multimapped = contig_aggregation[contig_aggregation["origin_gene_nunique"] != 1][ + "origin_orthogroup_nunique" +].value_counts() multi_og_contigs = number_ogs_per_multimapped[number_ogs_per_multimapped.index != 1].sum() single_og_contigs = number_ogs_per_multimapped[number_ogs_per_multimapped.index == 1].sum() @@ -88,20 +81,77 @@ unmpapped_contig_lengths_multi.append(record_length) -result_dict = {"total_contigs": len_all_contigs_ids, "total_bases": sum(mapped_contig_lengths) + sum(unmpapped_contig_lengths_single) + sum(unmpapped_contig_lengths_multi), - "mapped_contigs": len(mapped_ids), "chimeric_mapped_contigs": len(chimeric_mapped_ids), "length_filtered_contigs": len(length_filtered_ids), - "unmapped_contigs": len(unmapped_contigs), "multi_mapped_contigs": len(multi_mapped_contigs), "multi_mapped_contigs_single_og": - single_og_contigs, "multi_mapped_contigs_multi_og": multi_og_contigs, "single_mapped_contigs": len(single_mapped_contigs), - "mapped_contig_bases": sum(mapped_contig_lengths), "unmapped_contig_bases": sum(unmpapped_contig_lengths_single) + sum(unmpapped_contig_lengths_multi),} +result_dict = { + "total_contigs": len_all_contigs_ids, + "total_bases": sum(mapped_contig_lengths) + + sum(unmpapped_contig_lengths_single) + + sum(unmpapped_contig_lengths_multi), + "mapped_contigs": len(mapped_ids), + "length_filtered_contigs": len(length_filtered_ids), + "unmapped_contigs": len(unmapped_contigs), + "multi_mapped_contigs": len(multi_mapped_contigs), + "multi_mapped_contigs_single_og": single_og_contigs, + "multi_mapped_contigs_multi_og": multi_og_contigs, + "single_mapped_contigs": len(single_mapped_contigs), + "mapped_contig_bases": sum(mapped_contig_lengths), + "unmapped_contig_bases": sum(unmpapped_contig_lengths_single) + sum(unmpapped_contig_lengths_multi), +} pd.DataFrame.from_dict(result_dict, orient="index", columns=[prefix]).to_csv(f"{prefix}_scores.tsv", sep="\t") -pd.DataFrame({prefix: mapped_contig_lengths + unmpapped_contig_lengths_single + unmpapped_contig_lengths_multi}).to_csv(f"{prefix}_all_contig_lengths.tsv", sep="\t", index=False) +pd.DataFrame({prefix: mapped_contig_lengths + unmpapped_contig_lengths_single + unmpapped_contig_lengths_multi}).to_csv( + f"{prefix}_all_contig_lengths.tsv", sep="\t", index=False +) pd.DataFrame({prefix: mapped_contig_lengths}).to_csv(f"{prefix}_mapped_contig_lengths.tsv", sep="\t", index=False) -pd.DataFrame({prefix: unmpapped_contig_lengths_single}).to_csv(f"{prefix}_unmapped_contig_lengths_single.tsv", sep="\t", index=False) -pd.DataFrame({prefix: unmpapped_contig_lengths_multi}).to_csv(f"{prefix}_unmapped_contig_lengths_multi.tsv", sep="\t", index=False) +pd.DataFrame({prefix: unmpapped_contig_lengths_single}).to_csv( + f"{prefix}_unmapped_contig_lengths_single.tsv", sep="\t", index=False +) +pd.DataFrame({prefix: unmpapped_contig_lengths_multi}).to_csv( + f"{prefix}_unmapped_contig_lengths_multi.tsv", sep="\t", index=False +) -pd.DataFrame({prefix: pd.concat([pd.Series(mapped_contig_lengths + unmpapped_contig_lengths_single + unmpapped_contig_lengths_multi).describe(), pd.Series({"total_length": sum(mapped_contig_lengths + unmpapped_contig_lengths_single + unmpapped_contig_lengths_multi)})])}).to_csv(f"{prefix}_all_contig_lengths_summary.tsv", sep="\t") -pd.DataFrame({prefix: pd.concat([pd.Series(mapped_contig_lengths).describe(), pd.Series({"total_length": sum(mapped_contig_lengths)})])}).to_csv(f"{prefix}_mapped_contig_lengths_summary.tsv", sep="\t") -pd.DataFrame({prefix: pd.concat([pd.Series(unmpapped_contig_lengths_multi).describe(), pd.Series({"total_length": sum(unmpapped_contig_lengths_multi)})])}).to_csv(f"{prefix}_unmapped_contig_lengths_multi_summary.tsv", sep="\t") -pd.DataFrame({prefix: pd.concat([pd.Series(unmpapped_contig_lengths_single).describe(), pd.Series({"total_length": sum(unmpapped_contig_lengths_single)})])}).to_csv(f"{prefix}_unmapped_contig_lengths_single_summary.tsv", sep="\t") \ No newline at end of file +pd.DataFrame( + { + prefix: pd.concat( + [ + pd.Series( + mapped_contig_lengths + unmpapped_contig_lengths_single + unmpapped_contig_lengths_multi + ).describe(), + pd.Series( + { + "total_length": sum( + mapped_contig_lengths + unmpapped_contig_lengths_single + unmpapped_contig_lengths_multi + ) + } + ), + ] + ) + } +).to_csv(f"{prefix}_all_contig_lengths_summary.tsv", sep="\t") +pd.DataFrame( + { + prefix: pd.concat( + [pd.Series(mapped_contig_lengths).describe(), pd.Series({"total_length": sum(mapped_contig_lengths)})] + ) + } +).to_csv(f"{prefix}_mapped_contig_lengths_summary.tsv", sep="\t") +pd.DataFrame( + { + prefix: pd.concat( + [ + pd.Series(unmpapped_contig_lengths_multi).describe(), + pd.Series({"total_length": sum(unmpapped_contig_lengths_multi)}), + ] + ) + } +).to_csv(f"{prefix}_unmapped_contig_lengths_multi_summary.tsv", sep="\t") +pd.DataFrame( + { + prefix: pd.concat( + [ + pd.Series(unmpapped_contig_lengths_single).describe(), + pd.Series({"total_length": sum(unmpapped_contig_lengths_single)}), + ] + ) + } +).to_csv(f"{prefix}_unmapped_contig_lengths_single_summary.tsv", sep="\t") diff --git a/bin/generate_blocks.py b/bin/generate_blocks.py index ef05393..87f37d8 100755 --- a/bin/generate_blocks.py +++ b/bin/generate_blocks.py @@ -18,18 +18,18 @@ bed_df = pd.read_csv(bed_file, sep="\t", header=None) -#the read name contains the read numbered, so I need to truncate to see if it is the orginal cds -bed_df['truncated_red_name'] = bed_df[3].str.split('_', expand=True)[[0,1]].agg('_'.join, axis=1) -#we deduplicated the sequences, but we need to include the reads of sequences that were deduplicated, so map to all the original reads -bed_df["mapped_group"] = bed_df['truncated_red_name'].map(lambda x: dict_val.get(x, [x])) +# the read name contains the read numbered, so I need to truncate to see if it is the orginal cds +bed_df["truncated_red_name"] = bed_df[3].str.split("_", expand=True)[[0, 1]].agg("_".join, axis=1) +# we deduplicated the sequences, but we need to include the reads of sequences that were deduplicated, so map to all the original reads +bed_df["mapped_group"] = bed_df["truncated_red_name"].map(lambda x: dict_val.get(x, [x])) bed_df = bed_df.explode("mapped_group") -#filter reads that are not the original cds, i.e. incorrectly mapped -bed_df = bed_df[bed_df[0]==bed_df["mapped_group"]] +# filter reads that are not the original cds, i.e. incorrectly mapped +bed_df = bed_df[bed_df[0] == bed_df["mapped_group"]] blocks = [] for group_name, group in bed_df.groupby(0): - #sort by start position + # sort by start position sorted_group = group.sort_values(1) new_cds = True block_index = 0 @@ -49,13 +49,17 @@ current_end = row[2] block_index += 1 fragment_count = 1 - #add the last block + # add the last block blocks.append((group_name, block_index, current_start, current_end, fragment_count)) blocks_df = pd.DataFrame(blocks, columns=["cds", "block_index", "start", "end", "fragment_count"]) blocks_df.to_csv(f"{prefix}_blocks.tsv", sep="\t", index=False) -cds_to_block_dict = blocks_df.groupby("cds").apply(lambda group: group[["block_index", "start", "end"]].values.tolist(), include_groups=False).to_dict() +cds_to_block_dict = ( + blocks_df.groupby("cds") + .apply(lambda group: group[["block_index", "start", "end"]].values.tolist(), include_groups=False) + .to_dict() +) in_data_set = 0 not_in_data_set = 0 @@ -65,7 +69,7 @@ for record in SeqIO.parse(ref_transcriptome, "fasta"): if record.id in cds_to_block_dict: for block in cds_to_block_dict[record.id]: - block_seq = record.seq[block[1]:(block[2]+1)] + block_seq = record.seq[block[1] : (block[2] + 1)] block_name = f"{record.id}_block{block[0]}" block_records.append(SeqIO.SeqRecord(block_seq, block_name, description="")) diff --git a/bin/make_histogramms.py b/bin/make_histogramms.py index 4319df6..33b8ead 100755 --- a/bin/make_histogramms.py +++ b/bin/make_histogramms.py @@ -51,4 +51,3 @@ plt.tight_layout() plt.savefig(f"{prefix}_log_histograms.png", dpi=900) plt.close() - diff --git a/bin/merge_assembler_og_counts.py b/bin/merge_assembler_og_counts.py index f515d09..9aff171 100755 --- a/bin/merge_assembler_og_counts.py +++ b/bin/merge_assembler_og_counts.py @@ -5,11 +5,10 @@ import re import sys - -assembler_count_dir_path = sys.argv[1] # mergequantsffiles -assembler_mapping_dir = sys.argv[2] # minimap2filter #{assembler}_mapping_dedup.json" +assembler_count_dir_path = sys.argv[1] # mergequantsffiles +assembler_mapping_dir = sys.argv[2] # minimap2filter #{assembler}_mapping_dedup.json" gene_summary_file = sys.argv[3] -assembler = sys.argv[4] # "rnaspades" +assembler = sys.argv[4] # "rnaspades" assembler_counts = pl.read_csv(assembler_count_dir_path, separator="\t") @@ -22,20 +21,27 @@ assembler_samples = [col for col in assembler_counts.columns if re.search(r"_sample_", col)] -# now i need to map the assemblers contig -> gene than gene -> og than aggregate -assembler_with_og = assembler_counts.with_columns( - pl.col("Name").cast(pl.Utf8), -).with_columns( - pl.col("Name").map_elements(lambda x: assembler_mapping_dict.get(str(x), None), return_dtype=pl.Utf8).alias("gene_block_name"), -).with_columns( - pl.col("gene_block_name").str.replace(r"_block\d+$", "", literal=False).alias("gene_name") -).with_columns( - pl.col("gene_name").map_elements(lambda x: gene_to_og.get(str(x), None), return_dtype=pl.Utf8).alias("orthogroups") -).group_by("orthogroups").agg((pl.col(assembler_samples).sum())) +# now i need to map the assemblers contig -> gene than gene -> og than aggregate +assembler_with_og = ( + assembler_counts.with_columns( + pl.col("Name").cast(pl.Utf8), + ) + .with_columns( + pl.col("Name") + .map_elements(lambda x: assembler_mapping_dict.get(str(x), None), return_dtype=pl.Utf8) + .alias("gene_block_name"), + ) + .with_columns(pl.col("gene_block_name").str.replace(r"_block\d+$", "", literal=False).alias("gene_name")) + .with_columns( + pl.col("gene_name") + .map_elements(lambda x: gene_to_og.get(str(x), None), return_dtype=pl.Utf8) + .alias("orthogroups") + ) + .group_by("orthogroups") + .agg((pl.col(assembler_samples).sum())) +) # prepare count file for dge tools assembler_with_og.rename({"orthogroups": "Name"}) assembler_count_file = f"{assembler}_merged_orthogroups.tsv" assembler_with_og.write_csv(assembler_count_file, separator="\t") - - diff --git a/bin/merge_bowtie2_logs.py b/bin/merge_bowtie2_logs.py index f1d68bb..8714bca 100755 --- a/bin/merge_bowtie2_logs.py +++ b/bin/merge_bowtie2_logs.py @@ -15,11 +15,13 @@ def merge_bowtie2_logs(log_files, column_name): log_dfs = pd.merge(log_dfs, log_df, left_index=True, right_index=True) merged_df = log_dfs.mean(axis=1) - pd.DataFrame(merged_df, columns=[f"{column_name}_mean"]).to_csv(f"{column_name}_mean_logs.tsv", sep="\t", index=True) + pd.DataFrame(merged_df, columns=[f"{column_name}_mean"]).to_csv( + f"{column_name}_mean_logs.tsv", sep="\t", index=True + ) log_dfs.to_csv(f"{column_name}_merged_logs.tsv", sep="\t", index=True) len_args = len(sys.argv) -log_files = sys.argv[1:len_args-1] -col_name = sys.argv[len_args-1] +log_files = sys.argv[1 : len_args - 1] +col_name = sys.argv[len_args - 1] merge_bowtie2_logs(log_files, col_name) diff --git a/bin/merge_dataframes.py b/bin/merge_dataframes.py index bf2b754..1d27d2f 100755 --- a/bin/merge_dataframes.py +++ b/bin/merge_dataframes.py @@ -16,6 +16,7 @@ def merge_dfs(df_files): dfs = dfs.fillna(0) dfs.to_csv(sys.stdout, sep="\t", index=True) + len_args = len(sys.argv) df_files = sys.argv[1:len_args] merge_dfs(df_files) diff --git a/bin/merge_og_counts.py b/bin/merge_og_counts.py index 0a5d657..64b0592 100755 --- a/bin/merge_og_counts.py +++ b/bin/merge_og_counts.py @@ -14,9 +14,7 @@ samples = [col for col in mt_genes.columns if re.search(r"_sample_", col)] # aggregate by orthogroups -mt_ogs_aggregated = mt_genes.group_by("orthogroup").agg( - (pl.col(samples).sum()) -) +mt_ogs_aggregated = mt_genes.group_by("orthogroup").agg((pl.col(samples).sum())) mt_ogs_aggregated = mt_ogs_aggregated.rename({"orthogroup": "Name"}) args = "group_1 10 group_2 10 \t".split(" ") diff --git a/bin/minimap2_classification.py b/bin/minimap2_classification.py index c0317f3..c0d1a79 100755 --- a/bin/minimap2_classification.py +++ b/bin/minimap2_classification.py @@ -13,7 +13,7 @@ def deduplicate_by_key(df, group_key, seed): max_cols = ["Mapping quality", "Number of matching bases in the mapping"] new_max_cols = ["Max Mapping quality", "Max Number of matching bases in the mapping"] - minimap_mapping[new_max_cols] = (minimap_mapping.groupby(group)[max_cols].transform("max")) + minimap_mapping[new_max_cols] = minimap_mapping.groupby(group)[max_cols].transform("max") min_col = ["Query sequence length"] new_min_col = ["Min Query sequence length"] @@ -21,33 +21,34 @@ def deduplicate_by_key(df, group_key, seed): minimap_mapping[new_min_col] = minimap_mapping.groupby(group)[min_col].transform("min") # calculate a cumulative score to include hits with less then optimal matches - minimap_mapping["cumulative_score"] = 5 * (minimap_mapping[max_cols[0]] == minimap_mapping[new_max_cols[0]]) + 3 * (minimap_mapping[max_cols[1]] == minimap_mapping[new_max_cols[1]]) + 1* (minimap_mapping[min_col[0]] == minimap_mapping[new_min_col[0]]) + minimap_mapping["cumulative_score"] = ( + 5 * (minimap_mapping[max_cols[0]] == minimap_mapping[new_max_cols[0]]) + + 3 * (minimap_mapping[max_cols[1]] == minimap_mapping[new_max_cols[1]]) + + 1 * (minimap_mapping[min_col[0]] == minimap_mapping[new_min_col[0]]) + ) minimap_mapping["max_cumulative_score"] = minimap_mapping.groupby(group)["cumulative_score"].transform("max") - minimap_mapping["optimal_hit"] = (minimap_mapping["cumulative_score"] == minimap_mapping["max_cumulative_score"]) + minimap_mapping["optimal_hit"] = minimap_mapping["cumulative_score"] == minimap_mapping["max_cumulative_score"] - cooptimal_bools = (minimap_mapping[minimap_mapping["optimal_hit"]].groupby(group).transform("size") > 1) + cooptimal_bools = minimap_mapping[minimap_mapping["optimal_hit"]].groupby(group).transform("size") > 1 minimap_mapping.loc[cooptimal_bools.index, "cooptimal_hits"] = cooptimal_bools - minimap_mapping["cooptimal_hits"] = ( - minimap_mapping["cooptimal_hits"] - .astype("boolean") - .fillna(False) - .astype(bool) - ) - cooptimal_hits_selection_idx = minimap_mapping[minimap_mapping["cooptimal_hits"]].groupby(group).sample(n=1, random_state=seed).index + minimap_mapping["cooptimal_hits"] = minimap_mapping["cooptimal_hits"].astype("boolean").fillna(False).astype(bool) + cooptimal_hits_selection_idx = ( + minimap_mapping[minimap_mapping["cooptimal_hits"]].groupby(group).sample(n=1, random_state=seed).index + ) minimap_mapping.loc[cooptimal_hits_selection_idx, "selected_cooptimal_hit"] = True minimap_mapping["selected_cooptimal_hit"] = ( - minimap_mapping["selected_cooptimal_hit"] - .astype("boolean") - .fillna(False) - .astype(bool) - ) - selection_idx = minimap_mapping[minimap_mapping["optimal_hit"] & (~minimap_mapping["cooptimal_hits"] | minimap_mapping["selected_cooptimal_hit"])].index + minimap_mapping["selected_cooptimal_hit"].astype("boolean").fillna(False).astype(bool) + ) + selection_idx = minimap_mapping[ + minimap_mapping["optimal_hit"] + & (~minimap_mapping["cooptimal_hits"] | minimap_mapping["selected_cooptimal_hit"]) + ].index return df.loc[selection_idx] @@ -76,10 +77,11 @@ def count_ns_in_contigs(contigs_path): minimap_path = sys.argv[1] -contig_path = sys.argv[2] -gene_summary_path = sys.argv[3] -seed = int(sys.argv[4]) -prefix = sys.argv[5] +minimap_overlap_path = sys.argv[2] +contig_path = sys.argv[3] +gene_summary_path = sys.argv[4] +seed = int(sys.argv[5]) +prefix = sys.argv[6] mapping_col_names = [ "contig", # Query sequence name @@ -93,49 +95,53 @@ def count_ns_in_contigs(contigs_path): "Target end coordinate on the original strand", "Number of matching bases in the mapping", "Number bases, including gaps, in the mapping", - "Mapping quality", "1", "2", "3", "4", "5", "6" + "Mapping quality", + "1", + "2", + "3", + "4", + "5", + "6", ] json_cols = ["contig", "block_id"] -minimap_mapping = pd.read_csv( - minimap_path, sep="\t", header=None, names=mapping_col_names, dtype={'contig': str} -) +minimap_mapping = pd.read_csv(minimap_path, sep="\t", header=None, names=mapping_col_names, dtype={"contig": str}) + +minimap_overlap_df = pd.read_csv(minimap_overlap_path, sep="\t", dtype={"contig": str}) gene_summary = pd.read_csv(gene_summary_path) gene_to_og = gene_summary.set_index("gene_name")["orthogroup"].to_dict() ns_per_contig_df = count_ns_in_contigs(contig_path) -minimap_mapping["correctly_mapped_bases"] = minimap_mapping["Number of matching bases in the mapping"] / minimap_mapping["Target sequence length"] +minimap_mapping["correctly_mapped_bases"] = ( + minimap_mapping["Number of matching bases in the mapping"] / minimap_mapping["Target sequence length"] +) # -> "Number of matching bases in the mapping"/ "Target sequence length" (blocks are targets) # Rationale: Number of matching bases/ length of block -> block bases should be reconstructed -minimap_mapping["covered_reference"] = minimap_mapping["Number bases, including gaps, in the mapping"] / minimap_mapping["Target sequence length"] +minimap_mapping["covered_reference"] = ( + minimap_mapping["Number bases, including gaps, in the mapping"] / minimap_mapping["Target sequence length"] +) pre_filter_mappings = minimap_mapping.shape[0] minimap_mapping = minimap_mapping[minimap_mapping["correctly_mapped_bases"] >= 0.95] minimap_mapping = minimap_mapping[ - (minimap_mapping["covered_reference"] <= 1.05) & - (minimap_mapping["covered_reference"] >= 0.95) + (minimap_mapping["covered_reference"] <= 1.05) & (minimap_mapping["covered_reference"] >= 0.95) ] +minimap_mapping = minimap_mapping[~minimap_mapping["contig"].isin(minimap_overlap_df["contig"])] + minimap_mapping = minimap_mapping.reset_index(drop=True) -minimap_mapping["block_size"] = ( - minimap_mapping - .groupby(["block_id"]) - .transform("size") -) +minimap_mapping["block_size"] = minimap_mapping.groupby(["block_id"]).transform("size") -minimap_mapping["contig_size"] = ( - minimap_mapping - .groupby(["contig"]) - .transform("size") -) +minimap_mapping["contig_size"] = minimap_mapping.groupby(["contig"]).transform("size") minimap_mapping = minimap_mapping.sort_values( - by=["Mapping quality", "Number of matching bases in the mapping", "Query sequence length"], ascending=[False, False, True] + by=["Mapping quality", "Number of matching bases in the mapping", "Query sequence length"], + ascending=[False, False, True], ) minimap_mapping_chosen = minimap_mapping.copy() @@ -148,142 +154,188 @@ def count_ns_in_contigs(contigs_path): minimap_mapping.loc[minimap_mapping_chosen.index, "assigned"] = True -minimap_mapping["assigned"] = ( - minimap_mapping["assigned"] - .astype("boolean") - .fillna(False) - .astype(bool) -) +minimap_mapping["assigned"] = minimap_mapping["assigned"].astype("boolean").fillna(False).astype(bool) minimap_mapping["is_duplicated"] = (~(minimap_mapping["assigned"])) & (minimap_mapping["block_size"] > 1) -minimap_mapping["is_joined"] = (minimap_mapping["contig_size"] > 1) +minimap_mapping["is_joined"] = minimap_mapping["contig_size"] > 1 -minimap_mapping["blocks_aggregated"] = ( - minimap_mapping - .groupby("contig")["block_id"] - .transform(lambda x: [list(x)] * len(x)) +minimap_mapping["blocks_aggregated"] = minimap_mapping.groupby("contig")["block_id"].transform( + lambda x: [list(x)] * len(x) ) minimap_mapping["gene"] = minimap_mapping["block_id"].str.replace(r"^((?:[^_]+_){1}[^_]+).*", r"\1", regex=True) -minimap_mapping["genes_aggregated"] = ( - minimap_mapping - .groupby("contig")["gene"] - .transform(lambda x: [list(x)] * len(x)) -) +minimap_mapping["genes_aggregated"] = minimap_mapping.groupby("contig")["gene"].transform(lambda x: [list(x)] * len(x)) minimap_mapping["is_one_gene"] = minimap_mapping["genes_aggregated"].apply(lambda x: len(set(x)) == 1) minimap_mapping["og"] = minimap_mapping["gene"].map(gene_to_og) -minimap_mapping["ogs_aggregated"] = ( - minimap_mapping - .groupby("contig")["og"] - .transform(lambda x: [list(x)] * len(x)) -) +minimap_mapping["ogs_aggregated"] = minimap_mapping.groupby("contig")["og"].transform(lambda x: [list(x)] * len(x)) minimap_mapping["is_one_og"] = minimap_mapping["ogs_aggregated"].apply(lambda x: len(set(x)) == 1) -minimap_mapping["is_recovered"] = (minimap_mapping["assigned"] & (~minimap_mapping["is_joined"])) +minimap_mapping["is_recovered"] = minimap_mapping["assigned"] & (~minimap_mapping["is_joined"]) minimap_mapping["is_block_recovered"] = minimap_mapping.groupby("block_id")["is_recovered"].transform("max") -block_is_super_recovered = minimap_mapping[minimap_mapping["assigned"] & (minimap_mapping["is_joined"]) & (minimap_mapping["is_one_gene"])]["blocks_aggregated"].explode() +block_is_super_recovered = minimap_mapping[ + minimap_mapping["assigned"] & (minimap_mapping["is_joined"]) & (minimap_mapping["is_one_gene"]) +]["blocks_aggregated"].explode() minimap_mapping["is_block_super_recovered"] = minimap_mapping["block_id"].isin(block_is_super_recovered) -minimap_mapping["is_block_super_recovered_or_recovered"] = minimap_mapping["is_block_recovered"] | minimap_mapping["is_block_super_recovered"] +minimap_mapping["is_block_super_recovered_or_recovered"] = ( + minimap_mapping["is_block_recovered"] | minimap_mapping["is_block_super_recovered"] +) # calculate blocks but subtract the blocks already assigned -super_recovered_amount_of_blocks = minimap_mapping[minimap_mapping["assigned"] & (minimap_mapping["is_joined"]) & (minimap_mapping["is_one_gene"])]["contig_size"].sum() - (minimap_mapping["assigned"] & (minimap_mapping["is_joined"]) & ~(minimap_mapping["is_one_gene"]) & (minimap_mapping["is_one_og"]) & minimap_mapping["is_block_recovered"]).sum() +super_recovered_amount_of_blocks = ( + minimap_mapping[minimap_mapping["assigned"] & (minimap_mapping["is_joined"]) & (minimap_mapping["is_one_gene"])][ + "contig_size" + ].sum() + - ( + minimap_mapping["assigned"] + & (minimap_mapping["is_joined"]) + & ~(minimap_mapping["is_one_gene"]) + & (minimap_mapping["is_one_og"]) + & minimap_mapping["is_block_recovered"] + ).sum() +) -og_recovered_amount_of_blocks = minimap_mapping[minimap_mapping["assigned"] & (minimap_mapping["is_joined"]) & ~(minimap_mapping["is_one_gene"]) & (minimap_mapping["is_one_og"])]["contig_size"].sum() - (minimap_mapping["assigned"] & (minimap_mapping["is_joined"]) & ~(minimap_mapping["is_one_gene"]) & minimap_mapping["is_one_og"] & minimap_mapping["is_block_super_recovered_or_recovered"]).sum() +og_recovered_amount_of_blocks = ( + minimap_mapping[ + minimap_mapping["assigned"] + & (minimap_mapping["is_joined"]) + & ~(minimap_mapping["is_one_gene"]) + & (minimap_mapping["is_one_og"]) + ]["contig_size"].sum() + - ( + minimap_mapping["assigned"] + & (minimap_mapping["is_joined"]) + & ~(minimap_mapping["is_one_gene"]) + & minimap_mapping["is_one_og"] + & minimap_mapping["is_block_super_recovered_or_recovered"] + ).sum() +) -dedupped_rest = minimap_mapping[~minimap_mapping["contig"].isin(minimap_mapping[minimap_mapping["assigned"]]["contig"])].groupby("contig").head(1) +dedupped_rest = ( + minimap_mapping[~minimap_mapping["contig"].isin(minimap_mapping[minimap_mapping["assigned"]]["contig"])] + .groupby("contig") + .head(1) +) duplicated_non_chimeric_sum = (dedupped_rest["is_duplicated"] & ~dedupped_rest["is_joined"]).sum() duplicated_chimeric_sum = (dedupped_rest["is_duplicated"] & dedupped_rest["is_joined"]).sum() -recovered = ( - minimap_mapping[(minimap_mapping["assigned"] & (~minimap_mapping["is_joined"])) | (minimap_mapping["assigned"] & (minimap_mapping["is_joined"]) & (minimap_mapping["is_one_gene"]))] -) +recovered = minimap_mapping[ + (minimap_mapping["assigned"] & (~minimap_mapping["is_joined"])) + | (minimap_mapping["assigned"] & (minimap_mapping["is_joined"]) & (minimap_mapping["is_one_gene"])) +] # add final contig classification -og_recovered = ( - minimap_mapping[(minimap_mapping["assigned"] & (minimap_mapping["is_joined"]) & ~(minimap_mapping["is_one_gene"]) & (minimap_mapping["is_one_og"]))] -) +og_recovered = minimap_mapping[ + ( + minimap_mapping["assigned"] + & (minimap_mapping["is_joined"]) + & ~(minimap_mapping["is_one_gene"]) + & (minimap_mapping["is_one_og"]) + ) +] -chimeric = ( - minimap_mapping[minimap_mapping["assigned"] & (minimap_mapping["is_joined"]) & ~(minimap_mapping["is_one_gene"]) & ~(minimap_mapping["is_one_og"])] -) +chimeric = minimap_mapping[ + minimap_mapping["assigned"] + & (minimap_mapping["is_joined"]) + & ~(minimap_mapping["is_one_gene"]) + & ~(minimap_mapping["is_one_og"]) +] -duplicated_non_chimeric = ( - dedupped_rest[(dedupped_rest["is_duplicated"] & ~dedupped_rest["is_joined"])] -) +duplicated_non_chimeric = dedupped_rest[(dedupped_rest["is_duplicated"] & ~dedupped_rest["is_joined"])] -duplicated_chimeric = ( - dedupped_rest[(dedupped_rest["is_duplicated"] & dedupped_rest["is_joined"])] -) +duplicated_chimeric = dedupped_rest[(dedupped_rest["is_duplicated"] & dedupped_rest["is_joined"])] -minimap_mapping.loc[ - (minimap_mapping["assigned"] & (~minimap_mapping["is_joined"])), "category" -] = "minimap2_single_recovered" +minimap_mapping.loc[(minimap_mapping["assigned"] & (~minimap_mapping["is_joined"])), "category"] = ( + "minimap2_single_recovered" +) minimap_mapping.loc[ (minimap_mapping["assigned"] & (minimap_mapping["is_joined"]) & (minimap_mapping["is_one_gene"])), "category" -] = "minimap2_merged_recovered" +] = "minimap2_merged_recovered" -#super_recovered_amount_of_blocks = minimap_mapping[ -# minimap_mapping["assigned"] & (minimap_mapping["is_joined"]) & (minimap_mapping["is_one_gene"])]["contig_size"].sum() - +# super_recovered_amount_of_blocks = minimap_mapping[ +# minimap_mapping["assigned"] & (minimap_mapping["is_joined"]) & (minimap_mapping["is_one_gene"])]["contig_size"].sum() - # (minimap_mapping["assigned"] & (minimap_mapping["is_joined"]) & ~(minimap_mapping["is_one_gene"]) & (minimap_mapping["is_one_og"]) & minimap_mapping["is_block_recovered"]).sum() minimap_mapping.loc[ - (minimap_mapping["assigned"] & (minimap_mapping["is_joined"]) & ~(minimap_mapping["is_one_gene"]) & (minimap_mapping["is_one_og"])), "category" -] = "minimap2_orthologous_recovered" + ( + minimap_mapping["assigned"] + & (minimap_mapping["is_joined"]) + & ~(minimap_mapping["is_one_gene"]) + & (minimap_mapping["is_one_og"]) + ), + "category", +] = "minimap2_orthologous_recovered" minimap_mapping.loc[ - (minimap_mapping["assigned"] & (minimap_mapping["is_joined"]) & ~(minimap_mapping["is_one_gene"]) & ~(minimap_mapping["is_one_og"])), "category" + ( + minimap_mapping["assigned"] + & (minimap_mapping["is_joined"]) + & ~(minimap_mapping["is_one_gene"]) + & ~(minimap_mapping["is_one_og"]) + ), + "category", ] = "minimap2_chimeric" -minimap_mapping.loc[ - (dedupped_rest[dedupped_rest["is_duplicated"] & ~dedupped_rest["is_joined"]]).index, "category" -] = "minimap2_duplicated_non_chimeric" +minimap_mapping.loc[(dedupped_rest[dedupped_rest["is_duplicated"] & ~dedupped_rest["is_joined"]]).index, "category"] = ( + "minimap2_duplicated_non_chimeric" +) -minimap_mapping.loc[ - (dedupped_rest[dedupped_rest["is_duplicated"] & dedupped_rest["is_joined"]]).index, "category" -] = "minimap2_duplicated_chimeric" +minimap_mapping.loc[(dedupped_rest[dedupped_rest["is_duplicated"] & dedupped_rest["is_joined"]]).index, "category"] = ( + "minimap2_duplicated_chimeric" +) minimap_mapping.to_csv("debug.tsv", sep="\t") -minimap_mapping.dropna( - subset=["category"], inplace=True -) +minimap_mapping.dropna(subset=["category"], inplace=True) + +minimap_mapping = pd.concat([minimap_mapping, minimap_overlap_df], ignore_index=True) # calculate how many Ns are in the assigned contigs -contigs_n_df = pd.merge( - minimap_mapping, - ns_per_contig_df, - how='left', - right_on='contig_id', - left_on='contig' -) +contigs_n_df = pd.merge(minimap_mapping, ns_per_contig_df, how="left", right_on="contig_id", left_on="contig") -recovered_with_n = (minimap_mapping["category"].isin(["minimap2_single_recovered", "minimap2_merged_recovered", "minimap2_orthologous_recovered"]) & (contigs_n_df["total_Ns"] > 0)).sum() +recovered_with_n = ( + minimap_mapping["category"].isin( + ["minimap2_single_recovered", "minimap2_merged_recovered", "minimap2_orthologous_recovered"] + ) + & (contigs_n_df["total_Ns"] > 0) +).sum() chimeric_with_n = ((minimap_mapping["category"] == "minimap2_chimeric") & (contigs_n_df["total_Ns"] > 0)).sum() -minimap_mapping[["contig", "block_id", "category", "blocks_aggregated"]].to_csv(f"{prefix}_minimap2_categories.tsv", sep="\t") +minimap_mapping[["contig", "block_id", "category", "blocks_aggregated"]].to_csv( + f"{prefix}_minimap2_categories.tsv", sep="\t" +) save_dict(recovered, f"{prefix}_recovered", json_cols) categories = minimap_mapping["category"].value_counts() -scores = pd.Series([super_recovered_amount_of_blocks, og_recovered_amount_of_blocks, recovered_with_n, chimeric_with_n], index=["minimap2_merged_recovered_blocks", "minimap2_orthologous_recovered_blocks", "minimap2_recovered_with_n", "minimap2_chimeric_with_n"]) +scores = pd.Series( + [super_recovered_amount_of_blocks, og_recovered_amount_of_blocks, recovered_with_n, chimeric_with_n], + index=[ + "minimap2_merged_recovered_blocks", + "minimap2_orthologous_recovered_blocks", + "minimap2_recovered_with_n", + "minimap2_chimeric_with_n", + ], +) -categories.add(scores, fill_value=0).sort_index().to_csv(f"{prefix}_minimap2_category_counts.tsv", sep="\t", header=[prefix]) +categories.add(scores, fill_value=0).sort_index().to_csv( + f"{prefix}_minimap2_category_counts.tsv", sep="\t", header=[prefix] +) ns_per_contig_df.to_csv(f"{prefix}_contig_n_stats.tsv", sep="\t", index=False) diff --git a/bin/minimap2_overlap_selection.py b/bin/minimap2_overlap_selection.py index 3c40e4f..bb58fb2 100755 --- a/bin/minimap2_overlap_selection.py +++ b/bin/minimap2_overlap_selection.py @@ -11,12 +11,11 @@ def deduplicate_by_key(df, group_key, seed): return df minimap_mapping = df.copy() group = [group_key] - print(df) max_cols = ["Mapping quality", "Number of matching bases in the mapping"] new_max_cols = ["Max Mapping quality", "Max Number of matching bases in the mapping"] - minimap_mapping[new_max_cols] = (minimap_mapping.groupby(group)[max_cols].transform("max")) + minimap_mapping[new_max_cols] = minimap_mapping.groupby(group)[max_cols].transform("max") min_col = ["Query sequence length"] new_min_col = ["Min Query sequence length"] @@ -24,71 +23,47 @@ def deduplicate_by_key(df, group_key, seed): minimap_mapping[new_min_col] = minimap_mapping.groupby(group)[min_col].transform("min") # calculate a cumulative score to include hits with less then optimal matches - minimap_mapping["cumulative_score"] = 5 * (minimap_mapping[max_cols[0]] == minimap_mapping[new_max_cols[0]]) + 3 * (minimap_mapping[max_cols[1]] == minimap_mapping[new_max_cols[1]]) + 1* (minimap_mapping[min_col[0]] == minimap_mapping[new_min_col[0]]) + minimap_mapping["cumulative_score"] = ( + 5 * (minimap_mapping[max_cols[0]] == minimap_mapping[new_max_cols[0]]) + + 3 * (minimap_mapping[max_cols[1]] == minimap_mapping[new_max_cols[1]]) + + 1 * (minimap_mapping[min_col[0]] == minimap_mapping[new_min_col[0]]) + ) minimap_mapping["max_cumulative_score"] = minimap_mapping.groupby(group)["cumulative_score"].transform("max") - minimap_mapping["optimal_hit"] = (minimap_mapping["cumulative_score"] == minimap_mapping["max_cumulative_score"]) + minimap_mapping["optimal_hit"] = minimap_mapping["cumulative_score"] == minimap_mapping["max_cumulative_score"] - cooptimal_bools = (minimap_mapping[minimap_mapping["optimal_hit"]].groupby(group).transform("size") > 1) + cooptimal_bools = minimap_mapping[minimap_mapping["optimal_hit"]].groupby(group).transform("size") > 1 if sum(cooptimal_bools): minimap_mapping.loc[cooptimal_bools.index, "cooptimal_hits"] = cooptimal_bools minimap_mapping["cooptimal_hits"] = ( - minimap_mapping["cooptimal_hits"] - .astype("boolean") - .fillna(False) - .astype(bool) - ) - cooptimal_hits_selection_idx = minimap_mapping[minimap_mapping["cooptimal_hits"]].groupby(group).sample(n=1, random_state=seed).index + minimap_mapping["cooptimal_hits"].astype("boolean").fillna(False).astype(bool) + ) + cooptimal_hits_selection_idx = ( + minimap_mapping[minimap_mapping["cooptimal_hits"]].groupby(group).sample(n=1, random_state=seed).index + ) minimap_mapping.loc[cooptimal_hits_selection_idx, "selected_cooptimal_hit"] = True minimap_mapping["selected_cooptimal_hit"] = ( - minimap_mapping["selected_cooptimal_hit"] - .astype("boolean") - .fillna(False) - .astype(bool) - ) - selection_idx = minimap_mapping[minimap_mapping["optimal_hit"] & (~minimap_mapping["cooptimal_hits"] | minimap_mapping["selected_cooptimal_hit"])].index + minimap_mapping["selected_cooptimal_hit"].astype("boolean").fillna(False).astype(bool) + ) + selection_idx = minimap_mapping[ + minimap_mapping["optimal_hit"] + & (~minimap_mapping["cooptimal_hits"] | minimap_mapping["selected_cooptimal_hit"]) + ].index return df.loc[selection_idx] - + else: - optimal_selection = minimap_mapping[minimap_mapping["optimal_hit"]].index + optimal_selection = minimap_mapping[minimap_mapping["optimal_hit"]].index return df.loc[optimal_selection] -def save_dict(filtered_df, output_file_name, cols): - key_val_dict = filtered_df[json_cols].set_index(json_cols[0])[json_cols[1]].to_dict() - with open(f"{output_file_name}.json", "w") as f: - json.dump(key_val_dict, f) - - -def count_ns_in_contigs(contigs_path): - all_entries = [] - - for record in SeqIO.parse(contigs_path, "fasta"): - total_ns = record.seq.count("N") - max_run = 0 - current = 0 - for c in record.seq: - if c == "N": - current += 1 - max_run = max(max_run, current) - else: - current = 0 - all_entries.append([record.id, total_ns, max_run]) - - return pd.DataFrame(all_entries, columns=["contig_id", "total_Ns", "max_N_run"]) - - minimap_path = sys.argv[1] -contig_path = sys.argv[2] -seed = int(sys.argv[3]) -prefix = sys.argv[4] - - +seed = int(sys.argv[2]) +prefix = sys.argv[3] mapping_col_names = [ "contig", # Query sequence name @@ -102,22 +77,27 @@ def count_ns_in_contigs(contigs_path): "Target end coordinate on the original strand", "Number of matching bases in the mapping", "Number bases, including gaps, in the mapping", - "Mapping quality", "1", "2", "3", "4", "5", "6" + "Mapping quality", + "1", + "2", + "3", + "4", + "5", + "6", ] json_cols = ["contig", "block_id"] -minimap_mapping = pd.read_csv( - minimap_path, sep="\t", header=None, names=mapping_col_names, dtype={'contig': str} -) - +minimap_mapping = pd.read_csv(minimap_path, sep="\t", header=None, names=mapping_col_names, dtype={"contig": str}) -ns_per_contig_df = count_ns_in_contigs(contig_path) - -minimap_mapping["correctly_mapped_bases"] = minimap_mapping["Number of matching bases in the mapping"] / minimap_mapping["Target sequence length"] +minimap_mapping["correctly_mapped_bases"] = ( + minimap_mapping["Number of matching bases in the mapping"] / minimap_mapping["Target sequence length"] +) # -> "Number of matching bases in the mapping"/ "Target sequence length" (blocks are targets) # Rationale: Number of matching bases/ length of block -> block bases should be reconstructed -minimap_mapping["covered_reference"] = minimap_mapping["Number bases, including gaps, in the mapping"] / minimap_mapping["Target sequence length"] +minimap_mapping["covered_reference"] = ( + minimap_mapping["Number bases, including gaps, in the mapping"] / minimap_mapping["Target sequence length"] +) pre_filter_mappings = minimap_mapping.shape[0] @@ -125,28 +105,20 @@ def count_ns_in_contigs(contigs_path): minimap_mapping = minimap_mapping[minimap_mapping["correctly_mapped_bases"] >= 0.95] minimap_mapping = minimap_mapping[ - (minimap_mapping["covered_reference"] <= 1.05) & - (minimap_mapping["covered_reference"] >= 0.95) + (minimap_mapping["covered_reference"] <= 1.05) & (minimap_mapping["covered_reference"] >= 0.95) ] minimap_mapping.to_csv("filtered_df.tsv", sep="\t") minimap_mapping = minimap_mapping.reset_index(drop=True) -minimap_mapping["block_size"] = ( - minimap_mapping - .groupby(["block_id"]) - .transform("size") -) +minimap_mapping["block_size"] = minimap_mapping.groupby(["block_id"]).transform("size") -minimap_mapping["contig_size"] = ( - minimap_mapping - .groupby(["contig"]) - .transform("size") -) +minimap_mapping["contig_size"] = minimap_mapping.groupby(["contig"]).transform("size") minimap_mapping = minimap_mapping.sort_values( - by=["Mapping quality", "Number of matching bases in the mapping", "Query sequence length"], ascending=[False, False, True] + by=["Mapping quality", "Number of matching bases in the mapping", "Query sequence length"], + ascending=[False, False, True], ) minimap_mapping_chosen = minimap_mapping.copy() @@ -157,18 +129,4 @@ def count_ns_in_contigs(contigs_path): minimap_mapping_chosen["category"] = "overlap_block" -contigs_n_df = pd.merge( - minimap_mapping_chosen, - ns_per_contig_df, - how='left', - right_on='contig_id', - left_on='contig' -) - -contigs_n_df[["contig", "block_id", "category", "total_Ns", "max_N_run"]].to_csv(f"{prefix}_minimap2_overlap_blocks.tsv", sep="\t") - -save_dict(minimap_mapping_chosen, f"{prefix}_overlap_recovered", json_cols) - -recovered_with_n = (contigs_n_df["total_Ns"] > 0).sum() - -pd.Series([recovered_with_n], index=["minimap2_overlap_with_n"]).to_csv(f"{prefix}_overlap_n_counts.tsv", sep="\t", header=[prefix]) +minimap_mapping_chosen[["contig", "block_id", "category"]].to_csv(f"{prefix}_mapped_overlap_blocks.tsv", sep="\t") diff --git a/bin/parse_bowtie2_logs.py b/bin/parse_bowtie2_logs.py index 20e098a..f9d7840 100755 --- a/bin/parse_bowtie2_logs.py +++ b/bin/parse_bowtie2_logs.py @@ -4,6 +4,7 @@ import pandas as pd import sys + def parse_bowtie2_log(filepath): stats = { "total_reads": None, @@ -19,7 +20,7 @@ def parse_bowtie2_log(filepath): "overall_alignment_rate": None, } - with open(filepath, 'r') as file: + with open(filepath, "r") as file: for line in file: line = line.strip() if match := re.match(r"^(\d+) reads; of these:", line): @@ -49,10 +50,11 @@ def parse_bowtie2_log(filepath): return stats + log_file = sys.argv[1] colname = sys.argv[2] logs_dict = parse_bowtie2_log(log_file) -logs_df = pd.DataFrame.from_dict(logs_dict, orient='index') +logs_df = pd.DataFrame.from_dict(logs_dict, orient="index") logs_df.columns = [colname] -logs_df.to_csv(sys.stdout, sep='\t') \ No newline at end of file +logs_df.to_csv(sys.stdout, sep="\t") diff --git a/bin/sort_dataframe.py b/bin/sort_dataframe.py index ae3edd8..0f1562b 100755 --- a/bin/sort_dataframe.py +++ b/bin/sort_dataframe.py @@ -7,7 +7,7 @@ axis = int(sys.argv[2]) prefix = sys.argv[3] -df = pd.read_csv(df_file, sep='\t', index_col=0) +df = pd.read_csv(df_file, sep="\t", index_col=0) df = df.sort_index(axis=axis) diff --git a/bin/summarize_mapping_stats.py b/bin/summarize_mapping_stats.py index 53c1380..dd405e2 100755 --- a/bin/summarize_mapping_stats.py +++ b/bin/summarize_mapping_stats.py @@ -19,12 +19,17 @@ number_of_reads_with_duplicated_mapping = exploded_genes[exploded_genes.duplicated()].nunique() -summary_genes = pd.DataFrame([summary_df["origin_gene_nunique"].describe(), summary_df["orthogroup_nunique"].describe()]).T +summary_genes = pd.DataFrame( + [summary_df["origin_gene_nunique"].describe(), summary_df["orthogroup_nunique"].describe()] +).T summary_genes.columns = [f"{assembler_name}_{col}" for col in summary_genes.columns] -summary_mapping_stats = pd.DataFrame([number_of_genes_with_one_read_mapped, number_of_reads_with_duplicated_mapping], columns=[f"{assembler_name}_mapping_stats"], - index=["number_of_genes_with_one_read_mapped", "number_of_reads_with_duplicated_read_mapping"]) +summary_mapping_stats = pd.DataFrame( + [number_of_genes_with_one_read_mapped, number_of_reads_with_duplicated_mapping], + columns=[f"{assembler_name}_mapping_stats"], + index=["number_of_genes_with_one_read_mapped", "number_of_reads_with_duplicated_read_mapping"], +) summary_mapping_stats.to_csv(f"{assembler_name}_mapping_stats.tsv", sep="\t") summary_genes.to_csv(f"{assembler_name}_contig_gene_stats.tsv", sep="\t") diff --git a/example/dataset.config b/example/dataset.config index af72d5e..c1b477f 100644 --- a/example/dataset.config +++ b/example/dataset.config @@ -2,7 +2,7 @@ params { // A nf-core sample sheet (csv) build from the Marbel reads https://nf-co.re/demo/1.0.1/docs/usage/ or see the example file: reads.csv reads = "" // Multiple contig files supplied with path *.fa - contigs = "" + contigs = "PATH/*.fa" // Files generated by marbel in summary dir reference_cds = "summary/metatranscriptome_reference.fasta" diff --git a/example/resources.config b/example/resources.config index ba17740..e755bf3 100644 --- a/example/resources.config +++ b/example/resources.config @@ -70,4 +70,8 @@ process { executor = 'local' ext.args = '123456' } + + withName: MINIMAP2_MAP_OVERLAP { + ext.suffix = "-overlap" + } } \ No newline at end of file diff --git a/modules/local/minimap2/map/main.nf b/modules/local/minimap2/map/main.nf index b7b9154..33e8517 100644 --- a/modules/local/minimap2/map/main.nf +++ b/modules/local/minimap2/map/main.nf @@ -14,16 +14,16 @@ process MINIMAP2_MAP { output: - tuple val(meta), path("${meta.id}_mapping.tsv") , emit: mapping - path "versions.yml" , emit: versions + tuple val(meta), path("${prefix}_mapping${suffix}.tsv") , emit: mapping + path "versions.yml" , emit: versions when: task.ext.when == null || task.ext.when script: def args = task.ext.args ?: '' - def prefix = task.ext.prefix ?: "${meta.id}" - + prefix = task.ext.prefix ?: "${meta.id}" + suffix = task.ext.suffix ?: '' """ @@ -31,7 +31,7 @@ process MINIMAP2_MAP { $args \\ -t $task.cpus \\ ${reference} \\ - ${reads} > ${prefix}_mapping.tsv + ${reads} > ${prefix}_mapping${suffix}.tsv cat <<-END_VERSIONS > versions.yml diff --git a/modules/local/scripts/calculate_final_scores/main.nf b/modules/local/scripts/calculate_final_scores/main.nf deleted file mode 100644 index e70e8ae..0000000 --- a/modules/local/scripts/calculate_final_scores/main.nf +++ /dev/null @@ -1,31 +0,0 @@ -process CALCULATEFINALSCORES { - tag "$meta.id" - label "process_medium" - //conda "${moduleDir}/environment.yml" - //TODO: create a custom container with pandas and jq - container "quay.io/tensulin/orb_toolchain:1.0" - - input: - tuple val(meta), path(merged_scores), path(blocks_df), path(chim_blocks_df) - - output: - tuple val(meta), path("${prefix}_orb_scores.tsv") , emit: final_scores - path "versions.yml" , emit: versions - - when: - task.ext.when == null || task.ext.when - - script: - def args = task.ext.args ?: '' - prefix = task.ext.prefix ?: "${meta.id}" - - """ - calculate_ratios.py ${merged_scores} ${blocks_df} ${chim_blocks_df} > ${prefix}_orb_scores.tsv - - cat <<-END_VERSIONS > versions.yml - "${task.process}": - python: "\$(python3 --version | sed 's/Python //')" - pandas: "\$(python3 -c 'import pandas as pd; print(pd.__version__)')" - END_VERSIONS - """ -} diff --git a/modules/local/scripts/categorize_contigs/main.nf b/modules/local/scripts/categorize_contigs/main.nf index 96921c4..d153656 100644 --- a/modules/local/scripts/categorize_contigs/main.nf +++ b/modules/local/scripts/categorize_contigs/main.nf @@ -5,11 +5,11 @@ process CATEGORIZECONTIGS { container "quay.io/tensulin/orb_toolchain:1.0" input: - tuple val(meta), path(contig_ids), path(length_filtered_contig_ids), path(minimap2_categories), path(mapped_chim_scores), path(assembler_mapping), path(gene_summary), path(contigs_fasta) + tuple val(meta), path(contig_ids), path(length_filtered_contig_ids), path(minimap2_categories), path(assembler_mapping), path(gene_summary), path(contigs_fasta) output: tuple val(meta), path("${prefix}_contigs_categorised.tsv") , emit: contig_categorisation - path "versions.yml" , emit: versions + path "versions.yml" , emit: versions when: @@ -20,8 +20,7 @@ process CATEGORIZECONTIGS { prefix = task.ext.prefix ?: "${meta.id}" """ - categorize_contigs.py ${contig_ids} ${minimap2_categories} ${mapped_chim_scores} ${length_filtered_contig_ids} \\ - ${assembler_mapping} ${gene_summary} ${prefix} + categorize_contigs.py ${contig_ids} ${minimap2_categories} ${length_filtered_contig_ids} ${assembler_mapping} ${gene_summary} ${prefix} cat <<-END_VERSIONS > versions.yml "${task.process}": diff --git a/modules/local/scripts/extract_mapped_ids/main.nf b/modules/local/scripts/extract_mapped_ids/main.nf index dfcf2fd..7774731 100644 --- a/modules/local/scripts/extract_mapped_ids/main.nf +++ b/modules/local/scripts/extract_mapped_ids/main.nf @@ -5,7 +5,7 @@ process EXTRACTMAPPEDIDS { container "quay.io/biocontainers/jq:1.5--4" input: - tuple val(meta), path(mapped_contigs), path(mapped_contigs_chim) + tuple val(meta), path(mapped_contigs) output: tuple val(meta), path("${prefix}_mapped_ids.txt") , emit: mapped_ids @@ -19,7 +19,7 @@ process EXTRACTMAPPEDIDS { prefix = task.ext.prefix ?: "${meta.id}" """ - jq -r 'keys[]' ${mapped_contigs} ${mapped_contigs_chim} > ${prefix}_mapped_ids.txt + jq -r 'keys[]' ${mapped_contigs} > ${prefix}_mapped_ids.txt cat <<-END_VERSIONS > versions.yml "${task.process}": diff --git a/modules/local/scripts/gather_results/main.nf b/modules/local/scripts/gather_results/main.nf index f876d35..81f661e 100644 --- a/modules/local/scripts/gather_results/main.nf +++ b/modules/local/scripts/gather_results/main.nf @@ -6,7 +6,7 @@ process GATHERRESULTS { container "quay.io/tensulin/orb_toolchain:1.0" input: - tuple val(meta), path(contig_ids), path(length_filtered_contig_ids), path(minimap2_categories), path(mapped_chim_scores), path(assembler_mapping), path(gene_summary), path(contigs_fasta) + tuple val(meta), path(contig_ids), path(length_filtered_contig_ids), path(minimap2_categories), path(assembler_mapping), path(gene_summary), path(contigs_fasta) output: tuple val(meta), path("${prefix}_scores.tsv") , emit: scores @@ -31,8 +31,8 @@ process GATHERRESULTS { prefix = task.ext.prefix ?: "${meta.id}" """ - gather_results.py ${contig_ids} ${minimap2_categories} ${mapped_chim_scores} ${length_filtered_contig_ids} \\ - ${assembler_mapping} ${gene_summary} ${contigs_fasta} ${prefix} + gather_results.py ${contig_ids} ${minimap2_categories} ${length_filtered_contig_ids} ${assembler_mapping} \\ + ${gene_summary} ${contigs_fasta} ${prefix} cat <<-END_VERSIONS > versions.yml "${task.process}": diff --git a/modules/local/scripts/minimap2classification/main.nf b/modules/local/scripts/minimap2classification/main.nf index 80434e9..b35b2bb 100644 --- a/modules/local/scripts/minimap2classification/main.nf +++ b/modules/local/scripts/minimap2classification/main.nf @@ -5,7 +5,7 @@ process MINIMAP2CLASSIFICATION { container "quay.io/tensulin/orb_toolchain:1.0" input: - tuple val(meta), path(mapping), path(contigs), path(gene_summary) + tuple val(meta), path(mapping), path(overlap_selection), path(contigs), path(gene_summary) output: tuple val(meta), path("${prefix}_recovered.json") , emit: map @@ -23,7 +23,7 @@ process MINIMAP2CLASSIFICATION { prefix = task.ext.prefix ?: "${meta.id}" """ - minimap2_classification.py ${mapping} ${contigs} ${gene_summary} ${args} ${prefix} + minimap2_classification.py ${mapping} ${overlap_selection} ${contigs} ${gene_summary} ${args} ${prefix} cat <<-END_VERSIONS > versions.yml "${task.process}": diff --git a/modules/local/scripts/minimap2overlapselection/main.nf b/modules/local/scripts/minimap2overlapselection/main.nf index 44618bc..f9f689a 100644 --- a/modules/local/scripts/minimap2overlapselection/main.nf +++ b/modules/local/scripts/minimap2overlapselection/main.nf @@ -5,12 +5,10 @@ process MINIMAP2OVERLAPSELECTION { container "quay.io/tensulin/orb_toolchain:1.0" input: - tuple val(meta), path(mapping), path(contigs) + tuple val(meta), path(mapping) output: - tuple val(meta), path("${prefix}_overlap_recovered.json") , emit: map - tuple val(meta), path("${prefix}_minimap2_overlap_blocks.tsv"), emit: categories - tuple val(meta), path("${prefix}_overlap_n_counts.tsv"), emit: n_counts + tuple val(meta), path("${prefix}_mapped_overlap_blocks.tsv"), emit: overlap_blocks path "versions.yml" , emit: versions when: @@ -22,7 +20,7 @@ process MINIMAP2OVERLAPSELECTION { prefix = task.ext.prefix ?: "${meta.id}" """ - minimap2_overlap_selection.py ${mapping} ${contigs} ${args} ${prefix} + minimap2_overlap_selection.py ${mapping} ${args} ${prefix} cat <<-END_VERSIONS > versions.yml "${task.process}": diff --git a/modules/nf-core/custom/dumpsoftwareversions/templates/dumpsoftwareversions.py b/modules/nf-core/custom/dumpsoftwareversions/templates/dumpsoftwareversions.py index e55b8d4..0ae5918 100755 --- a/modules/nf-core/custom/dumpsoftwareversions/templates/dumpsoftwareversions.py +++ b/modules/nf-core/custom/dumpsoftwareversions/templates/dumpsoftwareversions.py @@ -3,7 +3,6 @@ """Provide functions to merge multiple versions.yml files.""" - import platform from textwrap import dedent @@ -12,9 +11,7 @@ def _make_versions_html(versions): """Generate a tabular HTML output of all versions for MultiQC.""" - html = [ - dedent( - """\\ + html = [dedent("""\\