607 微生物源追踪（SourceTracker2 / FEAST）¶

                  ┌─ 皮肤 (15%)
                  │
Sink样本（婴儿肠道）├─ 母乳 (40%)
                  │
                  ├─ 口腔 (20%)
                  │
                  └─ 未知 (25%)

SourceTracker2的核心思想：
1. 假设sink样本是多个source的混合
2. 用Dirichlet-Multinomial模型描述混合过程
3. 用Gibbs采样（MCMC）迭代估计每个source的贡献比例
4. 同时估计"unknown"（未知源）的比例

数学表达（简化版）：
Sink = p1×Source1 + p2×Source2 + ... + pk×Sourcek + pu×Unknown
其中 p1+p2+...+pk+pu = 1（所有比例之和为1）

# ========== 安装SourceTracker2 ==========
# SourceTracker2是Python包

# 方法1：pip安装
pip install sourcetracker  # pip安装

# 方法2：conda安装
conda install -c bioconda sourcetracker  # conda安装

# 方法3：从GitHub安装（最新版）
git clone https://github.com/biota/sourcetracker2.git  # 克隆仓库
cd sourcetracker2  # 进入目录
pip install .      # 安装

# ========== 准备输入数据 ==========
# SourceTracker2需要两个输入文件：
# 1. OTU/ASV丰度表（BIOM格式或TSV格式）
# 2. 样本元数据（标注source/sink分组）

import pandas as pd  # 导入pandas

# 准备元数据文件（mapping file）
# 关键列：SampleID, SourceSink, Env
metadata = pd.DataFrame({
    "SampleID": [
        "skin1", "skin2", "skin3",         # 皮肤source样本
        "oral1", "oral2", "oral3",          # 口腔source样本
        "gut1", "gut2", "gut3",             # 肠道source样本
        "infant1", "infant2", "infant3"     # 婴儿sink样本
    ],
    "SourceSink": [
        "source", "source", "source",      # source样本标记为"source"
        "source", "source", "source",
        "source", "source", "source",
        "sink", "sink", "sink"              # sink样本标记为"sink"
    ],
    "Env": [
        "Skin", "Skin", "Skin",            # source的环境类型
        "Oral", "Oral", "Oral",
        "Gut", "Gut", "Gut",
        "Infant", "Infant", "Infant"       # sink的环境标签（不参与计算）
    ]
})

# 保存元数据
metadata.to_csv("mapping_file.txt", sep="\t", index=False)  # TSV格式保存
print("元数据准备完毕")
print(metadata)

# ========== 命令行运行SourceTracker2 ==========

# 基本用法
sourcetracker2 gibbs \
    -i otu_table.biom \              # 输入：OTU/ASV丰度表（BIOM格式）
    -m mapping_file.txt \            # 输入：元数据文件
    -o st2_results/ \                # 输出目录
    --source_sink_column SourceSink \  # 元数据中标注source/sink的列名
    --source_column_value source \   # source的标签值
    --sink_column_value sink \       # sink的标签值
    --source_category_column Env \   # 元数据中标注环境类型的列名
    --jobs 8                         # 并行8个作业

# 高级参数
sourcetracker2 gibbs \
    -i otu_table.biom \
    -m mapping_file.txt \
    -o st2_results_advanced/ \
    --source_sink_column SourceSink \
    --source_column_value source \
    --sink_column_value sink \
    --source_category_column Env \
    --burnin 100 \                   # 烧入期100次（预热迭代，丢弃结果）
    --draws_per_restart 10 \         # 每次重启后抽取10个样本
    --restarts 5 \                   # 重启5次（提高结果稳定性）
    --alpha1 0.001 \                 # source的Dirichlet先验参数
    --alpha2 0.1 \                   # source环境间的Dirichlet先验
    --beta 10 \                      # sink的Dirichlet先验
    --jobs 16                        # 并行16个作业

# 输出文件解读：
# mixing_proportions.txt — 核心结果！每个sink中各source的比例
# mixing_proportions_stds.txt — 比例的标准差

# ========== R语言：SourceTracker2结果可视化 ==========
library(ggplot2)    # 可视化
library(tidyr)      # 数据整理
library(dplyr)      # 数据操作

