群体遗传学分析（选择信号 / 群体结构 / Fst / PCA / ADMIXTURE）¶

# 基本质控
plink2 --vcf all_populations.vcf.gz \
  --maf 0.05 \
  --geno 0.1 \
  --mind 0.1 \
  --hwe 1e-6 \
  --max-alleles 2 \
  --make-bed \
  --out popgen_qc

# 计算样本间亲缘关系
plink2 --bfile popgen_qc \
  --king-cutoff 0.0884 \
  --out king_results

# 去除近亲
plink2 --bfile popgen_qc \
  --remove king_results.king.cutoff.out.id \
  --make-bed \
  --out popgen_unrelated

echo "质控后SNP数: $(wc -l < popgen_unrelated.bim)"
echo "质控后样本数: $(wc -l < popgen_unrelated.fam)"

# 严格LD修剪（r²<0.1）
plink2 --bfile popgen_unrelated \
  --indep-pairwise 50 10 0.1 \
  --out ld_prune

# 排除长程LD区域
cat > long_range_ld_hg38.txt << 'EOF'
6 25000000 35000000 MHC
8 7000000 13000000 INV
EOF

plink2 --bfile popgen_unrelated \
  --extract ld_prune.prune.in \
  --exclude range long_range_ld_hg38.txt \
  --make-bed \
  --out popgen_pruned

echo "LD修剪后SNP数: $(wc -l < popgen_pruned.bim)"

# PLINK2 PCA
plink2 --bfile popgen_pruned \
  --pca 20 \
  --out pca_results

# smartPCA（更多统计选项）
# 准备par文件
cat > smartpca.par << 'EOF'
genotypename: popgen_pruned.bed
snpname: popgen_pruned.bim
indivname: popgen_pruned.fam
evecoutname: smartpca.evec
evaloutname: smartpca.eval
numoutevec: 20
numthreads: 16
EOF

smartpca -p smartpca.par

# R可视化PCA
library(ggplot2)
library(ggrepel)

# 读取PCA结果
pca <- read.table("pca_results.eigenvec", header=FALSE)
eigenval <- read.table("pca_results.eigenval", header=FALSE)

# 添加群体信息
pop_info <- read.table("population_labels.txt", header=TRUE)
pca_df <- merge(pca, pop_info, by.x="V2", by.y="IID")

# 计算方差解释比例
pct_var <- eigenval$V1 / sum(eigenval$V1) * 100

# 绘制PCA
ggplot(pca_df, aes(x=V3, y=V4, color=Population)) +
  geom_point(size=2, alpha=0.7) +
  labs(
    x = paste0("PC1 (", round(pct_var[1], 1), "%)"),
    y = paste0("PC2 (", round(pct_var[2], 1), "%)"),
    title = "Population Structure (PCA)"
  ) +
  theme_minimal(base_size=14) +
  scale_color_brewer(palette="Set2")

ggsave("pca_plot.pdf", width=9, height=7)

# 对K=2到K=10运行ADMIXTURE
for K in $(seq 2 10); do
  admixture --cv popgen_pruned.bed ${K} -j16 --seed=42 \
    | tee admixture_K${K}.log
done

# 提取CV error
grep "CV error" admixture_K*.log | \
  awk '{print NR+1, $NF}' > cv_errors.txt

# 可视化CV error选择最优K

# ADMIXTURE结果可视化
library(ggplot2)
library(tidyr)
library(dplyr)

# 读取最优K的Q矩阵（假设K=4最优）
K <- 4
q_matrix <- read.table(paste0("popgen_pruned.", K, ".Q"))
fam <- read.table("popgen_pruned.fam")
pop_info <- read.table("population_labels.txt", header=TRUE)

# 整理数据
q_df <- cbind(IID=fam$V2, q_matrix)
q_df <- merge(q_df, pop_info, by="IID")
colnames(q_df)[2:(K+1)] <- paste0("Ancestry_", 1:K)

# 转为长格式
q_long <- q_df %>%
  pivot_longer(cols=starts_with("Ancestry_"), 
               names_to="Ancestry", values_to="Proportion") %>%
  arrange(Population, Ancestry)

