导入和转换数据

本页目录

本页准备
把已有表格转换成需要的列
BED 和 PAF：坐标只转换一次
GenBank：保留记录身份
可选的区间对象
完成后检查

将已有表格、BED、PAF 或 GenBank 转成 ggchord 输入，保留记录身份，并且只转换一次坐标约定。

本页准备

在新的 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)

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)

把已有表格转换成需要的列

这里使用随包的实际序列、连接和基因，只重命名列来演示角色映射；没有额外生成示意序列或比对。

代码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)

转换返回普通数据框，并记录 ggchord_conversion 信息。标识原样保留；缺少完整序列长度时应补充可靠来源，不能把最大注释端点当作全长。核酸与蛋白比对的最小列和单位见导入比对结果。

BED 和 PAF：坐标只转换一次

BED 例子把已知 pUC19 注释写为最小 BED3，再读回；PAF 直接读取公开教程中的真实 minimap2 输出，保留原始可选标签和值。显示筛选保留 mapq ≥30 的记录，绘图只接收六个端点角色。

代码# 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

从真实 minimap2 PAF 结果绘制的两条 EMALE 序列 · 点击放大查看

BED 的零起点半开区间由读取器转成一基包含坐标；不要再手工加一次 1。 BED3 没有链列，读回时链为未知；需要链方向时应使用含 strand 的原始格式。 PAF 的负相对方向表现为递减的 subject 端点，这是合法反向命中，不能一律排序去掉方向。比对质量和 pident 来自生产工具报告的字段，不在绘图时重新计算碱基比对。

GenBank：保留记录身份

代码# 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")
}

这个文件是 NCBI 的原始 L09137.2 记录，名称是 pUC19c；不要把它强行改成随包 pUC19。返回值包括 seq_data、features$annotation 与 report，分段位置及记录/注释来源保留。这里需要所选 Python 已安装 Biopython；函数不自动安装依赖。 unsupported = "report" 将不支持的位置记入报告，不把远端或复杂位置猜成区间。

可选的区间对象

代码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))
}

代码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)
  )
}

项目	结果
GenBank 可选读取器	依赖未安装，已跳过
可选区间对象适配器	依赖未安装，已跳过
可选GInteractions适配器	依赖未安装，已跳过

这里选择 pUC19 中原本未知链的注释；转换后的 * 仍表示未知。完整长度来自已知序列长度和 Seqinfo。GInteractions 可以转换成连接，但仍需明确两端序列及长度。缺少可选依赖时结果会标记跳过，不把源码存在当作运行通过。详见数据格式说明。当前不直接解析 SAM/BAM 的 CIGAR；需先由相应工具正确解析，再通过表格角色映射转换，不能仅靠改列名处理剪接。

完成后检查

先查看返回列与转换报告，再选绘图列。PAF 反向命中保留端点方向；GenBank 记录身份不能因为展示需要而改成另一条随包质粒。

继续：导入核酸和蛋白比对 →