614_群体遗传学PCA_ADMIXTURE¶

# 安装PLINK（数据处理+PCA）
conda install -c bioconda plink  # PLINK 1.9
conda install -c bioconda plink2  # PLINK 2.0（更快）

# 安装ADMIXTURE
conda install -c bioconda admixture  # ADMIXTURE

# 安装smartpca（EIGENSOFT的PCA工具）
conda install -c bioconda eigensoft  # 可选

# 将VCF文件转换为PLINK的bed/bim/fam格式
plink --vcf input.vcf.gz \    # 输入VCF
    --make-bed \               # 输出二进制格式
    --out dataset \            # 输出前缀
    --allow-extra-chr          # 允许非标准染色体名

# === 第1步：过滤个体（样本水平）===
plink --bfile dataset \
    --mind 0.1 \               # 去除缺失率>10%的个体
    --make-bed \
    --out dataset_qc1

# === 第2步：过滤位点（SNP水平）===
plink --bfile dataset_qc1 \
    --geno 0.05 \              # 去除缺失率>5%的SNP
    --maf 0.05 \               # 去除MAF<5%的SNP
    --hwe 1e-6 \               # 去除严重偏离HWE的SNP
    --make-bed \
    --out dataset_qc2

# === 第3步：LD修剪（PCA和ADMIXTURE必须做！）===
# 先标记高LD的SNP
plink --bfile dataset_qc2 \
    --indep-pairwise 50 5 0.2 \  # 窗口50SNP，步长5，r²>0.2则去除
    --out ld_pruned

# 然后只保留独立的SNP
plink --bfile dataset_qc2 \
    --extract ld_pruned.prune.in \  # 保留这些SNP
    --make-bed \
    --out dataset_final

echo "QC前SNP数: $(wc -l < dataset.bim)"     # QC前
echo "QC后SNP数: $(wc -l < dataset_final.bim)"  # QC后

# 计算个体间的IBS（相同等位基因比例）
plink --bfile dataset_final \
    --genome \                 # 计算亲缘系数
    --out relatedness

# 去除一级亲属（PI_HAT > 0.25）
awk '$10 > 0.25 {print $1, $2}' relatedness.genome > related_individuals.txt
# 手动选择要去除的个体，保留表型信息更完整的那个

# === 运行PCA ===
plink --bfile dataset_final \
    --pca 20 \                 # 计算前20个主成分
    --out pca_result           # 输出前缀

# 输出文件：
# pca_result.eigenval — 特征值（每个PC解释的变异量）
# pca_result.eigenvec — 特征向量（每个个体在各PC上的坐标）

library(ggplot2)  # 加载ggplot2

# === 读取PCA结果 ===
eigenvec <- read.table("pca_result.eigenvec", header = FALSE)  # 读取坐标
eigenval <- read.table("pca_result.eigenval", header = FALSE)  # 读取特征值

# 设置列名
colnames(eigenvec) <- c("FID", "IID", paste0("PC", 1:20))  # 命名列

# 读取群体信息
pop_info <- read.table("population_info.txt", header = TRUE)  # 群体信息
# 格式：FID IID Population
data <- merge(eigenvec, pop_info, by = c("FID", "IID"))  # 合并

# 计算各PC解释的方差比例
pve <- eigenval$V1 / sum(eigenval$V1) * 100  # 方差解释比例(%)

# === 画PCA散点图 ===
ggplot(data, aes(x = PC1, y = PC2, color = Population)) +
    geom_point(size = 2, alpha = 0.7) +  # 散点
    labs(
        title = "Population Structure (PCA)",
        x = paste0("PC1 (", round(pve[1], 1), "%)"),  # 加方差比例
        y = paste0("PC2 (", round(pve[2], 1), "%)"),
        color = "Population"
    ) +
    theme_minimal() +  # 简洁主题
    theme(legend.position = "right")  # 图例在右边

ggsave("pca_plot.pdf", width = 10, height = 7)  # 保存

# === 画特征值碎石图（确定有多少个有意义的PC）===
scree_df <- data.frame(PC = 1:20, PVE = pve[1:20])
ggplot(scree_df, aes(x = PC, y = PVE)) +
    geom_bar(stat = "identity", fill = "steelblue") +  # 柱状图
    geom_line(color = "red") +  # 折线
    labs(x = "Principal Component", y = "Variance Explained (%)") +
    theme_minimal()

ggsave("scree_plot.pdf", width = 8, height = 5)  # 保存

