细胞色素 c：序列比对与系统发育分析

Author

Guangchuang Yu
School of Basic Medical Sciences, Southern Medical University

Published

October 6, 2026

数据文件包含 8 个物种的细胞色素 c 蛋白序列。完成以下分析：

Human–Tuna 的 identity、score 和比对长度；
局部比对结果与多序列比对；
保守列、距离矩阵和 NJ 树；
bootstrap 支持度与 WAG 模型 ML 树比较。

数据文件

FASTA 标题保留 NCBI、UniProt 或 PDB accession。

环境准备

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
install.packages(c("ape", "phangorn", "ggplot2", "dplyr", "yulab.utils"))
BiocManager::install(c("Biostrings", "pwalign", "DECIPHER", "ggtree", "tidytree", "treeio"))

主线需要 Biostrings、pwalign、DECIPHER、ape、phangorn 和 ggtree。muscle、rentrez 和网络下载不是必需依赖；目录中的 FASTA 和预计算 MSA 支持离线运行。

1. 读取真实序列

cytc_file <- "cytc_8species.fasta"
cytc <- readAAStringSet(cytc_file)
names(cytc) <- sub("_.*$", "", names(cytc))
cytc

AAStringSet object of length 8:
    width seq                                               names               
[1]   105 MGDVEKGKKIFIMKCSQCHTVEK...FVGIKKKEERADLIAYLKKATNE Human
[2]   105 MGDVEKGKKIFIMKCSQCHTVEK...FVGIKKKEERADLIAYLKKATNE Chimpanzee
[3]   105 MGDVEKGKKIFVQKCAQCHTVEK...FAGIKKKGERADLIAYLKKATNE Mouse
[4]   105 MGDIEKGKKIFVQKCSQCHTVEK...FAGIKKKSERVDLIAYLKDATSK Chicken
[5]   103 GDVAKGKKTFVQKCAQCHTVENG...IFAGIKKKGERQDLVAYLKSATS Tuna
[6]   108 MGVPAGDVEKGKKLFVQRCAQCH...IFAGLKKPNERGDLIAYLKSATK Fruitfly
[7]   108 AEFKAGSAKKGATLFKTRCLQCH...AFGGLKKEKDRNDLITYLKKATE Yeast
[8]   112 ASFSEAPPGNPDAGAKIFKTKCA...FPGLKKPQDRADLIAYLKKATSS Wheat

width(cytc)

[1] 105 105 105 105 103 108 108 112

先核对序列名称、长度和 FASTA 标题中的 accession。蛋白序列长度不同并不一定表示数据错误，后续的比对会显式处理插入和缺失。

2. 人和金枪鱼：80.6% 到底是什么量？

80.6% 是全长序列局部比对中的 identity：对齐位置上完全相同的氨基酸比例。它不是“同源性的百分比”，也不能单独证明同源。

pload(pwalign)

── loading the package: pwalign ──

human_tuna <- pairwiseAlignment(
  cytc[["Human"]], cytc[["Tuna"]],
  substitutionMatrix = "BLOSUM62",
  gapOpening = 10, gapExtension = 0.5,
  type = "global"
)

data.frame(
  score = score(human_tuna),
  identity_percent = round(pid(human_tuna), 1),
  aligned_length = nchar(as.character(alignedPattern(human_tuna)))
)

  score identity_percent aligned_length
1   436             80.6            105

cat("Human: ", as.character(alignedPattern(human_tuna)), "\\n", sep = "")

Human: MGDVEKGKKIFIMKCSQCHTVEKGGKHKTGPNLHGLFGRKTGQAPGYSYTAANKNKGIIWGEDTLMEYLENPKKYIPGTKMIFVGIKKKEERADLIAYLKKATNE\n

cat("Tuna : ", as.character(alignedSubject(human_tuna)), "\\n", sep = "")

Tuna : -GDVAKGKKTFVQKCAQCHTVENGGKHKVGPNLWGLFGRKTGQAEGYSYTDANKSKGIVWNNDTLMEYLENPKKYIPGTKMIFAGIKKKGERQDLVAYLKSATS-\n

