系统发育分析与进化树构建¶

白话类比：
想象你有 100 个人的 DNA，想知道谁和谁是亲戚。
做法是：比对每个人 DNA 的差异位点，差异越少 = 关系越近。
最后画一张"家族树"，把所有人的亲缘关系展示出来。

系统发育做的是同样的事，只不过对象换成了：
- 不同物种（比如人、猩猩、猴子是怎么进化来的）
- 不同细菌（T2D 患者肠道里的菌谁和谁是近亲）
- 不同基因家族（同一个基因在不同物种里怎么演变的）

         ┌── 物种A
    ┌────┤                    ← 这个分叉点 = A和B的共同祖先
    │    └── 物种B
────┤                         ← 这个分叉点 = (A,B)和(C,D)的共同祖先
    │    ┌── 物种C
    └────┤
         └── 物种D

读法规则：
1. 分支越短 = 进化差异越小 = 关系越近
2. A和B 共享一个最近共同祖先（比 C、D 近）
3. 分支上的数字 = Bootstrap 支持值（越高越可信）

白话：
先算出所有序列两两之间的"距离分数"（差异有多大），
填成一张距离矩阵表，然后按"最近的先合并"的策略把树拼出来。
就像聚类分析——先把最像的两个聚一起，再往外扩。

白话：
给定一个进化模型（比如"A突变成T的概率是多少"），
尝试所有可能的树形结构，算每棵树"生成当前序列数据的概率"，
挑概率最高的那棵树作为最终结果。

类比：你看到一组考试分数，想猜班级的平均分和标准差。
ML 的做法是：试不同的平均分/标准差组合，
哪个组合最可能产生你看到的这组分数，就选哪个。

白话：
和 ML 类似，但加入了"先验信念"。
比如你觉得"哺乳动物大概 6500 万年前分化"，
贝叶斯法会把这个先验知识和序列数据结合起来，
用 MCMC 采样跑很多轮，最后给出一个"后验概率分布"。

类比：你猜明天下雨的概率。
ML 纯看天气数据，贝叶斯还会考虑"这个季节一般多雨"这个先验。

原始序列 → 多序列比对 → 修剪 → 模型选择 → 建树 → 评估 → 可视化
 (FASTA)   (MAFFT)   (trimAl)  (ModelFinder) (IQ-TREE) (Bootstrap) (iTOL/ggtree)

# ====== 从 NCBI 下载 16S 序列 ======
# 方法一：用 Entrez Direct（NCBI 命令行工具）
# 搜索 Bacteroides 属的 16S rRNA 基因序列，下载前 50 条
esearch -db nucleotide -query "Bacteroides[Organism] AND 16S ribosomal RNA[Title]" | \
  efetch -format fasta | \       # 下载 FASTA 格式
  head -500 > bacteroides_16s.fasta  # 保存前 ~50 条序列

# 方法二：从 SILVA 数据库下载高质量 16S 序列
# 访问 https://www.arb-silva.de/ 手动下载
# SILVA 的优势：序列经过质控，分类信息更准确

# 方法三：从自己的项目中提取（比如 T2D 项目的 Kraken2 分类结果）
# 提取特定物种的序列
seqkit grep -r -p "Bacteroides" all_contigs.fasta > target_sequences.fasta

# ====== 安装 MAFFT ======
conda install -c bioconda mafft  # conda 安装

# ====== 运行比对 ======
# --auto：自动选择最合适的算法（根据序列数量和长度）
# --thread 4：用 4 个 CPU 核心并行计算
mafft --auto --thread 4 input_sequences.fasta > aligned.fasta

# ====== 不同精度的比对策略 ======
# L-INS-i：最精确，适合 <200 条序列（推荐发文章用）
mafft --localpair --maxiterate 1000 --thread 4 input.fasta > aligned.fasta

# FFT-NS-2：快速模式，适合 >2000 条序列
mafft --retree 2 --thread 4 input.fasta > aligned.fasta