# 结构图
ggplot(q_long, aes(x=IID, y=Proportion, fill=Ancestry)) +
  geom_bar(stat="identity", width=1) +
  facet_grid(~Population, scales="free_x", space="free_x") +
  scale_fill_brewer(palette="Set2") +
  labs(y="Ancestry Proportion", title=paste0("ADMIXTURE K=", K)) +
  theme_minimal() +
  theme(
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    panel.spacing = unit(0.1, "lines")
  )

ggsave("admixture_barplot.pdf", width=12, height=4)

# CV error图
cv <- read.table("cv_errors.txt")
ggplot(cv, aes(x=V1, y=V2)) +
  geom_line() + geom_point() +
  labs(x="K", y="Cross-validation error", title="ADMIXTURE CV Error") +
  scale_x_continuous(breaks=2:10) +
  theme_minimal()

# SNPhylo流水线
snphylo.sh -v filtered.vcf \
  -l 0.1 \
  -m 0.1 \
  -M 0.05 \
  -b 1000 \
  -t 16

# 或手动方式：PLINK距离矩阵 → NJ树
plink2 --bfile popgen_pruned \
  --distance square 1-ibs \
  --out dist_matrix

# 用R构建NJ树

library(ape)
library(ggtree)

# 读取距离矩阵
dist_mat <- read.table("dist_matrix.dist", header=FALSE, skip=1)
ids <- read.table("dist_matrix.dist.id")
rownames(dist_mat) <- ids$V2
colnames(dist_mat) <- ids$V2

# NJ树
nj_tree <- nj(as.dist(dist_mat))
nj_tree <- root(nj_tree, outgroup="outgroup_sample", resolve.root=TRUE)

# 可视化
ggtree(nj_tree, layout="rectangular") +
  geom_tiplab(size=2) +
  theme_tree2()

# VCFtools计算全局Fst
vcftools --vcf filtered.vcf \
  --weir-fst-pop population_A.txt \
  --weir-fst-pop population_B.txt \
  --out fst_A_vs_B

# 滑窗Fst
vcftools --vcf filtered.vcf \
  --weir-fst-pop population_A.txt \
  --weir-fst-pop population_B.txt \
  --fst-window-size 50000 \
  --fst-window-step 25000 \
  --out fst_windowed_A_vs_B

# 多群体Fst（PLINK2）
plink2 --bfile popgen_qc \
  --fst POP \
  --pheno population_assignment.txt \
  --out pairwise_fst

# 计算滑窗π
vcftools --vcf pop_A.vcf \
  --window-pi 50000 \
  --window-pi-step 25000 \
  --out pi_pop_A

vcftools --vcf pop_B.vcf \
  --window-pi 50000 \
  --window-pi-step 25000 \
  --out pi_pop_B

# VCFtools计算Tajima's D
vcftools --vcf pop_A.vcf \
  --TajimaD 50000 \
  --out tajima_pop_A

# 使用selscan计算iHS
selscan --ihs \
  --vcf pop_A_phased.vcf.gz \
  --map genetic_map.map \
  --out ihs_pop_A \
  --threads 16

# 标准化iHS
norm --ihs --files ihs_pop_A.ihs.out --bins 100

# 计算XP-EHH
selscan --xpehh \
  --vcf pop_A_phased.vcf.gz \
  --vcf-ref pop_B_phased.vcf.gz \
  --map genetic_map.map \
  --out xpehh_A_vs_B \
  --threads 16

norm --xpehh --files xpehh_A_vs_B.xpehh.out --bins 100

# PSMC分析流程
# 1. 生成diploid consensus序列
samtools mpileup -C50 -uf ref.fa sample.bam | \
  bcftools call -c | \
  vcfutils.pl vcf2fq -d 10 -D 100 | \
  gzip > diploid.fq.gz

# 2. 创建PSMC输入
fq2psmcfa -q20 diploid.fq.gz > diploid.psmcfa

# 3. 运行PSMC
psmc -N25 -t15 -r5 -p "4+25*2+4+6" -o output.psmc diploid.psmcfa

# 4. Bootstrap（100次）
seq 100 | xargs -I {} -P 16 \
  psmc -N25 -t15 -r5 -b -p "4+25*2+4+6" \
  -o round-{}.psmc diploid.psmcfa

# 5. 绘图（需要指定突变率和世代时间）
psmc_plot.pl -u 1.4e-08 -g 25 -p output output.psmc

# SMC++分析
# 1. 转换VCF为SMC格式
smc++ vcf2smc filtered.vcf.gz smc_data/chr1.smc.gz chr1 Pop1:sample1,sample2,sample3