练习 1：同时报告 score、identity_percent 和 aligned_length。只写“80.6%”而不写指标名称，别人无法判断你报告的是什么。

3. 把人序列作为查询：局部比对结果

local_hits <- lapply(setdiff(names(cytc), "Human"), function(sp) {
  hit <- pairwiseAlignment(
    cytc[["Human"]], cytc[[sp]],
    substitutionMatrix = "BLOSUM62",
    gapOpening = 10, gapExtension = 0.5,
    type = "local"
  )
  data.frame(
    species = sp,
    score = as.numeric(score(hit)),
    identity_percent = round(as.numeric(pid(hit)), 1),
    aligned_length = nchar(as.character(alignedPattern(hit)))
  )
}) |> bind_rows() |> arrange(desc(score))
local_hits

     species score identity_percent aligned_length
1 Chimpanzee   565            100.0            105
2      Mouse   522             91.4            105
3    Chicken   507             87.6            105
4       Tuna   457             80.6            103
5   Fruitfly   447             78.4            102
6      Wheat   385             67.0            103
7      Yeast   373             64.7            102

局部比对分数和 identity 随演化距离总体下降，但它们仍然是“线索”：打分矩阵、缺口惩罚、比对范围都会影响数值。下一步需要把多个物种放在同一个坐标系里比较。

4. 多序列比对与保守位点

cytc_msa.fasta 是由 DECIPHER::AlignSeqs() 对同一组序列计算得到的预计算结果。优先读取该文件；删除该文件后可运行备用分支重新计算。

msa_file <- "cytc_msa.fasta"

if (file.exists(msa_file)) {
  message("使用目录中附带的预计算 MSA：", msa_file)
  cytc_msa <- readAAStringSet(msa_file)
} else {
  pload(DECIPHER)
  cytc_msa <- AlignSeqs(cytc, verbose = FALSE)
  writeXStringSet(cytc_msa, filepath = msa_file)
}

使用目录中附带的预计算 MSA：cytc_msa.fasta

cytc_msa

AAStringSet object of length 8:
    width seq                                               names               
[1]   112 -------MGDVEKGKKIFIMKCS...FVGIKKKEERADLIAYLKKATNE Human
[2]   112 -------MGDVEKGKKIFIMKCS...FVGIKKKEERADLIAYLKKATNE Chimpanzee
[3]   112 -------MGDVEKGKKIFVQKCA...FAGIKKKGERADLIAYLKKATNE Mouse
[4]   112 -------MGDIEKGKKIFVQKCS...FAGIKKKSERVDLIAYLKDATSK Chicken
[5]   112 --------GDVAKGKKTFVQKCA...FAGIKKKGERQDLVAYLKSATS- Tuna
[6]   112 ---MGVPAGDVEKGKKLFVQRCA...FAGLKKPNERGDLIAYLKSATK- Fruitfly
[7]   112 ---AEFKAGSAKKGATLFKTRCL...FGGLKKEKDRNDLITYLKKATE- Yeast
[8]   112 ASFSEAPPGNPDAGAKIFKTKCA...FPGLKKPQDRADLIAYLKKATSS Wheat

msa_mat <- as.matrix(cytc_msa)
conserved <- apply(msa_mat, 2, function(col) {
  col <- col[col != "-"]
  length(col) > 0 && length(unique(col)) == 1
})

data.frame(
  aligned_columns = ncol(msa_mat),
  conserved_columns = sum(conserved),
  conserved_percent = round(100 * mean(conserved), 1)
)

  aligned_columns conserved_columns conserved_percent
1             112                54              48.2

练习 2：比较“完全一致的列”和“含有保守替换的列”。完全一致不是唯一重要证据，功能相关位点的保守替换同样需要结合结构和生物学解释。

5. 距离矩阵与 NJ 树