# 查看比对结果的基本统计
# 比对后的序列应该等长（因为插入了 gap）
grep -c ">" aligned.fasta   # 统计序列数量

# ====== 安装 trimAl ======
conda install -c bioconda trimal  # conda 安装

# ====== 自动修剪（推荐） ======
# -automated1：自动选择最优修剪策略（基于比对特征）
# 适合大多数系统发育分析场景
trimal -in aligned.fasta -out trimmed.fasta -automated1

# ====== 手动设置阈值 ======
# -gt 0.8：只保留 80% 以上序列都有碱基（非 gap）的列
# -cons 60：只保留保守度 >60% 的列
trimal -in aligned.fasta -out trimmed.fasta -gt 0.8 -cons 60

# ====== 用 Gblocks（另一个修剪工具） ======
# -t=d 表示 DNA 序列（-t=p 表示蛋白质）
# -b4=5 表示允许连续非保守列最多 5 个
Gblocks aligned.fasta -t=d -b4=5 -b5=h

# ====== 查看修剪前后的变化 ======
# 修剪前的列数
grep -v ">" aligned.fasta | head -1 | wc -c
# 修剪后的列数（应该变少了，去掉了噪音位点）
grep -v ">" trimmed.fasta | head -1 | wc -c

# ====== 什么是进化模型？ ======
# 进化模型描述碱基/氨基酸替换的概率。
# 比如 GTR 模型：A→T、A→C、A→G 各有不同的替换概率。
# 选错模型 → 建出来的树不准确。

# ====== 用 IQ-TREE 的 ModelFinder 自动选模型 ======
# -s：输入比对文件
# -m MF：只做模型选择，不建树
# 结果会输出到 trimmed.fasta.iqtree 文件中
iqtree2 -s trimmed.fasta -m MF -T AUTO

# 查看结果（找 Best-fit model 那行）
grep "Best-fit model" trimmed.fasta.iqtree
# 典型输出：Best-fit model: GTR+F+I+G4 chosen according to BIC

# ====== 常见 DNA 替换模型 ======
# JC（Jukes-Cantor）：最简单，所有碱基替换概率相同
# K2P（Kimura 2-parameter）：区分转换和颠换
# HKY：在 K2P 基础上允许碱基频率不等
# GTR（General Time Reversible）：最复杂，6 个独立替换率
# +I：允许一部分位点不变（invariant sites）
# +G4：位点间速率异质性（4 个 gamma 类别）

# ====== 安装 IQ-TREE 2 ======
conda install -c bioconda iqtree  # conda 安装

# ====== 标准建树命令（推荐，一步完成模型选择+建树+Bootstrap） ======
# -s：输入修剪后的比对文件
# -m MFP：先用 ModelFinder 选最优模型，然后用该模型建树
# -B 1000：做 1000 次 ultrafast bootstrap 评估分支可信度
# -alrt 1000：同时做 SH-aLRT 检验（双重验证）
# -T AUTO：自动选择最优线程数
iqtree2 -s trimmed.fasta -m MFP -B 1000 -alrt 1000 -T AUTO

# ====== 指定模型建树 ======
# 如果已经知道最佳模型（比如 Step 4 选出来的）
iqtree2 -s trimmed.fasta -m GTR+F+I+G4 -B 1000 -T AUTO

# ====== 输出文件说明 ======
# trimmed.fasta.treefile  → Newick 格式的树文件（用于可视化）
# trimmed.fasta.iqtree    → 详细分析报告（模型参数、日志似然值等）
# trimmed.fasta.log       → 运行日志
# trimmed.fasta.contree   → Consensus tree（含 bootstrap 值）
# trimmed.fasta.splits.nex → Bootstrap 分裂信息

# ====== 用 RAxML-NG 建树（替代方案，v2.0.1） ======
# --all：模型选择+建树+Bootstrap 一步完成
# --model：指定替换模型
# --bs-trees 100：做 100 次标准 bootstrap
# --threads auto：自动线程
raxml-ng --all --msa trimmed.fasta --model GTR+G --bs-trees 100 --threads auto

