1.2 生成与导入比对数据

本页目录

从FASTA与GFF3开始
用BLAST产生自己的结果
已有其他格式时

先有序列与分析结果，再由ggchord读取；没有比对表时不应为了连线伪造同源关系。

从FASTA与GFF3开始

解压教学输入包，使当前目录包含 examples/。FASTA提供完整长度，GFF3提供区间和链；保留原始记录，在显示表中选择必要字段。

准备代码 · 在本页运行一次

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

show_alignments 保留32条原始HSP，show_links 明确筛为12条；show_genes 保留全部266条CDS，show_features 按长度至少500 bp、序列内注释顺序分四组各选一个，并补充有来源的负链首尾注释用于展示。筛选不是重跑比对，也不证明功能同源。

坐标排查时可用tidyr将四个端点列转成长表，原始连接表保持原样：

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

用BLAST产生自己的结果

以下命令是分析者在安装BLAST后对自己的文件执行的生产步骤，网站构建不运行它们。

代码blastn -query query.fna -subject subject.fna -outfmt 6 -out result.tsv

代码library(ggchord)
blast_source <- read_blast(file = "result.tsv")

无自定义字段列表的 -outfmt 6 为标准12列表，7另带注释；完全不写outfmt是其他格式。query/subject ID由读取器映射到qaccver/saccver。参见NCBI字段说明和核酸/蛋白比对教程。随站历史HSP只是格式转换，没有冒充本轮新BLAST结果。

已有其他格式时

来源	入口	注意
已读表格	as_ggchord_seq/feature/link	显式mapping映射角色，保留原表
BED、PAF	read_bed、read_paf	读取器只转换一次零基半开坐标
GenBank	read_genbank	需要可用Biopython；保留记录与分段来源
GRanges、Seqinfo	as_ggchord_feature、as_ggchord_seq	完整长度来自Seqinfo，未知链不转成正链
GInteractions	as_ggchord_link、as_ggchord_seq	第一锚点为query；仍需可靠完整长度

这些适配器对应用户issue #2提出的需求。当前代码已有入口；issue中的历史“开发版未经测试”不等于现阶段所有组合均已验证。本地缺失可选包的路径会明确跳过。格式专题给出可选对象的完整代码和限制。

继续：1.3 从骨架逐层画起。