这里先用原始 p-distance，即对齐后不同氨基酸位点的比例。模型校正是后续扩展，不把“观察差异”直接当成真实演化距离。

pload(phangorn)

── loading the package: phangorn ──

phy_aln <- as.phyDat(cytc_msa, type = "AA")
p_dist <- dist.hamming(phy_aln, ratio = TRUE)
round(as.matrix(p_dist), 3)

           Human Chimpanzee Mouse Chicken  Tuna Fruitfly Yeast Wheat
Human      0.000      0.000 0.080   0.116 0.179    0.214 0.339 0.321
Chimpanzee 0.000      0.000 0.080   0.116 0.179    0.214 0.339 0.321
Mouse      0.080      0.080 0.000   0.071 0.134    0.143 0.339 0.321
Chicken    0.116      0.116 0.071   0.000 0.143    0.170 0.339 0.348
Tuna       0.179      0.179 0.134   0.143 0.000    0.170 0.348 0.357
Fruitfly   0.214      0.214 0.143   0.170 0.170    0.000 0.366 0.375
Yeast      0.339      0.339 0.339   0.339 0.348    0.366 0.000 0.366
Wheat      0.321      0.321 0.321   0.348 0.357    0.375 0.366 0.000

tree_nj <- NJ(p_dist)
tree_nj <- root(tree_nj, outgroup = "Yeast", resolve.root = TRUE)

ggtree(tree_nj) +
  geom_tiplab() +
  ggtitle("NJ tree from amino-acid p-distance")

练习 3：先不看物种名称，只看拓扑，写出人、黑猩猩、小鼠和鸡的大致关系；再打开名称检查你的判断。树的拓扑是关系假设，不是自动生成的生物学结论。

6. bootstrap：树对重采样是否稳定？

set.seed(20261005)
bs_trees <- bootstrap.phyDat(
  phy_aln,
  FUN = function(x) NJ(dist.hamming(x, ratio = TRUE)),
  bs = 50
)
tree_nj_bs <- plotBS(tree_nj, bs_trees, p = 50, type = "none")
tree_nj_bs$node.label

[1] 1.00 0.46 0.54 1.00 0.52 1.00   NA

ggtree(tree_nj_bs) +
  geom_tiplab() +
  geom_text2(aes(subset = !isTip, label = label), hjust = -0.2) +
  ggtitle("NJ tree with bootstrap support (50 replicates)")

bootstrap 重采样次数为 50。低支持度分支应标记为不确定，不应解释成确定的物种关系。

7. 扩展：用蛋白替换模型做 ML 树

fit0 <- pml(tree_nj, phy_aln, model = "WAG")
fit_ml <- optim.pml(
  fit0, optNni = TRUE, optBf = TRUE,
  optQ = FALSE, optInv = FALSE, optGamma = FALSE,
  rearrangement = "NNI", control = pml.control(trace = 0)
)

Warning: I unrooted the tree

tree_ml <- root(fit_ml$tree, outgroup = "Yeast", resolve.root = TRUE)

cat("ML log-likelihood:", round(fit_ml$logLik, 2), "\\n")

ML log-likelihood: -873.09 \n

ggtree(tree_ml) +
  geom_tiplab() +
  ggtitle("ML tree under the WAG amino-acid model")

练习 4：把 NJ 和 ML 的拓扑并排比较：哪些分支一致，哪些分支改变？改变本身不是“谁对谁错”的证明，而是提醒我们：距离定义、替换模型和数据支持度都会影响树。

8. 实验报告要求

提交一页以内的结果报告，至少包含：

Human–Tuna 的 score、identity 和比对长度；
8 物种 MSA 的长度与保守列比例；
NJ 树及至少一个内部节点的 bootstrap 支持度；
NJ 与 ML 的一处相同点和一处差异；
一段边界说明：为什么“identity 高”仍不能单独等同于“同源已证实”。

本实验的判断链是：identity 是观测量，MSA 提供可比坐标，树提出关系假设，bootstrap 和模型比较帮助我们表达不确定性。