Add option for genome information annotations (#50)

* Add option --annotate-genomes to annotate genome summary info

- to the right of the genomes in the ribbon plot
- Genome size + number of chromosomes

* Update new parameter name to --annotate-genome-info

* Fix for pylint

* Refactor parameter section of ntsynt_viz.smk

* Update README.md

* Update README.md

* Rename genome info annotation parameter in smk, R script

- Also update .pylintrc

* Update passing of name_conversion in ntsynt_viz.py

* Remove name conversion from a test in CI

* Update Example 2 ubuntu orders

* Update CI

* Update expected orders TSV

Author: lcoombe

README.md

@@ -66,8 +66,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] [--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]
+                     [--no-arrow] [--ribbon_adjust RIBBON_ADJUST] [-f] [-n] [-v]
 
 Visualizing multi-genome synteny
 
@@ -114,8 +114,10 @@ output arguments:
                         Output format of plot [png]
   --scale SCALE         Length of scale bar in bases [100e6]
   --height HEIGHT       Height of plot in cm [20]
-  --width WIDTH         Width of plot in cm [50]
+  --width WIDTH         Width of plot in cm [55]
   --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
   --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]

azure-pipelines.yml

@@ -53,15 +53,17 @@ 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  --name_conversion great-apes.name-conversions.tsv  --prefix great-apes_ribbon-plots_no-tree --ribbon_adjust 0.15 --scale 1e9 --target-genome Homo_sapiens --format pdf --force
-      diff -q great-apes_ribbon-plots_no-tree_est-distances.order.tsv expected/great-apes_ribbon-plots_no-tree_est-distances.order.tsv
+      ntsynt_viz.py --blocks great-apes.ntSynt.synteny_blocks.tsv --fais fais.tsv --prefix great-apes_ribbon-plots_no-tree --ribbon_adjust 0.15 --scale 1e9 --target-genome GCA_000001405.15_GRCh38_genomic.chr-only.nums.fa  --format pdf --force
+      cat great-apes_ribbon-plots_no-tree_est-distances.order.tsv
+      diff -q great-apes_ribbon-plots_no-tree_est-distances.order.tsv expected/great-apes_ribbon-plots_no-tree_est-distances_no-names.order.tsv
     displayName: Run Example 2 test (PDF output)
 
   - script: |
       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  --name_conversion great-apes.name-conversions.tsv  --prefix great-apes_ribbon-plots_no-tree --ribbon_adjust 0.15 --scale 2e9 --target-genome Pan_paniscus --format svg --force
+      ntsynt_viz.py --blocks great-apes.ntSynt.synteny_blocks.tsv --fais fais.tsv  --name_conversion great-apes.name-conversions.tsv  --prefix great-apes_ribbon-plots_no-tree \
+      --ribbon_adjust 0.15 --scale 2e9 --target-genome Pan_paniscus --format svg --annotate-genome-info --force
     displayName: Run Example 3 test (SVG output)
 
   - script: |

bin/.pylintrc

@@ -6,7 +6,7 @@ max-attributes = 10
 max-line-length = 120
 max-branches = 20 
 max-public-methods=50
-max-statements = 70 
+max-statements = 90 
 max-module-lines = 1500
 max-bool-expr = 7
 max-positional-arguments = 10

bin/ntsynt_viz.py

@@ -93,9 +93,13 @@ def main():
     output_group.add_argument("--scale", help="Length of scale bar in bases [100e6]", required=False, type=float,
                         default=100e6)
     output_group.add_argument("--height", help="Height of plot in cm [20]", required=False, type=int, default=20)
-    output_group.add_argument("--width", help="Width of plot in cm [50]", required=False, type=int, default=50)
+    output_group.add_argument("--width", help="Width of plot in cm [55]", required=False, type=int, default=55)
     output_group.add_argument("--dpi", help="Resolution of output plot - png output only [300]",
                               required=False, type=int, default=300)
+    output_group.add_argument("--annotate-genome-info",
+                              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; "
                         "do not draw arrows indicating reverse-complementation",
                         action="store_true")
@@ -128,7 +132,6 @@ def main():
             f"--config " \
             f"prefix={args.prefix} " \
             f"blocks={args.blocks} " \
-            f"name_conversion={args.name_conversion} " \
             f"fai='{args.fais}' " \
             f"ribbon_ratio={args.ribbon_adjust} " \
             f"scale={args.scale} " \