# 读取SourceTracker2的混合比例结果
proportions <- read.delim("st2_results/mixing_proportions.txt",  # 读取结果
                           row.names = 1,                         # 第一列为行名
                           check.names = FALSE)                   # 不修改列名

# 转换为长格式（方便ggplot绑定）
props_long <- proportions %>%
  tibble::rownames_to_column("Sink_Sample") %>%  # 行名转列
  pivot_longer(cols = -Sink_Sample,               # 转为长格式
               names_to = "Source",                # source列名
               values_to = "Proportion")           # 比例值列名

# 堆叠柱状图
ggplot(props_long, aes(x = Sink_Sample, y = Proportion, fill = Source)) +  # 绑定数据
  geom_bar(stat = "identity", width = 0.8) +   # 堆叠柱状图
  scale_y_continuous(labels = scales::percent,  # Y轴显示百分比
                     limits = c(0, 1)) +
  labs(title = "Microbial Source Tracking Results",  # 标题
       x = "Sink Samples",                           # X轴
       y = "Source Proportion",                       # Y轴
       fill = "Source Environment") +                 # 图例
  theme_minimal() +                                   # 简约主题
  theme(axis.text.x = element_text(angle = 45, hjust = 1))  # X轴标签旋转

# 饼图（单个样本）
single_sample <- props_long %>% filter(Sink_Sample == "infant1")  # 选一个样本
ggplot(single_sample, aes(x = "", y = Proportion, fill = Source)) +
  geom_bar(stat = "identity", width = 1) +  # 柱状图
  coord_polar("y") +                         # 转为饼图
  labs(title = "Source Contributions to infant1") +
  theme_void()                                # 无背景

FEAST (Fast Expectation-maximization for Microbial Source Tracking)

核心区别于SourceTracker2：
- SourceTracker2用Gibbs采样（慢但理论性强）
- FEAST用EM算法（快约300倍！）

FEAST的EM算法流程：
E步：基于当前估计的混合比例，计算每个OTU来自各source的概率
M步：基于E步的概率，更新混合比例
重复E-M步直到收敛

FEAST优势：
1. 速度快300倍（大数据集时优势明显）
2. 支持多个sink同时分析
3. 自动估计unknown比例

# ========== R语言：安装和运行FEAST ==========

# 安装FEAST
# install.packages("devtools")  # 如果没有devtools
devtools::install_github("cozygene/FEAST")  # 从GitHub安装

library(FEAST)  # 加载FEAST

# 准备输入数据
# 1. OTU丰度矩阵（行=样本，列=OTU/ASV）
# 2. 元数据（必须包含Env, SourceSink, id列）

# 读取OTU表
otu_table <- read.delim("otu_table.tsv",   # 读取OTU表
                         row.names = 1,      # 第一列为行名
                         check.names = FALSE)

# 准备元数据（FEAST格式）
metadata <- data.frame(
  Env = c(rep("Skin", 3), rep("Oral", 3),     # 环境类型
          rep("Gut", 3), rep("Infant", 3)),
  SourceSink = c(rep("Source", 9),              # source/sink标签
                 rep("Sink", 3)),
  id = c(rep(1, 9),                             # sink id（所有source共享id=1表示都对应sink 1）
         1, 2, 3),                              # 每个sink有独立id
  row.names = rownames(otu_table)               # 行名=样本名
)

# 注意：FEAST的元数据格式与SourceTracker2不同！
# id列：每个sink样本有唯一id，对应的source样本共享同一个id

# 运行FEAST
feast_result <- FEAST(
  C = as.matrix(otu_table),     # OTU计数矩阵（必须是矩阵格式）
  metadata = metadata,          # 元数据
  different_sources_flag = 0,   # 0=所有sink共享同一组source
  dir_path = "feast_output/",   # 输出目录
  outfile = "feast_result"      # 输出文件前缀
)

# 查看结果
print(feast_result)
# 输出为比例矩阵：行=sink样本，列=source环境+Unknown

# ========== FEAST：每个sink有不同的source组合 ==========

