群体遗传结构分析 — ADMIXTURE¶

VCF文件 → PLINK格式转换 → LD pruning → ADMIXTURE运行(K=1~10) → CV error选最佳K → 可视化

# 安装 ADMIXTURE（下载二进制文件）
wget https://dalexander.github.io/admixture/binaries/admixture_linux-1.3.0.tar.gz  # 下载压缩包
tar -xzf admixture_linux-1.3.0.tar.gz  # 解压
chmod +x admixture  # 添加执行权限
export PATH=$PATH:$(pwd)  # 添加到环境变量

# 安装 PLINK（用于格式转换和LD pruning）
conda install -c bioconda plink  # 通过conda安装plink

# 将 VCF 转换为 PLINK 的 bed/bim/fam 格式
plink --vcf input.vcf.gz \       # 输入VCF文件
      --make-bed \                # 生成bed格式
      --allow-extra-chr \         # 允许非标准染色体名
      --set-missing-var-ids @:# \ # 给没有ID的SNP设置ID（染色体:位置）
      --out mydata                # 输出文件前缀

# 输出文件说明：
# mydata.bed — 基因型的二进制文件（核心数据）
# mydata.bim — 每个SNP的信息（染色体、位置、等位基因）
# mydata.fam — 每个样本的信息（家系、性别、表型）

# 为什么要 LD pruning？
# ADMIXTURE 假设 SNP 之间是独立的，连锁的 SNP 会让结果偏差

# 第一步：标记要保留/删除的 SNP
plink --bfile mydata \            # 输入bed文件
      --indep-pairwise 50 10 0.1 \ # 窗口50个SNP，步长10，r²阈值0.1
      --out mydata_pruned          # 输出前缀

# 第二步：只保留独立的 SNP
plink --bfile mydata \            # 输入原始bed文件
      --extract mydata_pruned.prune.in \ # 只保留通过筛选的SNP
      --make-bed \                 # 生成新的bed文件
      --out mydata_ldpruned        # 输出前缀

# 检查保留了多少 SNP
wc -l mydata_ldpruned.bim  # 统计SNP数量
# 一般保留几万到几十万个SNP就够了

# 对 K=1 到 K=10 分别运行 ADMIXTURE
for K in $(seq 1 10); do
    admixture --cv \              # 启用交叉验证（计算CV error）
              mydata_ldpruned.bed \ # 输入文件
              $K \                 # 当前K值
              -j4 \                # 使用4个线程
              | tee log_K${K}.out  # 同时输出到屏幕和日志文件
done

# 输出文件说明（以K=3为例）：
# mydata_ldpruned.3.Q — Q矩阵，每行是一个个体，3列分别是3个祖先群体的比例
# mydata_ldpruned.3.P — P矩阵，每行是一个SNP在各祖先群体中的等位基因频率

# 提取所有K值的CV error
grep -h "CV" log_K*.out | awk '{print $3,$4}' > cv_errors.txt
# 输出格式：(K=1): 0.54321
# 输出格式：(K=2): 0.43210
# ...

# 用R画CV error曲线
cat << 'EOF' > plot_cv.R
# 读取CV error数据
cv <- read.table("cv_errors.txt")                # 读取文件
cv$K <- as.numeric(gsub("[^0-9]", "", cv$V1))    # 提取K值
cv$error <- as.numeric(cv$V2)                     # 提取CV error

# 画折线图
pdf("cv_error_plot.pdf", width=8, height=5)       # 保存为PDF
plot(cv$K, cv$error,                              # X轴K值，Y轴CV error
     type="b",                                     # 画线+点
     xlab="K (Number of ancestral populations)",   # X轴标签
     ylab="Cross-validation error",                # Y轴标签
     main="ADMIXTURE CV Error",                    # 标题
     pch=16, col="steelblue")                      # 点的样式和颜色
abline(v=cv$K[which.min(cv$error)],               # 在最小CV error处画竖线
       col="red", lty=2)                           # 红色虚线
dev.off()                                          # 关闭PDF设备
EOF
Rscript plot_cv.R  # 执行R脚本

# 最小 CV error 对应的 K 就是"最佳K"
# 但注意：如果多个K的CV error很接近，应该报告一个范围而不是单个K

# 方法1：用R画经典的ADMIXTURE柱状图
cat << 'EOF' > plot_admixture.R
# 读取最佳K的Q矩阵（假设K=3）
Q <- read.table("mydata_ldpruned.3.Q")    # 每行一个个体，每列一个祖先群体
fam <- read.table("mydata_ldpruned.fam")  # 读取样本信息