@@ -139,6 +142,8 @@ def main():
             f"width={args.width} " \
             f"min_seq_length={args.seq_length} "
 
+    if args.name_conversion:
+        cmd += f"name_conversion={args.name_conversion} "
     if args.centromeres:
         cmd += f"centromeres={args.centromeres} "
     if args.normalize:
@@ -155,6 +160,8 @@ def main():
         cmd += f"dpi={args.dpi} "
     if args.order:
         cmd += f"order={args.order} "
+    if args.annotate_genome_info:
+        cmd += "annotate_genome_info=True "
     if args.tree:
         cmd += f"tree={args.tree} "
         target = "gggenomes_ribbon_plot_tree"

bin/ntsynt_viz.smk

@@ -5,30 +5,30 @@ import os
 
 
 # Parameters
-synteny_blocks = config["blocks"] if "blocks" in config else None
-name_conversion = config["name_conversion"] if "name_conversion" in config else None
-name_conversion = None if name_conversion == "None" else name_conversion
-prefix = config["prefix"] if "prefix" in config else "ntSynt_gggenomes"
-indel_threshold = config["indel_threshold"] if "indel_threshold" in config else 50000
-min_len = config["min_length"] if "min_length" in config else 50000
-fais = config["fai"] if "fai" in config else None
-ribbon_ratio = config["ribbon_ratio"] if "ribbon_ratio" in config else 0.1
-cladogram_ratio = config["cladogram_adjust"] if "cladogram_adjust" in config else 0.1
-scale = config["scale"] if "scale" in config else 1e9
-centromeres = config["centromeres"] if "centromeres" in config else None
-normalize = config["normalize"] if "normalize" in config else False
+synteny_blocks = config.get("blocks", None)
+name_conversion = config.get("name_conversion", None)
+prefix = config.get("prefix", "ntSynt_gggenomes")
+indel_threshold = config.get("indel_threshold", 50000)
+min_len = config.get("min_length", 50000)
+fais = config.get("fai", None)
+ribbon_ratio = config.get("ribbon_ratio", 0.1)
+cladogram_ratio = config.get("cladogram_adjust", 0.1)
+scale = config.get("scale", 1e9)
+centromeres = config.get("centromeres", None)
+normalize = config.get("normalize", False)
 blocks_no_suffix = os.path.basename(synteny_blocks).removesuffix(".tsv")
-tree = config["tree"] if "tree" in config else None
-format_img = config["format"] if "format" in config else "png"
-plot_height = config["height"] if "height" in config else 20
-plot_width = config["width"] if "width" in config else 50
-no_arrow = config["no_arrow"] if "no_arrow" in config else False
-target_genome = config["target_genome"] if "target_genome" in config else False
-haplotypes = config["haplotypes"] if "haplotypes" in config else []
-keep = " ".join(config["keep"]) if "keep" in config else []
-min_seq_length = config["min_seq_length"] if "min_seq_length" in config else 500000
-res = config["dpi"] if "dpi" in config else 300
+tree = config.get("tree", None)
+format_img = config.get("format", "png")
+plot_height = config.get("height", 20)
+plot_width = config.get("width", 50)
+no_arrow = config.get("no_arrow", False)
+target_genome = config.get("target_genome", False)
+haplotypes = config.get("haplotypes", [])
+keep = " ".join(config.get("keep", []))
+min_seq_length = config.get("min_seq_length", 500000)
+res = config.get("dpi", 300)
 orders = config.get("order", [])
+annotate_genome_info = config.get("annotate_genome_info", False)
 
 def sort_fais(fai_list, name_conversion, orders):
     "Based on the name conversion TSV, sort the FAIs based on orders"
@@ -271,10 +271,11 @@ rule ribbon_plot:
         arrow = "--no-arrow" if no_arrow else "",
         haplotypes = f"--haplotypes {rules.nudges.output.nudges}" if haplotypes else "",
         resolution = f"--dpi {res}" if format_img == "png" else "",
+        annotate_genome_info = "--annotate-genome-info" if annotate_genome_info else "",
     shell:
         "ntsynt_viz_plot_synteny_blocks_ribbon_plot.R -s {input.sequences} -l {input.links} -p {params.prefix} --ratio {params.ratio}" 
         " --scale {params.scale} -c {input.colour_feats} --format {params.out_img_format} --height {params.height} --width {params.width}"
