Project protein intervals to genomic coordinates

On this page

Prepare this page
Start with source-backed CDS
Segments, phase and partial regions
BLASTp query spans and two-ended links
Check your result

Project known protein intervals to genomic positions through real CDS, retaining phase, strand and segmentation explicitly.

Prepare this page

In a fresh R session, run the code below, then follow this page in order. Colours and layout are optional display choices; inputs select only the fields needed.

Preparation code — run once on this page

Codelibrary(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

Start with source-backed CDS

Run the identical-CDS comparison first. It uses real bundled plasmid annotations and DNA. Complete CDS use explicit phase 0 and translation table 11.

Code# 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

This example selects residues 2–20 for display; it does not claim a detected protein domain. Projection returns genomic intervals with strand/segment context. Actual domains require exact parent CDS IDs and amino-acid ranges; annotation names or maximum feature endpoints cannot supply missing evidence.

Segments, phase and partial regions

build_ggchord_cds() constructs a segment map; project_ggchord_features() then splits amino-acid intervals through it. Crossing a boundary can return several rows; multiplying coordinates by three is insufficient. Inspect policies/reports for partial CDS, phase and truncation. DNA matching, interval projection and actual translation/alignment are distinct operations.

BLASTp query spans and two-ended links

Code# 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

This ungapped case relies on verified identical coding DNA. Real BLASTp HSPs may contain gaps: query spans can still be projected, but paired links need the subject CDS and justified ungapped correspondence. See real BLASTp import/projection. This version does not reinterpret aa as bp or add native protein-coordinate plotting.

Check your result

Check parent CDS IDs, projection reports and returned segments. The selected residue range is a display example, not a detected domain.

Continue: Import nucleotide and protein alignments →