# IQ-TREE 的 ultrafast bootstrap（UFBoot）已经在 Step 5 中完成（-B 1000）
# 结果直接写在 .treefile 中的分支上

# ====== 解读 Bootstrap 值 ======
# UFBoot ≥ 95：强支持（分支可信）
# UFBoot 70-94：中等支持
# UFBoot < 70：弱支持（分支不可靠，可能需要更多数据）
#
# SH-aLRT ≥ 80：强支持
# SH-aLRT < 80：弱支持
#
# 论文中通常报告两个值：SH-aLRT/UFBoot，比如 95/99

# ====== 传统 Bootstrap（更慢但更保守） ======
# -b 100：做 100 次传统 bootstrap（比 UFBoot 慢 10-40 倍）
iqtree2 -s trimmed.fasta -m GTR+F+I+G4 -b 100 -T AUTO

# 1. 打开浏览器访问 https://itol.embl.de/
# 2. 上传 .treefile 文件
# 3. 在线编辑：
#    - 切换树的布局（矩形/圆形/无根）
#    - 给分支上色（按物种分类）
#    - 添加 bootstrap 值显示
#    - 添加热图、柱状图等注释
# 4. 导出为 PDF/SVG/PNG

# iTOL 注释文件格式示例（添加颜色标签）
# 创建 colors.txt，上传到 iTOL
# 格式：序列名称, 颜色, 标签
cat > itol_colors.txt << 'EOF'
DATASET_COLORSTRIP
SEPARATOR TAB
DATASET_LABEL   Phylum
COLOR   #ff0000

DATA
Bacteroides_fragilis    #e41a1c Bacteroidetes
Escherichia_coli    #377eb8 Proteobacteria
Lactobacillus_acidophilus   #4daf4a Firmicutes
EOF

# ====== 安装 ggtree ======
# BiocManager 安装（R 包管理器）
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("ggtree")  # 安装 ggtree

# ====== 基础绘图 ======
library(ggtree)      # 加载 ggtree 包
library(treeio)      # 加载树文件读取包

# 读取 IQ-TREE 输出的树文件
tree <- read.tree("trimmed.fasta.treefile")  # 读取 Newick 格式

# 基础进化树
ggtree(tree) +                               # 画树的主体
  geom_tiplab(size = 3) +                    # 添加叶节点标签（物种名）
  geom_nodelab(aes(label = label),           # 添加 Bootstrap 值
               size = 2, hjust = -0.2) +
  theme_tree2()                               # 添加坐标轴（显示进化距离）

# ====== 进阶：圆形树 + 分类着色 ======
# 准备分类信息
group_info <- data.frame(
  label = tree$tip.label,                     # 物种名
  phylum = c("Bacteroidetes", "Firmicutes",   # 对应的门分类
             "Proteobacteria", "Firmicutes")  # ...根据实际数据填
)

# 画圆形树并按门着色
ggtree(tree, layout = "circular") %<+% group_info +  # 圆形布局 + 合并分类信息
  geom_tiplab(aes(color = phylum), size = 2) +        # 标签按门着色
  geom_tippoint(aes(color = phylum), size = 1.5) +    # 叶节点圆点着色
  scale_color_brewer(palette = "Set1") +               # 使用 Set1 配色方案
  theme(legend.position = "right")                     # 图例放右边

# 保存图片
ggsave("phylogenetic_tree.pdf", width = 10, height = 8)  # 保存为 PDF

16S rRNA 基因是微生物分类的"金标准"基因：
1. 几乎所有细菌/古菌都有（通用性好）
2. 约 1500bp，包含 9 个可变区（V1-V9），提供分类信息
3. 保守区用于设计通用引物，可变区用于区分物种
4. SILVA、Greengenes、RDP 等数据库收录了几百万条 16S 序列