-        " {params.centromeres} {params.arrow} {params.haplotypes} --colour_indices {input.colour_seqs} {params.resolution}"
+        " {params.centromeres} {params.arrow} {params.haplotypes} --colour_indices {input.colour_seqs} {params.resolution} {params.annotate_genome_info}"
 
 rule ribbon_plot_tree:
     input: 
@@ -297,7 +298,8 @@ rule ribbon_plot_tree:
         arrow = "--no-arrow" if no_arrow else "",
         haplotypes = f"--haplotypes {rules.nudges.output.nudges}" if haplotypes else "",
         resolution = f"--dpi {res}" if format_img == "png" else "",
+        annotate_genome_info = "--annotate-genome-info" if annotate_genome_info else "",
     shell:
         "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}"
+        " --order {input.orders} {params.centromeres} {params.arrow} {params.haplotypes} --colour_indices {input.colour_seqs} {params.resolution} {params.annotate_genome_info}"

bin/ntsynt_viz_plot_synteny_blocks_ribbon_plot.R

@@ -38,13 +38,20 @@ parser$add_argument("--ratio",
                                  "Increase if the labels are cut-off,",
                                  "decrease to decrease space between ribbon plot and cladogram"),
                     default = 0.1, required = FALSE, type = "double")
+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")
 parser$add_argument("-p", "--prefix",
                     help = "Output prefix for PNG image (default synteny_gggenomes_plot)", required = FALSE,
                     default = "synteny_gggenomes_plot")
 parser$add_argument("--format", help = "Output format for image (png, pdf or svg)", required = FALSE,
                     default = "png", choices = c("png", "pdf", "svg"))
 parser$add_argument("--order", help = "TSV file with desired order of tip labels (only used if --tree specified).", required = FALSE)
 parser$add_argument("--dpi", help = "Output plot resolution - for PNG only (default 300)", default = 300, required = FALSE, type = "integer")
+parser$add_argument("--annotate-genome-info", help = "Add annotations about number of sequences and genome size to the right of each bin",
+                    action = "store_true")
 
 args <- parser$parse_args()
 
@@ -124,6 +131,32 @@ get_y_coord <- function(haplotypes, bin_id, y, end=FALSE) {
   }
 }
 
+format_genome_size <- function(bp) {
+  dplyr::case_when(
+    bp >= 1e9 ~ paste0(sprintf("%.1f", bp / 1e9), " Gbp"),
+    bp >= 1e6 ~ paste0(sprintf("%.1f", bp / 1e6), " Mbp"),
+    bp >= 1e3 ~ paste0(sprintf("%.1f", bp / 1e3), " kbp"),
+    TRUE      ~ paste0(bp, " bp")
+  )
+}
+
+# Return dataframe with bin annotations
+get_bin_annotations <- function(plot){
+  bin_stats <- plot %>%
+  pull_seqs() %>%                          # or seqs(gggenomes_obj)
+  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
+  ) %>%
+  mutate(label = paste0(" ", 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!
 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]]
@@ -185,6 +218,14 @@ make_plot <- function(links, sequences, painting, colours_df, add_scale_bar = FA
             axis.ticks.x = element_blank())
   }
 
+  if (args$annotate_genome_info) {
+    bin_stats <- get_bin_annotations(plot)
+    
+    plot <- plot + geom_text(data = bin_stats, aes(x = x_right, y = y, label = label, hjust = 0),
+                             size = 4) +
+                  expand_limits(x = max(bin_stats$x_right) + (xmax * (args$right_ratio)))
+  }
+
   return(plot)
 
 }

tests/expected/great-apes_ribbon-plots_no-tree_est-distances_no-names.order.tsv

@@ -0,0 +1,6 @@
+GCA_000001405.15_GRCh38_genomic.chr-only.nums.fa
+GCF_002880755.1_Clint_PTRv2_genomic.chr-only.rename.fa
+GCF_013052645.1_Mhudiblu_PPA_v0_genomic.chr-only.rename.fa
+GCF_008122165.1_Kamilah_GGO_v0_genomic.chr-only.rename.fa
+GCF_028885625.2_NHGRI_mPonPyg2-v2.0_pri_genomic.chr.renamed.fa
+GCA_021498465.1_ASM2149846v1_genomic.renamed.chr.fa