导入核酸和蛋白比对

本页目录

本页准备
默认表格列就够用
文件布局与绘图最小列是两回事
BLASTp 可以读取，但坐标单位不同
来源与限制
完成后检查

使用 BLAST 标准表格输出，无需早期自定义 17 列。区分“文件有哪些列”“绘图需要哪些列”及核酸/蛋白单位。

本页准备

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

代码library(ggchord)
library(dplyr)

默认表格列就够用

不需要早期教程中的 17 列自定义字段。使用 BLAST 标准表格输出即可：

代码blastn -query query.fna -subject subject.fna -outfmt 6 -out result.tsv
blastp -query proteins.faa -subject reference.faa -outfmt 7 -out result.tsv

6 是无注释表格，7 加上注释头；不再列出自定义字段时，两者使用默认 12 列。完全省略 -outfmt 得到的是 pairwise 文本，当前表格读取器不解析该格式。默认列为 query ID、subject ID、pident、length、mismatch、gapopen、qstart、qend、 sstart、send、evalue、bitscore。读取器将两个 ID 的角色命名为 qaccver、saccver，保留字符串原值，并不要求它们一定是 accession version。 NCBI 格式说明。

文件布局与绘图最小列是两回事

文件里有 12 列，不表示绘图必须接收 12 列。运行本页准备代码加载包后，使用下面的真实记录：

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

blast_source 保留全部统计和来源信息；blast_links 只有六个端点角色。 pident 只有在按相似度配色/筛选时才需要，length 只有按比对长度筛选/加权时才需要。 qlen、slen 用于相应覆盖率计算，不能为了基础出图强制要求它们。无表头的自定义列须明确 col_names，读取器不猜测任意六列的含义。也可对任意已读入的表使用 as_ggchord_link(x, mapping = ...)，见格式转换。

仓库 examples/blastn-standard/ 的文件是历史原始 HSP 的前 12 列转换，保留行和值，不是本轮重新计算的 BLAST 结果。旧 .o7 文件保持原样，转换脚本和 hash 在来源清单中。

BLASTp 可以读取，但坐标单位不同

下面是真实 BLASTP 2.9.0+ 结果与对应查询 CDS 注释，来自公开 EMALE 教程。先保留完整输入，再选择有完整 CDS 映射、足够长度和较小 evalue 的查询 HSP。

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

项目	结果
原始 BLASTp 比对记录数	2948
具有完整注释的查询 CDS 数	24
明确筛选后保留的查询 HSP 数	3
所选查询起点（aa）	50
所选查询终点（aa）	289
投影后的查询范围行数	1

BLASTp 的 qstart/qend/sstart/send 是氨基酸坐标。这里将查询蛋白的 HSP 范围转换到已知 CDS 的基因组坐标，保留负链和 phase；有缺口的 HSP 仍可投影查询范围，但这只是范围，不能解释为逐个核苷酸的配对。原始弱命中、gapopen 和统计仍在原表中。

要得到两端的核苷酸连接，还须提供 subject 的 CDS 映射。不能从蛋白 accession、最大端点或同名注释猜出来。project_ggchord_links() 的 ungapped = TRUE 必须有证据，不能为了让默认 BLASTp 文件通过而随意打开。投影教程使用两条编码 DNA 完全一致的真实质粒 CDS，展示可证明无缺口的双端例子。

本版不直接绘制以 aa 为轴的蛋白骨架。blastx、tblastn 的两侧单位也可能不同，必须逐侧核对和转换；读取成功不代表可以直接用 bp 坐标绘图。

来源与限制

EMALE 原始输出来自 gggenomes 的公开数据与教程。本项目只复用数据记录，用 ggchord 自行绘图；没有复制上游图或绘图代码。下载、原始成员、行号范围、SHA256 和筛选规则见教学输入包的 sources.json。表格解析不重新运行比对，也不验证命中的生物学功能。

完成后检查

基础连接表只需六列。BLASTp 读取后仍是 aa 坐标，必须通过有证据的 CDS 映射转换后才能用于本版 bp 绘图。

继续：将蛋白坐标转换为基因组坐标 →