群体间基因流检测 — TreeMix¶

# conda 安装
conda install -c bioconda treemix  # 安装TreeMix

# 或源码编译
# 需要依赖：GSL, Boost
sudo apt-get install libgsl-dev libboost-all-dev  # 安装依赖
wget https://bitbucket.org/nygcresearch/treemix/downloads/treemix-1.13.tar.gz
tar -xzf treemix-1.13.tar.gz
cd treemix-1.13
./configure && make && sudo make install

# TreeMix 需要特殊的等位基因频率输入格式
# 每行一个 SNP，每列一个群体，格式为 "ref_count,alt_count"

# 从 VCF 生成 TreeMix 输入
# 方法1：用 plink + 自定义脚本

# 先生成 PLINK 频率文件
plink --vcf input.vcf.gz \
      --freq \                     # 计算等位基因频率
      --family \                    # 按群体（FID）分组
      --out allele_freq             # 输出前缀

# 方法2：用 stacks 的 populations 模块（如果用的是RAD数据）
# populations -P stacks_dir -M popmap.txt --treemix

# 方法3：直接从 VCF 用 Python 脚本转换
cat << 'PYEOF' > vcf_to_treemix.py
import gzip
import sys
from collections import defaultdict

# 读取群体分配文件（样本名 群体名）
pop_map = {}                                        # 样本->群体映射
with open(sys.argv[2]) as f:                        # 打开群体文件
    for line in f:
        sample, pop = line.strip().split()          # 解析样本和群体
        pop_map[sample] = pop

# 读取VCF并计算频率
with gzip.open(sys.argv[1], 'rt') as f:            # 打开VCF
    for line in f:
        if line.startswith('##'):                   # 跳过注释行
            continue
        if line.startswith('#CHROM'):               # 表头行
            samples = line.strip().split('\t')[9:]  # 获取样本名
            pops = sorted(set(pop_map.values()))    # 获取群体列表
            print(' '.join(pops))                   # 打印群体名表头
            continue
        fields = line.strip().split('\t')
        counts = defaultdict(lambda: [0, 0])        # 每个群体的ref/alt计数
        for i, sample in enumerate(samples):
            gt = fields[9+i].split(':')[0]          # 获取基因型
            pop = pop_map.get(sample, None)
            if pop is None: continue
            alleles = gt.replace('|', '/').split('/')
            for a in alleles:
                if a == '0': counts[pop][0] += 1     # ref计数
                elif a == '1': counts[pop][1] += 1   # alt计数
        # 输出格式：ref,alt  ref,alt  ...
        out = []
        for pop in pops:
            out.append(f"{counts[pop][0]},{counts[pop][1]}")
        print(' '.join(out))
PYEOF

python3 vcf_to_treemix.py input.vcf.gz popmap.txt > treemix_input.txt
gzip treemix_input.txt  # 压缩（TreeMix支持gz输入）

# 运行 m=0 到 m=5 的模型（0到5条迁移边）
for m in $(seq 0 5); do
    treemix -i treemix_input.txt.gz \  # 输入文件
            -m $m \                      # 迁移边数量
            -o treemix_m${m} \          # 输出前缀
            -root outgroup_pop \         # 指定外群（重要！）
            -bootstrap \                 # 做bootstrap
            -k 500 \                     # SNP block大小（考虑LD）
            -seed 12345                  # 随机种子（可重复）
done

# 输出文件说明：
# treemix_m0.treeout.gz   — 树拓扑（Newick格式）
# treemix_m0.vertices.gz  — 节点信息
# treemix_m0.edges.gz     — 边信息
# treemix_m0.covse.gz     — 协方差矩阵标准误
# treemix_m0.modelcov.gz  — 模型拟合的协方差矩阵
# treemix_m0.llik         — 对数似然值

# ========== R 脚本：用 OptM 确定最佳 m ==========
# install.packages("OptM")
library(OptM)                                       # 加载OptM包

# 读取所有TreeMix运行的似然值
optm_result <- optM("treemix_m",                   # 输出文件前缀
                    method="Evanno")                # 使用Evanno方法

# 画图
pdf("OptM_result.pdf", width=8, height=6)
plot_optM(optm_result)                              # 画Δm图
dev.off()

# Δm 最高点对应的 m 就是最佳迁移边数量
# 类似于 ADMIXTURE 的 CV error 选 K

# ========== R 脚本：画 TreeMix 图 ==========
# TreeMix 自带R脚本
source("treemix-1.13/src/plotting_funcs.R")        # 加载画图函数

# 画迁移树（最常用的图）
pdf("treemix_m3.pdf", width=10, height=8)
plot_tree("treemix_m3")                             # 画带迁移边的树
# 迁移边的颜色表示迁移比例（从黄色到红色）
# 箭头方向表示基因流方向
dev.off()

# 画残差矩阵（检查模型拟合）
pdf("treemix_m3_residuals.pdf", width=8, height=8)
plot_resid("treemix_m3",                            # 画残差图
           pop_order="pop_order.txt")               # 群体排列顺序
dev.off()
# 正残差（蓝色）= 实际比模型更相似 → 可能缺少迁移边
# 负残差（红色）= 实际比模型更不同 → 模型可能过拟合

# D-统计量是检测基因流的另一种方法
# 比 TreeMix 更简单直接，常用于验证

# 用 ADMIXTOOLS 计算 D-统计量
# 安装
conda install -c bioconda admixtools  # 安装ADMIXTOOLS

# 准备参数文件
# D(P1, P2; P3, O) 检验 P3 是否与 P1 或 P2 存在基因流
# D > 0 → P3 与 P2 基因流
# D < 0 → P3 与 P1 基因流
# D ≈ 0 → 无基因流

cat << 'EOF' > d_stat_params.txt
genotypename: mydata.geno
snpname: mydata.snp
indivname: mydata.ind
popfilename: pop_combinations.txt
EOF

# pop_combinations.txt 格式：
# P1  P2  P3  Outgroup
echo "PopA PopB PopC Outgroup" > pop_combinations.txt

# 运行 D-统计量
qpDstat -p d_stat_params.txt > d_stat_results.txt

# 结果解读：
# |Z-score| > 3 → 基因流显著

# === TreeMix 速查 ===

# 运行（m条迁移边）
treemix -i input.txt.gz -m 3 -o out_m3 -root outgroup -bootstrap -k 500

# 批量运行 m=0~5
for m in $(seq 0 5); do treemix -i in.txt.gz -m $m -o out_m${m} -root OG -k 500; done

# 选择最佳m（R）
library(OptM); optM("out_m")

# 可视化（R）
source("plotting_funcs.R")
plot_tree("out_m3")
plot_resid("out_m3", pop_order="order.txt")

# D-统计量验证
qpDstat -p params.txt > d_results.txt
# |Z| > 3 = 显著基因流