跳到正文
0.16.6 文档
全站目录 · 按任务查找

1.5 筛选与整理连接

本页目录

减少连接拥挤,同时保留每一步的来源:筛选、去重、合并和聚束有不同含义。

本页准备

在新的 R 会话中先运行下方代码,再按页面顺序继续。配色和排版是可选的显示设置,输入只选择需要的列。 涉及文件时,先解压教学输入,使工作目录包含 examples/。

准备代码 · 在本页运行一次
代码
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

从原始连接表开始

先用真实记录工作流创建 seq_real、link_real 和 real_filtered,保留原表及筛选报告。

代码
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
项目 结果
筛选后的比对数 8
重叠去重后数量 8
合并后的显示区间数 8
聚束输出行数 28

理解每步转换

筛选保留满足已声明阈值的观测;互惠重叠去重选择已有代表观测,不累加相似度; 合并把邻近兼容区间组合为显示外包范围;聚束按邻近 query/subject 位置形成汇总。 逐步检查 $report 和源行成员,特别注意冲突的额外元数据。 合并或聚束范围不是重新计算的核苷酸比对,估计长度保留标记,不凭空生成相似度。

密度是显式统计

原始比对记录的区间密度
原始比对记录的区间密度 · 点击放大查看

stat_ribbon_density() 映射计算后的密度,stat_ribbon_bundle() 在 build 时执行 明确聚合,普通区间连接几何层仍逐行绘制。计数、已报告长度和相似度加权有各自单位, 更不透明的连接不一定表示更强的独立同源证据。没有方向证据的点锚,不能因端点 相等就获得聚合方向。

选择区间区间连接或点连接

代码
# 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
从已有 HSP 区间计算中点的点连接
从已有 HSP 区间计算中点的点连接 · 点击放大查看

点位置由随包真实 HSP 的中点向下取整得到,颜色使用已报告的 pident。line 使用 qpos/spos,而非区间端点,箭头由 query 指向 subject。点连接默认避障为 none,区间连接为 uniform。 全局弦图旋转属于坐标系,直线模式选对与相邻行默认见三种模式。

点连接的颜色条明确指定 available_aes = "link_colour"。普通 ggplot2 颜色条和 包装器默认值不接受全部扩展角色;自定义图例时应明确选择对应美学角色。

完成后检查

先检查各步骤 report 和成员行,再读图。聚合带的范围不是新比对;点连接是已有 HSP 中点的摘要。