vcf2treemix.sh $FILE.vcf.gz $FILE.clust
.clust 格式是三列，第一列和第二列相同均为样品名称，第三列是样品所属的群体信息
利用生成的频率文件进行下一步
迁移0-10次，bootstrap 计算1-10次
for i in {0..10}; do for j in {1..10}; do echo treemix -i populations_11pops.treemix.gz -root DK,HK -k 500 -m ${i} -bootstrap -o migration_${i}_bootstrap_${j} ">" ${i}.${j}.out ; done; done > treemix.sh
运行过程中很多命令会跑不出来，不清楚什么原因，如果批量运行时候，最终生成非常多文件，可以使用
grep -L . *llik 命令查看有那些命令没有成功运行，可以单独再次运行。
grep -L . *llik | awk -F . '{print $1"."$2"."$3}'> failed.txt
grep -f failed.txt treemix.sh > treemix_failed.sh
可以放到一条命令执行：
grep -L . *llik | awk -F .ll '{print $1}' > tmp.file && grep -f tmp.file treemix.sh > treemix_failed.sh
library(OptM)
setwd("Treemix生成的所有文件")
source("treemix-1.13/src/plotting_funcs.R")
evanno <- optM(folder ="treemix生成的文件", method = "Evanno",tsv="evanno.txt")
plot_optM(evanno,method = "Evanno",pdf = "evanno.pdf")
par(mfrow=c(3,4))
getwd()
for(edge in 0:10){
# 如果plot_tree函数提供了边距参数，你也可以在此处调整
plot_tree(cex=0.9, paste0("sample.",edge,".1"))
title(paste(edge,"edges"))
}
for(edge in 0:10){
plot_resid(stem=paste0("sample.",edge,".1"),pop_order="../../popmap.list.txt")
title(paste(edge,"edges"))
}
# 保存树图
pdf("trees.pdf", width=12, height=8) # 适合横向排列的5个图
par(mar=c(4, 4, 2, 1))
par(mfrow=c(2,3))
for(edge in 1:5){
plot_tree(cex=0.8, paste0("sample.4",".",edge))
title(paste(edge,"edges"))
}
dev.off()
# 保存残差图
pdf("residuals.pdf", width=12, height=8)
par(mar=c(4, 4, 2, 1))
par(mfrow=c(2,3))
for(edge in 1:5){
plot_resid(stem=paste0("sample.4",".",edge), pop_order="../../popmap.list.txt")
}
dev.off()
2025.6.26更新
大家可以尝试
运行过程中，使用ps -ef | grep treemix时刻关注后台命令的运行时间情况，如果一个命令超过了5分钟，基本上运行过程中出现了nan，推荐修改block为较大值比如1000以上，另外可以在相同文件下，单独把该命令重新运行，出现正常结果（即.llik文件有明确输出，不为空），则利用kill -9 杀死卡死的命令，则剩余命令即可继续运行。