常见应用：
- 新物种发现：测 16S 序列，和已知物种建树，看它在树上的位置
- 菌群结构分析：把样本中所有 OTU/ASV 的代表序列建树
- 分类验证：Kraken2 分类结果不确定时，用 16S 系统发育确认

# ====== 从 QIIME2 的 ASV 代表序列开始 ======
# 假设你已经用 QIIME2 的 DADA2 得到了 rep-seqs.fasta

# Step 1: 加入参考序列（外群+已知物种）
# 从 SILVA 下载几条已知物种的 16S 序列作为参考
cat rep-seqs.fasta reference_16s.fasta > combined.fasta

# Step 2: 比对（16S 推荐用 MAFFT 的 L-INS-i）
mafft --localpair --maxiterate 1000 --thread 4 combined.fasta > aligned_16s.fasta

# Step 3: 修剪
trimal -in aligned_16s.fasta -out trimmed_16s.fasta -automated1

# Step 4: 建树
iqtree2 -s trimmed_16s.fasta -m MFP -B 1000 -alrt 1000 -T AUTO

# Step 5: 可视化（上传到 iTOL 或用 ggtree）

该 T2D 项目用 Kraken2 做了物种分类，得到了差异菌种列表。
系统发育分析可以在以下环节"锦上添花"：

1. 差异菌种的进化关系
   Kraken2 告诉你 T2D 患者中 Bacteroides fragilis 富集，
   你可以把所有差异菌种的 16S 序列建树，看看：
   - 富集的菌是不是集中在某个进化分支？（系统发育聚类 = phylogenetic clustering）
   - 这暗示"整个进化支系"与疾病相关，而不仅是单个物种

2. UniFrac 距离（系统发育多样性）
   普通的 Bray-Curtis 距离只看"物种有没有、多不多"
   UniFrac 还考虑了物种之间的进化关系：
   - 两个样本都有 Bacteroides，但一个是 B. fragilis，一个是 B. thetaiotaomicron
   - Bray-Curtis 会说它们完全不同（不同物种）
   - UniFrac 会说它们比较像（同一个属，进化距离近）

3. 功能推断的辅助
   进化关系近的细菌通常有相似的代谢功能（系统发育保守性）
   如果你发现 T2D 富集菌在树上聚在一起，
   可以推测它们可能共享某些影响血糖代谢的功能通路

面试话术模板：
"该 T2D 项目用 Kraken2 做物种分类，发现了若干差异菌种。
 为了进一步理解这些差异菌的进化关系，我提取了它们的 16S
 代表序列，用 MAFFT 比对、trimAl 修剪、IQ-TREE 建最大
 似然树。结果发现 T2D 富集菌集中在 Bacteroidetes 的
 某个分支上，提示该进化支系可能整体与疾病相关。
 同时在 Beta 多样性分析中，我对比了 Bray-Curtis 和
 weighted UniFrac 距离，后者考虑了物种间的系统发育关系，
 分组效果更明显。"

参考回答：
"三种方法的核心区别在于信息利用程度和统计框架：

NJ（距离法）：先把序列间的差异压缩成一个距离数字，再按最近邻
策略合并。速度最快但丢失了原始序列信息，适合快速预览。

ML（最大似然法）：在给定替换模型下，搜索能让数据出现概率最大
的那棵树。利用了完整序列信息，是发文章的首选。所用的 IQ-TREE
就是 ML 方法，它还内置了 ModelFinder 自动选模型和 ultrafast
bootstrap。

贝叶斯法：在 ML 基础上加入先验分布，用 MCMC 采样得到后验
概率。适合分子钟和分化时间估算，但计算量最大。

我在实际项目中优先用 IQ-TREE 的 ML 方法，需要分化时间
估算时才考虑 BEAST2 的贝叶斯方法。"

参考回答：
"Bootstrap 是通过对比对矩阵的列进行有放回抽样，重新建树
很多次（比如 1000 次），看某个分支在多少比例的重抽样树中
出现。本质上是评估'这个分支是不是因为少数几个位点才出现的'。

