From 2606c955e3b31113d844a0315921e6d4addd3015 Mon Sep 17 00:00:00 2001 From: Stefan Janssen Date: Tue, 24 Feb 2026 11:28:46 +0100 Subject: [PATCH 01/13] add ability to report class size in recovery plot --- plotting/plot_include_orb.py | 26 +++++++++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/plotting/plot_include_orb.py b/plotting/plot_include_orb.py index 14c2980..83dfa72 100644 --- a/plotting/plot_include_orb.py +++ b/plotting/plot_include_orb.py @@ -18,7 +18,7 @@ from statannotations.Annotator import Annotator -def plot_recovery(fp_orb_basedir, settings, num_columns:int=3, verbose=True, detail_view=False): +def plot_recovery(fp_orb_basedir, settings, num_columns:int=3, verbose=True, detail_view=False, report_percent=True): """Plots contig recovery. Parameters @@ -34,6 +34,10 @@ def plot_recovery(fp_orb_basedir, settings, num_columns:int=3, verbose=True, det Maximal number of panels in a row. verbose : boolean Report progress on sys.stderr + report_percent : boolean + If true, legend items also mention the fraction of contigs classified as + such across all environments / all assemblers. This gives a feeling for + importance of a class. Returns ------- @@ -107,6 +111,26 @@ def plot_recovery(fp_orb_basedir, settings, num_columns:int=3, verbose=True, det sorted_handels.append(handles[labels.index(label)]) handles = sorted_handels labels = sorted_labels + + if report_percent: + # print percent of contigs in each class in legend + # compute ratio classes across all environments and assembler (except missed blocks) + denominators_noneutral = data_recovery[[c['label'] + for _, c in settings['contig_classes'].items() + if c['class'] != 'neutral' + if c['label'] in data_recovery.columns]].sum() + percent = (denominators_noneutral / denominators_noneutral.sum()) + # compute ratio of missed blocks: neutral / (good + neutral) blocks + denominators_rhs = data_recovery[[c['label'] + for _, c in settings['contig_classes'].items() + if c['class'] != 'bad' + if c['label'] in data_recovery.columns]].sum() + percent = pd.concat([percent, (denominators_rhs[[c['label'] + for _, c in settings['contig_classes'].items() + if c['class'] == 'neutral']] / denominators_rhs.sum())]).apply(lambda x: '%.2f%%' % (x * 100)) + # updates legend labels + labels = ['%s (%s)' % (l, percent.loc[l]) if l in percent.index else l for l in labels] + ax_good.legend(handles, labels, ncol=8, bbox_to_anchor=(-0.1, -0.20)) else: ax_good.legend().remove() From 9df98bbc2a22f6a5daab5516fd8c493e604a753a Mon Sep 17 00:00:00 2001 From: Timo Wentong Lin Date: Thu, 26 Feb 2026 16:51:53 +0100 Subject: [PATCH 02/13] update overlap, remove overlap, improve categorisation, update config --- bin/calculate_ratios.py | 22 ------ bin/categorize_contigs.py | 34 +++----- bin/gather_results.py | 28 +++---- bin/minimap2_classification.py | 19 +++-- bin/minimap2_overlap_selection.py | 50 +----------- example/dataset.config | 2 +- example/resources.config | 4 + modules/local/minimap2/map/main.nf | 10 +-- .../scripts/calculate_final_scores/main.nf | 31 -------- .../local/scripts/categorize_contigs/main.nf | 7 +- .../local/scripts/extract_mapped_ids/main.nf | 4 +- modules/local/scripts/gather_results/main.nf | 6 +- .../scripts/minimap2classification/main.nf | 4 +- .../scripts/minimap2overlapselection/main.nf | 8 +- plotting/plot_include_orb.py | 3 +- plotting/style.yaml | 21 +++-- test_nextflow.config | 14 ++-- workflows/refbasedassemblereval.nf | 79 +++++-------------- 18 files changed, 101 insertions(+), 245 deletions(-) delete mode 100755 bin/calculate_ratios.py delete mode 100644 modules/local/scripts/calculate_final_scores/main.nf 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/categorize_contigs.py b/bin/categorize_contigs.py index a06680e..7189818 100755 --- a/bin/categorize_contigs.py +++ b/bin/categorize_contigs.py @@ -27,15 +27,13 @@ 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"]) @@ -45,11 +43,6 @@ except (pl.exceptions.NoDataError, OSError): 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() except (pl.exceptions.NoDataError, OSError): @@ -94,17 +87,18 @@ 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( + +unmapped_contigs = unassigned_contigs.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)) +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")) @@ -118,10 +112,6 @@ .alias("category") ]) -contigs_with_cat_no_mapped = contigs_with_cat.filter( - pl.col("category") != "mapped_contigs" -) - print("before concat") result = pl.concat( diff --git a/bin/gather_results.py b/bin/gather_results.py index d2481fa..76aaf9a 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,18 +30,15 @@ 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())] @@ -89,7 +81,7 @@ 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), + "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),} diff --git a/bin/minimap2_classification.py b/bin/minimap2_classification.py index c0317f3..e22871f 100755 --- a/bin/minimap2_classification.py +++ b/bin/minimap2_classification.py @@ -76,10 +76,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 @@ -99,7 +100,11 @@ def count_ns_in_contigs(contigs_path): json_cols = ["contig", "block_id"] minimap_mapping = pd.read_csv( - minimap_path, sep="\t", header=None, names=mapping_col_names, dtype={'contig': str} + 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) @@ -120,6 +125,8 @@ def count_ns_in_contigs(contigs_path): (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"] = ( @@ -264,6 +271,8 @@ def count_ns_in_contigs(contigs_path): 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, diff --git a/bin/minimap2_overlap_selection.py b/bin/minimap2_overlap_selection.py index 3c40e4f..e2f54e0 100755 --- a/bin/minimap2_overlap_selection.py +++ b/bin/minimap2_overlap_selection.py @@ -11,7 +11,6 @@ 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"] @@ -59,36 +58,10 @@ def deduplicate_by_key(df, group_key, seed): 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 @@ -111,9 +84,6 @@ def count_ns_in_contigs(contigs_path): 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"] # -> "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 @@ -157,18 +127,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") \ No newline at end of file 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/plotting/plot_include_orb.py b/plotting/plot_include_orb.py index 83dfa72..bfe7ce0 100644 --- a/plotting/plot_include_orb.py +++ b/plotting/plot_include_orb.py @@ -131,7 +131,7 @@ def plot_recovery(fp_orb_basedir, settings, num_columns:int=3, verbose=True, det # updates legend labels labels = ['%s (%s)' % (l, percent.loc[l]) if l in percent.index else l for l in labels] - ax_good.legend(handles, labels, ncol=8, bbox_to_anchor=(-0.1, -0.20)) + ax_good.legend(handles, labels, ncol=8, bbox_to_anchor=(0.9, -0.20)) #-0.1 else: ax_good.legend().remove() ax_bad.legend().remove() @@ -329,7 +329,6 @@ def plot_DEgenes(fp_orb_basedir, settings, forOrthogroups=False, num_columns:int shortening_dict = settings["de_image"]["orthogroup_assembler_shortening"] else: shortening_dict = settings["de_image"]["gene_assembler_shortening"] - for ex_env, ex_ass in shortening_dict.items(): if environment == ex_env: pdata = DEfeatures.loc[environment, [ass for ass in DEfeatures.loc[environment, :].index.levels[0] if ass != ex_ass], :] diff --git a/plotting/style.yaml b/plotting/style.yaml index e7cc816..f6a7349 100644 --- a/plotting/style.yaml +++ b/plotting/style.yaml @@ -12,13 +12,13 @@ contig_classes: goodblockcount: True singlerecovery: True simplelabel: 'recovered' - chimeric_mapped_contigs: + overlap_block: class: good color: 'green' label: 'overlapping annotation' goodblockcount: True singlerecovery: False - simplelabel: 'recovered' + simplelabel: 'overlap_block' multi_mapped_contigs_single_og: class: good color: 'lightgreen' @@ -118,22 +118,22 @@ simplelabels: class: bad color: '#F00F18' label: 'chimera' - legend_position: 6 + legend_position: 7 incomplete: class: bad color: '#FECB47' label: 'incomplete' - legend_position: 4 + legend_position: 5 hallucinated: class: bad color: '#B84303' label: 'hallucinated' - legend_position: 5 + legend_position: 6 length_filtered: class: bad color: '#FEEAA2' label: 'length filtered' - legend_position: 3 + legend_position: 4 recovered: class: good color: 'darkgreen' @@ -143,12 +143,17 @@ simplelabels: class: good color: 'lightgreen' label: 'orthologous recovered' + legend_position: 2 + overlap_block: + class: good + color: 'green' + label: 'overlapping annotation' legend_position: 1 missed: class: neutral color: '#dddddd' label: 'missed blocks' - legend_position: 2 + legend_position: 3 # rename ugly technical names with pretty tool / environment labels labels: @@ -184,4 +189,4 @@ de_image: freshwater: Trinity seawater: Trinity orthogroup_assembler_shortening: - healhty_gut: rnaSPAdes \ No newline at end of file + healthy_gut: rnaSPAdes \ No newline at end of file diff --git a/test_nextflow.config b/test_nextflow.config index a01ba53..fdc53e0 100644 --- a/test_nextflow.config +++ b/test_nextflow.config @@ -1,8 +1,4 @@ params { - //holds ncbi accession with genomic information - base_gtf_dir = "/vol/jlab/tlin/all_project/in_silico_dataset/reference_information/references/ncbi_dataset/data" - //translates the edgar names back to accessions - translation_df = "/homes/tlin/Projects/refbasedassemblereval/input_files/full_translation_df.csv" read_length = 126 calour_min_reads = 1000 } @@ -16,7 +12,7 @@ process { time = '20h' errorStrategy = { task.attempt <= 2 ? 'retry' : 'finish' } maxRetries = 2 - clusterOptions = '--exclude=rs700-bcf-5.intra' //rs700-bcf-1.intra,rs700-bcf-4.intra + clusterOptions = '' // --exclude=rs700-bcf-5.intra //rs700-bcf-1.intra,rs700-bcf-4.intra beforeScript = """ var=/var/scratch/\$USER/process_tmp_\$(date +%s%N)\$RANDOM @@ -54,6 +50,10 @@ process { ext.args = '1234' //change for different seed } + withName: MINIMAP2_MAP_OVERLAP { + ext.suffix = "-overlap" + } + withName: MINIMAP2OVERLAPSELECTION { time = '5h' ext.args = '1234' //change for different seed @@ -82,7 +82,6 @@ process { ext.args = "-m ${params.read_length + 1}" } - withName: GATHERRESULTS { memory = '380G' } @@ -150,9 +149,6 @@ process { } } -podman { - runOptions = '--group-add keep-groups --cgroup-manager=cgroupfs --volume=/homes/tlin/Projects/orb/bin/:/container/bin/' -} apptainer { envWhitelist = "http_proxy,https_proxy,ftp_proxy" diff --git a/workflows/refbasedassemblereval.nf b/workflows/refbasedassemblereval.nf index 0077936..4d4c87b 100644 --- a/workflows/refbasedassemblereval.nf +++ b/workflows/refbasedassemblereval.nf @@ -80,9 +80,8 @@ include { SEQKIT_SEQ } from '../modules/nf-core/seqkit/seq/main include { SEQKIT_GREP } from '../modules/nf-core/seqkit/grep/main' include { GATHERRESULTS } from '../modules/local/scripts/gather_results/main' include { CATEGORIZECONTIGS } from '../modules/local/scripts/categorize_contigs/main' -include { EXTRACTMAPPEDVALUES; EXTRACTMAPPEDVALUES as EXTRACTMAPPEDVALUES_OVERLAP } from '../modules/local/jq/extract_values/main' +include { EXTRACTMAPPEDVALUES } from '../modules/local/jq/extract_values/main' include { EXTRACTIDS } from '../modules/local/scripts/extract_ids/main' -include { CALCULATEFINALSCORES } from '../modules/local/scripts/calculate_final_scores/main' include { CALCULATEREFERENCEINFO } from '../modules/local/scripts/calculate_reference_info/main' include { MAKEHISTOGRAMS } from '../modules/local/scripts/make_histograms/main' include { SORTDATAFRAME } from '../modules/local/scripts/sort_dataframe/main' @@ -166,18 +165,6 @@ workflow ASSEMBLEREVAL { ch_versions = ch_versions.mix(MINIMAP2_MAP.out.versions) - MINIMAP2_MAP.out.mapping.join( - assembler_contigs - ).combine( - gene_summary - ).set { mapping_contigs_gene_summary } - - MINIMAP2CLASSIFICATION( - mapping_contigs_gene_summary - ) - - ch_versions = ch_versions.mix(MINIMAP2_MAP.out.versions) - assembler_contigs .combine(overlap_blocks_ch) .multiMap { it -> @@ -192,24 +179,26 @@ workflow ASSEMBLEREVAL { ch_versions = ch_versions.mix(MINIMAP2_MAP_OVERLAP.out.versions) - MINIMAP2_MAP_OVERLAP.out.mapping.join( - assembler_contigs.map { - [["id":it[0]["id"]], it[1]] - } - ).set { overlap_mapping_contigs } - - MINIMAP2OVERLAPSELECTION ( - overlap_mapping_contigs + MINIMAP2OVERLAPSELECTION( + MINIMAP2_MAP_OVERLAP.out.mapping ) ch_versions = ch_versions.mix(MINIMAP2OVERLAPSELECTION.out.versions) - MINIMAP2CLASSIFICATION.out.map.join( - MINIMAP2OVERLAPSELECTION.out.map - ).set { all_mapping } + MINIMAP2_MAP.out.mapping.join( + MINIMAP2OVERLAPSELECTION.out.overlap_blocks + ).join( + assembler_contigs + ).combine( + gene_summary + ).set { overlap_mapping_contig_gene_summary } + + MINIMAP2CLASSIFICATION( + overlap_mapping_contig_gene_summary + ) - EXTRACTMAPPEDIDS( - all_mapping + EXTRACTMAPPEDIDS( + MINIMAP2CLASSIFICATION.out.map ) ch_versions = ch_versions.mix(EXTRACTMAPPEDIDS.out.versions) @@ -296,17 +285,11 @@ workflow ASSEMBLEREVAL { ch_versions = ch_versions.mix(EXTRACTMAPPEDVALUES.out.versions) - EXTRACTMAPPEDVALUES_OVERLAP ( - MINIMAP2OVERLAPSELECTION.out.map - ) - - ch_versions = ch_versions.mix(EXTRACTMAPPEDVALUES_OVERLAP.out.versions) - assembler_contigs.join( SEQKIT_GREP.out.filter ).join( SEQKIT_SEQ.out.fastx - ). set { contigs_filtered } + ).set { contigs_filtered } EXTRACTIDS( contigs_filtered @@ -316,10 +299,6 @@ workflow ASSEMBLEREVAL { EXTRACTIDS.out.contig_ids.join( MINIMAP2CLASSIFICATION.out.categories - ).join( - EXTRACTMAPPEDVALUES_OVERLAP.out.values.map { - [["id":it[0]["id"]], it[1]] - } ).join( STACKDATAFRAMES.out.stacked_dfs ).combine( @@ -342,9 +321,8 @@ workflow ASSEMBLEREVAL { GATHERRESULTS.out.scores.mix( MINIMAP2CLASSIFICATION.out.category_counts - ).mix( - MINIMAP2OVERLAPSELECTION.out.n_counts - ).groupTuple().set { joined_scores } + ).groupTuple() + .set { joined_scores } STACKSCORES( joined_scores @@ -362,23 +340,6 @@ workflow ASSEMBLEREVAL { ) ch_versions = ch_versions.mix(MERGEDATAFRAMESMAPPING.out.versions) - - MERGEDATAFRAMESMAPPING.out.merged_dfs.combine( - blocks_tsv_ch.map { - it[1] - } - ) - .combine( - overlap_blocks_tsv_ch.map { - it[1] - } - ).set { all_scores_with_block } - - CALCULATEFINALSCORES( - all_scores_with_block - ) - - ch_versions = ch_versions.mix(CALCULATEFINALSCORES.out.versions) CALCULATEREFERENCEINFO( reference_cds_val, @@ -458,7 +419,7 @@ workflow ASSEMBLEREVAL { ch_versions = ch_versions.mix(MAKEHISTOGRAMS.out.versions) SORTDATAFRAME( - CALCULATEFINALSCORES.out.final_scores, + MERGEDATAFRAMESMAPPING.out.merged_dfs, Channel.of(1) ) From 83065510baad1462b6af054b37e62d826904a171 Mon Sep 17 00:00:00 2001 From: Stefan Janssen Date: Tue, 10 Mar 2026 12:51:59 +0100 Subject: [PATCH 03/13] different sorting of environments + flexible labels for recovery plot --- plotting/style.yaml | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/plotting/style.yaml b/plotting/style.yaml index f6a7349..8d4404b 100644 --- a/plotting/style.yaml +++ b/plotting/style.yaml @@ -172,13 +172,17 @@ labels: trinity: 'Trinity' environments: - oak: 'oak' - seawater: 'seawater' - healthy_gut: 'healthy gut' + oak: 'oak' moss: 'moss' + seawater: 'seawater' freshwater: 'freshwater' + healthy_gut: 'healthy gut' diseased_gut: 'diseased gut' + recovery_plot: + bad: 'artifacts' + good: 'reasonably good recovered' + read_length: 126 environment_for_legend_display: "healthy_gut" From 47fb581bb973df76ae0e849af84d6825bf42d505 Mon Sep 17 00:00:00 2001 From: Stefan Janssen Date: Tue, 10 Mar 2026 12:52:30 +0100 Subject: [PATCH 04/13] adding status bar to see progress of long running job --- plotting/orthogroups_identity.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/plotting/orthogroups_identity.py b/plotting/orthogroups_identity.py index 97f9a87..d516282 100644 --- a/plotting/orthogroups_identity.py +++ b/plotting/orthogroups_identity.py @@ -5,6 +5,7 @@ import json import os import pandas as pd +from tqdm import tqdm def edlib_identity(seq1, seq2): @@ -45,7 +46,7 @@ def calculate_og_seq_ident(summary_df, sequences_file, env_name, clear_cache=Fal if os.path.exists(f"data_cache/{og_identity_file}"): return pd.read_csv(f"data_cache/{og_identity_file}") - for og, genes in summary.groupby("orthogroup")["gene_name"]: + for og, genes in tqdm(summary.groupby("orthogroup")["gene_name"]): mean_id = compute_edlib_pairwise_identities(genes, sequence_db) og_names.append(og) mean_ids.append(mean_id) From 6e5c8b94fba52c9f11dfc59ecf902af5cebf35a7 Mon Sep 17 00:00:00 2001 From: Stefan Janssen Date: Tue, 10 Mar 2026 12:53:22 +0100 Subject: [PATCH 05/13] adding detonate data compilation and adding verbose msgs to report progres --- plotting/compile_data_orb.py | 42 +++++++++++++++++++++++++++++++++++- 1 file changed, 41 insertions(+), 1 deletion(-) diff --git a/plotting/compile_data_orb.py b/plotting/compile_data_orb.py index 38d3e67..ad6b336 100644 --- a/plotting/compile_data_orb.py +++ b/plotting/compile_data_orb.py @@ -1,3 +1,4 @@ +import sys from glob import glob from os.path import join, basename, relpath, sep import yaml @@ -486,6 +487,26 @@ def getdata_rnaquast(fp_orb_basedir:str, fp_quast_basedir:str, settings, verbose return pd.concat(data).set_index(['environment', 'assembler', 'metric']) +def getdata_detonate(fp_orb_basedir, fp_detonate_basedir:str, settings, verbose=True) -> pd.DataFrame: + # which environments to plot and in which order + environments = get_environments(fp_orb_basedir, settings) + + data = [] + for environment in tqdm(environments, disable=not verbose, desc='Compile Detonate data'): + detonate = pd.read_csv(join(fp_detonate_basedir, '%s_detonate_scores_merged.tsv' % environment), sep="\t", index_col=0) + detonate.index.name = 'metric' + detonate.columns = list(map(lambda x: settings['labels']['assemblers'].get(x, x), detonate.columns)) + detonate = detonate.stack().reset_index().rename(columns={'level_1': 'assembler', 0: 'score'}) + detonate['environment'] = environment + data.append(detonate) + + data = pd.concat(data).set_index(['environment', 'assembler', 'metric']) + # compute score_KC + data = (data.loc[:, :, 'weighted_kmer_recall'] - data.loc[:, :, 'inverse_compression_rate']).rename(columns={'score': 'score_{KC}'}) + + return data + + def getdata_block_recovery(fp_orb_basedir:str, fp_marbel_basedir:str, sequence_file:str, settings, verbose=True): # which environments to plot and in which order @@ -494,15 +515,31 @@ def getdata_block_recovery(fp_orb_basedir:str, fp_marbel_basedir:str, sequence_f recovered_blocks = dict() #fix_ass_labels(ax, dimension='Y') - for environment in tqdm(environments, disable=not verbose, desc='Compiling data for DE orthogroup plot'): + if verbose: + print("Compiling data for DE orthogroup plot:", file=sys.stderr) + for i, environment in enumerate(environments): #, disable=not verbose, desc=''): + if verbose: + print(" %i/%i: %s" % (i+1, len(environments), environment), file=sys.stderr) # obtain orthogroup information about genes + if verbose: + print(" a) read blocks.bed as pandas.DataFrame ...", end="", file=sys.stderr) blocks = pd.read_csv(join(fp_marbel_basedir, '%s_microbiome' % environment, "summary", "blocks.bed"), sep="\t", header=None, names=["gene", "start", "end", "block_name", "read_count"]) + if verbose: + print("done.", file=sys.stderr) blocks['block_length'] = blocks['end'] - blocks['start'] blocks['Coverage'] = blocks['read_count'] * settings["read_length"] / blocks['block_length'] + if verbose: + print(" b) read gene_summary.csv as pandas.DataFrame ...", end="", file=sys.stderr) genes_to_og = pd.read_csv(join(fp_marbel_basedir, '%s_microbiome' % environment, "summary", "gene_summary.csv"), sep=",").set_index('gene_name')["orthogroup"].to_dict() + if verbose: + print("done.", file=sys.stderr) blocks["orthogroup"] = blocks['gene'].map(genes_to_og) all_block_assignments = pd.DataFrame() + if verbose: + print(" c) calculating OG sequence identity:", file=sys.stderr) og_means = calculate_og_seq_ident(join(fp_marbel_basedir, '%s_microbiome' % environment, "summary", "gene_summary.csv"), sequence_file, environment) + if verbose: + print(" done.", file=sys.stderr) og_means.set_index("og_name", inplace=True) og_means_dict = og_means["mean_identity"].to_dict() blocks["mean_identity"] = blocks["orthogroup"].map(og_means_dict) @@ -521,6 +558,9 @@ def getdata_block_recovery(fp_orb_basedir:str, fp_marbel_basedir:str, sequence_f recovered_blocks[environment] = all_block_assignments + if verbose: + print("data compilation completed.", file=sys.stderr) + return recovered_blocks From c3da178f8ab1ff736c104dbca71dd2845d0faa0f Mon Sep 17 00:00:00 2001 From: Stefan Janssen Date: Tue, 10 Mar 2026 12:53:52 +0100 Subject: [PATCH 06/13] reformatting almost all plots for two col layout (except heatmap) --- plotting/plot_include_orb.py | 240 ++++++++++++++++++++++++----------- 1 file changed, 169 insertions(+), 71 deletions(-) diff --git a/plotting/plot_include_orb.py b/plotting/plot_include_orb.py index bfe7ce0..cbecbf5 100644 --- a/plotting/plot_include_orb.py +++ b/plotting/plot_include_orb.py @@ -10,7 +10,7 @@ from matplotlib_venn import venn2 import colorsys import seaborn as sns -from compile_data_orb import get_environments, getdata_DEgenes_mod, getdata_recovery, getdata_gene_recovery, getdata_runtime_memory, getdata_DEvennOrtho, getdata_DEorthogroups_mod, get_recovered_contigs, getdata_rnaquast, getdata_block_recovery, get_contigs, filter_for_assembler_with_ns +from compile_data_orb import get_environments, getdata_DEgenes_mod, getdata_recovery, getdata_gene_recovery, getdata_runtime_memory, getdata_DEvennOrtho, getdata_DEorthogroups_mod, get_recovered_contigs, getdata_rnaquast, getdata_block_recovery, get_contigs, filter_for_assembler_with_ns, getdata_detonate from tqdm import tqdm from scipy.cluster.hierarchy import linkage, leaves_list, dendrogram from scipy.spatial.distance import pdist, squareform @@ -18,6 +18,15 @@ from statannotations.Annotator import Annotator +def legendentries_col2rowordering(handles, labels, num_rows): + entries = list(zip(handles, labels)) + num_columns = int(np.ceil(len(entries) / num_rows)) + reordered_entries = [] + for i in range(num_columns): + reordered_entries.extend(entries[i::num_columns]) + return zip(*reordered_entries) + + def plot_recovery(fp_orb_basedir, settings, num_columns:int=3, verbose=True, detail_view=False, report_percent=True): """Plots contig recovery. @@ -53,9 +62,9 @@ def plot_recovery(fp_orb_basedir, settings, num_columns:int=3, verbose=True, det fig, axes = plt.subplots( int(np.ceil(len(environments) / num_columns)), - num_columns * 2, figsize=(2 * num_columns * 4, np.ceil(len(environments) / num_columns) * 5), - gridspec_kw={"wspace": 0.31, "hspace": 0.3}) - + num_columns * 2, figsize=(2 * num_columns * 4, np.ceil(len(environments) / num_columns) * 3), + gridspec_kw={"wspace": 0.31, "hspace": 0.45}) + if detail_view: view_class = "contig_classes" else: @@ -68,21 +77,24 @@ def plot_recovery(fp_orb_basedir, settings, num_columns:int=3, verbose=True, det ax_bad = axes[i // num_columns, (i % num_columns) * 2] ax_bad.invert_xaxis() orb.loc[:, [c['label'] for _, c in settings[view_class].items() if c['class'] == 'bad']].plot(kind='barh', stacked=True, ax=ax_bad, color={c['label']: c['color'] for _, c in settings['contig_classes'].items()}) + ax_bad.set_xlabel("number contigs") ax_top_bad = ax_bad.twiny() ax_top_bad.xaxis.set_label_position('top') ax_top_bad.set_xticks([]) - ax_top_bad.set_xlabel("weak") - ax_bad.xaxis.set_label_coords(1, -0.08) - ax_bad.text(-0.5, 1.05, chr(97+i), transform=ax_bad.transAxes, fontsize=16, fontweight='bold',) - + ax_top_bad.set_xlabel(settings['labels']['recovery_plot']['bad']) + ax_bad.xaxis.set_label_coords(1, -0.23) + ax_bad.text(-0.5, 1.00, chr(97+i), transform=ax_bad.transAxes, fontsize=16, fontweight='bold',) + ax_bad.tick_params(axis='x', labelrotation=20) + # good contigs ax_good = axes[i // num_columns, (i % num_columns) * 2 + 1] orb.loc[:, [c['label'] for _, c in settings[view_class].items() if (c['class'] == 'good') or (c['class'] == 'neutral')]].plot(kind='barh', stacked=True, ax=ax_good, color={c['label']: c['color'] for _, c in settings['contig_classes'].items()}) ax_good.set_yticks([]) - ax_good.set_xlabel("robust") + ax_good.set_xlabel(settings['labels']['recovery_plot']['good']) ax_good.xaxis.set_label_position('top') - + ax_good.tick_params(axis='x', labelrotation=20) + # concat right (=good) axis directly adjacent to left (=bad) axis ax_good.set_position([ ax_bad.get_position().x1, @@ -90,7 +102,11 @@ def plot_recovery(fp_orb_basedir, settings, num_columns:int=3, verbose=True, det ax_good.get_position().width, ax_good.get_position().height ]) - ax_bad.set_title(settings['labels']['environments'].get(environment, environment), loc='right', horizontalalignment='center') + # ax_bad.set_title(settings['labels']['environments'].get(environment, environment), loc='right', horizontalalignment='center') + ax_bad.text(0.06, 0.96, settings['labels']['environments'].get(environment, environment), + transform=ax_bad.transAxes, + fontsize=12, fontweight='bold', + va='top') # one legend for all panels if i+1 == len(environments): @@ -131,10 +147,17 @@ def plot_recovery(fp_orb_basedir, settings, num_columns:int=3, verbose=True, det # updates legend labels labels = ['%s (%s)' % (l, percent.loc[l]) if l in percent.index else l for l in labels] - ax_good.legend(handles, labels, ncol=8, bbox_to_anchor=(0.9, -0.20)) #-0.1 + # reorder entries in legend + xoffset = 0 + if detail_view is False: + handles, labels = legendentries_col2rowordering(handles, labels, 2) + else: + xoffset = 0.5 + ax_good.legend(handles, labels, ncol=4, bbox_to_anchor=(0.5 + xoffset, -0.3)) else: ax_good.legend().remove() ax_bad.legend().remove() + #TODO: print DF for check with orb base dir return fig @@ -150,15 +173,16 @@ def _color_desaturate(color, saturation=0.75): h, s, v = colorsys.rgb_to_hsv(*mcolors.to_rgb(color)) return colorsys.hsv_to_rgb(h, s, v*saturation) -def plot_gene_recovery(fp_orb_basedir, settings, num_columns:int=3, verbose=True): +def plot_gene_recovery(fp_orb_basedir, settings, num_columns:int=3, verbose=True, recovered_genes=None): # which environments to plot and in which order environments = get_environments(fp_orb_basedir, settings) # load data - recovered_genes = getdata_gene_recovery(fp_orb_basedir, settings) + if recovered_genes is None: + recovered_genes = getdata_gene_recovery(fp_orb_basedir, settings) fig, axes = plt.subplots( int(np.ceil(len(environments) / num_columns)), - num_columns, figsize=(num_columns * 5, np.ceil(len(environments) / num_columns) * 3), + num_columns, figsize=(num_columns * 5, np.ceil(len(environments) / num_columns) * 2.5), gridspec_kw={"wspace": 0.6, "hspace": 0.4} ) @@ -174,7 +198,11 @@ def plot_gene_recovery(fp_orb_basedir, settings, num_columns:int=3, verbose=True ax.axvline(x=recovered_genes.loc['core', environment], color=palette['core']) ax.set_ylabel("") ax.set_xlabel("number recovered genes") - ax.set_title(settings['labels']['environments'].get(environment, environment)) + #ax.set_title(settings['labels']['environments'].get(environment, environment)) + ax.text(0.95, 0.8, settings['labels']['environments'].get(environment, environment), + transform=ax.transAxes, + fontsize=12, fontweight='bold', + va='center', ha='right') ax.set_xscale('log') if i+1 == len(environments): @@ -182,10 +210,10 @@ def plot_gene_recovery(fp_orb_basedir, settings, num_columns:int=3, verbose=True mpatches.Patch(color=palette[assemblers[0]], label='exclusively recovered'), mpatches.Patch(color=_color_desaturate(palette['core'], 0.75), label='recovered by all'), mpatches.Patch(color=_color_desaturate(palette['shared'], 0.75), label='recovered by at least two')], - bbox_to_anchor=(-0.1, -0.25), ncols=3) + bbox_to_anchor=(0.5, -0.35), ncols=3) # panel labels - ax.text(-0.43, 1.05, chr(97+i), transform=ax.transAxes, fontsize=16, fontweight='bold',) + ax.text(-0.48, 0.95, chr(97+i), transform=ax.transAxes, fontsize=16, fontweight='bold',) return fig @@ -193,7 +221,7 @@ def plot_gene_recovery(fp_orb_basedir, settings, num_columns:int=3, verbose=True def plotTimeMemory(fp_caviar_basedir:str, settings, verbose=True): timemem = getdata_runtime_memory(fp_caviar_basedir, settings) - fig, axes = plt.subplots(1, 2, figsize=(10, 4), gridspec_kw={"wspace": 0.5}) + fig, axes = plt.subplots(1, 2, figsize=(10, 2.5), gridspec_kw={"wspace": 0.5}) for ax, (pType, factor, label, title) in zip(axes, [ ('CPU time (seconds)', 3600, 'CPU hours', 'Runtime'), @@ -210,11 +238,11 @@ def plotTimeMemory(fp_caviar_basedir:str, settings, verbose=True): ax.legend().remove() else: ax.axvline(x=64, linestyle='-.', color='gray', zorder=-1, label="64 GB laptop") - ax.legend(bbox_to_anchor=(1.1, -0.15), ncols=7) + ax.legend(bbox_to_anchor=(1.1, -0.25), ncols=7) # panel labels for i, ax in enumerate(axes): - ax.text(-0.43, 1.05, chr(97+i), transform=ax.transAxes, fontsize=16, fontweight='bold',) + ax.text(-0.43, 1.00, chr(97+i), transform=ax.transAxes, fontsize=16, fontweight='bold',) return fig @@ -287,8 +315,8 @@ def plot_DEgenes(fp_orb_basedir, settings, forOrthogroups=False, num_columns:int fig, axes = plt.subplots( int(np.ceil(len(environments) / num_columns)), - num_columns, figsize=(num_columns * 7, np.ceil(len(environments) / num_columns) * 4), - gridspec_kw={"hspace": 0.31, "wspace": 0.5}) + num_columns, figsize=(num_columns * 7, np.ceil(len(environments) / num_columns) * 2.5), + gridspec_kw={"hspace": 0.5, "wspace": 0.5}) BARWIDTH=0.8 palette = {'True Positive': '#238cc3', @@ -347,12 +375,12 @@ def plot_DEgenes(fp_orb_basedir, settings, forOrthogroups=False, num_columns:int ax.legend(handles=[mpatches.Patch(color=palette[cls], label=cls) for cls in palette.keys()] + \ [Line2D([0], [0], color=palette['True Positive'], lw=2, label='Positive')], #title='Category', - bbox_to_anchor=(-0.4, -0.25), + bbox_to_anchor=(0.3, -0.35), ncols=4) return fig -def plot_DEvennOrtho(fp_orb_basedir:str, fp_marbel_basedir:str, settings, num_columns:int=3, verbose=True): +def plot_DEvennOrtho(fp_orb_basedir:str, fp_marbel_basedir:str, settings, verbose=True): # which environments to plot and in which order environments = get_environments(fp_orb_basedir, settings) @@ -364,36 +392,32 @@ def plot_DEvennOrtho(fp_orb_basedir:str, fp_marbel_basedir:str, settings, num_co len(environments), figsize=(len(environments) * 3, len(environments)/2 * 3), #gridspec_kw={"wspace": 0.6, "hspace": 0.4} ) - palette = {'DE orthogroups': 'blue', - 'DE genes contained in orthogroups': 'orange'} + palette = {'DE orthologous groups': 'blue', + 'DE genes contained in orthologous groups': 'orange'} labels = ["", ""] - reordered_environments = [] - for start in range(num_columns): - reordered_environments.extend(environments[start::num_columns]) - - for i, environment in tqdm(enumerate(reordered_environments), disable=not verbose, desc='Drawing panels for DE Venn diagrams'): + for i, environment in tqdm(enumerate(environments), disable=not verbose, desc='Drawing panels for DE Venn diagrams'): ax = axes[0][i] venn2([set(truth[environment][truth[environment]['DEorthogroup']]['orthogroup'].values), set(truth[environment][truth[environment]['DEgene']]['orthogroup'].values)], labels, - ax=ax, set_colors=(palette['DE orthogroups'], palette['DE genes contained in orthogroups'])) + ax=ax, set_colors=(palette['DE orthologous groups'], palette['DE genes contained in orthologous groups'])) ax.set_title(settings['labels']['environments'].get(environment, environment)) - ax.set_ylabel("all orthogroups") + ax.set_ylabel(r'$\mathbf{all}$'+ "\northologous groups") ax = axes[1][i] venn2([set(truth[environment][truth[environment]['DEorthogroup'] & (truth[environment]['OGsize'] == 'multi_gene_OG')]['orthogroup'].values), set(truth[environment][truth[environment]['DEgene'] & (truth[environment]['OGsize'] == 'multi_gene_OG')]['orthogroup'].values)], labels, - ax=ax, set_colors=(palette['DE orthogroups'], palette['DE genes contained in orthogroups'])) - ax.set_ylabel("only multi gene orthogroups") + ax=ax, set_colors=(palette['DE orthologous groups'], palette['DE genes contained in orthologous groups'])) + ax.set_ylabel("only " + r'$\mathbf{multi}$' + " genes\northologous groups") ax = axes[2][i] venn2([set(truth[environment][truth[environment]['DEorthogroup'] & (truth[environment]['OGsize'] == 'single_gene_OG')]['orthogroup'].values), set(truth[environment][truth[environment]['DEgene'] & (truth[environment]['OGsize'] == 'single_gene_OG')]['orthogroup'].values)], labels, - ax=ax, set_colors=(palette['DE orthogroups'], palette['DE genes contained in orthogroups'])) - ax.set_ylabel("only single gene orthogroups") + ax=ax, set_colors=(palette['DE orthologous groups'], palette['DE genes contained in orthologous groups'])) + ax.set_ylabel("only " + r'$\mathbf{single}$' + " gene\northologous groups") if i == 0: for row in range(len(axes)): @@ -410,6 +434,8 @@ def plot_DEvennOrtho(fp_orb_basedir:str, fp_marbel_basedir:str, settings, num_co def plot_heatmap(fp_orb_basedir:str, settings, num_columns:int=3, verbose=True): # which environments to plot and in which order environments = get_environments(fp_orb_basedir, settings) + # re-order such that paired-environments are on top of each other + environments = environments[::2] + environments[1::2] # re-order environments such that two additional "environments" # are spiked in for the color map and the combined environment def _spikein(environments, num_cols=3, spikeelements=['colormap', 'all six environments']): @@ -421,7 +447,7 @@ def _spikein(environments, num_cols=3, spikeelements=['colormap', 'all six envir for chunk in chunks[len(spikeelements):]: reordered.extend(chunk) return reordered - ext_environments = _spikein(environments) + ext_environments = _spikein(environments, num_cols=num_columns) num_columns += 1 # load data @@ -515,8 +541,8 @@ def plot_rnaquast(fp_orb_basedir:str, fp_quast_basedir:str, settings, num_column fig, axes = plt.subplots( int(np.ceil(len(environments) / num_columns)), - num_columns * 2, figsize=(2 * num_columns * 4, np.ceil(len(environments) / num_columns) * 5), - gridspec_kw={"wspace": 0.31, "hspace": 0.3}) + num_columns * 2, figsize=(2 * num_columns * 4, np.ceil(len(environments) / num_columns) * 2.5), + gridspec_kw={"wspace": 0.31, "hspace": 0.5}) for i, environment in tqdm(enumerate(environments), disable=not verbose, desc='Draw RNAquast panels'): order = list(data_recovery[data_recovery['environment'] == environment].sort_values(by='recovery_rank').index) @@ -528,19 +554,19 @@ def plot_rnaquast(fp_orb_basedir:str, fp_quast_basedir:str, settings, num_column ax_bad.invert_xaxis() ax_bad.set_ylabel("") ax_bad.set_xlabel("number contigs") - ax_bad.set_title(settings['labels']['environments'].get(environment, environment), loc='right', horizontalalignment='center') + #ax_bad.set_title(settings['labels']['environments'].get(environment, environment), loc='right', horizontalalignment='center') ax_top_bad = ax_bad.twiny() ax_top_bad.xaxis.set_label_position('top') ax_top_bad.set_xticks([]) - ax_top_bad.set_xlabel("weak") - ax_bad.xaxis.set_label_coords(1, -0.08) + ax_top_bad.set_xlabel(settings['labels']['recovery_plot']['bad']) + ax_bad.xaxis.set_label_coords(1, -0.18) ax_bad.text(-0.5, 1.05, chr(97+i), transform=ax_bad.transAxes, fontsize=16, fontweight='bold',) ax_good = axes[i // num_columns, (i % num_columns) * 2 + 1] sns.barplot(data=quast.loc[environment, :, '95%-assembled isoforms'], x='score', y='assembler', ax=ax_good, order=order, - color=settings['contig_classes']['mapped_contigs']['color']) + color=settings['contig_classes']['minimap2_single_recovered']['color']) ax_good.set_yticks([]) - ax_good.set_xlabel("robust") + ax_good.set_xlabel(settings['labels']['recovery_plot']['good']) ax_good.xaxis.set_label_position('top') ax_good.set_ylabel("") @@ -555,30 +581,88 @@ def plot_rnaquast(fp_orb_basedir:str, fp_quast_basedir:str, settings, num_column if i+1 == len(environments): ax_good.legend(handles=[ mpatches.Patch(color=settings['contig_classes']['multi_mapped_contigs_multi_og']['color'], label='Misassemblies'), - mpatches.Patch(color=settings['contig_classes']['mapped_contigs']['color'], label='95%-assembled isoforms')], + mpatches.Patch(color=settings['contig_classes']['minimap2_single_recovered']['color'], label='95%-assembled isoforms')], ncol=2, bbox_to_anchor=(-1.8, -0.20)) + + ax_bad.text(0.06, 0.96, settings['labels']['environments'].get(environment, environment), + transform=ax_bad.transAxes, + fontsize=12, fontweight='bold', + va='top') + + return fig + + +def plot_rnaquast_singlescore(fp_orb_basedir:str, fp_quast_basedir:str, settings, field:str='Database coverage', num_columns:int=3, verbose=True): + # which environments to plot and in which order + environments = get_environments(fp_orb_basedir, settings) + + # load data + quast = getdata_rnaquast(fp_orb_basedir, fp_quast_basedir, settings) + data_recovery = getdata_recovery(fp_orb_basedir, settings, verbose) + + fig, axes = plt.subplots( + int(np.ceil(len(environments) / num_columns)), + num_columns, figsize=(num_columns * 5, np.ceil(len(environments) / num_columns) * 2.5), + gridspec_kw={"wspace": 0.7, "hspace": 0.6}) + + for i, environment in tqdm(enumerate(environments), disable=not verbose, desc='Draw RNAquast panels'): + order = list(data_recovery[data_recovery['environment'] == environment].sort_values(by='recovery_rank').index) + ax = axes[i // num_columns, i % num_columns] + sns.barplot(data=quast.loc[environment, :, field], x='score', y='assembler', ax=ax, order=order) + ax.set_ylabel("") + ax.set_xlabel(field) + ax.set_title(settings['labels']['environments'].get(environment, environment)) + ax.text(-0.5, 1.05, chr(97+i), transform=ax.transAxes, fontsize=16, fontweight='bold',) return fig -def plot_block_correlation(fp_orb_basedir, fp_marbel_basedir, sequence_file, settings, field, verbose=True, test = 'Mann-Whitney-gt'): +def plot_detonate(fp_orb_basedir:str, fp_detonate_basedir:str, settings, num_columns:int=3, verbose=True): + # which environments to plot and in which order + environments = get_environments(fp_orb_basedir, settings) + + # load data + detonate = getdata_detonate(fp_orb_basedir, fp_detonate_basedir, settings) + data_recovery = getdata_recovery(fp_orb_basedir, settings, verbose) + + fig, axes = plt.subplots( + int(np.ceil(len(environments) / num_columns)), + num_columns, figsize=(num_columns * 5, np.ceil(len(environments) / num_columns) * 2.5), + gridspec_kw={"wspace": 0.7, "hspace": 0.6}) + + for i, environment in tqdm(enumerate(environments), disable=not verbose, desc='Draw RNAquast panels'): + order = list(data_recovery[data_recovery['environment'] == environment].sort_values(by='recovery_rank').index) + ax = axes[i // num_columns, i % num_columns] + sns.barplot(data=detonate.loc[environment, :], x='score_{KC}', y='assembler', ax=ax, order=order) + ax.set_ylabel("") + ax.set_xlabel(r'$%s$' % ax.get_xlabel()) + #ax.set_xlabel(field) + ax.set_title(settings['labels']['environments'].get(environment, environment)) + ax.text(-0.5, 1.05, chr(97+i), transform=ax.transAxes, fontsize=16, fontweight='bold',) + + return fig + +def plot_block_correlation(fp_orb_basedir, fp_marbel_basedir, sequence_file, settings, field, verbose=True, test = 'Mann-Whitney-gt', data_block_recovery=None, num_columns:int=3): # field: block_length, read_mean_count, mean_identity, Coverage environments = get_environments(fp_orb_basedir, settings) - fig, axes = plt.subplots(2, int(len(environments)/2), - figsize=(len(environments)/2*4, 2*4), - gridspec_kw={"wspace": 0.81, "hspace": 0.4} - ) + fig, axes = plt.subplots( + int(np.ceil(len(environments) / num_columns)), + num_columns, + figsize=(int(np.ceil(len(environments) / num_columns)) * 4, num_columns * 6), + gridspec_kw={"wspace": 0.6, "hspace": 0.4}) + axes = axes.flatten() hue_order = ['recovered', 'missed'] - data_block_recovery = getdata_block_recovery(fp_orb_basedir, fp_marbel_basedir, sequence_file, settings, verbose=verbose) + if data_block_recovery is None: + data_block_recovery = getdata_block_recovery(fp_orb_basedir, fp_marbel_basedir, sequence_file, settings, verbose=verbose) assemblers = settings['labels']['assemblers'] # 'Mann-Whitney', 't-test_ind', 'Wilcoxon' logscale = False for i, environment in tqdm(enumerate(environments), disable=not verbose, desc=f'Drawing panels for {field} plot'): ax = axes[i] plotdata = data_block_recovery[environment] - print(plotdata.columns) + # assembler ordered_assembler = sorted(list(plotdata["assembler"].unique())) if field == 'block_length': @@ -587,23 +671,27 @@ def plot_block_correlation(fp_orb_basedir, fp_marbel_basedir, sequence_file, set sns.boxplot(data=plotdata, orient='h', y='assembler', x=field, hue='category', ax=ax, palette={'recovered': 'lightgreen', 'missed': '#dddddd'}, showfliers=False, hue_order=hue_order, order=ordered_assembler) - if logscale: ax.set_xscale('log') if environment != settings["environment_for_legend_display"]: ax.legend().remove() else: - ax.legend(bbox_to_anchor=(1.1, 0.95), title='blocks') + ax.legend(bbox_to_anchor=(1.6, -0.2), title='blocks', ncols=2) #fix_ass_labels(ax, dimension='Y') ax.set_ylabel("") ax.set_xlabel({'read_mean_count': 'coverage: mean read number', 'block_length': 'block length (bps)', 'mean_identity': 'sequence similarity'}.get(field, field)) + ax.set_title(settings['labels']['environments'].get(environment, environment)) + ax.text(-0.43, 1.05, chr(97+i), transform=ax.transAxes, fontsize=16, fontweight='bold',) - annotator = Annotator(ax, [((assembler, 'recovered'), (assembler, 'missed')) for assembler in ordered_assembler], - data=plotdata, x=field, y='assembler', hue='category', orient='h', hue_order=hue_order, order=ordered_assembler) - annotator.configure(test=test, text_format='star', loc='inside', comparisons_correction="fdr_bh", correction_format="default") - testres = annotator.apply_and_annotate() + if test is not None: + annotator = Annotator(ax, [((assembler, 'recovered'), (assembler, 'missed')) for assembler in ordered_assembler], + data=plotdata, x=field, y='assembler', hue='category', orient='h', hue_order=hue_order, order=ordered_assembler, verbose=0) + annotator.configure(test=test, text_format='star', loc='inside', comparisons_correction="fdr_bh", correction_format="default") + testres = annotator.apply_and_annotate() + else: + testres = pd.DataFrame(columns=[1]) # gather mean values with expections of those assembler that yield no sig. difference except_assembler = [ @@ -614,27 +702,28 @@ def plot_block_correlation(fp_orb_basedir, fp_marbel_basedir, sequence_file, set except_assembler = [] mean_field_diff = data_block_recovery[environment].sort_values(['assembler']).set_index(['assembler']).loc[[a for a in ordered_assembler if a not in except_assembler], :].groupby('category')[field].median() - ax.set_title(settings['labels']['environments'].get(environment, environment)) ax_top = ax.twiny() ax_top.set_xticks([]) ax_top.set_xlabel('⌀ recovered ~ %.2f, ⌀ missed ~ %.2f' % (mean_field_diff.loc['recovered'], mean_field_diff['missed'])) - ax.text(-0.43, 1.05, chr(97+i), transform=ax.transAxes, fontsize=16, fontweight='bold',) if (field == 'block_length') and (logscale is False): ax.set_xticks([0, 100, 200, 300, 400, 500, 1000, 2000, 3000], labels=["0", "", "", "", "", "500", "1000", "2000", "3000"]) return fig -def plot_nruns(fp_orb_basedir, settings, verbose=True): +def plot_nruns(fp_orb_basedir, settings, num_columns:int=3, verbose=True, contigs=None): sns.set_style("ticks") environments = get_environments(fp_orb_basedir, settings) - fig, axes = plt.subplots(2, len(environments), - figsize=(len(environments)*4, 2*2), - gridspec_kw={"wspace": 0.31, "hspace": 0.9} - ) + fig, axes = plt.subplots( + int(np.ceil(len(environments) / num_columns)), + num_columns * 2, + figsize=(num_columns * 6, int(np.ceil(len(environments) / num_columns)) * 1.5), + gridspec_kw={"wspace": 0.51, "hspace": 1.3}) + axes = axes.flatten() - contigs = filter_for_assembler_with_ns(fp_orb_basedir, settings) + if contigs is None: + contigs = filter_for_assembler_with_ns(fp_orb_basedir, settings) for i, environment in tqdm(enumerate(environments), disable=not verbose, desc=f'Drawing panels for Ns plot'): n_run_data = contigs[environment] @@ -644,12 +733,13 @@ def plot_nruns(fp_orb_basedir, settings, verbose=True): ax_bad.set_xlim((0, plot_data[0].max()*1.1)) ax_bad.invert_xaxis() ax_bad.set_ylabel("") - ax_bad.set_title(settings['labels']['environments'].get(environment, environment), loc='right', horizontalalignment='center') + #ax_bad.set_title(settings['labels']['environments'].get(environment, environment), loc='right', horizontalalignment='center') + # add panel label ax_bad.text(-0.5, 1.05, chr(97+i), transform=ax_bad.transAxes, fontsize=16, fontweight='bold',) ax_bad.set_xlabel("number contigs") - ax_bad.xaxis.set_label_coords(1, -0.3) + ax_bad.xaxis.set_label_coords(1, -0.4) ax_top_bad = ax_bad.twiny() ax_top_bad.set_xticks([]) ax_top_bad.set_xlabel("weak recovery") @@ -669,10 +759,18 @@ def plot_nruns(fp_orb_basedir, settings, verbose=True): ax_good.get_position().width, ax_good.get_position().height ]) + + ax_good.text(0.95, 0.90, settings['labels']['environments'].get(environment, environment), + transform=ax_good.transAxes, + fontsize=12, fontweight='bold', + va='top', ha='right') # create one joined legend if environment == settings["environment_for_legend_display"]: - ax_good.legend(bbox_to_anchor=(1.8, 1.0)) + handles, lables = ax_good.get_legend_handles_labels() + ax_good.legend( + handles, [l.replace('read length', 'read length (%ibp)' % settings['read_length']) for l in lables], + bbox_to_anchor=(2.8, -0.9), ncol=2) else: ax_good.legend().remove() ax_bad.legend().remove() From 595945b54c3a13f04086c66a6583297529901e9c Mon Sep 17 00:00:00 2001 From: Timo Wentong Lin Date: Tue, 10 Mar 2026 17:43:06 +0100 Subject: [PATCH 07/13] update gitignore, style.yaml --- .gitignore | 3 ++- plotting/style.yaml | 8 ++++---- 2 files changed, 6 insertions(+), 5 deletions(-) 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/plotting/style.yaml b/plotting/style.yaml index 8d4404b..5a006a0 100644 --- a/plotting/style.yaml +++ b/plotting/style.yaml @@ -32,28 +32,28 @@ contig_classes: label: 'minimap2 merged recovered' goodblockcount: False singlerecovery: True - simplelabel: 'none' + simplelabel: 'recovered' minimap2_merged_recovered_blocks: class: good_block color: 'aquamarine' label: 'minimap2 merged recovered blocks' goodblockcount: True singlerecovery: False - simplelabel: 'recovered' + simplelabel: 'none' minimap2_orthologous_recovered: class: good color: 'lime' label: 'minimap2 orthologous recovered' goodblockcount: False singlerecovery: False - simplelabel: 'none' + simplelabel: 'orthologous_recovered' minimap2_orthologous_recovered_blocks: class: good_block color: 'lime' label: 'minimap2 orthologous recovered blocks' goodblockcount: True singlerecovery: False - simplelabel: 'orthologous_recovered' + simplelabel: 'none' missed: class: neutral color: '#dddddd' From 67f2e94d3f220f6d3f013ed5ea28c43d469ce135 Mon Sep 17 00:00:00 2001 From: Stefan Janssen Date: Wed, 18 Mar 2026 23:23:55 +0100 Subject: [PATCH 08/13] split chimera into mosaic and diffuse --- plotting/style.yaml | 31 ++++++++++++++++++------------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/plotting/style.yaml b/plotting/style.yaml index 5a006a0..62dd3d1 100644 --- a/plotting/style.yaml +++ b/plotting/style.yaml @@ -64,10 +64,10 @@ contig_classes: minimap2_chimeric: class: bad color: '#c32aff' - label: 'minimap2 chimera' + label: 'mosaic chimera' goodblockcount: False singlerecovery: False - simplelabel: 'chimera' + simplelabel: 'minimap2_chimeric' minimap2_duplicated_non_chimeric: class: bad color: '#d9c1db' @@ -93,10 +93,10 @@ contig_classes: multi_mapped_contigs_multi_og: class: bad color: '#F00F18' - label: 'chimera' + label: 'diffuse chimera' goodblockcount: False singlerecovery: False - simplelabel: 'chimera' + simplelabel: 'multi_mapped_contigs_multi_og' single_mapped_contigs: class: bad color: '#FECB47' @@ -114,21 +114,26 @@ contig_classes: # oprtion: legend position count from zero simplelabels: - chimera: + minimap2_chimeric: class: bad - color: '#F00F18' - label: 'chimera' + color: '#FECB47' + label: 'mosaic chimera' + legend_position: 8 + multi_mapped_contigs_multi_og: + class: bad + color: '#c32aff' + label: 'diffuse chimera' legend_position: 7 incomplete: class: bad - color: '#FECB47' + color: '#F00F18' label: 'incomplete' - legend_position: 5 + legend_position: 6 hallucinated: class: bad color: '#B84303' label: 'hallucinated' - legend_position: 6 + legend_position: 5 length_filtered: class: bad color: '#FEEAA2' @@ -143,12 +148,12 @@ simplelabels: class: good color: 'lightgreen' label: 'orthologous recovered' - legend_position: 2 + legend_position: 1 overlap_block: class: good color: 'green' label: 'overlapping annotation' - legend_position: 1 + legend_position: 2 missed: class: neutral color: '#dddddd' @@ -181,7 +186,7 @@ labels: recovery_plot: bad: 'artifacts' - good: 'reasonably good recovered' + good: 'reasonably well recovered' read_length: 126 From d25daa199f04865efe145135e892bf34f1cd63d8 Mon Sep 17 00:00:00 2001 From: Stefan Janssen Date: Wed, 18 Mar 2026 23:24:23 +0100 Subject: [PATCH 09/13] add ability to rank DE by F1 and add functions to collect genomic positions of hits --- plotting/compile_data_orb.py | 110 +++++++++++++++++++++++++++++++---- 1 file changed, 99 insertions(+), 11 deletions(-) diff --git a/plotting/compile_data_orb.py b/plotting/compile_data_orb.py index ad6b336..72b6b98 100644 --- a/plotting/compile_data_orb.py +++ b/plotting/compile_data_orb.py @@ -272,7 +272,7 @@ def getdata_DEgenes(fp_orb_basedir:str, settings, verbose=True): return confusion, recovered_contigs -def getdata_DEgenes_mod(fp_orb_basedir: str, settings, verbose=True): +def getdata_DEgenes_mod(fp_orb_basedir: str, settings, verbose=True, sortF1=False): conversion_map = { "deseq2_DE_TN": "True Negative", "deseq2_DE_FP": "False Positive", @@ -312,11 +312,16 @@ def getdata_DEgenes_mod(fp_orb_basedir: str, settings, verbose=True): # add rank information to return dataframe for environment in confusion.index.levels[0]: - - order = list(reversed( - pd.pivot_table(data=confusion.loc[environment, :], index='assembler', columns='class', values='num_genes', aggfunc="sum").sort_values( - by= ['True Positive', 'False Positive', 'False Negative'], - ascending=[False, True, True]).index)) + cfg = pd.pivot_table(data=confusion.loc[environment, :], index='assembler', columns='class', values='num_genes', aggfunc="sum") + # compute F1 score: https://en.wikipedia.org/wiki/F-score#Definition + cfg['F1'] = 2*cfg['True Positive'] / (2*cfg['True Positive'] + cfg['False Positive'] + cfg['False Negative']) + if sortF1: + order = list(cfg.sort_values(by='F1', ascending=True).index) + else: + order = list(reversed( + cfg.sort_values( + by= ['True Positive', 'False Positive', 'False Negative'], + ascending=[False, True, True]).index)) confusion.loc[confusion.loc[environment, order, :].index, 'rank'] = [ rank for rank in list(reversed(range(1, len(order) + 1))) @@ -389,7 +394,7 @@ def getdata_DEorthogroups(fp_orb_basedir:str, fp_marbel_basedir:str, fp_ogtruth_ return confusion, recovered_contigs -def getdata_DEorthogroups_mod(fp_orb_basedir: str, settings, verbose=True): +def getdata_DEorthogroups_mod(fp_orb_basedir: str, settings, verbose=True, sortF1=False): conversion_map = { "deseq2_og_DE_TN": "True Negative", "deseq2_og_DE_FP": "False Positive", @@ -426,10 +431,16 @@ def getdata_DEorthogroups_mod(fp_orb_basedir: str, settings, verbose=True): # add rank information to return dataframe for environment in confusion.index.levels[0]: - order = list(reversed( - pd.pivot_table(data=confusion.loc[environment, :], index='assembler', columns='class', values='num_genes', aggfunc="sum").sort_values( - by= ['True Positive', 'False Positive', 'False Negative'], - ascending=[False, True, True]).index)) + cfg = pd.pivot_table(data=confusion.loc[environment, :], index='assembler', columns='class', values='num_genes', aggfunc="sum") + # compute F1 score: https://en.wikipedia.org/wiki/F-score#Definition + cfg['F1'] = 2*cfg['True Positive'] / (2*cfg['True Positive'] + cfg['False Positive'] + cfg['False Negative']) + if sortF1: + order = list(cfg.sort_values(by='F1', ascending=True).index) + else: + order = list(reversed( + cfg.sort_values( + by= ['True Positive', 'False Positive', 'False Negative'], + ascending=[False, True, True]).index)) confusion.loc[confusion.loc[environment, order, :].index, 'rank'] = [ rank for rank in list(reversed(range(1, len(order) + 1))) @@ -611,3 +622,80 @@ def filter_for_assembler_with_ns(fp_orb_basedir:str, settings, verbose=True): contigs[environment] = n_run_data return contigs + + +def get_genome_accession_map(fp_genomes:str): + """As genome accessions and assembly accessions differ, we here obtain a 1to1 map.""" + # get a mapping from refseqAccession to assemblyAccession + accessionmap = [] + for fp_report in tqdm(sorted(glob(join(fp_genomes, '**/*.jsonl'))), 'get reference genome accessions'): + accessions = pd.read_json(fp_report, lines=True) + accessionmap.append(accessions) + accessionmap = pd.concat(accessionmap) + accessionmap['refseqacc_noversion'] = accessionmap['refseqAccession'].apply(lambda x: x.split('.')[0]) + accessionmap = accessionmap.set_index('refseqacc_noversion').rename(columns={'refseqAccession': 'refseqAccession_versioned'}) + + return accessionmap + + +def get_blockpos_ingenes(fp_genomes:str, fp_marbel_basedir:str, env:str, accessionmap): + """Get relative start / end positions for each block in its according gene.""" + + # load marbel's gene summary + summary = pd.read_csv(join(fp_marbel_basedir, '%s_microbiome' % env, 'summary', 'gene_summary.csv'), sep=",") + # add refseqAccession + summary['anchor_refseqAccession'] = summary['origin_species'].apply(lambda x: '_'.join(x.split('_')[-2:])) + summary = summary.merge(accessionmap[['assemblyAccession', 'refseqAccession_versioned']], left_on='anchor_refseqAccession', right_index=True, how='left') + # obtain genomic positions of genes + annotations = [] + for acc_assembly in tqdm(summary['assemblyAccession'].unique(), 'collect block-coordinates and species annotations for %s' % env): + annot = pd.read_csv(join(fp_genomes, '%s/genomic.gff' % acc_assembly), sep="\t", comment='#', header=None) + # subset to CDS + annot = annot[annot[2] == 'CDS'] + # get locus tag name from attributes field + annot['locus_tag'] = annot[8].apply(lambda x: [entry.split('=')[-1] for entry in x.split(';') if entry.split('=')[0] == 'locus_tag'][0]) + with open(join(fp_genomes, '%s/genomic.gbff' % acc_assembly), 'r') as f: + for line in f.readlines(): + if line.startswith(' ORGANISM'): + annot['organism'] = line.split(' ORGANISM ')[-1].strip() + break + annotations.append(annot) + annotations = pd.concat(annotations).set_index([0, 'locus_tag']) + + # add genomic positions for genes + summary = summary.merge(annotations[[3,4,6, 'organism']].reset_index(), left_on='gene_name', right_on='locus_tag', how='left').rename(columns={0: 'assembly_accession', 3: 'gene_genomic_start', 4: 'gene_genomic_stop', 6: 'strand'}) + #summary['assembly_length'] = accessionmap.reset_index().set_index('refseqAccession_versioned').loc[summary['assembly_accession'].values, 'length'].values + _map = accessionmap.reset_index().set_index('refseqAccession_versioned')['length'].to_dict() + summary['assembly_length'] = summary['assembly_accession'].apply(lambda x: _map.get(x, np.nan)) + + return summary + + +def get_genomic_hit_positons(fp_orb_basedir:str, fp_marbel_basedir:str, env:str, assembler:str, summary): + mapping = pd.read_csv('%s/%s/minimap2/%s_mapping.tsv' % (fp_orb_basedir, env, assembler), sep="\t", header=None).rename( + columns={0: 'Query sequence name', + 1: "Query sequence length", + 2: "Query start coordinate", + 3: "Query end coordinate", + 5: 'Target sequence name', + 7: 'Target start coordinate on the original strand', + 8: 'Target end coordinate on the original strand'}) + # obtain gene name from block name + mapping['gene_name'] = mapping['Target sequence name'].apply(lambda x: x.split('_block')[0]) + + # obtain relative positions of blocks regarding a gene + blocks = pd.read_csv(join(fp_marbel_basedir, '%s_microbiome' % env, 'summary', 'blocks.bed'), sep="\t", header=None, names=['gene_name', 'start', 'stop', 'block', 'unknown']) + + # add relative positions of blocks within genes + mapping = mapping.merge(blocks.set_index('block')[['start', 'stop']].rename(columns={'start': 'rel_block_start', 'stop': 'rel_block_stop'}), left_on='Target sequence name', right_index=True, how='left') + # add genomic gene positions + mapping = mapping.merge(summary.set_index('gene_name')[['gene_genomic_start', 'gene_genomic_stop', 'strand', 'assembly_accession', 'assembly_length', 'assemblyAccession', 'organism']], left_on='gene_name', right_index=True, how='left') + # compute genomic block positions + mapping['block_genomic_start'] = mapping['gene_genomic_start'] + mapping['rel_block_start'] + mapping['block_genomic_stop'] = mapping['block_genomic_start'] + (mapping['rel_block_stop'] - mapping['rel_block_start'] + 1) + + # compute genomic HIT positions + mapping['hit_start'] = mapping['block_genomic_start'] + mapping['Target start coordinate on the original strand'] + mapping['hit_stop'] = mapping['block_genomic_start'] + mapping['Target end coordinate on the original strand'] + + return mapping \ No newline at end of file From 5577544a1247ac875e0fbcc689b8dc8e49799b79 Mon Sep 17 00:00:00 2001 From: Stefan Janssen Date: Wed, 18 Mar 2026 23:25:44 +0100 Subject: [PATCH 10/13] as with compile_data_orb.py: f1 ranking and hit plotting --- plotting/plot_include_orb.py | 120 ++++++++++++++++++++++++++++++----- 1 file changed, 104 insertions(+), 16 deletions(-) diff --git a/plotting/plot_include_orb.py b/plotting/plot_include_orb.py index cbecbf5..4a183c5 100644 --- a/plotting/plot_include_orb.py +++ b/plotting/plot_include_orb.py @@ -8,6 +8,7 @@ from matplotlib import patches as mpatches from matplotlib.lines import Line2D from matplotlib_venn import venn2 +import matplotlib.ticker as ticker import colorsys import seaborn as sns from compile_data_orb import get_environments, getdata_DEgenes_mod, getdata_recovery, getdata_gene_recovery, getdata_runtime_memory, getdata_DEvennOrtho, getdata_DEorthogroups_mod, get_recovered_contigs, getdata_rnaquast, getdata_block_recovery, get_contigs, filter_for_assembler_with_ns, getdata_detonate @@ -18,9 +19,11 @@ from statannotations.Annotator import Annotator -def legendentries_col2rowordering(handles, labels, num_rows): +def legendentries_col2rowordering(handles, labels, num_rows, dummiesatpositions=[]): entries = list(zip(handles, labels)) num_columns = int(np.ceil(len(entries) / num_rows)) + for pos in dummiesatpositions: + entries.insert(pos, (Line2D([], [], linestyle="none"), "")) reordered_entries = [] for i in range(num_columns): reordered_entries.extend(entries[i::num_columns]) @@ -150,10 +153,10 @@ def plot_recovery(fp_orb_basedir, settings, num_columns:int=3, verbose=True, det # reorder entries in legend xoffset = 0 if detail_view is False: - handles, labels = legendentries_col2rowordering(handles, labels, 2) + handles, labels = legendentries_col2rowordering(handles, labels, 2, dummiesatpositions=[4]) else: - xoffset = 0.5 - ax_good.legend(handles, labels, ncol=4, bbox_to_anchor=(0.5 + xoffset, -0.3)) + xoffset = 0.2 + ax_good.legend(handles, labels, ncol=5, bbox_to_anchor=(1 + xoffset, -0.3)) else: ax_good.legend().remove() ax_bad.legend().remove() @@ -230,15 +233,20 @@ def plotTimeMemory(fp_caviar_basedir:str, settings, verbose=True): plotdata['unit'] = plotdata['value'] / factor plotdata = plotdata.sort_values(by='unit') sns.boxplot(data=plotdata, x='unit', y='assembler', orient='h', ax=ax, color='lightgray') - sns.stripplot(data=plotdata, x='unit', y='assembler', orient='h', ax=ax, hue='environment', hue_order=settings['labels']['environments']) + sns.stripplot(data=plotdata, x='unit', y='assembler', orient='h', ax=ax, hue='environment', hue_order=list(settings['labels']['environments'].values())) ax.set_xlabel(label) ax.set_ylabel("") ax.set_title(title) if ax != axes[-1]: + ax.axvline(x=4*24, linestyle=':', color='gray', zorder=-1, label="4 CPU days") + rt_handles, rt_labels = ax.get_legend_handles_labels() ax.legend().remove() else: ax.axvline(x=64, linestyle='-.', color='gray', zorder=-1, label="64 GB laptop") - ax.legend(bbox_to_anchor=(1.1, -0.25), ncols=7) + mem_handles, mem_labels = ax.get_legend_handles_labels() + mem_handles.insert(len(mem_handles)-1, rt_handles[-1]) + mem_labels.insert(len(mem_labels)-1, rt_labels[-1]) + ax.legend(mem_handles, mem_labels, bbox_to_anchor=(0.5, -0.25), ncols=4) # panel labels for i, ax in enumerate(axes): @@ -301,17 +309,17 @@ def _create_color(row): return cmp -def plot_DEgenes(fp_orb_basedir, settings, forOrthogroups=False, num_columns:int=3, verbose=True): +def plot_DEgenes(fp_orb_basedir, settings, forOrthogroups=False, num_columns:int=3, verbose=True, sortF1=False): # which environments to plot and in which order environments = get_environments(fp_orb_basedir, settings) # load data - DEfeatures = getdata_DEgenes_mod(fp_orb_basedir, settings, verbose) + DEfeatures = getdata_DEgenes_mod(fp_orb_basedir, settings, verbose, sortF1=sortF1) if forOrthogroups is False: rank_data = getdata_recovery(fp_orb_basedir, settings, verbose=verbose) else: rank_data = DEfeatures.loc[:, :, 'True Positive']['rank'].reset_index().rename(columns={'rank': 'recovery_rank'}).set_index('assembler').copy() - DEfeatures = getdata_DEorthogroups_mod(fp_orb_basedir, settings) + DEfeatures = getdata_DEorthogroups_mod(fp_orb_basedir, settings, sortF1=sortF1) fig, axes = plt.subplots( int(np.ceil(len(environments) / num_columns)), @@ -347,7 +355,10 @@ def plot_DEgenes(fp_orb_basedir, settings, forOrthogroups=False, num_columns:int tick_label.set_color(color) ax.set_title(settings['labels']['environments'].get(environment, environment)) - ax.set_xlabel('number genes') + if forOrthogroups: + ax.set_xlabel('number orthologous groups') + else: + ax.set_xlabel('number genes') num_positives = DEfeatures.loc[environment, :].reset_index().set_index('truth').loc[True].groupby('assembler')['num_genes'].sum().iloc[0] ax.axvline(x=num_positives, color=palette['True Positive'], label='Positive') @@ -582,7 +593,7 @@ def plot_rnaquast(fp_orb_basedir:str, fp_quast_basedir:str, settings, num_column ax_good.legend(handles=[ mpatches.Patch(color=settings['contig_classes']['multi_mapped_contigs_multi_og']['color'], label='Misassemblies'), mpatches.Patch(color=settings['contig_classes']['minimap2_single_recovered']['color'], label='95%-assembled isoforms')], - ncol=2, bbox_to_anchor=(-1.8, -0.20)) + ncol=2, bbox_to_anchor=(-0.5, -0.30)) ax_bad.text(0.06, 0.96, settings['labels']['environments'].get(environment, environment), transform=ax_bad.transAxes, @@ -642,7 +653,8 @@ def plot_detonate(fp_orb_basedir:str, fp_detonate_basedir:str, settings, num_col return fig -def plot_block_correlation(fp_orb_basedir, fp_marbel_basedir, sequence_file, settings, field, verbose=True, test = 'Mann-Whitney-gt', data_block_recovery=None, num_columns:int=3): + +def plot_block_correlation(fp_orb_basedir, fp_marbel_basedir, sequence_file, settings, field, verbose=True, test = 'Mann-Whitney-gt', data_block_recovery=None, num_columns:int=3, logscale=False): # field: block_length, read_mean_count, mean_identity, Coverage environments = get_environments(fp_orb_basedir, settings) @@ -658,13 +670,13 @@ def plot_block_correlation(fp_orb_basedir, fp_marbel_basedir, sequence_file, set data_block_recovery = getdata_block_recovery(fp_orb_basedir, fp_marbel_basedir, sequence_file, settings, verbose=verbose) assemblers = settings['labels']['assemblers'] # 'Mann-Whitney', 't-test_ind', 'Wilcoxon' - logscale = False + for i, environment in tqdm(enumerate(environments), disable=not verbose, desc=f'Drawing panels for {field} plot'): ax = axes[i] plotdata = data_block_recovery[environment] # assembler - ordered_assembler = sorted(list(plotdata["assembler"].unique())) + ordered_assembler = sorted(list(plotdata["assembler"].unique()), key=lambda x: x.lower()) if field == 'block_length': plotdata = plotdata[plotdata[field] > settings["read_length"]] @@ -742,14 +754,14 @@ def plot_nruns(fp_orb_basedir, settings, num_columns:int=3, verbose=True, contig ax_bad.xaxis.set_label_coords(1, -0.4) ax_top_bad = ax_bad.twiny() ax_top_bad.set_xticks([]) - ax_top_bad.set_xlabel("weak recovery") + ax_top_bad.set_xlabel(settings['labels']['recovery_plot']['bad']) ax_good = axes[i*2+1] sns.barplot(data=plot_data[plot_data['class'] != 'missed'], y='assembler', x=0, hue='max_n_longer_read_length', ax=ax_good) ax_good.tick_params(left = False) ax_good.set_yticks([]) ax_good.set_ylabel("") - ax_good.set_xlabel("robust recovery") + ax_good.set_xlabel(settings['labels']['recovery_plot']['good']) ax_good.xaxis.set_label_position('top') ax_good.set_xlim((0, plot_data[0].max()*1.1)) @@ -775,3 +787,79 @@ def plot_nruns(fp_orb_basedir, settings, num_columns:int=3, verbose=True, contig ax_good.legend().remove() ax_bad.legend().remove() return fig + + +def plot_hit_genomes_mapping(mapping, ax=None): + assert mapping['Query sequence name'].unique().shape[0] == 1 + CONTIG_GAP = 100000 + hitYpos = -2 + + plotinfo = mapping.copy() + + # one line per organism + for y, (taxon, g) in enumerate(plotinfo.groupby('assemblyAccession')): + plotinfo.loc[g.index, 'ypos'] = y + # one segment per genome assembly contig + x_pos = 0 + for i, (contig, g2) in enumerate(g.sort_values(by='assembly_length').groupby('assembly_accession')): + for field in ['gene_genomic_start', 'gene_genomic_stop', + 'block_genomic_start', 'block_genomic_stop', + 'hit_start', 'hit_stop']: + plotinfo.loc[g2.index, 'xpos_%s' % field] = x_pos + plotinfo.loc[g2.index, field] + plotinfo.loc[g2.index, 'xpos_assembly_contig_start'] = x_pos + plotinfo.loc[g2.index, 'xpos_assembly_contig_stop'] = x_pos + g2['assembly_length'].iloc[0] + x_pos += g2['assembly_length'].iloc[0] + CONTIG_GAP + + if ax is None: + fig, axes = plt.subplots(1, 1, figsize=(10, 3)) + ax = axes + else: + fig = None + + ylabels = [] + total_width = plotinfo['xpos_assembly_contig_stop'].max() + colors = sns.color_palette(n_colors=mapping.shape[0]) + color_index = 0 + for y, (taxon, g) in enumerate(plotinfo.groupby('assemblyAccession')): + # plot horizontal lines for the contigs of the genome assembly + ylabels.append((y, g['organism'].iloc[0] + "\n" + taxon)) + for _, asscontig in g.groupby('assembly_accession').head(1).iterrows(): + ax.plot([asscontig['xpos_assembly_contig_start'], asscontig['xpos_assembly_contig_stop']], + [asscontig['ypos'], asscontig['ypos']], color='black', linewidth=1, zorder=-1) + + # plot boxes (will be reduced to vertical lines in most cases) for genes in genomes + for gene, g_gene in g.sort_values(by='xpos_gene_genomic_start').groupby('gene_name').head(1).iterrows(): + width = max(g_gene['xpos_gene_genomic_stop'] - g_gene['xpos_gene_genomic_start'] + 1, total_width / 500) + ax.add_patch(mpatches.Rectangle((g_gene['xpos_gene_genomic_start'], g_gene['ypos'] - 0.2), width, 0.2*2)) + + for _, hit in g.sort_values(by='Query start coordinate').iterrows(): + hit_start, hit_stop = hit['Query start coordinate'], hit['Query end coordinate'] + if True: + hit_start = hit_start / hit['Query sequence length'] * total_width + hit_stop = hit_stop / hit['Query sequence length'] * total_width + ax.add_patch(mpatches.Polygon([[hit_start, hitYpos], [hit_stop, hitYpos], + [hit['xpos_hit_start'], hit['ypos']], [hit['xpos_hit_stop'], hit['ypos']]], + closed=True, facecolor=colors[color_index], edgecolor='k', linewidth=2, alpha=0.6)) + color_index += 1 + + # plot horizontal line for hit contig + ax.plot([0, total_width], [hitYpos, hitYpos], color='blue', linewidth=1, zorder=-1) + ylabels.append((hitYpos, plotinfo.iloc[0, :]['Query sequence name'])) + + ax.set_yticks([pos for pos, _ in ylabels], [label for _, label in ylabels]) + ax_top = ax.twiny() + + # give reasonable xticks for hit (but in genomic coordinates) + hit_xticks = ticker.MaxNLocator(nbins=7).tick_values(0, plotinfo.iloc[0, :]['Query sequence length']) + ax.set_xticks([x / hit['Query sequence length'] * total_width for x in hit_xticks], list(map(int, hit_xticks))) + ax.set_xlabel('contig position (bp)') + + genome_xticks = ticker.MaxNLocator(nbins=7).tick_values(0, total_width) + ax_top.set_xlim(ax.get_xlim()) + ax_top.set_xticks(genome_xticks, ['%.1f' % (x / 1000000) for x in genome_xticks]) + ax_top.set_xlabel('genomic position (Mbp)') + + if fig is None: + return ax + else: + return fig \ No newline at end of file From 02b7520ccfef292737d0ea33aeef0ac56c7ff2f4 Mon Sep 17 00:00:00 2001 From: Timo Wentong Lin Date: Tue, 31 Mar 2026 13:58:09 +0200 Subject: [PATCH 11/13] ci: remove nf-core related linting, remove run --- .github/workflows/ci.yml | 27 ---------------------- .github/workflows/linting.yml | 42 +---------------------------------- 2 files changed, 1 insertion(+), 68 deletions(-) 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..ede4165 100644 --- a/.github/workflows/linting.yml +++ b/.github/workflows/linting.yml @@ -65,44 +65,4 @@ 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 + allow-repeats: false \ No newline at end of file From 4d9b4c2444354707c4c034a3efd1a36641314041 Mon Sep 17 00:00:00 2001 From: Timo Wentong Lin Date: Tue, 31 Mar 2026 14:26:12 +0200 Subject: [PATCH 12/13] ci: fix prettier linting --- .github/ISSUE_TEMPLATE/bug_report.yml | 100 +++++++-------- .github/workflows/linting.yml | 2 +- .nf-core.yml | 46 +++---- CITATIONS.md | 32 +++-- README.md | 12 +- assets/schema_input.json | 2 +- nextflow_schema.json | 15 +-- plotting/style.yaml | 176 +++++++++++++------------- 8 files changed, 187 insertions(+), 198 deletions(-) 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/linting.yml b/.github/workflows/linting.yml index ede4165..0b6eb95 100644 --- a/.github/workflows/linting.yml +++ b/.github/workflows/linting.yml @@ -65,4 +65,4 @@ jobs: Thanks again for your contribution! repo-token: ${{ secrets.GITHUB_TOKEN }} - allow-repeats: false \ No newline at end of file + allow-repeats: false 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 Date: Tue, 31 Mar 2026 15:22:45 +0200 Subject: [PATCH 13/13] ci: black linting reformat --- bin/calculate_calour_confusion.py | 17 +- bin/calculate_dge_confusion.py | 22 +- bin/calculate_og_confusion.py | 17 +- bin/calculate_reference_info.py | 12 +- bin/calour_calculation.py | 6 +- bin/categorize_contigs.py | 113 +- bin/check_samplesheet.py | 8 +- bin/check_samplesheet_reads.py | 8 +- bin/count_genes_and_ogs_per_contig.py | 6 +- bin/gather_results.py | 90 +- bin/generate_blocks.py | 24 +- bin/make_histogramms.py | 1 - bin/merge_assembler_og_counts.py | 38 +- bin/merge_bowtie2_logs.py | 8 +- bin/merge_dataframes.py | 1 + bin/merge_og_counts.py | 4 +- bin/minimap2_classification.py | 263 ++-- bin/minimap2_overlap_selection.py | 80 +- bin/parse_bowtie2_logs.py | 8 +- bin/sort_dataframe.py | 2 +- bin/summarize_mapping_stats.py | 11 +- .../templates/dumpsoftwareversions.py | 17 +- plotting/compile_data_orb.py | 787 ++++++---- plotting/orthogroups_identity.py | 12 +- plotting/plot_include_orb.py | 1271 +++++++++++------ plotting/plot_orb_figures.py | 8 +- 26 files changed, 1845 insertions(+), 989 deletions(-) diff --git a/bin/calculate_calour_confusion.py b/bin/calculate_calour_confusion.py index 8ad82ea..f7abf2f 100755 --- a/bin/calculate_calour_confusion.py +++ b/bin/calculate_calour_confusion.py @@ -13,15 +13,15 @@ row_prefix = "calour_" -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") with open(map_file, "r") as file: contig_to_gene = json.load(file) assembler_dge["Name"] = assembler_dge["Name"].astype(str) - + assembler_dge["original_cds"] = assembler_dge["Name"].map(contig_to_gene) assembler_dge["gene_name"] = assembler_dge["original_cds"].str.replace(r"_block\d+$", "", regex=True) @@ -37,6 +37,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_dge_confusion.py b/bin/calculate_dge_confusion.py index 3efc10a..2c99731 100755 --- a/bin/calculate_dge_confusion.py +++ b/bin/calculate_dge_confusion.py @@ -15,8 +15,8 @@ log_name = sys.argv[6] row_prefix = sys.argv[7] -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") with open(map_file, "r") as file: @@ -36,7 +36,12 @@ if sum(group["de"]) > 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_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 7189818..34d02d4 100755 --- a/bin/categorize_contigs.py +++ b/bin/categorize_contigs.py @@ -23,17 +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" -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[ + 3 +] # "/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[ + 4 +] # "/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[ + 5 +] # "/vol/jlab/tlin/marbel_benchmarking_integration/benchmarking_sets_all_sparse_fixed_libsize/moss_microbiome/summary/gene_summary.csv" -prefix = sys.argv[6] # "rnaspades" +prefix = sys.argv[6] # "rnaspades" all_contigs_ids = pl.read_csv(all_contig_ids_fl, has_header=False, new_columns=["contigs"]) @@ -44,78 +54,81 @@ 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) -) +unassigned_contigs = all_contigs_ids.filter(~pl.col("contigs").is_in(mapped_ids)) assembler_contigs = assembler_mapping["contigs"].to_list() -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") -]) +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 76aaf9a..66be7af 100755 --- a/bin/gather_results.py +++ b/bin/gather_results.py @@ -44,19 +44,20 @@ 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() @@ -80,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), "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 e22871f..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] @@ -94,55 +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} -) +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() @@ -155,144 +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 e2f54e0..bb58fb2 100755 --- a/bin/minimap2_overlap_selection.py +++ b/bin/minimap2_overlap_selection.py @@ -15,7 +15,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"] @@ -23,42 +23,44 @@ 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] - minimap_path = sys.argv[1] seed = int(sys.argv[2]) prefix = sys.argv[3] @@ -75,19 +77,27 @@ def deduplicate_by_key(df, group_key, seed): "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_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] @@ -95,28 +105,20 @@ def deduplicate_by_key(df, group_key, seed): 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() @@ -127,4 +129,4 @@ def deduplicate_by_key(df, group_key, seed): minimap_mapping_chosen["category"] = "overlap_block" -minimap_mapping_chosen[["contig", "block_id", "category"]].to_csv(f"{prefix}_mapped_overlap_blocks.tsv", sep="\t") \ No newline at end of file +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/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("""\\