# 准备元数据（每个sink独立的source组）
metadata_multi <- data.frame(
  Env = c("Skin", "Skin", "Oral", "Oral", "Gut", "Gut",  # sink1的source
          "Skin", "Skin", "Oral", "Oral", "Gut", "Gut",  # sink2的source
          "Infant", "Infant"),                             # sink样本
  SourceSink = c(rep("Source", 12), rep("Sink", 2)),
  id = c(rep(1, 6), rep(2, 6), 1, 2),  # id=1对应sink1，id=2对应sink2
  row.names = rownames(otu_table)
)

# 运行FEAST（different_sources_flag=1表示每个sink有独立的source组）
feast_multi <- FEAST(
  C = as.matrix(otu_table),
  metadata = metadata_multi,
  different_sources_flag = 1,     # 每个sink有不同的source组
  dir_path = "feast_multi_output/",
  outfile = "feast_multi"
)

SourceID-NMF：基于非负矩阵分解的源追踪

核心思想：
将OTU丰度矩阵分解为 W（源特征矩阵）× H（混合比例矩阵）
通过NMF的非负约束，自然满足比例非负的要求

优势：
- 更准确的定量（特别是source贡献差异大时）
- 不需要预设混合模型
- 对非贡献source的误分配更少

FastST：支持方向性推断的快速源追踪

独特功能：
1. 不需要预先指定source和sink（自动推断方向性！）
2. 计算速度比FEAST更快
3. 方向性推断准确率>90%

适用场景：
- 不确定哪个是source、哪个是sink时
- 研究微生物在不同环境之间的传播方向

# ========== 案例：FMT（粪菌移植）后菌群来源追踪 ==========
# 场景：T2D患者接受FMT后，追踪肠道菌群来自供体还是自身

# 准备数据
# source1: 供体粪便（FMT前）
# source2: 患者自身粪便（FMT前）
# sink:    患者粪便（FMT后不同时间点）

metadata_fmt <- data.frame(
  Env = c("Donor", "Donor", "Donor",        # 供体source
          "Patient_Pre", "Patient_Pre", "Patient_Pre",  # 患者FMT前source
          "FMT_1w", "FMT_4w", "FMT_12w"),   # FMT后1/4/12周的sink
  SourceSink = c(rep("Source", 6), rep("Sink", 3)),
  id = c(rep(1, 3), rep(1, 3), 1, 1, 1),    # 共享id=1
  row.names = c("donor1", "donor2", "donor3",
                "patient_pre1", "patient_pre2", "patient_pre3",
                "patient_1w", "patient_4w", "patient_12w")
)

# 运行FEAST
fmt_result <- FEAST(
  C = as.matrix(otu_table_fmt),
  metadata = metadata_fmt,
  different_sources_flag = 0,
  dir_path = "fmt_tracking/",
  outfile = "fmt_source_tracking"
)

# 预期结果解读：
# FMT后1周：供体菌群占比高（如60%供体，20%自身，20%未知）
# FMT后4周：供体菌群比例下降
# FMT后12周：部分回到自身菌群为主
# 这可以评估FMT的持久性效果

========== 源追踪工具对比 ==========
工具          年份  算法            速度     准确性  方向性推断
SourceTracker2 2011  Gibbs采样       慢      高      不支持
FEAST         2019  EM算法          快(300x) 高      不支持
STENSL        2022  FEAST+L1正则    快       更高    不支持
SourceID-NMF  2024  非负矩阵分解    中       更高    不支持
FastST        2025  新框架          最快     高      支持(>90%)

========== 选择建议 ==========
样本少(<50)、需要稳定结果  → SourceTracker2
样本多(>50)、需要速度      → FEAST
需要最高准确性             → SourceID-NMF
不确定source/sink方向      → FastST

========== 输入数据要求 ==========
1. OTU/ASV丰度表（计数值，非相对丰度）
2. 元数据文件（标注source/sink分组）
3. Source样本：每个环境类型建议 ≥3 个重复
4. OTU数：建议过滤到100-500个（太多太慢）
5. 聚合层面：种/属/科都可以，属水平最常用

========== 结果解读 ==========
比例 > 0.3 → 该source是主要来源
比例 0.1-0.3 → 有一定贡献
比例 < 0.1 → 贡献微小（可能是假阳性）
Unknown > 0.5 → source库不够全面

========== 关键R包 ==========
FEAST：devtools::install_github("cozygene/FEAST")
可视化：ggplot2, tidyr, dplyr
数据处理：phyloseq, biomformat