# 读取群体标签（需要自己准备）
# pop_info.txt 格式：样本ID  群体标签
pop <- read.table("pop_info.txt", header=FALSE)  # 读取群体信息

# 按群体排序（让同一群体的个体排在一起）
ord <- order(pop$V2)                       # 按群体名排序
Q_sorted <- Q[ord, ]                       # 排序Q矩阵
pop_sorted <- pop[ord, ]                   # 排序群体信息

# 画柱状图
pdf("admixture_K3.pdf", width=15, height=5)  # 创建宽PDF
barplot(t(as.matrix(Q_sorted)),               # 转置后画堆叠柱状图
        col=c("#E41A1C","#377EB8","#4DAF4A"), # 三种颜色对应三个祖先
        border=NA,                             # 不画柱子边框
        space=0,                               # 柱子之间无间距
        ylab="Ancestry proportion",            # Y轴标签
        main="ADMIXTURE K=3")                  # 标题

# 在底部标注群体名
# 用 tapply 计算每个群体的中间位置
mid <- tapply(1:nrow(Q_sorted), pop_sorted$V2, mean)  # 每个群体的中间位置
axis(1, at=mid, labels=names(mid), las=2, cex.axis=0.7)  # 添加X轴标签
dev.off()
EOF
Rscript plot_admixture.R  # 执行画图

# 方法2：用 pong 进行更专业的可视化（推荐）
# pong 可以自动对齐多次运行结果，避免颜色翻转问题
pip install pong  # 安装pong

# 准备pong需要的文件列表（filemap）
# 格式：K值  运行编号  Q文件路径
echo "3 1 mydata_ldpruned.3.Q" > pong_filemap.txt

# 运行pong
pong -m pong_filemap.txt \     # 输入filemap
     -i pop_info.txt \          # 群体标签文件
     -n pop_order.txt           # 群体排列顺序（可选）
# pong 会启动一个网页服务器，浏览器打开即可查看交互式结果

# ADMIXTURE 使用随机种子，每次运行结果可能略有不同
# 建议每个K值运行5-10次，检查结果稳定性

for K in $(seq 1 10); do
    for rep in $(seq 1 5); do
        admixture --cv \
                  -s ${rep}${K} \     # 设置不同的随机种子
                  mydata_ldpruned.bed \
                  $K \
                  -j4
        # 重命名输出文件，避免覆盖
        mv mydata_ldpruned.${K}.Q mydata_ldpruned.K${K}.rep${rep}.Q  # 重命名Q矩阵
        mv mydata_ldpruned.${K}.P mydata_ldpruned.K${K}.rep${rep}.P  # 重命名P矩阵
    done
done

# 如果你知道部分个体的群体归属（参考群体），可以用有监督模式
# 在 .fam 文件的第6列（表型列）标注群体编号
# 已知群体标为 1, 2, 3...；未知个体标为 -

admixture --supervised \          # 启用有监督模式
          mydata_ldpruned.bed \   # 输入文件
          3                        # K值等于已知群体数

# evaladmix 检查 ADMIXTURE 结果是否很好地拟合了数据
# 如果某些群体对之间的残差很大，说明模型拟合不好

# 安装 evaladmix
git clone https://github.com/GenisGE/evalAdmix.git  # 克隆仓库
cd evalAdmix && make                                  # 编译

# 运行评估
./evalAdmix -plink mydata_ldpruned \  # 输入PLINK文件前缀
            -fname mydata_ldpruned.3.P \  # P矩阵
            -qname mydata_ldpruned.3.Q \  # Q矩阵
            -o evaladmix_K3.corres        # 输出相关性残差矩阵

# === ADMIXTURE 完整流程速查 ===

# 1. VCF → PLINK
plink --vcf input.vcf.gz --make-bed --allow-extra-chr --out mydata

# 2. LD pruning
plink --bfile mydata --indep-pairwise 50 10 0.1 --out pruned
plink --bfile mydata --extract pruned.prune.in --make-bed --out data_ld

# 3. 运行 ADMIXTURE（K=1~10）
for K in $(seq 1 10); do admixture --cv data_ld.bed $K -j4 | tee log_K${K}.out; done

# 4. 提取 CV error
grep "CV" log_K*.out

# 5. 可视化
# R: barplot(t(Q), col=colors, border=NA, space=0)
# pong: pong -m filemap.txt -i pop_info.txt