# 2. 估计群体历史
smc++ estimate -o smc_output/ 1.4e-8 smc_data/*.smc.gz

# 3. 绘图
smc++ plot smc_plot.pdf smc_output/model.final.json

#!/bin/bash
# ====================================================
# 群体遗传学分析完整流程
# ====================================================

set -euo pipefail
THREADS=16
VCF="all_populations.vcf.gz"
POP_FILE="population_info.txt"  # IID Population
OUTDIR="popgen_results"
mkdir -p ${OUTDIR}/{qc,pca,admixture,fst,selection,history}

# ---- Step 1-2: 质控+LD修剪 ----
echo "=== 数据质控 ==="
plink2 --vcf ${VCF} --maf 0.05 --geno 0.1 --mind 0.1 --hwe 1e-6 \
  --max-alleles 2 --make-bed --out ${OUTDIR}/qc/popgen_qc

plink2 --bfile ${OUTDIR}/qc/popgen_qc --king-cutoff 0.0884 --out ${OUTDIR}/qc/king
plink2 --bfile ${OUTDIR}/qc/popgen_qc \
  --remove ${OUTDIR}/qc/king.king.cutoff.out.id \
  --make-bed --out ${OUTDIR}/qc/popgen_unrel

plink2 --bfile ${OUTDIR}/qc/popgen_unrel --indep-pairwise 50 10 0.1 --out ${OUTDIR}/qc/ld
plink2 --bfile ${OUTDIR}/qc/popgen_unrel \
  --extract ${OUTDIR}/qc/ld.prune.in \
  --make-bed --out ${OUTDIR}/qc/popgen_pruned

echo "QC后: $(wc -l < ${OUTDIR}/qc/popgen_unrel.bim) SNPs, $(wc -l < ${OUTDIR}/qc/popgen_unrel.fam) samples"
echo "LD修剪后: $(wc -l < ${OUTDIR}/qc/popgen_pruned.bim) 独立SNPs"

# ---- Step 3: PCA ----
echo "=== PCA ==="
plink2 --bfile ${OUTDIR}/qc/popgen_pruned --pca 20 --out ${OUTDIR}/pca/pca

# ---- Step 4: ADMIXTURE ----
echo "=== ADMIXTURE ==="
for K in $(seq 2 8); do
  admixture --cv ${OUTDIR}/qc/popgen_pruned.bed ${K} -j${THREADS} --seed=42 2>&1 \
    | tee ${OUTDIR}/admixture/log_K${K}.txt
  mv popgen_pruned.${K}.Q ${OUTDIR}/admixture/
  mv popgen_pruned.${K}.P ${OUTDIR}/admixture/
done
grep "CV error" ${OUTDIR}/admixture/log_K*.txt > ${OUTDIR}/admixture/cv_summary.txt

# ---- Step 5: Fst ----
echo "=== Fst ==="
# 提取各群体样本列表
awk '$2=="PopA" {print $1}' ${POP_FILE} > ${OUTDIR}/fst/pop_A.txt
awk '$2=="PopB" {print $1}' ${POP_FILE} > ${OUTDIR}/fst/pop_B.txt

vcftools --gzvcf ${VCF} \
  --weir-fst-pop ${OUTDIR}/fst/pop_A.txt \
  --weir-fst-pop ${OUTDIR}/fst/pop_B.txt \
  --fst-window-size 50000 --fst-window-step 25000 \
  --out ${OUTDIR}/fst/fst_A_vs_B

# ---- Step 6: π和Tajima's D ----
echo "=== π和Tajima's D ==="
vcftools --gzvcf ${VCF} --keep ${OUTDIR}/fst/pop_A.txt \
  --window-pi 50000 --window-pi-step 25000 \
  --out ${OUTDIR}/selection/pi_A

vcftools --gzvcf ${VCF} --keep ${OUTDIR}/fst/pop_A.txt \
  --TajimaD 50000 \
  --out ${OUTDIR}/selection/tajD_A

# ---- Step 7: iHS ----
echo "=== 选择性扫描 ==="
selscan --ihs --vcf pop_A_phased.vcf.gz \
  --map genetic_map.map --out ${OUTDIR}/selection/ihs_A --threads ${THREADS}
norm --ihs --files ${OUTDIR}/selection/ihs_A.ihs.out --bins 100

echo "=== 分析完成 ==="