# === 对不同K值运行ADMIXTURE ===
# K从1到10各跑一次，找最佳K值
for K in $(seq 1 10); do
    admixture \
        --cv \                    # 计算交叉验证误差（关键！）
        -j8 \                     # 8个线程
        dataset_final.bed \       # 输入文件（PLINK格式）
        ${K} \                    # K值
        2>&1 | tee log_K${K}.txt  # 同时输出到屏幕和文件
done

# === 提取CV error ===
grep -h "CV" log_K*.txt | awk '{print $3, $4}' | \
    sed 's/(K=//;s/)://' > cv_errors.txt  # 提取CV误差

# 查看结果
cat cv_errors.txt
# 选择CV error最低的K值

library(ggplot2)  # 加载ggplot2
library(reshape2)  # 数据重塑

# === 1. 画CV error图，确定最佳K ===
cv <- read.table("cv_errors.txt", header = FALSE)  # 读取CV误差
colnames(cv) <- c("K", "CV_error")  # 命名列

ggplot(cv, aes(x = K, y = CV_error)) +
    geom_line(color = "blue") +  # 折线
    geom_point(color = "red", size = 3) +  # 点
    scale_x_continuous(breaks = 1:10) +  # x轴刻度
    labs(
        title = "Cross-Validation Error",
        x = "K (Number of ancestral populations)",
        y = "CV Error"
    ) +
    theme_minimal()

ggsave("cv_error_plot.pdf", width = 8, height = 5)  # 保存

# === 2. 画ADMIXTURE柱状图（以最佳K=3为例）===
best_K <- 3  # 假设最佳K是3

# 读取Q矩阵
q_matrix <- read.table(
    paste0("dataset_final.", best_K, ".Q"),  # 读取Q文件
    header = FALSE
)
colnames(q_matrix) <- paste0("Ancestry_", 1:best_K)  # 命名列

# 添加个体信息
fam <- read.table("dataset_final.fam", header = FALSE)  # 读取fam文件
q_matrix$IID <- fam$V2  # 个体ID
q_matrix$Population <- pop_info$Population[match(fam$V2, pop_info$IID)]

# 按群体排序
q_matrix <- q_matrix[order(q_matrix$Population), ]
q_matrix$Index <- 1:nrow(q_matrix)  # 添加索引

# 长格式转换
q_long <- melt(q_matrix, 
    id.vars = c("IID", "Population", "Index"),  # 保留的列
    variable.name = "Ancestry",  # 变量名
    value.name = "Proportion"    # 值名
)

# 画ADMIXTURE柱状图
ggplot(q_long, aes(x = Index, y = Proportion, fill = Ancestry)) +
    geom_bar(stat = "identity", width = 1) +  # 堆叠柱状图
    facet_grid(~Population, scales = "free_x", space = "free_x") +  # 按群体分面
    labs(
        title = paste0("ADMIXTURE (K=", best_K, ")"),
        x = "Individuals",
        y = "Ancestry Proportion"
    ) +
    theme_minimal() +
    theme(
        axis.text.x = element_blank(),  # 隐藏x轴文字
        axis.ticks.x = element_blank(),  # 隐藏x轴刻度
        strip.text = element_text(angle = 45, hjust = 1)  # 群体标签斜45度
    ) +
    scale_fill_brewer(palette = "Set2")  # 配色方案

ggsave("admixture_plot.pdf", width = 14, height = 4)  # 保存

#!/bin/bash
# 群体结构分析完整流程

INPUT="raw_data"      # 原始PLINK文件前缀
OUTPUT="pop_structure" # 输出目录
THREADS=16            # 线程数

mkdir -p ${OUTPUT}

# === QC ===
plink --bfile ${INPUT} \
    --mind 0.1 --geno 0.05 --maf 0.05 --hwe 1e-6 \
    --make-bed --out ${OUTPUT}/qc

# === LD pruning ===
plink --bfile ${OUTPUT}/qc \
    --indep-pairwise 50 5 0.2 \
    --out ${OUTPUT}/ld

plink --bfile ${OUTPUT}/qc \
    --extract ${OUTPUT}/ld.prune.in \
    --make-bed --out ${OUTPUT}/final

# === PCA ===
plink --bfile ${OUTPUT}/final --pca 20 --out ${OUTPUT}/pca

# === ADMIXTURE (K=1~10) ===
cd ${OUTPUT}
for K in $(seq 1 10); do
    admixture --cv -j${THREADS} final.bed ${K} 2>&1 | tee log_K${K}.txt
done
grep -h "CV" log_K*.txt > cv_errors.txt
cd ..

echo "分析完成！使用R脚本进行可视化"