## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(eval = FALSE, collapse = TRUE, comment = "#>") ## ----example-1, eval=FALSE---------------------------------------------------- # install.packages(c("remotes", "dplyr", "purrr", "tidyr")) # remotes::install_github("DangJem/ggchord@v0.16.6") # library(ggchord) # packageVersion("ggchord") ## ----example-2, eval=FALSE---------------------------------------------------- # install.packages("pkgload") # pkgload::load_all(".") # packageVersion("ggchord") ## ----example-3, eval=FALSE---------------------------------------------------- # library(ggchord) # library(dplyr) # data(seq_data_example) # data(gene_data_example) # data(ribbon_data_example) # ids <- c("MT108731.1", "MT118296.1") # first_sequences <- seq_data_example |> # filter(accver %in% ids) |> # select(accver, length) # first_genes <- gene_data_example |> # filter(accver %in% ids) |> # transmute(accver, start, end, strand = source_strand) # first_links <- ribbon_data_example |> # filter(qaccver %in% ids, saccver %in% ids) |> # select(qaccver, saccver, qstart, qend, sstart, send) # # first_colours <- setNames(c("#E64B35", "#233E4A"), ids) # p_first <- ggchord( # seq_data = first_sequences, ribbon_data = first_links, # gene_data = first_genes # ) + # geom_link_ribbon(fill = "#82B6AD", alpha = .55) + # geom_seq(colour = "#233E4A", show.legend = FALSE) + # geom_gene(aes(gene_fill = accver), position = position_strand(.10), width = .045) + # geom_seq_label(size = 4) + # coord_collinear() + # scale_gene_fill_manual( # values = first_colours, name = "Sequence", # guide = guide_ggchord_legend(nrow = 1, theme = theme( # legend.title = element_text(size = 13, face = "bold"), legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), legend.margin = margin(l = 7, unit = "mm") # )) # ) + # scale_seq_position_continuous(n.breaks = 3, minor_breaks = NULL) + # theme_ggchord_publication(axis.text = element_text(size = 11)) + # theme( # legend.position = "bottom", # plot.background = element_rect(fill = "#FFFCF5", colour = NA), # panel.background = element_blank() # ) # p_first ## ----example-4, eval=FALSE---------------------------------------------------- # ggsave("first-comparison.png", p_first, width = 10, height = 5, units = "in", dpi = 144) ## ----example-5, eval=FALSE---------------------------------------------------- # library(ggchord) # library(dplyr) # data(seq_data_example) # data(ribbon_data_example) # data(gene_data_example) # seqs <- seq_data_example |> # select(accver, length) # links <- ribbon_data_example |> # select(qaccver, saccver, qstart, qend, sstart, send) # # The fixture preserves biological strands separately from its old display strands. # genes <- gene_data_example |> # transmute(accver, start, end, strand = source_strand) # # # Optional appearance choices; none of these columns or scales is required input. # library(ggchord) # seq_palette <- setNames( # c("#E64B35", "#233E4A", "#269B82", "#D59B17"), # c("MT108731.1", "MT118296.1", "OQ646790.1", "OR222515.1") # ) # gene_palette <- c("+" = "#CE1536", "-" = "#346B9D", "." = "#8291A6") # # library(ggchord) # feature_palette <- c( # CDS = "#993366", promoter = "#31849B", rep_origin = "#F1E51D", # primer_bind = "#A020F0", misc_feature = "#A6ACB3", protein_bind = "#31849B" # ) # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } # # library(dplyr) # data(plasmid_example_pUC19) # plasmid_seq <- plasmid_example_pUC19 |> # select(accver, length) # feature_source <- find_common_features(plasmid_example_pUC19) # features <- feature_source |> # select(accver, start, end, strand, type, anno) ## ----example-6, eval=FALSE---------------------------------------------------- # p_chord <- ggchord(seq_data = seqs, ribbon_data = links, gene_data = genes) + # geom_link_ribbon(aes(ribbon_fill = qaccver), alpha = .7) + # scale_ribbon_fill_manual(values = seq_palette, guide = "none") + # geom_seq(show.legend = FALSE) + # geom_gene(position = position_strand(.1), width = .04) + # geom_seq_label(size = 4, seq_label_radius = .84, seq_label_hjust = .5) + # coord_chord() + # scale_seq_colour_manual(values = seq_palette) + # scale_gene_fill_manual(values = gene_palette) + # doc_theme # p_chord ## ----example-7, eval=FALSE---------------------------------------------------- # anchors <- genes |> # group_by(accver) |> # slice_max(end - start, n = 1, with_ties = FALSE) |> # ungroup() |> # select(accver, start, end) # p_collinear <- ggchord(seq_data = seqs, ribbon_data = links, gene_data = genes) + # geom_link_ribbon(aes(ribbon_fill = qaccver), alpha = .7) + # scale_ribbon_fill_manual(values = seq_palette, guide = "none") + # geom_seq() + # geom_gene(position = position_strand(.1)) + # geom_seq_label(size = 4) + # coord_collinear(anchor = anchors, align = "center", orientation = "auto") + # doc_axis() + # scale_seq_colour_manual(values = seq_palette) + # scale_gene_fill_manual(values = gene_palette) + # doc_theme # p_collinear ## ----example-8, eval=FALSE---------------------------------------------------- # circle_features <- features |> # filter(type %in% c("CDS", "rep_origin")) # p_circular <- ggchord(seq_data = plasmid_seq) + # geom_seq() + # geom_feature(aes(feature_fill = type), data = circle_features, position = position_plasmid()) + # geom_feature_label_repel(data = circle_features, position = position_plasmid(), size = 4) + # geom_seq_center_label() + # scale_feature_fill_manual(values = feature_palette, breaks = unique(circle_features$type)) + # coord_circular() + # doc_theme # p_circular ## ----example-9, eval=FALSE---------------------------------------------------- # library(ggchord) # library(dplyr) # seq_data_example |> select(accver, length) |> slice_head(n = 3) # ribbon_data_example |> select(qaccver, saccver, qstart, qend, sstart, send) |> slice_head(n = 3) # gene_data_example |> transmute(accver, start, end, strand = source_strand, anno) |> slice_head(n = 3) ## ----example-10, eval=FALSE--------------------------------------------------- # # Optional appearance choices; none of these columns or scales is required input. # library(ggchord) # seq_palette <- setNames( # c("#E64B35", "#233E4A", "#269B82", "#D59B17"), # c("MT108731.1", "MT118296.1", "OQ646790.1", "OR222515.1") # ) # gene_palette <- c("+" = "#CE1536", "-" = "#346B9D", "." = "#8291A6") # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } ## ----example-11, eval=FALSE--------------------------------------------------- # library(ggchord) # library(dplyr) # library(purrr) # library(tidyr) # # Original repository records; extract the tutorial input bundle first. # show_sequences <- read_fasta_lengths(files = "examples/fasta/*.fna") |> # select(accver, length) # show_annotations <- read_gff3(files = "examples/gff3/*.gff3", unstranded = "keep") # show_genes <- show_annotations |> # filter(type == "CDS") |> # select(accver, start, end, strand, anno) # show_alignments <- read_blast(files = "examples/blastn-standard/*.tsv") # # At most three substantial observed HSPs per pair; no coordinates are invented. # show_links <- show_alignments |> # filter(length >= 900, pident >= 75) |> # group_by(qaccver, saccver) |> # slice_max(length, n = 3, with_ties = FALSE) |> # ungroup() # # Select up to four source CDS per sequence, spread across annotation order. # # No feature is moved or invented; keep complete source tables separately. # show_features <- show_genes |> # filter(end - start + 1 >= 500) |> # group_by(accver) |> # arrange(start, .by_group = TRUE) |> # mutate(display_bin = ntile(row_number(), 4), informative = !is.na(anno) & anno != "hypothetical protein") |> # group_by(accver, display_bin) |> # arrange(desc(informative), desc(end - start), .by_group = TRUE) |> # slice_head(n = 1) |> # ungroup() |> # select(accver, start, end, strand, anno) # # Also retain well-separated negative-strand CDS, if not already represented. # show_reverse <- show_genes |> # filter(strand == "-", end - start + 1 >= 500) |> # group_by(accver) |> # arrange(start, .by_group = TRUE) |> # slice(c(1, n())) |> # ungroup() # show_features <- bind_rows(show_features, show_reverse) |> # distinct(accver, start, end, strand, .keep_all = TRUE) |> # arrange(accver, start) # show_labels <- show_features # show_optimization <- optimize_ggchord_layout( # show_sequences, show_links, # optimize = "orientation", weight = "length" # ) ## ----example-12, eval=FALSE--------------------------------------------------- # link_endpoints <- show_links |> # select(qaccver, saccver, qstart, qend, sstart, send) |> # pivot_longer( # cols = c(qstart, qend, sstart, send), # names_to = c("side", "endpoint"), names_pattern = "([qs])(start|end)", # values_to = "position" # ) # head(link_endpoints) ## ----example-13, eval=FALSE--------------------------------------------------- # library(ggchord) # blast_source <- read_blast(file = "result.tsv") ## ----example-14, eval=FALSE--------------------------------------------------- # # Optional appearance choices; none of these columns or scales is required input. # library(ggchord) # seq_palette <- setNames( # c("#E64B35", "#233E4A", "#269B82", "#D59B17"), # c("MT108731.1", "MT118296.1", "OQ646790.1", "OR222515.1") # ) # gene_palette <- c("+" = "#CE1536", "-" = "#346B9D", "." = "#8291A6") # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } # # library(ggchord) # library(dplyr) # library(purrr) # library(tidyr) # # Original repository records; extract the tutorial input bundle first. # show_sequences <- read_fasta_lengths(files = "examples/fasta/*.fna") |> # select(accver, length) # show_annotations <- read_gff3(files = "examples/gff3/*.gff3", unstranded = "keep") # show_genes <- show_annotations |> # filter(type == "CDS") |> # select(accver, start, end, strand, anno) # show_alignments <- read_blast(files = "examples/blastn-standard/*.tsv") # # At most three substantial observed HSPs per pair; no coordinates are invented. # show_links <- show_alignments |> # filter(length >= 900, pident >= 75) |> # group_by(qaccver, saccver) |> # slice_max(length, n = 3, with_ties = FALSE) |> # ungroup() # # Select up to four source CDS per sequence, spread across annotation order. # # No feature is moved or invented; keep complete source tables separately. # show_features <- show_genes |> # filter(end - start + 1 >= 500) |> # group_by(accver) |> # arrange(start, .by_group = TRUE) |> # mutate(display_bin = ntile(row_number(), 4), informative = !is.na(anno) & anno != "hypothetical protein") |> # group_by(accver, display_bin) |> # arrange(desc(informative), desc(end - start), .by_group = TRUE) |> # slice_head(n = 1) |> # ungroup() |> # select(accver, start, end, strand, anno) # # Also retain well-separated negative-strand CDS, if not already represented. # show_reverse <- show_genes |> # filter(strand == "-", end - start + 1 >= 500) |> # group_by(accver) |> # arrange(start, .by_group = TRUE) |> # slice(c(1, n())) |> # ungroup() # show_features <- bind_rows(show_features, show_reverse) |> # distinct(accver, start, end, strand, .keep_all = TRUE) |> # arrange(accver, start) # show_labels <- show_features # show_optimization <- optimize_ggchord_layout( # show_sequences, show_links, # optimize = "orientation", weight = "length" # ) ## ----example-15, eval=FALSE--------------------------------------------------- # p_learn_seq <- ggchord(seq_data = show_sequences) + # geom_seq() + # geom_seq_label(size = 4) + # scale_seq_colour_manual( # values = seq_palette, name = "Sequence", # guide = sequence_guide # ) + # scale_seq_position_continuous(breaks = c(0, 20000, 40000), minor_breaks = NULL) + # coord_chord() + # doc_theme # p_learn_seq ## ----example-16, eval=FALSE--------------------------------------------------- # p_learn_links <- p_learn_seq + # geom_link_ribbon( # aes(ribbon_fill = pident), # data = show_links, # adjacent_only = FALSE, link_avoid = "smooth" # ) + # scale_ribbon_fill_gradientn( # colours = identity_palette, values = identity_values, limits = c(0, 100), # breaks = identity_breaks, name = "Identity (%)", guide = identity_guide # ) # p_learn_links ## ----example-17, eval=FALSE--------------------------------------------------- # p_learn_genes <- p_learn_links + # geom_gene(data = show_features, position = position_strand(0.08)) + # scale_gene_fill_manual(values = gene_palette, breaks = c("+", "-"), name = "Strand") + # guides(gene_fill = strand_guide) # p_learn_genes ## ----example-18, eval=FALSE--------------------------------------------------- # p_learn_labels <- p_learn_genes + # geom_gene_label_repel( # data = show_features, label_layout = "arc", label_wrap = 24, # label_max_lines = 3, position = position_strand(0.08), size = 3 # ) # p_learn_labels ## ----example-19, eval=FALSE--------------------------------------------------- # # Optional appearance choices; none of these columns or scales is required input. # library(ggchord) # seq_palette <- setNames( # c("#E64B35", "#233E4A", "#269B82", "#D59B17"), # c("MT108731.1", "MT118296.1", "OQ646790.1", "OR222515.1") # ) # gene_palette <- c("+" = "#CE1536", "-" = "#346B9D", "." = "#8291A6") # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } # # library(ggchord) # library(dplyr) # library(purrr) # library(tidyr) # # Original repository records; extract the tutorial input bundle first. # show_sequences <- read_fasta_lengths(files = "examples/fasta/*.fna") |> # select(accver, length) # show_annotations <- read_gff3(files = "examples/gff3/*.gff3", unstranded = "keep") # show_genes <- show_annotations |> # filter(type == "CDS") |> # select(accver, start, end, strand, anno) # show_alignments <- read_blast(files = "examples/blastn-standard/*.tsv") # # At most three substantial observed HSPs per pair; no coordinates are invented. # show_links <- show_alignments |> # filter(length >= 900, pident >= 75) |> # group_by(qaccver, saccver) |> # slice_max(length, n = 3, with_ties = FALSE) |> # ungroup() # # Select up to four source CDS per sequence, spread across annotation order. # # No feature is moved or invented; keep complete source tables separately. # show_features <- show_genes |> # filter(end - start + 1 >= 500) |> # group_by(accver) |> # arrange(start, .by_group = TRUE) |> # mutate(display_bin = ntile(row_number(), 4), informative = !is.na(anno) & anno != "hypothetical protein") |> # group_by(accver, display_bin) |> # arrange(desc(informative), desc(end - start), .by_group = TRUE) |> # slice_head(n = 1) |> # ungroup() |> # select(accver, start, end, strand, anno) # # Also retain well-separated negative-strand CDS, if not already represented. # show_reverse <- show_genes |> # filter(strand == "-", end - start + 1 >= 500) |> # group_by(accver) |> # arrange(start, .by_group = TRUE) |> # slice(c(1, n())) |> # ungroup() # show_features <- bind_rows(show_features, show_reverse) |> # distinct(accver, start, end, strand, .keep_all = TRUE) |> # arrange(accver, start) # show_labels <- show_features # show_optimization <- optimize_ggchord_layout( # show_sequences, show_links, # optimize = "orientation", weight = "length" # ) ## ----example-20, eval=FALSE--------------------------------------------------- # show_curvature <- c( # "MT108731.1" = 0, "MT118296.1" = 0.85, # "OQ646790.1" = 1.1, "OR222515.1" = -0.4 # ) # show_radius <- c( # "MT108731.1" = 1, "MT118296.1" = 1.05, # "OQ646790.1" = 0.95, "OR222515.1" = 1.08 # ) # p_show_chord <- ggchord(show_sequences, show_links, show_features) + # geom_link_ribbon(aes(ribbon_fill = pident), link_avoid = "smooth", alpha = 0.58, linewidth = 0.18) + # geom_seq( # seq_orientation = show_optimization$seq_orientation, # seq_curvature = show_curvature, seq_radius = show_radius, # seq_gap = 0.06, linewidth = 1.1 # ) + # geom_gene(position = position_strand(0.07), width = 0.035) + # geom_gene_label_repel( # data = show_labels, label_layout = "radial", label_wrap = 24, # size = 3.1, label_max_lines = 3, position = position_strand(0.07) # ) + # geom_seq_label(size = 4, seq_label_radius = 0.84, seq_label_hjust = 0.5) + # scale_seq_colour_manual( # values = seq_palette, name = "Sequence", # guide = sequence_guide # ) + # scale_gene_fill_manual(values = gene_palette, breaks = c("+", "-"), name = "Strand") + # scale_ribbon_fill_gradientn( # colours = identity_palette, values = identity_values, # limits = c(0, 100), breaks = identity_breaks, name = "Identity (%)", guide = identity_guide # ) + # guides(gene_fill = strand_guide) + # scale_seq_position_continuous(breaks = c(0, 20000, 40000), minor_breaks = NULL) + # coord_chord(rotation = 50) + # doc_theme + # labs( # title = "Four genomes, one connected view", # subtitle = "Signed curvature and radius · selected strand-aware CDS · local ribbon clearance" # ) # p_show_chord ## ----example-21, eval=FALSE--------------------------------------------------- # library(ggchord) # library(dplyr) # data(seq_data_example) # data(ribbon_data_example) # data(gene_data_example) # seqs <- seq_data_example |> # select(accver, length) # links <- ribbon_data_example |> # select(qaccver, saccver, qstart, qend, sstart, send) # # The fixture preserves biological strands separately from its old display strands. # genes <- gene_data_example |> # transmute(accver, start, end, strand = source_strand) # # # Optional appearance choices; none of these columns or scales is required input. # library(ggchord) # seq_palette <- setNames( # c("#E64B35", "#233E4A", "#269B82", "#D59B17"), # c("MT108731.1", "MT118296.1", "OQ646790.1", "OR222515.1") # ) # gene_palette <- c("+" = "#CE1536", "-" = "#346B9D", "." = "#8291A6") # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } # # library(dplyr) # sequence_source <- read_fasta_lengths(files = "examples/fasta/*.fna", source_file = TRUE) # seq_real <- sequence_source |> # select(accver, length) # annotation_source <- read_gff3(files = "examples/gff3/*.gff3", unstranded = "keep", source_file = TRUE) # feature_real <- annotation_source |> # select(accver, start, end, strand, type, anno) # alignment_source <- read_blast(files = "examples/blastn-standard/*.tsv", source_file = TRUE) # link_real <- alignment_source |> # select(qaccver, saccver, qstart, qend, sstart, send, pident, length) # real_validation <- validate_ggchord_data(seq_real, link_real) # real_validation # # real_prepared <- prepare_ggchord_data(seq_real, link_real, # features = list(annotation = feature_real) # ) # real_filtered <- filter_ggchord_ribbons(link_real, min_pident = 80, min_length = 500) # # Explicit display subset: the six longest CDS per sequence; raw annotations stay intact. # real_display_features <- feature_real |> # group_by(accver) |> # slice_max(end - start + 1, n = 6, with_ties = FALSE) |> # ungroup() # real_labels <- real_display_features |> # filter(anno != "hypothetical protein") |> # group_by(accver) |> # slice_head(n = 2) |> # ungroup() # real_filtered$report ## ----example-22, eval=FALSE--------------------------------------------------- # unique_links <- deduplicate_ggchord_ribbons(real_filtered$data, # by = "overlap", # min_reciprocal_overlap = .9 # ) # merged_links <- merge_ggchord_ribbons(unique_links$data, # max_gap = 100, # min_pident_difference = 2 # ) # bundled_links <- bundle_ggchord_ribbons(link_real, seq_real, bins = 30) # p_density <- ggchord(seq_real, link_real, validate = "none") + # geom_seq(show.legend = FALSE) + # geom_seq_label(size = 4, seq_label_radius = .84, seq_label_hjust = .5) + # stat_ribbon_density(aes(ribbon_alpha = after_stat(density)), fill = "#269B82") + # scale_ribbon_alpha_continuous( # name = "Relative density", range = c(.1, .6), # guide = guide_ggchord_legend(position = "bottom", nrow = 1, theme = doc_legend_theme) # ) + # scale_seq_colour_manual(values = seq_palette) + # coord_chord() + # doc_theme # p_bundle <- ggchord(seq_real, link_real, validate = "none") + # geom_seq() + # stat_ribbon_bundle(bins = 30, fill = "#269B82") + # coord_chord() + # doc_theme # p_density ## ----example-23, eval=FALSE--------------------------------------------------- # # Midpoints summarize existing HSPs for display; they are not new alignments. # points <- ribbon_data_example |> # transmute( # qaccver, saccver, # qpos = floor((qstart + qend) / 2), # spos = floor((sstart + send) / 2), # pident # ) # p_points <- ggchord(seq_data = seqs) + # geom_seq() + # geom_seq_label(size = 4) + # geom_link_line(aes(link_colour = pident), # data = points, # arrow = arrow(length = unit(2, "mm")), linewidth = .7 # ) + # scale_link_colour_gradientn( # colours = identity_palette, values = identity_values, # limits = c(0, 100), breaks = identity_breaks, # name = "Identity (%)", guide = guide_ggchord_colourbar( # available_aes = "link_colour", # position = "bottom", direction = "horizontal", # theme = doc_legend_theme + theme(legend.key.width = unit(65, "mm"), legend.key.height = unit(5, "mm")) # ) # ) + # scale_seq_colour_manual(values = seq_palette) + # coord_chord() + # doc_theme # p_points ## ----example-24, eval=FALSE--------------------------------------------------- # # Optional appearance choices; none of these columns or scales is required input. # library(ggchord) # seq_palette <- setNames( # c("#E64B35", "#233E4A", "#269B82", "#D59B17"), # c("MT108731.1", "MT118296.1", "OQ646790.1", "OR222515.1") # ) # gene_palette <- c("+" = "#CE1536", "-" = "#346B9D", "." = "#8291A6") # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } ## ----example-25, eval=FALSE--------------------------------------------------- # library(dplyr) # sequence_source <- read_fasta_lengths(files = "examples/fasta/*.fna", source_file = TRUE) # seq_real <- sequence_source |> # select(accver, length) # annotation_source <- read_gff3(files = "examples/gff3/*.gff3", unstranded = "keep", source_file = TRUE) # feature_real <- annotation_source |> # select(accver, start, end, strand, type, anno) # alignment_source <- read_blast(files = "examples/blastn-standard/*.tsv", source_file = TRUE) # link_real <- alignment_source |> # select(qaccver, saccver, qstart, qend, sstart, send, pident, length) # real_validation <- validate_ggchord_data(seq_real, link_real) # real_validation ## ----example-26, eval=FALSE--------------------------------------------------- # real_prepared <- prepare_ggchord_data(seq_real, link_real, # features = list(annotation = feature_real) # ) # real_filtered <- filter_ggchord_ribbons(link_real, min_pident = 80, min_length = 500) # # Explicit display subset: the six longest CDS per sequence; raw annotations stay intact. # real_display_features <- feature_real |> # group_by(accver) |> # slice_max(end - start + 1, n = 6, with_ties = FALSE) |> # ungroup() # real_labels <- real_display_features |> # filter(anno != "hypothetical protein") |> # group_by(accver) |> # slice_head(n = 2) |> # ungroup() # real_filtered$report ## ----example-27, eval=FALSE--------------------------------------------------- # p_real_chord <- ggchord(seq_data = seq_real, ribbon_data = real_filtered$data, validate = "none") + # geom_link_ribbon(aes(ribbon_fill = pident), alpha = .65) + # geom_seq(show.legend = FALSE) + # geom_feature(aes(feature_fill = strand), # data = real_display_features, # width = .045, position = position_strand(.08) # ) + # geom_feature_label_repel( # data = real_labels, label_layout = "radial", # label_wrap = 20, size = 4, position = position_strand(.08) # ) + # geom_seq_label(size = 4, seq_label_radius = .84, seq_label_hjust = .5) + # coord_chord() + # scale_seq_colour_manual(values = seq_palette) + # scale_feature_fill_manual(values = gene_palette, breaks = c("+", "-")) + # scale_ribbon_fill_gradientn( # colours = identity_palette, values = identity_values, # limits = c(0, 100), breaks = identity_breaks, name = "Identity (%)", guide = identity_guide # ) + # doc_theme # p_real_chord ## ----example-28, eval=FALSE--------------------------------------------------- # p_real_linear <- ggchord(seq_real, real_filtered$data, validate = "none") + # geom_link_ribbon(fill = "#6BAFA6", alpha = .45) + # geom_seq() + # geom_seq_label(size = 4) + # geom_feature(data = feature_real, feature_fill = "#CE1536", show.legend = FALSE, position = position_strand(.18)) + # geom_seq_scale(length = 10000) + # coord_collinear() + # doc_axis() + # scale_seq_colour_manual(values = seq_palette) + # doc_theme # p_real_linear ## ----example-29, eval=FALSE--------------------------------------------------- # ggsave("phage-comparison.png", p_real_chord, width = 10, height = 8, units = "in", dpi = 144) ## ----example-30, eval=FALSE--------------------------------------------------- # # Optional appearance choices; none of these columns or scales is required input. # library(ggchord) # seq_palette <- setNames( # c("#E64B35", "#233E4A", "#269B82", "#D59B17"), # c("MT108731.1", "MT118296.1", "OQ646790.1", "OR222515.1") # ) # gene_palette <- c("+" = "#CE1536", "-" = "#346B9D", "." = "#8291A6") # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } ## ----example-31, eval=FALSE--------------------------------------------------- # library(ggchord) # library(dplyr) # library(purrr) # library(tidyr) # # Original repository records; extract the tutorial input bundle first. # show_sequences <- read_fasta_lengths(files = "examples/fasta/*.fna") |> # select(accver, length) # show_annotations <- read_gff3(files = "examples/gff3/*.gff3", unstranded = "keep") # show_genes <- show_annotations |> # filter(type == "CDS") |> # select(accver, start, end, strand, anno) # show_alignments <- read_blast(files = "examples/blastn-standard/*.tsv") # # At most three substantial observed HSPs per pair; no coordinates are invented. # show_links <- show_alignments |> # filter(length >= 900, pident >= 75) |> # group_by(qaccver, saccver) |> # slice_max(length, n = 3, with_ties = FALSE) |> # ungroup() # # Select up to four source CDS per sequence, spread across annotation order. # # No feature is moved or invented; keep complete source tables separately. # show_features <- show_genes |> # filter(end - start + 1 >= 500) |> # group_by(accver) |> # arrange(start, .by_group = TRUE) |> # mutate(display_bin = ntile(row_number(), 4), informative = !is.na(anno) & anno != "hypothetical protein") |> # group_by(accver, display_bin) |> # arrange(desc(informative), desc(end - start), .by_group = TRUE) |> # slice_head(n = 1) |> # ungroup() |> # select(accver, start, end, strand, anno) # # Also retain well-separated negative-strand CDS, if not already represented. # show_reverse <- show_genes |> # filter(strand == "-", end - start + 1 >= 500) |> # group_by(accver) |> # arrange(start, .by_group = TRUE) |> # slice(c(1, n())) |> # ungroup() # show_features <- bind_rows(show_features, show_reverse) |> # distinct(accver, start, end, strand, .keep_all = TRUE) |> # arrange(accver, start) # show_labels <- show_features # show_optimization <- optimize_ggchord_layout( # show_sequences, show_links, # optimize = "orientation", weight = "length" # ) ## ----example-32, eval=FALSE--------------------------------------------------- # show_curvature <- c( # "MT108731.1" = 0, "MT118296.1" = 0.85, # "OQ646790.1" = 1.1, "OR222515.1" = -0.4 # ) # show_radius <- c( # "MT108731.1" = 1, "MT118296.1" = 1.05, # "OQ646790.1" = 0.95, "OR222515.1" = 1.08 # ) # p_show_chord <- ggchord(show_sequences, show_links, show_features) + # geom_link_ribbon(aes(ribbon_fill = pident), link_avoid = "smooth", alpha = 0.58, linewidth = 0.18) + # geom_seq( # seq_orientation = show_optimization$seq_orientation, # seq_curvature = show_curvature, seq_radius = show_radius, # seq_gap = 0.06, linewidth = 1.1 # ) + # geom_gene(position = position_strand(0.07), width = 0.035) + # geom_gene_label_repel( # data = show_labels, label_layout = "radial", label_wrap = 24, # size = 3.1, label_max_lines = 3, position = position_strand(0.07) # ) + # geom_seq_label(size = 4, seq_label_radius = 0.84, seq_label_hjust = 0.5) + # scale_seq_colour_manual( # values = seq_palette, name = "Sequence", # guide = sequence_guide # ) + # scale_gene_fill_manual(values = gene_palette, breaks = c("+", "-"), name = "Strand") + # scale_ribbon_fill_gradientn( # colours = identity_palette, values = identity_values, # limits = c(0, 100), breaks = identity_breaks, name = "Identity (%)", guide = identity_guide # ) + # guides(gene_fill = strand_guide) + # scale_seq_position_continuous(breaks = c(0, 20000, 40000), minor_breaks = NULL) + # coord_chord(rotation = 50) + # doc_theme + # labs( # title = "Four genomes, one connected view", # subtitle = "Signed curvature and radius · selected strand-aware CDS · local ribbon clearance" # ) # p_show_chord ## ----example-33, eval=FALSE--------------------------------------------------- # ggsave("four-genome-chord.png", p_show_chord, width = 12, height = 10, units = "in", dpi = 144) ## ----example-34, eval=FALSE--------------------------------------------------- # # Optional appearance choices; none of these columns or scales is required input. # library(ggchord) # seq_palette <- setNames( # c("#E64B35", "#233E4A", "#269B82", "#D59B17"), # c("MT108731.1", "MT118296.1", "OQ646790.1", "OR222515.1") # ) # gene_palette <- c("+" = "#CE1536", "-" = "#346B9D", "." = "#8291A6") # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } # # library(ggchord) # library(dplyr) # library(purrr) # library(tidyr) # # Original repository records; extract the tutorial input bundle first. # show_sequences <- read_fasta_lengths(files = "examples/fasta/*.fna") |> # select(accver, length) # show_annotations <- read_gff3(files = "examples/gff3/*.gff3", unstranded = "keep") # show_genes <- show_annotations |> # filter(type == "CDS") |> # select(accver, start, end, strand, anno) # show_alignments <- read_blast(files = "examples/blastn-standard/*.tsv") # # At most three substantial observed HSPs per pair; no coordinates are invented. # show_links <- show_alignments |> # filter(length >= 900, pident >= 75) |> # group_by(qaccver, saccver) |> # slice_max(length, n = 3, with_ties = FALSE) |> # ungroup() # # Select up to four source CDS per sequence, spread across annotation order. # # No feature is moved or invented; keep complete source tables separately. # show_features <- show_genes |> # filter(end - start + 1 >= 500) |> # group_by(accver) |> # arrange(start, .by_group = TRUE) |> # mutate(display_bin = ntile(row_number(), 4), informative = !is.na(anno) & anno != "hypothetical protein") |> # group_by(accver, display_bin) |> # arrange(desc(informative), desc(end - start), .by_group = TRUE) |> # slice_head(n = 1) |> # ungroup() |> # select(accver, start, end, strand, anno) # # Also retain well-separated negative-strand CDS, if not already represented. # show_reverse <- show_genes |> # filter(strand == "-", end - start + 1 >= 500) |> # group_by(accver) |> # arrange(start, .by_group = TRUE) |> # slice(c(1, n())) |> # ungroup() # show_features <- bind_rows(show_features, show_reverse) |> # distinct(accver, start, end, strand, .keep_all = TRUE) |> # arrange(accver, start) # show_labels <- show_features # show_optimization <- optimize_ggchord_layout( # show_sequences, show_links, # optimize = "orientation", weight = "length" # ) ## ----example-35, eval=FALSE--------------------------------------------------- # p_learn_seq <- ggchord(seq_data = show_sequences) + # geom_seq() + # geom_seq_label(size = 4) + # scale_seq_colour_manual( # values = seq_palette, name = "Sequence", # guide = sequence_guide # ) + # scale_seq_position_continuous(breaks = c(0, 20000, 40000), minor_breaks = NULL) + # coord_chord() + # doc_theme # p_learn_seq # # p_learn_links <- p_learn_seq + # geom_link_ribbon( # aes(ribbon_fill = pident), # data = show_links, # adjacent_only = FALSE, link_avoid = "smooth" # ) + # scale_ribbon_fill_gradientn( # colours = identity_palette, values = identity_values, limits = c(0, 100), # breaks = identity_breaks, name = "Identity (%)", guide = identity_guide # ) # p_learn_links # # p_learn_genes <- p_learn_links + # geom_gene(data = show_features, position = position_strand(0.08)) + # scale_gene_fill_manual(values = gene_palette, breaks = c("+", "-"), name = "Strand") + # guides(gene_fill = strand_guide) # p_learn_genes # # p_learn_labels <- p_learn_genes + # geom_gene_label_repel( # data = show_features, label_layout = "arc", label_wrap = 24, # label_max_lines = 3, position = position_strand(0.08), size = 3 # ) # p_learn_labels ## ----example-36, eval=FALSE--------------------------------------------------- # # Reuse exactly the same data and layers; replace only the coordinate system. # p_linear_seq <- p_learn_seq + coord_collinear() # p_linear_links <- p_learn_links + coord_collinear() # p_linear_labels <- p_learn_labels + coord_collinear(ratio = 0.65) # p_linear_labels ## ----example-37, eval=FALSE--------------------------------------------------- # # Optional appearance choices; none of these columns or scales is required input. # library(ggchord) # seq_palette <- setNames( # c("#E64B35", "#233E4A", "#269B82", "#D59B17"), # c("MT108731.1", "MT118296.1", "OQ646790.1", "OR222515.1") # ) # gene_palette <- c("+" = "#CE1536", "-" = "#346B9D", "." = "#8291A6") # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } # # library(ggchord) # library(dplyr) # library(purrr) # library(tidyr) # # Original repository records; extract the tutorial input bundle first. # show_sequences <- read_fasta_lengths(files = "examples/fasta/*.fna") |> # select(accver, length) # show_annotations <- read_gff3(files = "examples/gff3/*.gff3", unstranded = "keep") # show_genes <- show_annotations |> # filter(type == "CDS") |> # select(accver, start, end, strand, anno) # show_alignments <- read_blast(files = "examples/blastn-standard/*.tsv") # # At most three substantial observed HSPs per pair; no coordinates are invented. # show_links <- show_alignments |> # filter(length >= 900, pident >= 75) |> # group_by(qaccver, saccver) |> # slice_max(length, n = 3, with_ties = FALSE) |> # ungroup() # # Select up to four source CDS per sequence, spread across annotation order. # # No feature is moved or invented; keep complete source tables separately. # show_features <- show_genes |> # filter(end - start + 1 >= 500) |> # group_by(accver) |> # arrange(start, .by_group = TRUE) |> # mutate(display_bin = ntile(row_number(), 4), informative = !is.na(anno) & anno != "hypothetical protein") |> # group_by(accver, display_bin) |> # arrange(desc(informative), desc(end - start), .by_group = TRUE) |> # slice_head(n = 1) |> # ungroup() |> # select(accver, start, end, strand, anno) # # Also retain well-separated negative-strand CDS, if not already represented. # show_reverse <- show_genes |> # filter(strand == "-", end - start + 1 >= 500) |> # group_by(accver) |> # arrange(start, .by_group = TRUE) |> # slice(c(1, n())) |> # ungroup() # show_features <- bind_rows(show_features, show_reverse) |> # distinct(accver, start, end, strand, .keep_all = TRUE) |> # arrange(accver, start) # show_labels <- show_features # show_optimization <- optimize_ggchord_layout( # show_sequences, show_links, # optimize = "orientation", weight = "length" # ) ## ----example-38, eval=FALSE--------------------------------------------------- # # Omit two short, crowded candidates from this linear display only. # # Both source rows remain in show_features and in the chord display. # show_linear_features <- show_features |> # filter(!anno %in% c("tail assembly protein", "gp4 protein")) # # Separate adjacent and cross-row links, retaining original endpoint values. # show_pair_rows <- show_links |> # mutate( # query_row = match(qaccver, show_sequences$accver), # subject_row = match(saccver, show_sequences$accver), # span = abs(query_row - subject_row) # ) # show_adjacent <- show_pair_rows |> filter(span == 1) # # One longest observed HSP per non-neighbour pair keeps the long routes legible. # show_cross <- show_pair_rows |> # filter(span > 1) |> # group_by(qaccver, saccver) |> # slice_max(length, n = 1, with_ties = FALSE) |> # ungroup() # # Inspect a backbone-only layout using the same alignment and gap as the final view. # show_linear_coord <- coord_collinear(align = "center", orientation = "auto", ratio = 0.7) # show_routing_base <- ggchord(show_sequences, show_links) + # geom_seq(seq_gap = 0.42) + show_linear_coord # show_reference <- get_chord_layout(show_routing_base)$sequence_reference # # Shared controls sit just beyond the intervening backbones, rather than far away. # show_cross_controls <- map2(show_cross$query_row, show_cross$subject_row, function(query_row, subject_row) { # between <- seq(min(query_row, subject_row) + 1L, max(query_row, subject_row) - 1L) # middle_ids <- show_sequences$accver[between] # right_edge <- max(show_reference$x_offset[middle_ids] + show_reference$lens[middle_ids] * show_reference$x_scale) # middle_y <- mean(show_reference$y[show_sequences$accver[c(query_row, subject_row)]]) # c(right_edge + if (abs(query_row - subject_row) == 3) 0.85 else 0.65, middle_y) # }) # p_show_collinear <- ggchord(show_sequences, gene_data = show_linear_features) + # geom_link_ribbon( # aes(ribbon_fill = pident), # data = show_adjacent, # link_avoid = "smooth", alpha = 0.48, linewidth = 0.22 # ) + # geom_link_ribbon( # aes(ribbon_fill = pident), # data = show_cross, # adjacent_only = FALSE, link_ctrl_point = show_cross_controls, # link_avoid = "smooth", alpha = 0.28, linewidth = 0.3 # ) + # geom_seq(seq_gap = 0.42, linewidth = 1) + # geom_gene(position = position_strand(0.09), width = 0.045) + # geom_gene_label_repel( # data = show_linear_features |> arrange(desc(nchar(anno)), accver, start), # label_layout = "auto", label_wrap = 24, # label_side = "auto", position = position_strand(0.09), size = 3.3 # ) + # geom_seq_label(size = 4) + # geom_seq_scale(length = 10000, offset = unit(c(0, 10), "mm")) + # scale_seq_colour_manual( # values = seq_palette, name = "Sequence", # guide = sequence_guide # ) + # scale_gene_fill_manual(values = gene_palette, breaks = c("+", "-"), name = "Strand") + # scale_ribbon_fill_gradientn( # colours = identity_palette, values = identity_values, # limits = c(0, 100), breaks = identity_breaks, name = "Identity (%)", guide = identity_guide # ) + # guides(gene_fill = strand_guide) + # scale_seq_position_continuous(breaks = c(0, 20000, 40000, 60000, 80000), minor_breaks = NULL) + # show_linear_coord + # doc_theme + # labs( # title = "Sequence architecture, across four genomes", # subtitle = "Selected labelled CDS · strand-separated arrows · routed cross-row comparisons" # ) # p_show_collinear ## ----example-39, eval=FALSE--------------------------------------------------- # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } ## ----example-40, eval=FALSE--------------------------------------------------- # library(dplyr) # data(plasmid_example_pUC19) # data(plasmid_example_pETDuet_1) # homolog_source <- bind_rows(plasmid_example_pUC19, plasmid_example_pETDuet_1) # homolog_seq <- homolog_source |> # select(accver, length) # cluster_features <- homolog_source |> # find_common_features() |> # filter(type == "CDS", anno == "AmpR") |> # select(accver, start, end, strand, anno) |> # arrange(match(accver, homolog_source$accver)) |> # mutate(feature_id = paste(accver, start, end, sep = ":")) # # Membership is computed from exact strand-corrected CDS DNA equality, not gene names. # reverse_complement <- function(sequence) { # bases <- strsplit(sequence, "")[[1]] # chartr("ACGT", "TGCA", paste(rev(bases), collapse = "")) # } # coding_sequences <- cluster_features |> # left_join(select(homolog_source, accver, sequence), by = "accver") |> # rowwise() |> # mutate( # coding_dna = substr(sequence, start, end), # coding_dna = if (strand == "-") reverse_complement(coding_dna) else coding_dna # ) |> # ungroup() # cds_dna <- pull(coding_sequences, coding_dna) ## ----example-41, eval=FALSE--------------------------------------------------- # membership <- tibble( # feature_id = cluster_features$feature_id, # cluster_id = paste0("exact_CDS_", match(cds_dna, unique(cds_dna))) # ) # homologs <- link_ggchord_clusters(cluster_features, membership, # pairs = tibble(qaccver = homolog_seq$accver[1], saccver = homolog_seq$accver[2]) # ) # loci <- find_ggchord_loci(homolog_seq, cluster_features, # targets = cluster_features$feature_id[cluster_features$anno == "AmpR"][1], upstream = 100, downstream = 100 # ) ## ----example-42, eval=FALSE--------------------------------------------------- # # Compare equivalent local windows instead of shrinking one gene on full plasmids. # homolog_windows <- cluster_features |> # select(accver, start, end) # local_homologs <- focus_ggchord_data(homolog_seq, homologs$data, cluster_features, # loci = homolog_windows, expand = 200 # ) # p_homologs <- ggchord( # seq_data = local_homologs$seq_data, # ribbon_data = local_homologs$ribbon_data, gene_data = local_homologs$gene_data # ) + # geom_link_ribbon(fill = "#269B82", alpha = .65) + # geom_seq(colour = "#233E4A", show.legend = FALSE) + # geom_feature( # data = local_homologs$gene_data, feature_fill = "#F1D334", # position = position_strand(.10), width = .07 # ) + # geom_feature_label_repel(data = local_homologs$gene_data, size = 4, label_fit = "shrink") + # geom_seq_label(size = 4) + # geom_seq_break() + # coord_collinear() + # doc_axis("source") + # doc_theme # p_homologs ## ----example-43, eval=FALSE--------------------------------------------------- # # Exact CDS DNA equality supplies evidence for these selected ungapped residues. # matching <- homologs$data |> # slice_head(n = 1) # paired_coding <- cluster_features |> # mutate(phase = 0) # paired_cds <- build_ggchord_cds(paired_coding) # protein_link <- tibble( # qfeature_id = matching$qfeature_id, sfeature_id = matching$sfeature_id, # qstart = 1, qend = 20, sstart = 1, send = 20 # ) # projected_links <- project_ggchord_links(protein_link, paired_cds, ungapped = TRUE) # projected_links ## ----example-44, eval=FALSE--------------------------------------------------- # ggsave("ampR-neighbourhood.png", p_homologs, width = 10, height = 5, units = "in", dpi = 144) ## ----example-45, eval=FALSE--------------------------------------------------- # library(ggchord) # library(dplyr) # data(seq_data_example) # data(ribbon_data_example) # data(gene_data_example) # seqs <- seq_data_example |> # select(accver, length) # links <- ribbon_data_example |> # select(qaccver, saccver, qstart, qend, sstart, send) # # The fixture preserves biological strands separately from its old display strands. # genes <- gene_data_example |> # transmute(accver, start, end, strand = source_strand) # # # Optional appearance choices; none of these columns or scales is required input. # library(ggchord) # seq_palette <- setNames( # c("#E64B35", "#233E4A", "#269B82", "#D59B17"), # c("MT108731.1", "MT118296.1", "OQ646790.1", "OR222515.1") # ) # gene_palette <- c("+" = "#CE1536", "-" = "#346B9D", "." = "#8291A6") # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } ## ----example-46, eval=FALSE--------------------------------------------------- # focus_loci <- tibble( # accver = seqs$accver, # start = floor(seqs$length * .25), end = floor(seqs$length * .75) # ) # focused <- focus_ggchord_data(seqs, links, genes, loci = focus_loci) # p_focus <- ggchord( # seq_data = focused$seq_data, ribbon_data = focused$ribbon_data, # gene_data = focused$gene_data # ) + # geom_seq() + # geom_link_ribbon(fill = "#6BAFA6", alpha = .5) + # geom_gene() + # geom_seq_break() + # geom_seq_label(size = 4) + # scale_seq_colour_manual(values = seq_palette) + # scale_gene_fill_manual(values = gene_palette) + # doc_axis("source") + # coord_collinear() + # doc_theme # p_focus ## ----example-47, eval=FALSE--------------------------------------------------- # windows <- tibble( # accver = rep(seqs$accver[1], 2), # start = c(1, floor(seqs$length[1] * .65)), end = c(floor(seqs$length[1] * .4), seqs$length[1]), # locus_id = c("left_window", "right_window"), row_id = seqs$accver[1], orientation = c(1, -1) # ) # segmented <- segment_ggchord_data(slice_head(seqs, n = 1), windows = windows) # p_segments <- ggchord(seq_data = segmented$seq_data) + # geom_seq() + # coord_collinear(segments = segmented$segments, segment_gap = 2000) # p_segments ## ----example-48, eval=FALSE--------------------------------------------------- # source_position <- focus_loci$start[1] + 100 # source_mark <- map_ggchord_position(p_focus, seqs$accver[1], source_position, source = TRUE) # p_marked <- p_focus + # ggplot2::geom_point( # data = source_mark, aes(x, y), # inherit.aes = FALSE, colour = "#D59B17", size = 3 # ) # p_marked ## ----example-49, eval=FALSE--------------------------------------------------- # # Optional appearance choices; none of these columns or scales is required input. # library(ggchord) # seq_palette <- setNames( # c("#E64B35", "#233E4A", "#269B82", "#D59B17"), # c("MT108731.1", "MT118296.1", "OQ646790.1", "OR222515.1") # ) # gene_palette <- c("+" = "#CE1536", "-" = "#346B9D", "." = "#8291A6") # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } ## ----example-50, eval=FALSE--------------------------------------------------- # library(ggchord) # library(dplyr) # library(purrr) # library(tidyr) # # Original repository records; extract the tutorial input bundle first. # show_sequences <- read_fasta_lengths(files = "examples/fasta/*.fna") |> # select(accver, length) # show_annotations <- read_gff3(files = "examples/gff3/*.gff3", unstranded = "keep") # show_genes <- show_annotations |> # filter(type == "CDS") |> # select(accver, start, end, strand, anno) # show_alignments <- read_blast(files = "examples/blastn-standard/*.tsv") # # At most three substantial observed HSPs per pair; no coordinates are invented. # show_links <- show_alignments |> # filter(length >= 900, pident >= 75) |> # group_by(qaccver, saccver) |> # slice_max(length, n = 3, with_ties = FALSE) |> # ungroup() # # Select up to four source CDS per sequence, spread across annotation order. # # No feature is moved or invented; keep complete source tables separately. # show_features <- show_genes |> # filter(end - start + 1 >= 500) |> # group_by(accver) |> # arrange(start, .by_group = TRUE) |> # mutate(display_bin = ntile(row_number(), 4), informative = !is.na(anno) & anno != "hypothetical protein") |> # group_by(accver, display_bin) |> # arrange(desc(informative), desc(end - start), .by_group = TRUE) |> # slice_head(n = 1) |> # ungroup() |> # select(accver, start, end, strand, anno) # # Also retain well-separated negative-strand CDS, if not already represented. # show_reverse <- show_genes |> # filter(strand == "-", end - start + 1 >= 500) |> # group_by(accver) |> # arrange(start, .by_group = TRUE) |> # slice(c(1, n())) |> # ungroup() # show_features <- bind_rows(show_features, show_reverse) |> # distinct(accver, start, end, strand, .keep_all = TRUE) |> # arrange(accver, start) # show_labels <- show_features # show_optimization <- optimize_ggchord_layout( # show_sequences, show_links, # optimize = "orientation", weight = "length" # ) ## ----example-51, eval=FALSE--------------------------------------------------- # # Omit two short, crowded candidates from this linear display only. # # Both source rows remain in show_features and in the chord display. # show_linear_features <- show_features |> # filter(!anno %in% c("tail assembly protein", "gp4 protein")) # # Separate adjacent and cross-row links, retaining original endpoint values. # show_pair_rows <- show_links |> # mutate( # query_row = match(qaccver, show_sequences$accver), # subject_row = match(saccver, show_sequences$accver), # span = abs(query_row - subject_row) # ) # show_adjacent <- show_pair_rows |> filter(span == 1) # # One longest observed HSP per non-neighbour pair keeps the long routes legible. # show_cross <- show_pair_rows |> # filter(span > 1) |> # group_by(qaccver, saccver) |> # slice_max(length, n = 1, with_ties = FALSE) |> # ungroup() # # Inspect a backbone-only layout using the same alignment and gap as the final view. # show_linear_coord <- coord_collinear(align = "center", orientation = "auto", ratio = 0.7) # show_routing_base <- ggchord(show_sequences, show_links) + # geom_seq(seq_gap = 0.42) + show_linear_coord # show_reference <- get_chord_layout(show_routing_base)$sequence_reference # # Shared controls sit just beyond the intervening backbones, rather than far away. # show_cross_controls <- map2(show_cross$query_row, show_cross$subject_row, function(query_row, subject_row) { # between <- seq(min(query_row, subject_row) + 1L, max(query_row, subject_row) - 1L) # middle_ids <- show_sequences$accver[between] # right_edge <- max(show_reference$x_offset[middle_ids] + show_reference$lens[middle_ids] * show_reference$x_scale) # middle_y <- mean(show_reference$y[show_sequences$accver[c(query_row, subject_row)]]) # c(right_edge + if (abs(query_row - subject_row) == 3) 0.85 else 0.65, middle_y) # }) # p_show_collinear <- ggchord(show_sequences, gene_data = show_linear_features) + # geom_link_ribbon( # aes(ribbon_fill = pident), # data = show_adjacent, # link_avoid = "smooth", alpha = 0.48, linewidth = 0.22 # ) + # geom_link_ribbon( # aes(ribbon_fill = pident), # data = show_cross, # adjacent_only = FALSE, link_ctrl_point = show_cross_controls, # link_avoid = "smooth", alpha = 0.28, linewidth = 0.3 # ) + # geom_seq(seq_gap = 0.42, linewidth = 1) + # geom_gene(position = position_strand(0.09), width = 0.045) + # geom_gene_label_repel( # data = show_linear_features |> arrange(desc(nchar(anno)), accver, start), # label_layout = "auto", label_wrap = 24, # label_side = "auto", position = position_strand(0.09), size = 3.3 # ) + # geom_seq_label(size = 4) + # geom_seq_scale(length = 10000, offset = unit(c(0, 10), "mm")) + # scale_seq_colour_manual( # values = seq_palette, name = "Sequence", # guide = sequence_guide # ) + # scale_gene_fill_manual(values = gene_palette, breaks = c("+", "-"), name = "Strand") + # scale_ribbon_fill_gradientn( # colours = identity_palette, values = identity_values, # limits = c(0, 100), breaks = identity_breaks, name = "Identity (%)", guide = identity_guide # ) + # guides(gene_fill = strand_guide) + # scale_seq_position_continuous(breaks = c(0, 20000, 40000, 60000, 80000), minor_breaks = NULL) + # show_linear_coord + # doc_theme + # labs( # title = "Sequence architecture, across four genomes", # subtitle = "Selected labelled CDS · strand-separated arrows · routed cross-row comparisons" # ) # p_show_collinear ## ----example-52, eval=FALSE--------------------------------------------------- # ggsave("four-genome-collinear.png", p_show_collinear, width = 12, height = 8, units = "in", dpi = 144) ## ----example-53, eval=FALSE--------------------------------------------------- # library(ggchord) # library(dplyr) # plasmid_example_pUC19 |> select(accver, length) # find_common_features(plasmid_example_pUC19) |> select(accver, start, end, strand, type, anno) |> slice_head(n = 5) ## ----example-54, eval=FALSE--------------------------------------------------- # library(ggchord) # feature_palette <- c( # CDS = "#993366", promoter = "#31849B", rep_origin = "#F1E51D", # primer_bind = "#A020F0", misc_feature = "#A6ACB3", protein_bind = "#31849B" # ) # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } # # library(dplyr) # data(plasmid_example_pUC19) # plasmid_seq <- plasmid_example_pUC19 |> # select(accver, length) # plasmid_feature_source <- find_common_features(plasmid_example_pUC19) # # Dedicated primer layers already show bindings; avoid duplicate primer annotations. # plasmid_features <- plasmid_feature_source |> # filter(type != "primer_bind") |> # select(accver, start, end, strand, type, anno) # plasmid_primers <- find_primer_bindings(plasmid_example_pUC19, set = "universal") # plasmid_sites <- find_restriction_sites(plasmid_example_pUC19) # display_sites <- filter_restriction_sites(plasmid_sites, enzymes = c("EcoRI", "BamHI", "HindIII")) # primer_names <- primer_catalog(set = "universal") ## ----example-55, eval=FALSE--------------------------------------------------- # plasmid_position <- position_feature_stack(spacing = 0.085, base_position = position_plasmid()) # p_plasmid_backbone <- ggchord(seq_data = plasmid_seq) + # geom_seq(seq_style = "double", colour = "#252525") + # geom_seq_center_label() + # coord_circular(rotation = 90) + # theme_ggchord_plasmid() + # theme( # panel.background = element_blank(), # plot.background = element_rect(fill = doc_background, colour = NA) # ) # p_plasmid_features <- p_plasmid_backbone + # geom_feature_plasmid(data = plasmid_features, position = plasmid_position) # p_plasmid_text <- p_plasmid_features + # geom_feature_label_repel( # data = plasmid_features, position = plasmid_position, # external = TRUE, label_fit = "shrink", label_min_size = 2.5 # ) # p_plasmid_text ## ----example-56, eval=FALSE--------------------------------------------------- # plasmid_position <- position_feature_stack(spacing = .10, base_position = position_plasmid()) # p_plasmid <- ggchord(seq_data = plasmid_seq) + # geom_seq(seq_style = "double") + # geom_feature_plasmid(aes(feature_fill = type), data = plasmid_features, position = plasmid_position) + # geom_feature_label_repel( # data = plasmid_features, position = plasmid_position, # size = 4 # ) + # geom_primer(data = plasmid_primers, colour = "#A020F0") + # geom_primer_label_repel(data = plasmid_primers, colour = "#A020F0", max_overlaps = 0, size = 3) + # geom_restriction_site(data = display_sites, colour = "#29364B", label_size = 3, leader = "radial") + # geom_seq_center_label() + # scale_feature_fill_manual(values = feature_palette) + # scale_seq_position_continuous() + # coord_circular(rotation = 90) + # theme_ggchord_plasmid() + # theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank() # ) # p_plasmid ## ----example-57, eval=FALSE--------------------------------------------------- # library(ggchord) # library(dplyr) ## ----example-58, eval=FALSE--------------------------------------------------- # data(plasmid_example_pUC19) # target <- setNames(plasmid_example_pUC19$sequence, plasmid_example_pUC19$accver) # primers <- find_primer_bindings(target, set = "universal") # sites <- find_restriction_sites(target, enzymes = c("EcoRI", "BamHI", "HindIII"), methylation = "none") # recognition_seq <- plasmid_example_pUC19 |> # select(accver, length) # p_sites <- ggchord(seq_data = recognition_seq) + # geom_seq() + # geom_primer(data = primers, colour = "#A020F0") + # geom_restriction_site(data = sites, colour = "#29364B") + # coord_circular() # p_sites ## ----example-59, eval=FALSE--------------------------------------------------- # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } ## ----example-60, eval=FALSE--------------------------------------------------- # library(ggchord) # library(dplyr) # library(purrr) # # Thirteen unique records; the compatibility spelling pUC19c is not an extra plasmid. # catalog_objects <- c( # "plasmid_example_pSB1C3", "plasmid_example_pUC19", "plasmid_example_pBR322", # "plasmid_example_pBluescript_II_SK_plus", "plasmid_example_pET_28a_plus", # "plasmid_example_pETDuet_1", "plasmid_example_pDONR221", # "plasmid_example_pcDNA3_1_plus", "plasmid_example_pTRE_Tight_BI", # "plasmid_example_pTRIPZ", "plasmid_example_pCAMBIA1300", # "plasmid_example_pEarleyGate_201", "plasmid_example_pSpCas9_BB_2A_GFP_PX458" # ) # read_catalog_record <- function(object) { # target <- new.env(parent = emptyenv()) # data(list = object, package = "ggchord", envir = target) # target[[object]] # } # catalog_sequences <- map(catalog_objects, read_catalog_record) # catalog_features <- map(catalog_sequences, find_common_features) # catalog_summary <- map2(catalog_sequences, catalog_features, function(sequence, features) { # tibble(accver = sequence$accver, length = sequence$length, features = nrow(features)) # }) |> # bind_rows() ## ----example-61, eval=FALSE--------------------------------------------------- # make_catalog_plot <- function(sequence, features) { # placement <- position_feature_stack(spacing = 0.085, base_position = position_plasmid()) # ggchord(seq_data = sequence |> select(accver, length)) + # geom_seq(seq_style = "double", colour = "#252525") + # geom_feature_plasmid(data = features, position = placement) + # geom_feature_label_repel( # data = features, position = placement, external = FALSE, # label_fit = "shrink", label_min_size = 2.5, max_overlaps = 0 # ) + # geom_seq_center_label() + # coord_circular(rotation = 90, fit = "manual", xlim = c(-1.35, 1.35), ylim = c(-1.35, 1.35)) + # theme_ggchord_plasmid() + # theme( # panel.background = element_blank(), # plot.background = element_rect(fill = doc_background, colour = NA) # ) # } # catalog_plots <- map2(catalog_sequences, catalog_features, make_catalog_plot) # names(catalog_plots) <- catalog_summary$accver # # Save any record at the same size used in the gallery. # # ggsave("pBR322.png", catalog_plots[["pBR322"]], width = 10, height = 8, dpi = 144) ## ----example-62, eval=FALSE--------------------------------------------------- # library(ggchord) # feature_palette <- c( # CDS = "#993366", promoter = "#31849B", rep_origin = "#F1E51D", # primer_bind = "#A020F0", misc_feature = "#A6ACB3", protein_bind = "#31849B" # ) # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } ## ----example-63, eval=FALSE--------------------------------------------------- # library(dplyr) # data(plasmid_example_pUC19) # plasmid_seq <- plasmid_example_pUC19 |> # select(accver, length) # plasmid_feature_source <- find_common_features(plasmid_example_pUC19) # # Dedicated primer layers already show bindings; avoid duplicate primer annotations. # plasmid_features <- plasmid_feature_source |> # filter(type != "primer_bind") |> # select(accver, start, end, strand, type, anno) # plasmid_primers <- find_primer_bindings(plasmid_example_pUC19, set = "universal") # plasmid_sites <- find_restriction_sites(plasmid_example_pUC19) # display_sites <- filter_restriction_sites(plasmid_sites, enzymes = c("EcoRI", "BamHI", "HindIII")) # primer_names <- primer_catalog(set = "universal") ## ----example-64, eval=FALSE--------------------------------------------------- # plasmid_position <- position_feature_stack(spacing = .10, base_position = position_plasmid()) # p_plasmid <- ggchord(seq_data = plasmid_seq) + # geom_seq(seq_style = "double") + # geom_feature_plasmid(aes(feature_fill = type), data = plasmid_features, position = plasmid_position) + # geom_feature_label_repel( # data = plasmid_features, position = plasmid_position, # size = 4 # ) + # geom_primer(data = plasmid_primers, colour = "#A020F0") + # geom_primer_label_repel(data = plasmid_primers, colour = "#A020F0", max_overlaps = 0, size = 3) + # geom_restriction_site(data = display_sites, colour = "#29364B", label_size = 3, leader = "radial") + # geom_seq_center_label() + # scale_feature_fill_manual(values = feature_palette) + # scale_seq_position_continuous() + # coord_circular(rotation = 90) + # theme_ggchord_plasmid() + # theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank() # ) # p_plasmid ## ----example-65, eval=FALSE--------------------------------------------------- # ggsave("pUC19-map.png", p_plasmid, width = 10, height = 9, units = "in", dpi = 144) ## ----example-66, eval=FALSE--------------------------------------------------- # library(ggchord) # library(dplyr) # data(plasmid_example_pETDuet_1) # show_plasmid <- plasmid_example_pETDuet_1 |> # select(accver, length) # show_plasmid_features <- find_common_features(plasmid_example_pETDuet_1) # show_plasmid_primers <- find_primer_bindings(plasmid_example_pETDuet_1, set = "reference") # show_plasmid_all_sites <- find_restriction_sites(plasmid_example_pETDuet_1) # show_plasmid_sites <- filter_restriction_sites( # show_plasmid_all_sites, # set = "unique_6plus", parent_set = "commercial_all" # ) ## ----example-67, eval=FALSE--------------------------------------------------- # show_plasmid_position <- position_feature_stack( # spacing = 0.085, base_position = position_plasmid() # ) # p_show_plasmid <- ggchord(show_plasmid, validate = "none") + # geom_seq(seq_style = "double", colour = "#252525") + # geom_feature_plasmid(data = show_plasmid_features, position = show_plasmid_position) + # geom_feature_label_repel( # data = show_plasmid_features, position = show_plasmid_position, # external = TRUE, max_overlaps = 0, # label_fit = "shrink", label_min_size = 2.5 # ) + # geom_primer(data = show_plasmid_primers, colour = "#A020F0") + # geom_primer_label_repel( # data = show_plasmid_primers, colour = "#A020F0", max_overlaps = 0 # ) + # geom_restriction_site(data = show_plasmid_sites, colour = "#252525") + # geom_seq_center_label() + # coord_circular( # rotation = 90, fit = "manual", # xlim = c(-3.4, 4.6), ylim = c(-1.6, 2.4) # ) + # theme_ggchord_plasmid() + # theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank() # ) # p_show_plasmid ## ----example-68, eval=FALSE--------------------------------------------------- # ggsave("pETDuet-1-map.png", p_show_plasmid, width = 16.5, height = 8, units = "in", dpi = 144) ## ----example-69, eval=FALSE--------------------------------------------------- # library(ggchord) # library(dplyr) # data(seq_data_example) # data(ribbon_data_example) # data(gene_data_example) # seqs <- seq_data_example |> # select(accver, length) # links <- ribbon_data_example |> # select(qaccver, saccver, qstart, qend, sstart, send) # # The fixture preserves biological strands separately from its old display strands. # genes <- gene_data_example |> # transmute(accver, start, end, strand = source_strand) ## ----example-70, eval=FALSE--------------------------------------------------- # seq_files <- read_fasta_lengths(files = "examples/fasta/*.fna", source_file = TRUE) # link_files <- read_blast(files = "examples/blastn-standard/*.tsv", source_file = TRUE) # feature_files <- read_gff3(files = "examples/gff3/*.gff3", unstranded = "keep", source_file = TRUE) ## ----example-71, eval=FALSE--------------------------------------------------- # validation <- validate_ggchord_data(seqs, links, genes) # cleaned <- clean_ggchord_data(seqs, links, genes) # prepared <- prepare_ggchord_data(seqs, links, features = list(genes = genes)) ## ----example-72, eval=FALSE--------------------------------------------------- # # Deliberately corrupt copies to demonstrate diagnostics; original rows remain intact. # problem_links <- ribbon_data_example |> # slice_head(n = 2) |> # select(qaccver, saccver, qstart, qend, sstart, send, pident) |> # mutate( # qaccver = replace(qaccver, 2, "absent"), # pident = replace(pident, 1, 103) # ) # problem_report <- validate_ggchord_data(seqs, problem_links) # repaired <- clean_ggchord_data(seqs, problem_links, unknown_id = "drop", invalid_pident = "clip") ## ----example-73, eval=FALSE--------------------------------------------------- # # Each source record represents a separate complete phage genome. # contigs <- seqs |> # mutate(genome = accver) # genomes <- tibble(genome = seqs$accver) # hierarchy <- validate_ggchord_hierarchy(contigs, genomes, strict = TRUE) ## ----example-74, eval=FALSE--------------------------------------------------- # library(ggchord) # library(dplyr) # data(seq_data_example) # data(ribbon_data_example) # data(gene_data_example) # seqs <- seq_data_example |> # select(accver, length) # links <- ribbon_data_example |> # select(qaccver, saccver, qstart, qend, sstart, send) # # The fixture preserves biological strands separately from its old display strands. # genes <- gene_data_example |> # transmute(accver, start, end, strand = source_strand) # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } # # library(dplyr) # data(plasmid_example_pUC19) # plasmid_seq <- plasmid_example_pUC19 |> # select(accver, length) # feature_source <- find_common_features(plasmid_example_pUC19) # features <- feature_source |> # select(accver, start, end, strand, type, anno) ## ----example-75, eval=FALSE--------------------------------------------------- # renamed_sequences <- seqs |> # rename(chrom = accver, size = length) # sequence_input <- as_ggchord_seq(renamed_sequences, mapping = c(accver = "chrom", length = "size")) # renamed_features <- genes |> # rename(chrom = accver, left = start, right = end) # feature_input <- as_ggchord_feature(renamed_features, # mapping = c(accver = "chrom", start = "left", end = "right") # ) # link_input <- as_ggchord_link(links) ## ----example-76, eval=FALSE--------------------------------------------------- # # BED3 carries intervals only; the source annotations remain in features. # bed_path <- tempfile(fileext = ".bed") # bed_rows <- tibble(accver = features$accver, start0 = features$start - 1L, end = features$end) # write.table(bed_rows, bed_path, sep = "\t", quote = FALSE, row.names = FALSE, col.names = FALSE) # bed_features <- read_bed(bed_path) # unlink(bed_path) # # Unchanged PAF rows from the public EMALE tutorial, not a simulated alignment. # paf_source <- read_paf("examples/format-inputs/emales-pair.paf") # paf_selected <- paf_source |> # filter(!is.na(mapq), mapq >= 30) # paf_links <- paf_selected |> # select(qaccver, saccver, qstart, qend, sstart, send) # paf_seqs <- bind_rows( # transmute(paf_source, accver = qaccver, length = qlen), # transmute(paf_source, accver = saccver, length = slen) # ) |> # distinct() # p_paf <- ggchord(seq_data = paf_seqs, ribbon_data = paf_links) + # geom_link_ribbon(fill = "#6BAFA6", alpha = .5) + # geom_seq(colour = "#233E4A") + # geom_seq_label(size = 4) + # coord_collinear() + # doc_axis() + # doc_theme # p_paf ## ----example-77, eval=FALSE--------------------------------------------------- # # The original NCBI L09137.2 record identifies pUC19c; keep that identity. # genbank_result <- NULL # if (nzchar(Sys.which("python3")) && identical(suppressWarnings(system2("python3", # c("-c", shQuote("from Bio import SeqIO")), # stdout = FALSE, stderr = FALSE # )), 0L)) { # genbank_result <- read_genbank("examples/format-inputs/pUC19c-L09137.2.gb") # } ## ----example-78, eval=FALSE--------------------------------------------------- # range_inputs <- NULL # if (requireNamespace("GenomicRanges", quietly = TRUE)) { # unknown <- features |> # filter(strand == ".") |> # slice_head(n = 1) # ranges <- GenomicRanges::GRanges( # seqnames = unknown$accver, # ranges = IRanges::IRanges(unknown$start, unknown$end), strand = "*", # seqlengths = setNames(plasmid_seq$length, plasmid_seq$accver) # ) # range_inputs <- list(seq_data = as_ggchord_seq(ranges), features = as_ggchord_feature(ranges)) # } ## ----example-79, eval=FALSE--------------------------------------------------- # interaction_inputs <- NULL # if (requireNamespace("GenomicRanges", quietly = TRUE) && # requireNamespace("InteractionSet", quietly = TRUE)) { # source_links <- ribbon_data_example |> # slice_head(n = 2) # lengths <- setNames(seq_data_example$length, seq_data_example$accver) # query <- GenomicRanges::GRanges( # seqnames = source_links$qaccver, # ranges = IRanges::IRanges( # pmin(source_links$qstart, source_links$qend), # pmax(source_links$qstart, source_links$qend) # ), # strand = ifelse(source_links$qstart <= source_links$qend, "+", "-"), # seqlengths = lengths # ) # subject <- GenomicRanges::GRanges( # seqnames = source_links$saccver, # ranges = IRanges::IRanges( # pmin(source_links$sstart, source_links$send), # pmax(source_links$sstart, source_links$send) # ), # strand = ifelse(source_links$sstart <= source_links$send, "+", "-"), # seqlengths = lengths # ) # interactions <- InteractionSet::GInteractions(query, subject) # interaction_inputs <- list( # seq_data = as_ggchord_seq(interactions), # ribbon_data = as_ggchord_link(interactions) # ) # } ## ----example-80, eval=FALSE--------------------------------------------------- # library(ggchord) # library(dplyr) ## ----example-81, eval=FALSE--------------------------------------------------- # blast_source <- read_blast("examples/blastn-standard/MT108731.1__MT118296.1.tsv") # blast_links <- blast_source |> # select(qaccver, saccver, qstart, qend, sstart, send) # # A six-column table is sufficient for drawing; extra fields stay in blast_source. # custom_path <- tempfile(fileext = ".tsv") # write.table(blast_links, custom_path, sep = "\t", quote = FALSE, row.names = FALSE, col.names = FALSE) # custom_source <- read_blast(custom_path, format = "custom", col_names = names(blast_links)) # custom_links <- custom_source |> # select(all_of(names(blast_links))) # unlink(custom_path) ## ----example-82, eval=FALSE--------------------------------------------------- # protein_source <- read_blast("examples/format-inputs/emales-mavirus-blastp.tsv", source_file = TRUE) # emale_source <- read_gff3("examples/format-inputs/emales-BVI_023A.gff3", source_file = TRUE) # # The source declares complete Prodigal CDS and translation table 11. # emale_coding <- emale_source |> # filter(grepl("partial=00", attributes)) |> # select(feature_id, accver, start, end, strand, phase) # emale_cds <- build_ggchord_cds(emale_coding, genetic_code = 11) # # Keep strong HSPs whose query protein has a known complete CDS. # # Gaps are retained: projecting one query span does not assert an ungapped pair. # protein_kept <- protein_source |> # filter(qaccver %in% emale_coding$feature_id, evalue <= 1e-10, length >= 50) |> # arrange(desc(bitscore)) # selected_hsp <- protein_kept |> # slice_head(n = 1) # query_interval <- tibble( # feature_id = "query_HSP", parent_id = selected_hsp$qaccver, # start = selected_hsp$qstart, end = selected_hsp$qend # ) # projected_query <- project_ggchord_features(query_interval, emale_cds) # # Subject projection additionally requires its own CDS map; it is not guessed from an accession. ## ----example-83, eval=FALSE--------------------------------------------------- # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } # # library(dplyr) # data(plasmid_example_pUC19) # data(plasmid_example_pETDuet_1) # homolog_source <- bind_rows(plasmid_example_pUC19, plasmid_example_pETDuet_1) # homolog_seq <- homolog_source |> # select(accver, length) # cluster_features <- homolog_source |> # find_common_features() |> # filter(type == "CDS", anno == "AmpR") |> # select(accver, start, end, strand, anno) |> # arrange(match(accver, homolog_source$accver)) |> # mutate(feature_id = paste(accver, start, end, sep = ":")) # # Membership is computed from exact strand-corrected CDS DNA equality, not gene names. # reverse_complement <- function(sequence) { # bases <- strsplit(sequence, "")[[1]] # chartr("ACGT", "TGCA", paste(rev(bases), collapse = "")) # } # coding_sequences <- cluster_features |> # left_join(select(homolog_source, accver, sequence), by = "accver") |> # rowwise() |> # mutate( # coding_dna = substr(sequence, start, end), # coding_dna = if (strand == "-") reverse_complement(coding_dna) else coding_dna # ) |> # ungroup() # cds_dna <- pull(coding_sequences, coding_dna) # # membership <- tibble( # feature_id = cluster_features$feature_id, # cluster_id = paste0("exact_CDS_", match(cds_dna, unique(cds_dna))) # ) # homologs <- link_ggchord_clusters(cluster_features, membership, # pairs = tibble(qaccver = homolog_seq$accver[1], saccver = homolog_seq$accver[2]) # ) # loci <- find_ggchord_loci(homolog_seq, cluster_features, # targets = cluster_features$feature_id[cluster_features$anno == "AmpR"][1], upstream = 100, downstream = 100 # ) # # # Compare equivalent local windows instead of shrinking one gene on full plasmids. # homolog_windows <- cluster_features |> # select(accver, start, end) # local_homologs <- focus_ggchord_data(homolog_seq, homologs$data, cluster_features, # loci = homolog_windows, expand = 200 # ) # p_homologs <- ggchord( # seq_data = local_homologs$seq_data, # ribbon_data = local_homologs$ribbon_data, gene_data = local_homologs$gene_data # ) + # geom_link_ribbon(fill = "#269B82", alpha = .65) + # geom_seq(colour = "#233E4A", show.legend = FALSE) + # geom_feature( # data = local_homologs$gene_data, feature_fill = "#F1D334", # position = position_strand(.10), width = .07 # ) + # geom_feature_label_repel(data = local_homologs$gene_data, size = 4, label_fit = "shrink") + # geom_seq_label(size = 4) + # geom_seq_break() + # coord_collinear() + # doc_axis("source") + # doc_theme # p_homologs ## ----example-84, eval=FALSE--------------------------------------------------- # # Run the exact-CDS comparison first; these are complete source-backed CDS. # coding <- cluster_features |> # mutate(phase = 0) # cds <- build_ggchord_cds(coding, genetic_code = 11) # # This selects residues for a display example; it does not claim a detected domain. # domains <- tibble(feature_id = "selected_residues", parent_id = coding$feature_id[1], start = 2, end = 20) # genomic <- project_ggchord_features(domains, cds) # genomic ## ----example-85, eval=FALSE--------------------------------------------------- # # Exact CDS DNA equality supplies evidence for these selected ungapped residues. # matching <- homologs$data |> # slice_head(n = 1) # paired_coding <- cluster_features |> # mutate(phase = 0) # paired_cds <- build_ggchord_cds(paired_coding) # protein_link <- tibble( # qfeature_id = matching$qfeature_id, sfeature_id = matching$sfeature_id, # qstart = 1, qend = 20, sstart = 1, send = 20 # ) # projected_links <- project_ggchord_links(protein_link, paired_cds, ungapped = TRUE) # projected_links ## ----example-86, eval=FALSE--------------------------------------------------- # library(ggchord) # library(dplyr) # data(seq_data_example) # data(ribbon_data_example) # data(gene_data_example) # seqs <- seq_data_example |> # select(accver, length) # links <- ribbon_data_example |> # select(qaccver, saccver, qstart, qend, sstart, send) # # The fixture preserves biological strands separately from its old display strands. # genes <- gene_data_example |> # transmute(accver, start, end, strand = source_strand) # # # Optional appearance choices; none of these columns or scales is required input. # library(ggchord) # seq_palette <- setNames( # c("#E64B35", "#233E4A", "#269B82", "#D59B17"), # c("MT108731.1", "MT118296.1", "OQ646790.1", "OR222515.1") # ) # gene_palette <- c("+" = "#CE1536", "-" = "#346B9D", "." = "#8291A6") # # library(ggchord) # feature_palette <- c( # CDS = "#993366", promoter = "#31849B", rep_origin = "#F1E51D", # primer_bind = "#A020F0", misc_feature = "#A6ACB3", protein_bind = "#31849B" # ) # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } # # library(dplyr) # data(plasmid_example_pUC19) # plasmid_seq <- plasmid_example_pUC19 |> # select(accver, length) # feature_source <- find_common_features(plasmid_example_pUC19) # features <- feature_source |> # select(accver, start, end, strand, type, anno) ## ----example-87, eval=FALSE--------------------------------------------------- # stack_features <- features |> # select(accver, start, end, strand, type, anno) # stack_position <- position_feature_stack(spacing = .12) # feature_shapes <- setNames(ifelse(names(feature_palette) == "CDS", "arrow", "block"), names(feature_palette)) # p_stack <- ggchord(seq_data = plasmid_seq) + # geom_seq(show.legend = FALSE) + # geom_feature(aes(feature_fill = type, feature_shape = type), # data = stack_features, # width = .08, position = stack_position # ) + # geom_feature_label_repel(data = stack_features, position = stack_position, size = 3.8, label_fit = "shrink", label_min_size = 2.8) + # scale_feature_shape_manual(values = feature_shapes) + # scale_feature_fill_manual(values = feature_palette) + # coord_collinear() + # doc_axis() + # doc_theme # p_stack ## ----example-88, eval=FALSE--------------------------------------------------- # regions <- tibble( # accver = links$qaccver[1], start = min(links$qstart[1], links$qend[1]), # end = max(links$qstart[1], links$qend[1]), category = "alignment interval" # ) # p_decorated <- ggchord(seq_data = seqs, ribbon_data = links) + # geom_seq_region(aes(region_fill = category), data = regions, colour = NA) + # geom_link_ribbon(fill = "#6BAFA6", alpha = .4) + # geom_seq() + # geom_seq_label(size = 4) + # scale_region_fill_manual(values = c("alignment interval" = "#F4D35E")) + # scale_seq_colour_manual(values = seq_palette) + # geom_seq_scale(length = 10000) + # coord_collinear() + # doc_axis() + # doc_theme # p_decorated ## ----example-89, eval=FALSE--------------------------------------------------- # library(ggchord) # library(dplyr) # data(seq_data_example) # data(ribbon_data_example) # data(gene_data_example) # seqs <- seq_data_example |> # select(accver, length) # links <- ribbon_data_example |> # select(qaccver, saccver, qstart, qend, sstart, send) # # The fixture preserves biological strands separately from its old display strands. # genes <- gene_data_example |> # transmute(accver, start, end, strand = source_strand) # # # Optional appearance choices; none of these columns or scales is required input. # library(ggchord) # seq_palette <- setNames( # c("#E64B35", "#233E4A", "#269B82", "#D59B17"), # c("MT108731.1", "MT118296.1", "OQ646790.1", "OR222515.1") # ) # gene_palette <- c("+" = "#CE1536", "-" = "#346B9D", "." = "#8291A6") # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } ## ----example-90, eval=FALSE--------------------------------------------------- # label_data <- gene_data_example |> # filter(anno != "hypothetical protein") |> # transmute(accver, start, end, strand = source_strand, anno) # p_labels <- ggchord(seq_data = seqs, gene_data = genes) + # geom_seq() + # geom_gene() + # geom_gene_label_repel(data = label_data, label_layout = "auto", size = 4, label_wrap = 22) + # geom_seq_label(size = 4) + # coord_collinear() + # doc_axis() + # scale_seq_colour_manual(values = seq_palette) + # scale_gene_fill_manual(values = gene_palette) + # doc_theme # p_labels ## ----example-91, eval=FALSE--------------------------------------------------- # # Compare layouts with at most the first and last informative label per sequence. # comparison_labels <- label_data |> # group_by(accver) |> # filter(row_number() %in% c(1L, n())) |> # ungroup() # label_plots <- lapply(c("auto", "radial", "arc"), function(layout) { # ggchord(seq_data = seqs, gene_data = genes) + # geom_seq() + # geom_gene() + # geom_gene_label_repel(data = comparison_labels, label_layout = layout, size = 4, label_wrap = 22) + # geom_seq_label(size = 4) + # scale_seq_colour_manual(values = seq_palette, name = "Sequence", # guide = guide_ggchord_legend(nrow = 2, theme = doc_legend_theme)) + # scale_gene_fill_manual(values = gene_palette, breaks = c("+", "-"), name = "Strand", guide = strand_guide) + # coord_chord(xlim = c(-1.8, 1.8), ylim = c(-1.6, 1.6)) + # theme_ggchord_publication(axis = element_blank()) + # doc_legend_theme + # theme( # legend.position = "bottom", # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank() # ) + # labs(title = layout) # }) # names(label_plots) <- c("auto", "radial", "arc") # p_manual_labels <- ggchord(seq_data = seqs, gene_data = genes) + # geom_seq() + # geom_gene() + # geom_gene_label(data = label_data, label_rotation = 30) + # coord_collinear() # label_plots$auto ## ----example-92, eval=FALSE--------------------------------------------------- # library(ggchord) # feature_palette <- c( # CDS = "#993366", promoter = "#31849B", rep_origin = "#F1E51D", # primer_bind = "#A020F0", misc_feature = "#A6ACB3", protein_bind = "#31849B" # ) # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } # # library(dplyr) # data(plasmid_example_pUC19) # plasmid_seq <- plasmid_example_pUC19 |> # select(accver, length) # feature_source <- find_common_features(plasmid_example_pUC19) # features <- feature_source |> # select(accver, start, end, strand, type, anno) ## ----example-93, eval=FALSE--------------------------------------------------- # library(dplyr) # data(plasmid_example_pUC19) # track_seq <- plasmid_example_pUC19 |> # select(accver, length) # starts <- seq.int(1L, track_seq$length, by = 100L) # gc_fraction <- function(sequence) { # bases <- strsplit(sequence, "")[[1]] # mean(bases %in% c("G", "C")) # } # bins <- tibble( # accver = track_seq$accver, # start = starts, # end = pmin(starts + 99L, track_seq$length) # ) |> # rowwise() |> # mutate(gc = gc_fraction(substr(plasmid_example_pUC19$sequence, start, end))) |> # ungroup() # head(bins) ## ----example-94, eval=FALSE--------------------------------------------------- # p_tracks <- ggchord(seq_data = track_seq) + # geom_seq() + # geom_seq_center_label() + # geom_track_bar(aes(track_value = gc, track_fill = gc), # data = bins, # position = position_track("gc", width = unit(10, "mm")) # ) + # scale_track_value_continuous("gc", limits = c(0, 1), unit = "GC fraction") + # scale_track_fill_gradientn("gc", colours = c("#E4F3EA", "#65B79B", "#27775C")) + # coord_circular() + # doc_theme # p_tracks ## ----example-95, eval=FALSE--------------------------------------------------- # line_bins <- bins |> # transmute( # accver, # pos = floor((start + end) / 2), # gc_mean = as.numeric(stats::filter(gc, rep(1 / 3, 3), sides = 2)) # ) |> # filter(is.finite(gc_mean)) # categories <- features |> # select(accver, start, end, type) # p_multitrack <- ggchord(seq_data = track_seq) + # geom_seq() + # geom_seq_center_label() + # geom_track_bar(aes(track_value = gc, track_fill = gc), # data = bins, # position = position_track("gc", width = unit(10, "mm")) # ) + # geom_track_line(aes(track_value = gc_mean), # data = line_bins, # position = position_track("gc", width = unit(10, "mm")), colour = "#233E4A" # ) + # scale_track_value_continuous("gc", limits = c(0, 1), unit = "GC fraction") + # scale_track_fill_gradientn("gc", colours = c("#E4F3EA", "#65B79B", "#27775C")) + # geom_track_tile(aes(track_fill = type), # data = categories, # position = position_track("features", order = 2, width = unit(2, "mm")) # ) + # scale_track_fill_manual("features", values = feature_palette) + # coord_circular() + # doc_theme + # theme(legend.box = "vertical", legend.box.just = "left") # focused_bins <- focus_ggchord_track(bins, track_seq, # loci = tibble(accver = track_seq$accver, start = 500, end = 1500), type = "interval" # ) # p_multitrack ## ----example-96, eval=FALSE--------------------------------------------------- # ggsave("pUC19-tracks.png", p_multitrack, width = 9, height = 8, units = "in", dpi = 144) ## ----example-97, eval=FALSE--------------------------------------------------- # library(ggchord) # library(dplyr) # data(seq_data_example) # data(ribbon_data_example) # data(gene_data_example) # seqs <- seq_data_example |> # select(accver, length) # links <- ribbon_data_example |> # select(qaccver, saccver, qstart, qend, sstart, send) # # The fixture preserves biological strands separately from its old display strands. # genes <- gene_data_example |> # transmute(accver, start, end, strand = source_strand) # # # Optional appearance choices; none of these columns or scales is required input. # library(ggchord) # seq_palette <- setNames( # c("#E64B35", "#233E4A", "#269B82", "#D59B17"), # c("MT108731.1", "MT118296.1", "OQ646790.1", "OR222515.1") # ) # gene_palette <- c("+" = "#CE1536", "-" = "#346B9D", "." = "#8291A6") # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } ## ----example-98, eval=FALSE--------------------------------------------------- # p_styled <- ggchord(seq_data = seqs, ribbon_data = links, gene_data = genes) + # geom_seq() + # geom_link_ribbon(aes(ribbon_fill = qaccver), alpha = .5) + # geom_gene() + # geom_seq_label(size = 4) + # coord_collinear() + # doc_axis() + # scale_seq_colour_manual(values = seq_palette) + # scale_gene_fill_manual(values = gene_palette) + # scale_ribbon_fill_manual(values = seq_palette) + # guides(gene_fill = guide_ggchord_legend(nrow = 1, theme = doc_legend_theme)) + # doc_theme # theme_variants <- list(minimal = theme_ggchord_minimal(), dark = theme_ggchord_dark(), plasmid = theme_ggchord_plasmid()) # p_styled ## ----example-99, eval=FALSE--------------------------------------------------- # library(ggchord) # library(dplyr) # data(seq_data_example) # data(ribbon_data_example) # data(gene_data_example) # seqs <- seq_data_example |> # select(accver, length) # links <- ribbon_data_example |> # select(qaccver, saccver, qstart, qend, sstart, send) # # The fixture preserves biological strands separately from its old display strands. # genes <- gene_data_example |> # transmute(accver, start, end, strand = source_strand) # # # Optional appearance choices; none of these columns or scales is required input. # library(ggchord) # seq_palette <- setNames( # c("#E64B35", "#233E4A", "#269B82", "#D59B17"), # c("MT108731.1", "MT118296.1", "OQ646790.1", "OR222515.1") # ) # gene_palette <- c("+" = "#CE1536", "-" = "#346B9D", "." = "#8291A6") # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } ## ----example-100, eval=FALSE-------------------------------------------------- # # Sequence IDs remain stable when the data is reordered; positional lists do not. # radii <- setNames(c(1.1, 1, 1, 1), seqs$accver) # widths <- c("+" = .045, "-" = .065) # p_parameters <- ggchord(seq_data = seqs, gene_data = genes) + # geom_seq(seq_radius = radii) + # geom_gene( # width = widths, # position = position_strand(c("+" = .08, "-" = -.1)) # ) + # geom_seq_label(size = 4) + # scale_seq_colour_manual(values = seq_palette) + # scale_gene_fill_manual(values = gene_palette) + # coord_chord() + # doc_theme # p_parameters ## ----example-101, eval=FALSE-------------------------------------------------- # library(ggchord) # library(dplyr) # data(seq_data_example) # data(ribbon_data_example) # data(gene_data_example) # seqs <- seq_data_example |> # select(accver, length) # links <- ribbon_data_example |> # select(qaccver, saccver, qstart, qend, sstart, send) # # The fixture preserves biological strands separately from its old display strands. # genes <- gene_data_example |> # transmute(accver, start, end, strand = source_strand) # # # Optional appearance choices; none of these columns or scales is required input. # library(ggchord) # seq_palette <- setNames( # c("#E64B35", "#233E4A", "#269B82", "#D59B17"), # c("MT108731.1", "MT118296.1", "OQ646790.1", "OR222515.1") # ) # gene_palette <- c("+" = "#CE1536", "-" = "#346B9D", "." = "#8291A6") # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } ## ----example-102, eval=FALSE-------------------------------------------------- # identity_links <- ribbon_data_example |> # select(qaccver, saccver, qstart, qend, sstart, send, pident) # p_mapped <- ggchord(seq_data = seqs, ribbon_data = identity_links) + # geom_seq() + # geom_link_ribbon(aes(ribbon_fill = pident), alpha = .6) + # scale_ribbon_fill_gradientn( # colours = identity_palette, values = identity_values, # limits = c(0, 100), breaks = identity_breaks, name = "Identity (%)", guide = identity_guide # ) + # scale_seq_colour_manual(values = seq_palette) + # coord_chord(rotation = 20) + # labs(title = "Sequence comparison") + # doc_theme # p_mapped ## ----example-103, eval=FALSE-------------------------------------------------- # library(ggchord) # library(dplyr) # data(seq_data_example) # data(ribbon_data_example) # data(gene_data_example) # seqs <- seq_data_example |> # select(accver, length) # links <- ribbon_data_example |> # select(qaccver, saccver, qstart, qend, sstart, send) # # The fixture preserves biological strands separately from its old display strands. # genes <- gene_data_example |> # transmute(accver, start, end, strand = source_strand) # # # Optional appearance choices; none of these columns or scales is required input. # library(ggchord) # seq_palette <- setNames( # c("#E64B35", "#233E4A", "#269B82", "#D59B17"), # c("MT108731.1", "MT118296.1", "OQ646790.1", "OR222515.1") # ) # gene_palette <- c("+" = "#CE1536", "-" = "#346B9D", "." = "#8291A6") # # library(ggchord) # doc_background <- "#FFFCF5" # identity_palette <- c("#303B71", "#36537D", "#269B82", "#A5D52C", "#ECEB31") # identity_values <- c(0, 0.5, 0.8, 0.9, 1) # identity_breaks <- c(0, 50, 80, 90, 100) # doc_legend_theme <- theme( # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.margin = margin(l = 4, r = 4, unit = "mm") # ) # sequence_guide <- guide_ggchord_legend(nrow = 1, theme = doc_legend_theme) # strand_guide <- guide_ggchord_legend( # nrow = 1, override.aes = list(strand = c("+", "-")), theme = doc_legend_theme # ) # identity_guide <- guide_ggchord_colourbar( # position = "left", direction = "vertical", # theme = theme( # legend.key.height = unit(75, "mm"), legend.key.width = unit(5, "mm"), # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), legend.margin = margin(r = 7, unit = "mm") # ) # ) # doc_theme <- theme_ggchord_publication( # axis.text = element_text(size = 11), # axis.text.check.overlap = TRUE # ) + # theme( # legend.position = "bottom", # legend.title = element_text(size = 13, face = "bold"), # legend.text = element_text(size = 11), # legend.key.width = unit(9, "mm"), legend.key.height = unit(5, "mm"), # legend.spacing.x = unit(9, "mm"), legend.spacing.y = unit(7, "mm"), # legend.box.spacing = unit(6, "mm"), # plot.background = element_rect(fill = doc_background, colour = NA), # panel.background = element_blank(), # plot.title = element_text(face = "bold", size = 18), # plot.subtitle = element_text(size = 10, colour = "#5C625A"), # plot.margin = margin(14, 18, 14, 18) # ) # doc_axis <- function(coordinates = "relative") { # scale_seq_position_continuous( # coordinates = coordinates, n.breaks = 3, minor_breaks = NULL, # labels = function(x) format(x, scientific = FALSE, big.mark = ",", trim = TRUE) # ) # } # # anchors <- genes |> # group_by(accver) |> # slice_max(end - start, n = 1, with_ties = FALSE) |> # ungroup() |> # select(accver, start, end) # p_collinear <- ggchord(seq_data = seqs, ribbon_data = links, gene_data = genes) + # geom_link_ribbon(aes(ribbon_fill = qaccver), alpha = .7) + # scale_ribbon_fill_manual(values = seq_palette, guide = "none") + # geom_seq() + # geom_gene(position = position_strand(.1)) + # geom_seq_label(size = 4) + # coord_collinear(anchor = anchors, align = "center", orientation = "auto") + # doc_axis() + # scale_seq_colour_manual(values = seq_palette) + # scale_gene_fill_manual(values = gene_palette) + # doc_theme # p_collinear ## ----example-104, eval=FALSE-------------------------------------------------- # ggsave("comparison.png", p_collinear, width = 8, height = 5, units = "in", dpi = 144) ## ----example-105, eval=FALSE-------------------------------------------------- # view_ggchord(p_collinear, width = 8, height = 5, units = "in") # diagnosis <- diagnose_ggchord(p_collinear, width = 8, height = 5, units = "in") ## ----example-106, eval=FALSE-------------------------------------------------- # geometry <- export_ggchord_layout(p_collinear, include = "seq") # rows <- geometry$seq |> # distinct(layer_id, component, source_row) ## ----example-107, eval=FALSE-------------------------------------------------- # full_layout <- get_chord_layout(p_collinear) # layer_tables <- export_ggchord_layout(p_collinear, include = c("seq", "gene"), original_data = TRUE) # layout_choice <- optimize_ggchord_layout(seqs, links, optimize = c("order", "orientation"), weight = "count") ## ----example-108, eval=FALSE-------------------------------------------------- # library(ggchord) # library(dplyr) # data(seq_data_example) # data(ribbon_data_example) # data(gene_data_example) # seqs <- seq_data_example |> # select(accver, length) # links <- ribbon_data_example |> # select(qaccver, saccver, qstart, qend, sstart, send) # # The fixture preserves biological strands separately from its old display strands. # genes <- gene_data_example |> # transmute(accver, start, end, strand = source_strand) # # library(dplyr) # data(plasmid_example_pUC19) # plasmid_seq <- plasmid_example_pUC19 |> # select(accver, length) # feature_source <- find_common_features(plasmid_example_pUC19) # features <- feature_source |> # select(accver, start, end, strand, type, anno) ## ----example-109, eval=FALSE-------------------------------------------------- # p_migrated <- ggchord(plasmid_seq) + # geom_seq() + # geom_feature(data = features, shape = "block", width = .05) + # geom_feature_label(data = features, label_rotation = 0) + # coord_circular() # p_offsets <- ggchord(seqs, gene_data = genes) + # geom_seq() + # geom_gene(position = position_strand(c("+" = .1, "-" = -.2)))