IQ-TREE 的 ultrafast bootstrap（UFBoot）阈值：≥95 为强支持。
传统 bootstrap 阈值：≥70 为可接受，≥90 为强支持。
贝叶斯后验概率阈值：≥0.95 为强支持。

注意 UFBoot 的数值通常比传统 bootstrap 高，所以阈值也更严。
论文中我会同时报告 SH-aLRT 和 UFBoot 两个值做交叉验证。"

参考回答：
"不能跳过，这是建树质量的基础。

多序列比对：把不同长度的序列对齐到同一坐标系，让每一列
代表同源位点。不比对的话，你比较的位置根本不对应，建出来
的树没有生物学意义。

修剪：比对后有些列质量很差（大量 gap 或者随机变异太多），
这些噪音位点会误导建树。trimAl 的 -automated1 会根据
比对质量自动选择修剪策略。

经验法则：修剪掉 10-30% 的列是正常的，如果修剪超过 50%
说明比对质量本身有问题，需要回头检查原始序列。"

参考回答：
"进化模型描述了碱基替换的概率规则。最简单的 JC 模型假设
所有替换概率相等，GTR 模型有 6 个独立替换率参数。

选模型的原因：不同基因/物种的进化方式不同。比如线粒体基因
转换（transition）远多于颠换（transversion），用 JC 模型
就低估了这种偏向，建出来的树分支长度不准确。

实践中我用 IQ-TREE 的 ModelFinder（-m MFP），它会根据
BIC 准则自动从几百个模型中选最优的。常见结果是 GTR+F+I+G4，
即一般时间可逆模型 + 经验碱基频率 + 不变位点 + gamma
速率异质性。"

参考回答：
"两者都是 Beta 多样性距离，区别在于是否考虑进化关系。

Bray-Curtis：纯粹比较物种丰度的差异，不管物种之间的关系。
如果样本 A 有 Bacteroides fragilis，样本 B 有 Bacteroides
thetaiotaomicron，它会当成完全不同的物种。

UniFrac：先建一棵系统发育树，然后比较两个样本在树上的
共享分支比例。上面的例子中，两个 Bacteroides 在树上很近，
所以 UniFrac 距离会比 Bray-Curtis 小。

weighted UniFrac 还考虑丰度权重，unweighted 只看有没有。
我在 T2D 项目中两种都做了：Bray-Curtis 看组成差异，
weighted UniFrac 看进化相关的组成差异，后者在 T2D 和
健康组之间的分离效果更显著。"

┌─────────────┐   ┌──────────┐   ┌──────────┐   ┌────────────┐   ┌──────────┐   ┌──────────┐
│ 序列收集     │ → │ 比对      │ → │ 修剪      │ → │ 模型选择    │ → │ 建树      │ → │ 可视化    │
│ NCBI/SILVA  │   │ MAFFT    │   │ trimAl   │   │ ModelFinder│   │ IQ-TREE  │   │ iTOL     │
│ seqkit      │   │ MUSCLE   │   │ Gblocks  │   │ (内置)      │   │ RAxML-NG │   │ ggtree   │
└─────────────┘   └──────────┘   └──────────┘   └────────────┘   └──────────┘   └──────────┘

# 一键完成（最常用）
iqtree2 -s input.fasta -m MFP -B 1000 -alrt 1000 -T AUTO

# 只做模型选择
iqtree2 -s input.fasta -m MF -T AUTO

# 指定模型建树
iqtree2 -s input.fasta -m GTR+F+I+G4 -B 1000 -T AUTO

# 蛋白质序列建树
iqtree2 -s protein.fasta -m MFP -B 1000 -T AUTO  # ModelFinder 自动选蛋白质模型

# 分区模型（不同基因/codon位置用不同模型）
iqtree2 -s concat.fasta -p partition.nex -m MFP+MERGE -B 1000 -T AUTO

# 重新运行（覆盖之前的检查点）
iqtree2 -s input.fasta -m MFP -B 1000 -T AUTO -redo