diff --git a/README.md b/README.md index d3a5383..9d8d29d 100644 --- a/README.md +++ b/README.md @@ -12,20 +12,22 @@ 4. [Installation](#install) 5. [Usage](#usage) 6. [Examples](#example) -7. [Citing](#citing) -8. [License](#license) +7. [Interactive plots](#interactive-plots) +8. [Citing](#citing) +9. [License](#license) ## Credit Written by Lauren Coombe ## Description -ntSynt-viz is an easy-to-use framework for generating ribbon plots combined with chromosome painting to visualize multi-genome synteny blocks. The tool is set-up to accept synteny blocks formatted in the [ntSynt](https://github.com/BirolLab/ntSynt) style, but any multi-genome synteny block file that adheres to the simple, BED-like TSV format of ntSynt can be visualized using ntSynt-viz. +ntSynt-viz is an easy-to-use framework for generating static and interactive ribbon plots combined with chromosome painting to visualize multi-genome synteny blocks. The tool is set-up to accept synteny blocks formatted in the [ntSynt](https://github.com/BirolLab/ntSynt) style, but any multi-genome synteny block file that adheres to the simple, BED-like TSV format of ntSynt can be visualized using ntSynt-viz. This flexible framework implements numerous features, including: * Option to normalize the strands of input chromosomes based on a target assembly * Synteny-guided ordering of assemblies from top-to-bottom, based on an input tree structure or distance estimates from the synteny blocks * Sorting chromosomes right-to-left based on synteny to adjacent assemblies * Colouring both the ribbons and chromosomes based on the target (top) assembly chromosomes +* Static and interactive ribbon plot outputs These features ensure that the output ribbon plots (powered by [gggenomes](https://thackl.github.io/gggenomes/)) are as easily understandable and as information-rich as possible. @@ -48,6 +50,8 @@ These features ensure that the output ribbon plots (powered by [gggenomes](https * [stringr](https://stringr.tidyverse.org/) * [ggplot2](https://ggplot2.tidyverse.org) * [svglite](https://cran.r-project.org/web/packages/svglite/index.html) + * [ggiraph](https://davidgohel.github.io/ggiraph/) + * [biopython](https://biopython.org/) ### Installing ntSynt-viz using conda ``` @@ -66,8 +70,8 @@ export PATH=/path/to/ntsynt-viz/github/ntSynt-viz/bin:$PATH ``` usage: ntsynt_viz.py [-h] --blocks BLOCKS --fais FAIS [FAIS ...] [--name_conversion NAME_CONVERSION] [--tree TREE] [--target-genome TARGET_GENOME] [--normalize] [--indel INDEL] [--length LENGTH] [--seq_length SEQ_LENGTH] [--keep KEEP [KEEP ...]] [--centromeres CENTROMERES] [--haplotypes HAPLOTYPES] - [--order ORDER] [--prefix PREFIX] [--format {png,pdf,svg}] [--scale SCALE] [--height HEIGHT] [--width WIDTH] [--dpi DPI] [--annotate-genome-info] - [--no-arrow] [--ribbon_adjust RIBBON_ADJUST] [-f] [-n] [-v] + [--order ORDER] [--prefix PREFIX] [--format {png,pdf,svg}] [--scale SCALE] [--height HEIGHT] [--width WIDTH] [--dpi DPI] + [--annotate-genome-info] [--optimize-ordering] [--no-arrow] [--ribbon_adjust RIBBON_ADJUST] [-f] [-n] [-v] Visualizing multi-genome synteny @@ -98,7 +102,6 @@ main plot formatting arguments: File listing haplotype assembly names: TSV, maternal/paternal assembly file names separated by tabs. --order ORDER Optional file specifying the order of genomes in the ribbon plot. If supplied, will override synteny distance-based ordering. If --tree supplied, the ordering must be compatible with the phylogenetic tree. - --no-arrow Only used with --normalize; do not draw arrows indicating reverse-complementation block filtering arguments: --indel INDEL Indel size threshold (used in computing synteny-based distances) [50000] @@ -118,6 +121,8 @@ output arguments: --dpi DPI Resolution of output plot - png output only [300] --annotate-genome-info Add annotations about number of sequences and genome size to the right of each genome in the ribbon plot + --optimize-ordering Optimize tree-guided genome sorting using inversions. Only use with strictly bifurcating trees. + --no-arrow Only used with --normalize; do not draw arrows indicating reverse-complementation --ribbon_adjust RIBBON_ADJUST Ratio for adjusting spacing beside ribbon plot. Increase if ribbon plot labels are cut off, and decrease to reduce the white space to the left of the ribbon plot [0.1] @@ -147,6 +152,25 @@ For more information about the output files from ntSynt-viz, check out our [wiki ## Using ntSynt-viz with pangenome graphs or synteny blocks from tools other than ntSynt To visualize synteny information in pangenome graphs or from synteny block detection tools other than ntSynt, the synteny blocks simply need to be converted to the straightforward [ntSynt format](https://github.com/BirolLab/ntsynt?tab=readme-ov-file#output-files). For convenience, we also provide some scripts to do this conversion in the `conversion_scripts` directory. +## Interactive plots +As of v1.1.0, ntSynt-viz generates an HTML file to allow interactive exploration of the ribbon plot. + +Interactive elements include: +* Hover/click on a ribbon to highlight the entire synteny block and display block coordinates and orientations +* Hover/click on a chromosome to display ID and length +* Click on one or more chromosomes in the legend to highlight the correseponding synteny blocks + +Example ribbon plot: XXX + + + + ## Citing Thank you for your [![Stars](https://img.shields.io/github/stars/BirolLab/ntSynt-viz.svg)](https://github.com/BirolLab/ntSynt-viz/stargazers) and for using and promoting this free software! We hope that ntSynt-viz (& ntSynt) is useful to you and your research. diff --git a/azure-pipelines.yml b/azure-pipelines.yml index 26d4796..6b29bb4 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -28,7 +28,7 @@ jobs: - script: | source activate ntsynt_viz_ci - mamba install --yes -c conda-forge -c bioconda quicktree r-base bioconductor-treeio r-ggpubr bioconductor-ggtree r-phytools r-dplyr r-argparse r-scales r-stringr pylint intervaltree snakemake r-ggplot2 r-gggenomes r-svglite + mamba install --yes -c conda-forge -c bioconda quicktree r-base bioconductor-treeio r-ggpubr bioconductor-ggtree r-phytools r-dplyr r-argparse r-scales r-stringr pylint intervaltree snakemake r-ggplot2 r-gggenomes r-svglite r-ggiraph biopython displayName: Install dependencies - script: | @@ -46,7 +46,7 @@ jobs: source activate ntsynt_viz_ci export PATH=$(pwd)/bin:$PATH cd tests - ntsynt_viz.py --blocks great-apes.ntSynt.synteny_blocks.tsv --fais fais.tsv --tree great-apes.mt-tree.nwk --name_conversion great-apes.name-conversions.tsv --normalize --prefix great-apes_ribbon-plots --ribbon_adjust 0.14 --scale 1e9 + ntsynt_viz.py --blocks great-apes.ntSynt.synteny_blocks.tsv --fais fais.tsv --tree great-apes.mt-tree.nwk --name_conversion great-apes.name-conversions.tsv --normalize --prefix great-apes_ribbon-plots --ribbon_adjust 0.14 --scale 1e9 --target-genome Homo_sapiens displayName: Run Example 1 test (PNG output) - script: | @@ -76,6 +76,9 @@ jobs: - publish: tests/great-apes_ribbon-plots_ribbon-plot_tree.png artifact: Example1 + - publish: tests/great-apes_ribbon-plots_ribbon-plot_tree.html + artifact: Example1_HTML + - publish: tests/great-apes_ribbon-plots_no-tree_ribbon-plot.pdf artifact: Example2 @@ -85,6 +88,9 @@ jobs: - publish: tests/great-apes_ribbon-plots_tree_order_ribbon-plot.png artifact: Example4 + - publish: tests/great-apes_ribbon-plots_tree_order_ribbon-plot.html + artifact: Example4_HTML + - job: displayName: macOS-latest @@ -113,7 +119,7 @@ jobs: - script: | source activate ntsynt_viz_ci - mamba install --yes -c conda-forge -c bioconda quicktree r-base bioconductor-treeio r-ggpubr bioconductor-ggtree r-phytools r-dplyr r-argparse r-scales r-stringr pylint intervaltree snakemake r-ggplot2 r-gggenomes r-svglite + mamba install --yes -c conda-forge -c bioconda quicktree r-base bioconductor-treeio r-ggpubr bioconductor-ggtree r-phytools r-dplyr r-argparse r-scales r-stringr pylint intervaltree snakemake r-ggplot2 r-gggenomes r-svglite r-ggiraph biopython displayName: Install dependencies - script: | diff --git a/bin/ntsynt_viz.py b/bin/ntsynt_viz.py index a61dc0c..8820cc0 100755 --- a/bin/ntsynt_viz.py +++ b/bin/ntsynt_viz.py @@ -67,7 +67,7 @@ def main(): action="store_true") block_filter_group.add_argument("--indel", help="Indel size threshold [50000]", default=50000, type=int) block_filter_group.add_argument("--length", help="Minimum synteny block length [100000]", default=100000, type=int) - block_filter_group.add_argument("--seq_length", help="Minimum sequence length [500000]", default=500000, type=float) + block_filter_group.add_argument("--seq_length", help="Minimum sequence length [500000]", default=500000, type=int) block_filter_group.add_argument("--keep", help="List of genome_name:chromosome to show in visualization. " "All chromosomes with links to the specified chromosomes will also be shown.", nargs="+", required=False, type=str) @@ -100,7 +100,11 @@ def main(): help="Add annotations about number of sequences " "and genome size to the right of each genome in the ribbon plot", action="store_true") - main_formatting_group.add_argument("--no-arrow", help="Only used with --normalize; " + output_group.add_argument("--optimize-ordering", + help="Optimize tree-guided genome sorting using inversions. " + "Only use with strictly bifurcating trees.", + action="store_true") + output_group.add_argument("--no-arrow", help="Only used with --normalize; " "do not draw arrows indicating reverse-complementation", action="store_true") output_group.add_argument( @@ -127,6 +131,8 @@ def main(): if len(args.fais) == 1: args.fais = read_fai_files(args.fais[0]) + print(f"Running ntSynt-viz {NTSYNT_VIZ_VERSION}...", flush=True) + cmd = f"snakemake -s {base_dir}/ntsynt_viz.smk " \ f"--cores 2 " \ f"--config " \ @@ -162,6 +168,8 @@ def main(): cmd += f"order={args.order} " if args.annotate_genome_info: cmd += "annotate_genome_info=True " + if args.optimize_ordering: + cmd += "optimize_ordering=True " if args.tree: cmd += f"tree={args.tree} " target = "gggenomes_ribbon_plot_tree" diff --git a/bin/ntsynt_viz.smk b/bin/ntsynt_viz.smk index cb4a132..ba7cb03 100644 --- a/bin/ntsynt_viz.smk +++ b/bin/ntsynt_viz.smk @@ -30,6 +30,8 @@ res = config.get("dpi", 300) orders = config.get("order", []) annotate_genome_info = config.get("annotate_genome_info", False) +optimize_ordering = config.get("optimize_ordering", False) + def sort_fais(fai_list, name_conversion, orders): "Based on the name conversion TSV, sort the FAIs based on orders" # Read in name conversions @@ -91,16 +93,20 @@ output_files = { } rule all: - input: output_files[format_img] + input: output_files[format_img], + f"{prefix}_ribbon-plot.html", + f"{prefix}_ribbon-plot_LLM-info.md", -ouput_files_tree = { +output_files_tree = { "png": f"{prefix}_ribbon-plot_tree.png", "pdf": f"{prefix}_ribbon-plot_tree.pdf", "svg": f"{prefix}_ribbon-plot_tree.svg", } rule gggenomes_ribbon_plot_tree: - input: ouput_files_tree[format_img] + input: output_files_tree[format_img], + f"{prefix}_ribbon-plot_tree.html", + f"{prefix}_ribbon-plot_tree_LLM-info.md" rule renaming: input: @@ -143,15 +149,25 @@ rule make_nj_tree: rule cladogram: input: - rules.make_nj_tree.output + tree = rules.make_nj_tree.output, + blocks = rules.renaming.output.renamed_blocks output: orders_tmp = temp(f"{prefix}_est-distances.order_tmp.tsv"), nwk_tmp = temp(f"{prefix}_est-distances_tmp.cladogram.nwk") params: prefix = f"{prefix}_est-distances", target = f"--target {target_genome}" if target_genome else "" - shell: - "ntsynt_viz_distance_cladogram.R --nwk {input} -p {params.prefix} --update_nwk {params.target}" + run: + if optimize_ordering: + shell(""" + set -eux -o pipefail + ntsynt_viz_distance_cladogram.R --nwk {input.tree} -p {params.prefix} --update_nwk {params.target} + ntsynt_viz_optimize_tree_topology.py --blocks {input.blocks} --tree {output.nwk_tmp} --out-order {output.orders_tmp} \ + --out-tree tmp_tree.nwk {params.target} + mv tmp_tree.nwk {output.nwk_tmp} + """) + else: + shell("set -eux -o pipefail; ntsynt_viz_distance_cladogram.R --nwk {input.tree} -p {params.prefix} --update_nwk {params.target}") rule orders: input: @@ -185,7 +201,7 @@ rule sort_blocks: output: sorted_blocks = f"{blocks_no_suffix}.renamed.sorted.blocks.tsv", intermediate_blocks = temp(f"{blocks_no_suffix}.renamed.sorted-tmp.tsv") if normalize else [], - chrom_oris = temp(f"{blocks_no_suffix}.renamed.sorted.chrom-orientations.tsv") if normalize else [] + chrom_oris = f"{blocks_no_suffix}.renamed.sorted.chrom-orientations.tsv" if normalize else [] params: intermediate_blocks = f"{blocks_no_suffix}.renamed.sorted-tmp.tsv", name_conversion = f"-c {name_conversion}" if name_conversion else "", @@ -250,7 +266,7 @@ rule chrom_paint: input: links = rules.gggenomes_files.output.links output: colour_feats = f"{prefix}.chrom-paint-feats.tsv" shell: - '''cat {input.links} |perl -ne 'chomp; @a=split("\t"); if(!defined $ct){{print "block_id\tseq_id\tbin_id\tstart\tend\tcolour_block\n"; $ct=1}} else {{print "$a[0]\t$a[1]\t$a[2]\t$a[3]\t$a[4]\t$a[10]\n"; print "$a[0]\t$a[5]\t$a[6]\t$a[7]\t$a[8]\t$a[10]\n";}}' > {output.colour_feats}''' + '''cat {input.links} |perl -ne 'chomp; @a=split("\t"); if(!defined $ct){{print "block_id\tseq_id\tbin_id\tstart\tend\tcolour_block\n"; $ct=1}} else {{print "$a[0]\t$a[1]\t$a[2]\t$a[3]\t$a[4]\t$a[12]\n"; print "$a[0]\t$a[5]\t$a[6]\t$a[7]\t$a[8]\t$a[12]\n";}}' > {output.colour_feats}''' rule ribbon_plot: input: @@ -260,7 +276,9 @@ rule ribbon_plot: haplotypes = rules.nudges.output.nudges if haplotypes else [], colour_seqs = rules.chrom_sorting.output.colour_info output: - out_img = output_files[format_img] + out_img = output_files[format_img] if format_img != "png" else [], + out_png = output_files["png"], + out_html = f"{prefix}_ribbon-plot.html" params: prefix = f"{prefix}_ribbon-plot", ratio = ribbon_ratio, @@ -287,7 +305,9 @@ rule ribbon_plot_tree: haplotypes = rules.nudges.output.nudges if haplotypes else [], colour_seqs = rules.chrom_sorting.output.colour_info output: - out_img = ouput_files_tree[format_img] + out_img = output_files_tree[format_img] if format_img != "png" else [], + out_png = output_files_tree["png"], + out_html = f"{prefix}_ribbon-plot_tree.html" params: prefix = f"{prefix}_ribbon-plot_tree", ratio = ribbon_ratio, @@ -303,3 +323,39 @@ rule ribbon_plot_tree: "ntsynt_viz_plot_synteny_blocks_ribbon_plot.R -s {input.sequences} -l {input.links} -p {params.prefix} --tree {input.tree}" " --ratio {params.ratio} --scale {params.scale} -c {input.colour_feats} --format {params.out_img_format} --height {params.height} --width {params.width}" " --order {input.orders} {params.centromeres} {params.arrow} {params.haplotypes} --colour_indices {input.colour_seqs} {params.resolution} {params.annotate_genome_info}" + +rule llm_instructions: + input: + ribbon_png = rules.ribbon_plot.output.out_png, # png + chrom_oris = f"{blocks_no_suffix}.renamed.sorted.chrom-orientations.tsv" if normalize else [], + synteny_tsv = rules.sort_blocks.output.sorted_blocks, + chrom_lengths = rules.chrom_sorting.output.sorted_seqs, + output: + markdown = f"{prefix}_ribbon-plot_LLM-info.md" + params: + normalize_opt = lambda wc, input: ( + f"--normalize {input.chrom_oris}" if normalize else "" + ), + top_n = max(50, len(fais)*2) + shell: + """ + ntsynt_viz_generate_llm_markdown.py -o {output.markdown} --top-n {params.top_n} --image {input.ribbon_png} --lengths {input.chrom_lengths} {params.normalize_opt} {input.synteny_tsv} + """ + +rule llm_instructions_tree: + input: + ribbon_png = rules.ribbon_plot_tree.output.out_png, # png + chrom_oris = f"{blocks_no_suffix}.renamed.sorted.chrom-orientations.tsv" if normalize else [], + synteny_tsv = rules.sort_blocks.output.sorted_blocks, + chrom_lengths = rules.chrom_sorting.output.sorted_seqs, + output: + markdown = f"{prefix}_ribbon-plot_tree_LLM-info.md" + params: + normalize_opt = lambda wc, input: ( + f"--normalize {input.chrom_oris}" if normalize else "" + ), + top_n = max(50, len(fais)*2) + shell: + """ + ntsynt_viz_generate_llm_markdown.py -o {output.markdown} --top-n {params.top_n} --image {input.ribbon_png} --top-n 50 --lengths {input.chrom_lengths} {params.normalize_opt} {input.synteny_tsv} + """ diff --git a/bin/ntsynt_viz_format_blocks_gggenomes.py b/bin/ntsynt_viz_format_blocks_gggenomes.py index d5ca58f..218b27c 100755 --- a/bin/ntsynt_viz_format_blocks_gggenomes.py +++ b/bin/ntsynt_viz_format_blocks_gggenomes.py @@ -94,7 +94,7 @@ def make_links_file(synteny_file, prefix, valid_blocks_set, target_assembly): target_assembly_chrom = None with open(f"{prefix}.links.tsv", 'w', encoding="utf-8") as fout: fout.write("block_id\tseq_id\tbin_id\tstart\tend\t"\ - "seq_id2\tbin_id2\tstart2\tend2\tstrand\tcolour_block\n") + "seq_id2\tbin_id2\tstart2\tend2\tstrand\tstrand1\tstrand2\tcolour_block\n") with open(synteny_file, 'r', encoding="utf-8") as fin: for line in fin: line = line.strip().split("\t") @@ -104,7 +104,8 @@ def make_links_file(synteny_file, prefix, valid_blocks_set, target_assembly): start_curr, end_curr = curr_block.start, curr_block.end relative_ori = "-" if curr_block.strand != prev_line.strand else "+" out_line = f"{curr_block.id}\t{prev_line.chrom}\t{prev_line.genome}\t{start_prev}\t{end_prev}\t"\ - f"{curr_block.chrom}\t{curr_block.genome}\t{start_curr}\t{end_curr}\t{relative_ori}" + f"{curr_block.chrom}\t{curr_block.genome}\t{start_curr}\t{end_curr}\t{relative_ori}\t"\ + f"{prev_line.strand}\t{curr_block.strand}" lines_to_print.append(out_line) if prev_line is not None and prev_line.id != curr_block.id: diff --git a/bin/ntsynt_viz_generate_llm_markdown.py b/bin/ntsynt_viz_generate_llm_markdown.py new file mode 100755 index 0000000..85a9577 --- /dev/null +++ b/bin/ntsynt_viz_generate_llm_markdown.py @@ -0,0 +1,559 @@ +#!/usr/bin/env python3 +""" +Generate a markdown "context" file summarizing an ntSynt ribbon plot, +intended to be handed to an LLM (of the user's choice) to produce a +short manuscript-ready summary. Does NOT call an LLM itself. +""" +import argparse +import csv +from collections import defaultdict +from dataclasses import dataclass +from pathlib import Path + +LLM_INSTRUCTIONS = """ +## Instructions +#### Output format +The tables below describe a multi-genome synteny (ribbon plot) comparison, +for use in a scientific manuscript. Using ONLY the data provided, write +2-3 sentences in clear, formal, manuscript-appropriate language summarizing +overall structural trends: the overall degree of synteny coverage/density +(how much of each genome is captured in synteny blocks); whether the blocks that are present +are largely collinear or show substantial structural rearrangement; any +notable large-scale fusions/fissions and which genomes are involved; and (optionally) +any standout individual inversions or genomes with a disproportionate +share of inversion relative to the others. + +Coverage/density and collinearity are two separate, both-reportable +questions -- do not let one substitute for the other. A set of genomes +can have blocks that stay in consistent order (collinear) while still +covering only a small, fragmented fraction of each chromosome's length +(low coverage), or vice versa. If coverage is low, sparse, or highly +variable across genomes, state this explicitly as a primary trend -- +do not omit it in favor of describing rearrangement/order alone. + +#### Synteny coverage / density description rules +Use the coverage table below (per-genome total syntenic length vs. +assembly length vs. % covered) to ground any statement about coverage. +Describe the overall pattern (e.g., coverage is low/high across most +genomes, or coverage varies substantially by genome/clade) rather than +listing each genome's percentage. Only call out an individual genome's +coverage figure if it stands out as a clear outlier from the rest. + +#### Fusion / fission description rules +When describing a specific fusion or fission, always name the chromosome +accession(s) involved for both genomes (e.g., \"A. example1 XXX.1 corresponds +to YYY.1 and ZZZ.1 in A. example2\"), not just the genomes. +When a chromosome in one genome corresponds to two chromosomes in all other +genomes, describe it as a single event in the odd-one-out genome (fusion or +fission, whichever direction is more parsimonious) rather than as independent +events in each of the other genomes. + +Chromosome correspondence table: use it only to identify fusions/fissions. +Do not report the % syntenic length values directly. Only pairs contributing a non-trivial share of a chromosome's +syntenic content are informative for calling a fusion/fission -- the +% column distinguishes small/stray correspondences from real ones. + +General: do not restate raw numbers verbatim from any table -- synthesize a +trend instead, and do not enumerate every row. + +#### Companion image usage +Write the summary as a direct visual description of the companion image -- +coverage/density, overall complexity, and any specific fusions, fissions, +or inversions should be described as things seen in the plot, not as +things reported in a table. +The whole image should be considered, not just the top genome or a single chromosome. + +When assessing coverage/density visually, note whether ribbons run +near-continuously along the full length of each chromosome bar (high +coverage) or appear as thin, scattered segments separated by substantial +uncolored/gray gaps (low coverage). This is a distinct, reportable +feature in its own right, separate from whether ribbons cross +(inversion) or split (fusion/fission) -- a plot can show clean, +uncrossed ribbons that still cover only a small fraction of each bar, +and that sparseness should be described rather than treated as +"collinear" or left unmentioned. + +Use the tables only silently, in the background, to: +- confirm that a visual feature (including apparent low coverage) is + genuine and not a rendering artifact, +- obtain the correct chromosome/accession identifiers to name it, +- decide whether an event is large enough to be worth describing. + +Never reference the tables, the markdown file, or the act of cross-referencing +in the output. Do not use phrases like "consistent with," "table-derived," +"the data shows," or similar meta-commentary -- state observations directly +(e.g., "Chromosome 8 shows an inversion in Han1..." not "the table confirms +an inversion consistent with..."). + +If a genuine event documented in the tables is too small to be discernible in +the image at this scale, omit it from the summary rather than describing an +invisible feature -- unless it is among the most significant events overall, +in which case describe it in general terms (e.g., "a smaller-scale inversion +is also present on...") without claiming it's visible. + +If something appears visually prominent in the image, verify it against the +tables before naming a specific chromosome/accession, but write the +conclusion as a plain visual observation. + + +#### Multi-genome framing +Remember that the comparisons are multi-genome, meaning that every synteny +block contains coordinates from each genome -- there is no concept of +separate pairwise comparisons when interpreting the plot. +In addition, the synteny blocks were generated using +an alignment-free approach, so always avoid any mention of +alignments (aligned, unaligned, etc.). + +#### Reverse-complementation arrows are not structural variation +Arrows under a chromosome are already explained in the image legend. +Do not mention arrows in the summary at all -- not their presence, +count, density, distribution, which genomes or chromosomes show them, +or what they mean. Never use arrow count or frequency as evidence of +structural variation, rearrangement, or any difference in genome +quality/complexity between genomes, and never compare genomes' arrow +patterns against each other, even indirectly (e.g., "several +chromosomes required reorientation in X" is still off-limits). +Arrows are a rendering/orientation detail for the reader to consult +the legend on, not a summary-worthy feature. + +#### Phylogenetic tree rules +If a phylogenetic tree is provided, derive all statements about which +genomes are \"closely related,\" \"sister taxa,\" etc. strictly from the tree +topology -- never from the structural consistency or correspondence tables. +Use the tables only to say whether structural patterns are consistent or +inconsistent with the tree relationships, not to define them. + +#### General language rules +Avoid non-scientific, unfalsifiable, or aesthetic language throughout, +including superlatives (exceptional, amazing, remarkable) and subjective +visual descriptors (modest, striking, subtle, understated, overwhelming). + +#### Disclaimer +Always print a separate disclaimer that a human expert should review the +summary before including it in a presentation or publication. +""" + +SYNTENY_COLS = [ + "block_id", "genome", "chrom", "start", "end", + "strand", "num_minimizers", "reason", +] + +@dataclass +class SyntenyRow: + """Describes a row in the synteny TSV file""" + block_id: int + genome: str + chrom: str + start: int + end: int + strand: str + num_minimizers: int + reason: str + +def load_synteny_tsv(path): + """Read synteny TSV file""" + rows = [] + with open(path, newline="", encoding="utf-8") as f: + reader = csv.reader(f, delimiter="\t") + for fields in reader: + if not fields: + continue + row = SyntenyRow( + block_id=int(fields[0]), + genome=fields[1], + chrom=fields[2], + start=int(fields[3]), + end=int(fields[4]), + strand=fields[5], + num_minimizers=int(fields[6]), + reason=fields[7] + ) + rows.append(row) + return rows + + +def load_normalization_tsv(path): + """Load the information about chromosome normalization""" + with open(path, newline="", encoding="utf-8") as f: + lines = [line for line in f if not line.startswith("#")] + reader = csv.reader(lines, delimiter="\t") + next(reader) # header + rows = [] + for fields in reader: + if not fields: + continue + rows.append({ + "genome": fields[0], + "chromosome": fields[1], + "relative_orientation": fields[2], + }) + return rows + +def read_lengths_tsv(path): + """Read a TSV with headers bin_id,seq_id,length,relative_orientation and return a dict of lengths.""" + lengths = defaultdict(dict) + with open(path, newline="", encoding="utf-8") as f: + reader = csv.reader(f, delimiter="\t") + next(reader) # header + for fields in reader: + if not fields: + continue + bin_id = fields[0] + seq_id = fields[1] + length = int(fields[2]) + lengths[bin_id][seq_id] = length + return lengths + + +def get_genome_order(rows): + """Genome order is fixed and identical within every block; read it + off the first block.""" + first_block_id = rows[0].block_id + order = [] + for row in rows: + if row.block_id != first_block_id: + break + order.append(row.genome) + return order + + +def normalization_summary(norm_rows): + """Per-genome fraction of chromosomes reverse-complemented during + normalization.""" + counts = defaultdict(lambda: [0, 0]) # genome -> [flipped, total] + for row in norm_rows: + counts[row["genome"]][1] += 1 + if row["relative_orientation"] == "-": + counts[row["genome"]][0] += 1 + return {genome: tuple(v) for genome, v in counts.items()} + + +def chromosome_correspondence(rows, genome_order, chrom_lengths): + """For each genome pair and direction, syntenic bp and % of genome A + chromosome's total syntenic length, per (chrom_a, chrom_b).""" + results = [] + genomes = genome_order + + by_genome = defaultdict(dict) + for r in rows: + by_genome[r.genome][r.block_id] = r + + for ga in genomes: + for gb in genomes: + if ga == gb: + continue + a_rows = by_genome[ga] + b_rows = by_genome[gb] + + pair_totals = defaultdict(int) + for block_id, a_row in a_rows.items(): + b_row = b_rows.get(block_id) + if b_row is None: + continue + chrom_a = a_row.chrom + chrom_b = b_row.chrom + length_a = a_row.end - a_row.start + pair_totals[(chrom_a, chrom_b)] += length_a + + for (chrom_a, chrom_b), bp in pair_totals.items(): + pct = 100 * bp / chrom_lengths[ga][chrom_a] + results.append({ + "genome_a": ga, "chrom_a": chrom_a, + "genome_b": gb, "chrom_b": chrom_b, + "syntenic_bp": bp, "pct_of_chrom_a": round(pct, 1), + }) + return results + + +def coverage_summary(rows, chrom_lengths, genome_order): + """Per-genome total syntenic length (bp of that genome's own + chromosomes covered by at least one synteny block) vs. that genome's + total assembly length, expressed as a % covered. + + This is distinct from the chromosome correspondence table: it answers + "how much of this genome is in a synteny block at all," not "how do + two genomes' chromosomes map onto each other." + """ + syntenic_totals = defaultdict(int) + for r in rows: + syntenic_totals[r.genome] += (r.end - r.start) + + results = [] + for genome in genome_order: + assembly_length = sum(chrom_lengths.get(genome, {}).values()) + syntenic_length = syntenic_totals.get(genome, 0) + pct = ( + round(100 * syntenic_length / assembly_length, 1) + if assembly_length > 0 else None + ) + results.append({ + "genome": genome, + "assembly_length": assembly_length, + "syntenic_length": syntenic_length, + "pct_covered": pct, + }) + return results + + +def inversion_summary(rows, genome_order, top_n=5): + """Per-genome proportion of syntenic length that is inverted relative + to the first genome in the provided order (strand '-'), plus the + top_n largest individually inverted blocks.""" + first_genome = genome_order[0] + + totals = defaultdict(int) + inverted_totals = defaultdict(int) + for r in rows: + totals[r.genome] += (r.end - r.start) + if r.strand == "-": + inverted_totals[r.genome] += (r.end - r.start) + + per_genome = [] + for genome in genome_order: + total_len = totals.get(genome, 0) + inv_len = inverted_totals.get(genome, 0) + pct = round(100 * inv_len / total_len, 1) if total_len > 0 else None + per_genome.append({ + "genome": genome, + "inverted_length": inv_len, + "pct_inverted": pct, + }) + + inverted_blocks = [ + r for r in rows if r.strand == "-" and r.genome != first_genome + ] + inverted_blocks.sort(key=lambda r: (r.end - r.start), reverse=True) + top = inverted_blocks[:top_n] + + return per_genome, top, first_genome + + +def build_header(): + """Build the markdown header""" + return [ + "# Synteny Ribbon Plot Summary Context", + "", + "", + "", + ] + + +def build_about_section(first_genome): + """Build the about information section""" + about_text = f""" +## About this plot +This data describes a multi-genome synteny comparison. A **synteny +block** is a genomic region conserved between genomes; **coverage/ +density** refers to what fraction of a genome's total length is +captured in synteny blocks at all, independent of block order; an +**inversion** is a block in reversed orientation relative to +{first_genome}, the first-listed genome (chosen only for its position +in the input order, not as a biological reference); a **fusion/ +fission** is a case where regions from one chromosome in one genome +correspond to multiple chromosomes in another. + +If a companion image is provided: each chromosome in the top genome is assigned a distinct +color; ribbon width corresponds to syntenic block length; a +twisted/crossed ribbon indicates an inverted block; uncolored/gray +regions of a chromosome bar are not captured in any synteny block; +genomes are arranged in the order listed below. +Optionally, a phylogenetic tree is rendered to the left of the ribbon plot. +Arrows under chromosomes indicate reverse complementation during normalization (also optional). + """ + return [about_text] + +def build_genome_table(genome_order, norm_summary): + """Build genome table for markdown""" + lines = [ + "## Genomes Compared", + "", + "Listed in the order shown in the ribbon plot.", + "", + "| Order | Genome | Normalization (optional) |", + "|---|---|---|", + ] + + for i, genome in enumerate(genome_order, 1): + if genome not in norm_summary: + note = "—" + else: + flipped, total = norm_summary[genome] + note = ( + "—" + if flipped == 0 + else f"{flipped}/{total} chromosomes reverse-complemented during normalization" + ) + + lines.append(f"| {i} | {genome} | {note} |") + + lines.append("") + return lines + + +def build_coverage_table(coverage_results): + """Build synteny coverage/density table for markdown""" + lines = [ + "## Synteny Coverage Summary", + "", + "Per-genome total length captured in synteny blocks, relative to " + "assembly length. Distinct from the chromosome correspondence " + "table below: this measures how much of each genome participates " + "in any synteny block at all, not how chromosomes map onto each " + "other.", + "", + "| Genome | Assembly Length | Total Syntenic Length | % Covered |", + "|---|---|---|---|", + ] + + for row in coverage_results: + pct = "N/A" if row["pct_covered"] is None else f"{row['pct_covered']}%" + lines.append( + f"| {row['genome']} | " + f"{row['assembly_length']} | " + f"{row['syntenic_length']} | " + f"{pct} |" + ) + + lines.append("") + return lines + + +def build_chromosome_table(chrom_results): + """Build chromosome correspondence table for markdown""" + lines = [ + "## Chromosome Correspondence (fusions / fissions)", + "", + "Only pairs contributing a non-trivial share:", + "", + "| Genome A | Chr A | Genome B | Chr B | Syntenic bp | % of Chr A Syntenic Length |", + "|---|---|---|---|---|---|", + ] + + for row in chrom_results: + lines.append( + f"| {row['genome_a']} | " + f"{row['chrom_a']} | " + f"{row['genome_b']} | " + f"{row['chrom_b']} | " + f"{row['syntenic_bp']} | " + f"{row['pct_of_chrom_a']}% |" + ) + + lines.append("") + return lines + +def build_inversion_table(inv_per_genome, inv_top, first_genome): + """Build table summarizing inversions""" + lines = [ + "## Inversion Summary", + "", + f"Relative to {first_genome} (first-listed genome; see note above).", + "", + "**Per-genome inverted proportion:**\n", + "| Genome | Inverted Length | % of Syntenic Length Inverted |", + "|---|---|---|" + ] + for row in inv_per_genome: + pct = "N/A" if row["pct_inverted"] is None else f"{row['pct_inverted']}%" + lines.append(f"| {row['genome']} | {row['inverted_length']} | {pct} |") + + lines.extend( + ["", + f"**Largest individual inversions (top {len(inv_top)}):**\n", + "| Genome | Chromosome | Location | Length |", + "|---|---|---|---|", + ] + ) + for row in inv_top: + lines.append(f"| {row.genome} | {row.chrom} | {row.start}-{row.end} | {row.end - row.start} |") + lines.append("") + + return lines + +def build_companion_image(image_path): + """Build companion image section""" + lines = [ + "## Companion Image (optional)", + "", + ] + + if image_path: + lines.append(f"`{image_path}`") + else: + lines.append("Not provided.") + + lines.extend([ + "", + """ +If available, provide this alongside the markdown file to a +vision-capable LLM for an additional, purely qualitative visual +impression -- it is not required to produce a valid summary. + """, + "", + ]) + + return lines + + +def build_markdown(genome_order, norm_summary, coverage_results, chrom_results, + inv_per_genome, inv_top, first_genome, image_path=None,): + """Build the markdown file""" + sections = [ + build_header(), + [LLM_INSTRUCTIONS], + build_about_section(first_genome), + build_coverage_table(coverage_results), + build_genome_table(genome_order, norm_summary), + build_chromosome_table(chrom_results), + build_inversion_table(inv_per_genome, inv_top, first_genome), + build_companion_image(image_path), + ] + + lines = [] + + for section in sections: + lines.extend(section) + + return "\n".join(lines) + + +def main(): + """Generate a LLM markdown file for the ntSynt-viz ribbon plot""" + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("synteny_tsv") + parser.add_argument("--normalize", + help="Path to normalization TSV (optional; " + "only needed if chromosomes were reverse-complemented during normalization)") + parser.add_argument("-o", "--output", default="ribbon_plot_summary_context.md") + parser.add_argument("--image", default=None, help="Path to companion PNG/PDF") + parser.add_argument("--top-n", type=int, default=5) + parser.add_argument("--lengths", type=str, required=True, + help="Path to TSV with headers bin_id,seq_id,length,relative_orientation") + args = parser.parse_args() + + rows = load_synteny_tsv(args.synteny_tsv) + norm_rows = [] + if args.normalize: + norm_rows = load_normalization_tsv(args.normalize) + genome_order = get_genome_order(rows) + norm_summary = normalization_summary(norm_rows) + chrom_lengths = read_lengths_tsv(args.lengths) + + chrom_results = chromosome_correspondence(rows, genome_order, chrom_lengths) + coverage_results = coverage_summary(rows, chrom_lengths, genome_order) + inv_per_genome, inv_top, first_genome = inversion_summary(rows, genome_order, top_n=args.top_n) + + md = build_markdown(genome_order, norm_summary, coverage_results, chrom_results, + inv_per_genome, inv_top, first_genome, image_path=args.image) + + Path(args.output).write_text(md, encoding="utf-8") + print(f"Wrote {args.output}") + + +if __name__ == "__main__": + main() diff --git a/bin/ntsynt_viz_optimize_tree_topology.py b/bin/ntsynt_viz_optimize_tree_topology.py new file mode 100755 index 0000000..0156752 --- /dev/null +++ b/bin/ntsynt_viz_optimize_tree_topology.py @@ -0,0 +1,372 @@ +#!/usr/bin/env python3 +""" +optimize_tree_topology.py + +Given an ntSynt synteny blocks TSV and a Newick tree, find the tree rotation +that minimizes the total pairwise inversion length between adjacent genomes. + +Outputs: + - TSV: one genome name per line, in optimized top-to-bottom order + - Newick: tree with subtrees rotated to match the optimal ordering +""" + +import argparse +import io +import sys +from collections import defaultdict +from Bio import Phylo +from Bio.Phylo.BaseTree import Clade, Tree + + +# --------------------------------------------------------------------------- +# Column indices (0-based) for the synteny blocks TSV +# --------------------------------------------------------------------------- +COL_BLOCK_ID = 0 +COL_GENOME = 1 +COL_CHROM = 2 +COL_START = 3 +COL_END = 4 +COL_STRAND = 5 +COL_MINIMIZERS = 6 +COL_DISCONTINUITY = 7 + + +def parse_args() -> argparse.Namespace: + "Parse arguments" + parser = argparse.ArgumentParser( + description="Optimize phylogenetic tree rotation to minimize visible inversions " + "in ntSynt-viz ribbon plots." + ) + parser.add_argument( + "--blocks", required=True, + help="ntSynt synteny blocks TSV file." + ) + parser.add_argument( + "--tree", required=True, + help="Newick tree file. Leaf names must match genome file names in the blocks TSV." + ) + parser.add_argument( + "--out-order", required=True, + help="Output TSV: one genome name per line, optimized top-to-bottom order." + ) + parser.add_argument( + "--out-tree", required=True, + help="Output Newick file with subtrees rotated to match the optimal ordering." + ) + parser.add_argument( + "--target", default=None, + help="If specified, this genome is always placed first (top-most) in the ordering." + ) + return parser.parse_args() + + +# --------------------------------------------------------------------------- +# Biopython helpers +# +# --------------------------------------------------------------------------- + +def _is_leaf(clade: Clade) -> bool: + return clade.is_terminal() + +def _get_children(clade: Clade) -> list[Clade]: + return clade.clades + +def _get_leaf_names(clade: Clade) -> list[str]: + return [leaf.name for leaf in clade.get_terminals()] + +def load_tree(newick_path: str) -> Tree: + """ + Read a Newick tree with Biopython. If the tree is unrooted (or has a + trifurcating/polytomous root, as is common from tools like RAxML/IQ-TREE), + midpoint rooting is applied, which also resolves the root into a + strictly bifurcating split. + """ + with open(newick_path, 'r', encoding="utf-8") as fh: + newick_str = fh.read().strip() + + tree = Phylo.read(io.StringIO(newick_str), "newick") + + if len(tree.root.clades) != 2: + print( + " Tree root is not bifurcating (likely unrooted input); " + "applying midpoint rooting.", + file=sys.stderr, + ) + tree.root_at_midpoint() + + if not tree.is_bifurcating(): + print("WARNING: --optimize-ordering option is not compatible with polytomies") + + return tree + +# --------------------------------------------------------------------------- +# Step 1: Build pairwise inversion cost matrix +# --------------------------------------------------------------------------- + +def _tally_block( + entries: list[tuple[str, str, int]], + cost: dict[tuple[str, str], int], +) -> None: + """ + Given all rows for a single synteny block, tally inverted pairs + directly into `cost`. Called once per block, then entries is discarded. + """ + for i, el_i in enumerate(entries): + for j in range(i + 1, len(entries)): + genome_a, strand_a, length_a = el_i + genome_b, strand_b, length_b = entries[j] + if genome_a == genome_b: + continue + if strand_a != strand_b: + inversion_length = (length_a + length_b) // 2 + cost[(genome_a, genome_b)] += inversion_length + cost[(genome_b, genome_a)] += inversion_length + + +def build_cost_matrix( + blocks_path: str, +) -> dict[tuple[str, str], int]: + """ + Stream the ntSynt blocks TSV one block at a time, tallying the total + length of inverted synteny blocks for each genome pair. + + Rows are assumed to be sorted by block ID (consecutive rows with the + same ID belong to the same block). Only the current block's rows are + held in memory at any point. + + Returns: + cost[(genome_a, genome_b)] = total inverted block length (symmetric) + """ + cost: dict[tuple[str, str], int] = defaultdict(int) + current_block_id: str | None = None + current_entries: list[tuple[str, str, int]] = [] + + with open(blocks_path, 'r', encoding="utf-8") as fh: + for line in fh: + fields = line.strip().split("\t") + block_id = fields[COL_BLOCK_ID] + genome = fields[COL_GENOME] + start = int(fields[COL_START]) + end = int(fields[COL_END]) + strand = fields[COL_STRAND] + length = abs(end - start) + + if block_id != current_block_id: + if current_entries: + _tally_block(current_entries, cost) + current_block_id = block_id + current_entries = [] + + current_entries.append((genome, strand, length)) + + # Flush the final block + if current_entries: + _tally_block(current_entries, cost) + + return dict(cost) + + +def get_genomes_from_cost_matrix( + cost: dict[tuple[str, str], int] +) -> list[str]: + "Returns the genomes from the cost matrix" + genomes: set[str] = set() + for a, b in cost: + genomes.add(a) + genomes.add(b) + return sorted(genomes) + + +# --------------------------------------------------------------------------- +# Step 2: Exhaustive memoized tree rotation optimization +# --------------------------------------------------------------------------- + +# Each subtree result is a list of options: +# (internal_cost, topmost_leaf, bottommost_leaf, ordering_tuple) +SubtreeOption = tuple[int, str, str, tuple[str, ...]] + +def build_subtree_options( + node: Clade, + cost: dict[tuple[str, str], int], + target_genome: str | None = None, +) -> list[SubtreeOption]: + """ + Recursively compute all valid leaf orderings for the subtree rooted at + `node`, memoizing results so each clade is computed exactly once. + + If target_genome is specified, any option where it would not appear as + the top leaf is discarded at nodes whose subtree contains it. + + Returns a list of (total_internal_cost, top_leaf, bottom_leaf, ordering). + """ + if _is_leaf(node): + return [(0, node.name, node.name, (node.name,))] + + children = _get_children(node) + if len(children) != 2: + raise ValueError( + f"Tree must be strictly bifurcating; node '{node.name}' " + f"has {len(children)} children. The option --optimize-ordering is not compatible with polytomies." + ) + + top_options = build_subtree_options(children[0], cost, target_genome) + bottom_options = build_subtree_options(children[1], cost, target_genome) + + # Determine whether the target genome falls in the top or bottom subtree, + # so we know which side must come first at this node. + top_leaf_names = set(_get_leaf_names(children[0])) + bottom_leaf_names = set(_get_leaf_names(children[1])) + target_in_top = target_genome in top_leaf_names if target_genome else False + target_in_bottom = target_genome in bottom_leaf_names if target_genome else False + + results: list[SubtreeOption] = [] + + for top_cost, top_top, top_bottom, top_order in top_options: + for bottom_cost, bottom_top, bottom_bottom, bottom_order in bottom_options: + + # Option A: [top subtree | bottom subtree] + # Skip if target must be topmost but would not be. + if not (target_in_bottom and not target_in_top): + boundary_a = cost.get((top_bottom, bottom_top), 0) + results.append(( + top_cost + bottom_cost + boundary_a, + top_top, + bottom_bottom, + top_order + bottom_order, + )) + + # Option B: [bottom subtree | top subtree] + # Skip if target must be topmost but would not be. + if not (target_in_top and not target_in_bottom): + boundary_b = cost.get((bottom_bottom, top_top), 0) + results.append(( + top_cost + bottom_cost + boundary_b, + bottom_top, + top_bottom, + bottom_order + top_order, + )) + return results + + +def optimize_topology( + newick_path: str, + cost: dict[tuple[str, str], int], + target_genome: str | None = None, +) -> tuple[list[str], int, Tree]: + """ + Load the tree, run exhaustive memoized optimization, and return: + - optimal genome ordering (list of names, top to bottom) + - optimal total inversion cost + - Biopython Tree object rotated to match the optimal ordering + + If target_genome is specified, it is guaranteed to be first in the ordering. + The tree is midpoint-rooted if it is not already strictly bifurcating + at the root (e.g. unrooted input). + """ + tree = load_tree(newick_path) + root = tree.root + + if target_genome is not None: + leaf_names = set(_get_leaf_names(root)) + if target_genome not in leaf_names: + raise ValueError( + f"--target-genome '{target_genome}' not found in tree leaves: " + f"{sorted(leaf_names)}" + ) + + options = build_subtree_options(root, cost, target_genome) + for option in options: + print(option, file=sys.stderr) + best_cost, _, _, best_order_tuple = min(options, key=lambda x: int(x[0])) + best_order = list(best_order_tuple) + + _rotate_tree_to_order(root, best_order) + + return best_order, best_cost, tree + + +def _rotate_tree_to_order(node: Clade, target_order: list[str]) -> None: + """ + Rotate subtrees of `node` in-place so that the leaf order matches + target_order (top-to-bottom). + """ + if _is_leaf(node): + return + + children = _get_children(node) + top_leaves = set(_get_leaf_names(children[0])) + + # Find where the top child's leaves sit in target_order + top_positions = [i for i, g in enumerate(target_order) if g in top_leaves] + bottom_positions = [i for i, g in enumerate(target_order) if g not in top_leaves] + + # If the top child's leaves all come after the bottom child's, swap + if top_positions and bottom_positions and min(top_positions) > min(bottom_positions): + node.clades[0], node.clades[1] = node.clades[1], node.clades[0] + children = _get_children(node) + + for child in children: + _rotate_tree_to_order(child, target_order) + + +# --------------------------------------------------------------------------- +# Step 3: Output +# --------------------------------------------------------------------------- + +def write_order_tsv(order: list[str], out_path: str) -> None: + "Write the optimized genome order to a TSV file, one genome per line." + with open(out_path, "w", encoding="utf-8") as fh: + for genome in order: + fh.write(genome + "\n") + + +def write_newick(tree: Tree, out_path: str) -> None: + "Write the rotated tree to a Newick file." + Phylo.write(tree, out_path, "newick") + + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + +def main() -> None: + "Parse arguments, build cost matrix, optimize tree topology, and write outputs." + args = parse_args() + + # --- Build cost matrix --- + print("Building pairwise inversion cost matrix...", file=sys.stderr) + cost = build_cost_matrix(args.blocks) + genomes = get_genomes_from_cost_matrix(cost) + print(f" Found {len(genomes)} genomes.", file=sys.stderr) + + if not cost: + print( + "Warning: no pairwise inversions found. All costs are zero; " + "the original tree ordering will be preserved.", + file=sys.stderr, + ) + + # --- Optimize --- + print("Optimizing tree topology...", file=sys.stderr) + if args.target: + print(f" Target genome (always first): {args.target}", file=sys.stderr) + best_order, best_cost, rotated_tree = optimize_topology( + args.tree, cost, target_genome=args.target + ) + + print(f" Optimal total inversion cost: {best_cost:,} bp", file=sys.stderr) + print(" Optimal ordering:", file=sys.stderr) + for i, genome in enumerate(best_order): + print(f" {i + 1}. {genome}", file=sys.stderr) + + # --- Write outputs --- + write_order_tsv(best_order, args.out_order) + print(f"Order written to: {args.out_order}", file=sys.stderr) + + write_newick(rotated_tree, args.out_tree) + print(f"Rotated tree written to: {args.out_tree}", file=sys.stderr) + + +if __name__ == "__main__": + main() diff --git a/bin/ntsynt_viz_plot_synteny_blocks_ribbon_plot.R b/bin/ntsynt_viz_plot_synteny_blocks_ribbon_plot.R index 63210dd..cffe76e 100755 --- a/bin/ntsynt_viz_plot_synteny_blocks_ribbon_plot.R +++ b/bin/ntsynt_viz_plot_synteny_blocks_ribbon_plot.R @@ -11,7 +11,8 @@ suppressPackageStartupMessages({ library(stringr) library(gggenomes) library(tidyr) - library(svglite)}) + library(svglite) + library(ggiraph)}) # Example script for generating ntSynt synteny ribbon plots using gggenomes @@ -42,7 +43,7 @@ parser$add_argument("--right-ratio", help = paste("Ratio adjustment for space on the right side of the ribbon plot.", "Increase if the labels on the right are cut-off,", "decrease to decrease space between ribbon plot and right edge of the plot"), - default = 0.03, required = FALSE, type = "double") + default = 0.07, required = FALSE, type = "double") parser$add_argument("-p", "--prefix", help = "Output prefix for PNG image (default synteny_gggenomes_plot)", required = FALSE, default = "synteny_gggenomes_plot") @@ -55,6 +56,10 @@ parser$add_argument("--annotate-genome-info", help = "Add annotations about numb args <- parser$parse_args() +############################# +# Prepare data +############################# + # Read in and prepare sequences sequences <- read.csv(args$sequences, sep = "\t", header = TRUE) %>% mutate(relative_orientation = if_else(relative_orientation == "+", "", "\u2190")) @@ -79,6 +84,8 @@ links_ntsynt <- read.csv(args$links, sep = "\t", header = TRUE) %>% mutate(bin_id = str_replace_all(bin_id, "_", " "), bin_id2 = str_replace_all(bin_id2, "_", " ")) +target_genome <- links_ntsynt %>% head(1) %>% select(bin_id) %>% pull() +print(paste("Target genome:", target_genome)) links_ntsynt$seq_id <- factor(links_ntsynt$seq_id, levels = input_chrom_order) links_ntsynt <- links_ntsynt %>% arrange(factor(seq_id, levels = input_chrom_order)) @@ -109,6 +116,23 @@ painting <- read.csv(args$painting, sep = "\t", header = TRUE) %>% # Read in the data frame with info about colours to choose for sequences colours_df <- read.csv(args$colour_indices, sep = "\t", header = TRUE) +# Read in haplotypes, or set to FALSE +if (! is.null(args$haplotypes)) { + haplotypes <- read.csv(args$haplotypes, sep = "\t", header = TRUE) %>% mutate(bin_id = str_replace_all(bin_id, "_", " ")) +} else { + haplotypes <- FALSE +} + +if (! is.null(args$centromeres)) { + centromeres <- read.csv(args$centromeres, sep = "\t", header = TRUE) +} else { + centromeres <- FALSE +} + +############################# +# Helper functions +############################# + # Get the y coordinates for the features get_y_coord <- function(haplotypes, bin_id, y, end=FALSE) { if (typeof(haplotypes) == "logical") { @@ -143,21 +167,92 @@ format_genome_size <- function(bp) { # Return dataframe with bin annotations get_bin_annotations <- function(plot){ bin_stats <- plot %>% - pull_seqs() %>% # or seqs(gggenomes_obj) + pull_seqs() %>% group_by(bin_id) %>% summarise( n_chr = n(), genome_size = sum(length), - y = max(y), # matches gggenomes' y layout - x_right = max(xend) # rightmost coordinate + y = max(y), + x_right = max(xend) ) %>% - mutate(label = paste0(" ", format_genome_size(genome_size), "; n=", n_chr)) + mutate(label = paste0("\u00A0\u00A0\u00A0", format_genome_size(genome_size), "; n=", n_chr)) bin_stats <- bin_stats %>% mutate(x_right = max(bin_stats$x_right)) return(bin_stats) } -# Make the ribbon plot - these layers can be fully customized as needed! +# Prepare a summary string for each synteny block +get_block_coord_info <- function(p, max_genome_len, max_chrom_len) { + block_coords <- pull_links(p) %>% + group_by(block_id) %>% + summarise( + coords = { + # Each row contributes two genomes; collect all unique combinations + g1 <- tibble(genome = bin_id, chrom = seq_id, start = start, end = end, strand = strand1) + g2 <- tibble(genome = bin_id2, chrom = seq_id2, start = start2, end = end2, strand = strand2) + bind_rows(g1, g2) %>% + distinct() %>% + mutate(line = paste0(stringr::str_pad(genome, width = max_genome_len, side="right", pad='\u00A0'), + " ", + stringr::str_pad(chrom, width = max_chrom_len, side="right", pad='\u00A0'), + ": ", + format(start, big.mark=","), " – ", + format(end, big.mark=","), " bp (", strand, ")")) %>% + pull(line) %>% + paste(collapse = "\n") + }, + .groups = "drop" + ) + return(block_coords) +} + +# Get the information/data about synteny links for interactive layers with ggiraph +get_link_info <- function(p, block_coords) { + link_data <- pull_links(p) %>% + left_join(block_coords, by = "block_id") %>% + mutate( + tooltip = paste0("Block ID: ", block_id, "\n", coords), + group_id = row_number() + ) %>% + rowwise() %>% + mutate(poly = list(tibble( + px = c(x, xend, xmax, xmin), + py = c(y, y, yend, yend), + tooltip = tooltip, + data_id = block_id, + colour_block = colour_block + ))) %>% + ungroup() %>% + select(group_id, poly) %>% + unnest(poly) + return(link_data) +} + +# Build chromosome -> block_id mapping (target genome seq_id only, as these appear in legend) +build_js_map <- function(p, target_genome) { + chrom_block_map <- pull_links(p) %>% + filter(bin_id == target_genome) %>% + select(block_id, seq_id) %>% + distinct() %>% + group_by(seq_id) %>% + summarise(block_ids = list(unique(block_id)), .groups = "drop") %>% + rename(chrom = seq_id) + + js_map_entries <- chrom_block_map %>% + rowwise() %>% + mutate(entry = paste0( + '"', chrom, '": [', + paste0('"', unlist(block_ids), '"', collapse = ","), + ']' + )) %>% + pull(entry) %>% + paste(collapse = ",\n") + + js_map <- paste0("const chromBlockMap = {\n", js_map_entries, "\n};") + return(js_map) +} + +# Make the ribbon plot - these layers can be fully customized as needed make_plot <- function(links, sequences, painting, colours_df, add_scale_bar = FALSE, centromeres = FALSE, add_arrow = FALSE, haplotypes = FALSE) { target_genome <- (sequences %>% head(1) %>% select(bin_id))[[1]] sequences_filt <- unique((sequences %>% filter(bin_id == target_genome))$seq_id) @@ -184,7 +279,6 @@ make_plot <- function(links, sequences, painting, colours_df, add_scale_bar = FA geom_bin_label(aes(label = bin_id, y = get_y_coord(haplotypes, bin_id, .data$y)), size = 6, fontface = "italic") + # label each bin - #geom_seq_label(aes(label = seq_id), vjust = 1.1, size = 4) + # Can add seq labels if desired theme(axis.text = element_text(size = 25, face = "italic"), legend.position = "bottom", legend.text = element_text(size = 15, margin = margin(r=10, l=3))) + @@ -226,26 +320,78 @@ make_plot <- function(links, sequences, painting, colours_df, add_scale_bar = FA expand_limits(x = max(bin_stats$x_right) + (xmax * (args$right_ratio))) } - return(plot) - -} - -# Read in haplotypes, or set to FALSE -if (! is.null(args$haplotypes)) { - haplotypes <- read.csv(args$haplotypes, sep = "\t", header = TRUE) %>% mutate(bin_id = str_replace_all(bin_id, "_", " ")) -} else { - haplotypes <- FALSE + # Prepare interactive components + # Combine all unique genome/chrom combinations across the dataset + all_elements <- bind_rows( + pull_links(p) %>% select(genome = bin_id, chrom = seq_id, start, end), + pull_links(p) %>% select(genome = bin_id2, chrom = seq_id2, start = start2, end = end2) + ) %>% distinct() + + # Find the maximum string length of the genome IDs to know how much to pad for hover boxes + max_genome_len <- max(nchar(all_elements$genome), na.rm = TRUE) + max_chrom_len <- max(nchar(all_elements$chrom), na.rm = TRUE) + + block_coords <- get_block_coord_info(p, max_genome_len, max_chrom_len) + link_data <- get_link_info(p, block_coords) + + # Add ribbon hover interactivity + plot <- plot + + geom_polygon_interactive( + data = link_data, + aes( + x = px, + y = py, + group = group_id, + tooltip = tooltip, + data_id = data_id, + fill = colour_block + ), + alpha = 0, + hover_nearest = FALSE + ) + + # Prepare information about chromosomes for interactive layers + seq_data <- pull_seqs(p) %>% + mutate( + length_fmt = format(length, big.mark = ",", scientific = FALSE), + tooltip = paste0( + "Genome: ", bin_id, "\n", + "Chromosome: ", seq_id, "\n", + "Length: ", length_fmt, " bp" + ) + ) + # Add interactive layer for chromosomes + plot <- plot + + geom_segment_interactive( + data = seq_data, + aes( + x = x, + xend = xend, + y = y, + yend = y, + tooltip = tooltip, + data_id = seq_id + ), + linewidth = 3, # Increased area to hit for hover box + alpha = 0, + hover_nearest = FALSE + ) + + # Build map of chromosome -> block_id for interactive legend + js_map <- build_js_map(p, target_genome) + + return(list(plot = plot, js_map = js_map)) } -if (! is.null(args$centromeres)) { - centromeres <- read.csv(args$centromeres, sep = "\t", header = TRUE) -} else { - centromeres <- FALSE -} +############################# +# Prepare plots +############################# # Make the ribbon plot -synteny_plot <- make_plot(links_ntsynt, sequences, painting, colours_df, add_scale_bar = TRUE, centromeres = centromeres, - add_arrow = !args$no_arrow, haplotypes = haplotypes) +synteny_plot_tmp <- make_plot(links_ntsynt, sequences, painting, colours_df, add_scale_bar = TRUE, centromeres = centromeres, + add_arrow = !args$no_arrow, haplotypes = haplotypes) +synteny_plot <- synteny_plot_tmp$plot +js_map <- synteny_plot_tmp$js_map if (is.null(args$tree)) { @@ -301,19 +447,77 @@ if (any_rc && !args$no_arrow) { plots <- ggarrange(plots, note, ncol = 1, heights = c(10, 1)) } -# Save the ribbon plot +# Save static plot in requested format if (args$format == "pdf") { - ggsave(paste(args$prefix, ".pdf", sep = ""), plots, + ggsave(paste0(args$prefix, ".pdf"), plots, units = "cm", width = args$width, height = args$height, bg = "white") - cat(paste("Plot saved:", paste(args$prefix, ".pdf", sep = ""), "\n", sep = " ")) + cat(paste("Plot saved:", paste0(args$prefix, ".pdf"), "\n")) } else if (args$format == "svg") { - ggsave(paste(args$prefix, ".svg", sep = ""), plots, + ggsave(paste0(args$prefix, ".svg"), plots, units = "cm", width = args$width, height = args$height, bg = "white") - cat(paste("Plot saved:", paste(args$prefix, ".svg", sep = ""), "\n", sep = " ")) -} else { - png(paste(args$prefix, ".png", sep = ""), units = "cm", width = args$width, height = args$height, - res = args$dpi, bg = "white") - print(plots) - dev.off() - cat(paste("Plot saved:", paste(args$prefix, ".png", sep = ""), "\n", sep = " ")) -} + cat(paste("Plot saved:", paste0(args$prefix, ".svg"), "\n")) +} +png(paste0(args$prefix, ".png"), units = "cm", width = args$width, height = args$height, + res = args$dpi, bg = "white") +print(plots) +dev.off() +cat(paste("Plot saved:", paste0(args$prefix, ".png"), "\n")) + + + +# Prepare interactive HTML +script_dir <- dirname(normalizePath(commandArgs(trailingOnly = FALSE)[ + grep("--file=", commandArgs(trailingOnly = FALSE)) +] |> sub("--file=", "", x = _))) +js_template <- paste( + readLines(paste(script_dir, "/ntsynt_viz_ribbon-interactive.js", sep=""), warn = FALSE), + collapse = "\n" +) + +js_inject <- gsub( + "__CHROM_BLOCK_MAP__", + js_map, + js_template, + fixed = TRUE +) + +interactive_plot <- girafe( + ggobj = plots, + width_svg = args$width / 2.54, # Converting to inches + height_svg = args$height / 2.54, + options = list( + opts_hover(css = "stroke:darkgrey; stroke-width:1; fill-opacity:0.1; transition: all 0.1s ease;"), + opts_hover_inv(css = "opacity:0.2;"), + opts_zoom(max = 10), + opts_toolbar(pngname = args$prefix), + opts_sizing(rescale = TRUE, width = 1), + opts_tooltip( + css = paste( + "background: rgba(255,255,255,0.9);", + "padding: 10px;", + "border: 1px solid black;", + "border-radius: 4px;", + "font-family: monospace;", + "font-size: 16px;" + ) + ) + ) +) + +# Inject CSS to ensure the interactive plot fills a browser window +css_override <- paste( + "", + sep = "\n" +) + +html_file <- paste0(args$prefix, ".html") +htmlwidgets::saveWidget(interactive_plot, html_file, selfcontained = TRUE, title = args$prefix) +html_content <- readLines(html_file, warn = FALSE) +head_close <- which(grepl("", html_content)) +html_content <- append(html_content, css_override, after = head_close - 1) +body_close <- which(grepl("", html_content)) +html_content <- append(html_content, js_inject, after = body_close - 1) +writeLines(html_content, html_file) +cat(paste("Interactive HTML saved:", html_file, "\n")) \ No newline at end of file diff --git a/bin/ntsynt_viz_ribbon-interactive.js b/bin/ntsynt_viz_ribbon-interactive.js new file mode 100644 index 0000000..77d740d --- /dev/null +++ b/bin/ntsynt_viz_ribbon-interactive.js @@ -0,0 +1,433 @@ + \ No newline at end of file diff --git a/bin/sort_ntsynt_blocks.py b/bin/sort_ntsynt_blocks.py index 3b1a5b2..a3c015a 100755 --- a/bin/sort_ntsynt_blocks.py +++ b/bin/sort_ntsynt_blocks.py @@ -13,6 +13,31 @@ faidx_re = re.compile(r'^(\S+)\.fai$') +def normalize_strands(sorted_blocks): + "Re-orient strands so they are relative to the top (first) genome in sorted_blocks" + if not sorted_blocks: + return sorted_blocks + + top_strand = sorted_blocks[0].strand + + normalized_blocks = [] + for block in sorted_blocks: + if top_strand == "+": + new_strand = block.strand + else: + new_strand = "+" if block.strand == "-" else "-" + normalized_blocks.append(block._replace(strand=new_strand)) + + return normalized_blocks + +def print_sorted_block(synteny_block_lines, sorting_order): + "Sort a single synteny block's lines by the given order, normalize strands, and print" + sorted_blocks = sorted(synteny_block_lines, key=lambda x: sorting_order[x.genome]) + sorted_blocks = normalize_strands(sorted_blocks) + for asm_block in sorted_blocks: + print(*asm_block, sep="\t") + + def sort_blocks(synteny_blocks_fin, sorting_order): "Sort the assemblies within each synteny block" @@ -28,15 +53,13 @@ def sort_blocks(synteny_blocks_fin, sorting_order): block = SyntenyBlockComp(*line[:6]) if curr_block_id is not None and block.id != curr_block_id: # Print out the sorted blocks, they are all read in - for asm_block in sorted(synteny_block_lines, key=lambda x: sorting_order[x.genome]): - print(*asm_block, sep="\t") + print_sorted_block(synteny_block_lines, sorting_order) synteny_block_lines = [] synteny_block_lines.append(block) curr_block_id = block.id # Print out the sorted blocks, they are all read in (making sure we have last block) - for asm_block in sorted(synteny_block_lines, key=lambda x: sorting_order[x.genome]): - print(*asm_block, sep="\t") + print_sorted_block(synteny_block_lines, sorting_order) def main(): diff --git a/tests/expected/great-apes_ribbon-plots_ribbon-plot_tree.html b/tests/expected/great-apes_ribbon-plots_ribbon-plot_tree.html new file mode 100644 index 0000000..52cb155 --- /dev/null +++ b/tests/expected/great-apes_ribbon-plots_ribbon-plot_tree.html @@ -0,0 +1,1675 @@ + + + + +great-apes_ribbon-plots_ribbon-plot_tree + + + + + +

+ +

+ + + +
+
+
+ + + + +