原理参考：GWAS(8)-基于Nucleotide diversity (π) and Tajima’s D评价基因的受选择趋势
操作步骤：
计算两个群体之间的FST、tajimaD以及PI
版本：VCFtools - 0.1.16
#FST分析
vcftools \
--gzvcf 142samples_basefilter_only_snp_qc.vcf.gz \
--weir-fst-pop /home/147_Tg_reseq/142_ped_final_gwas_result/F0_sample.txt \
--weir-fst-pop /home/147_Tg_reseq/142_ped_final_gwas_result/CK_sample.txt \
--fst-window-size 50000 \
--fst-window-step 5000 \
--out 142rna_F0_vs_CK_fst_50k
# Pi (CK)
vcftools --gzvcf 142samples_basefilter_only_snp_qc.vcf.gz --keep /home/147_Tg_reseq/142_ped_final_gwas_result/CK_sample.txt --window-pi 50000 --window-pi-step 5000 --out 142rna_CK_pi_50k
# Pi (FO)
nohup vcftools --gzvcf 142samples_basefilter_only_snp_qc.vcf.gz --keep /home/147_Tg_reseq/142_ped_final_gwas_result/F0_sample.txt --window-pi 50000 --window-pi-step 5000 --out 142rna_FO_pi_50k &
# Tajima's D (CK)
vcftools --gzvcf 142samples_basefilter_only_snp_qc.vcf.gz --keep /home/147_Tg_reseq/142_ped_final_gwas_result/CK_sample.txt --TajimaD 50000 --out tajima_142rna_CK
# Tajima's D (FO)
nohup vcftools --gzvcf 142samples_basefilter_only_snp_qc.vcf.gz --keep /home/147_Tg_reseq/142_ped_final_gwas_result/F0_sample.txt --TajimaD 50000 --out tajima_142rna_FO &对基因区间绘图
setwd("D:/1、博士期间/9、博士-实验以及项目/12、根腐病GWAS/基因型数据分析/142重测序家系填补/全部SNP的FST.TAJIMAD.PI分析结果")
library(ggplot2)
library(tidyr)
library(dplyr)
library(stringr)
# ===================== 1. 读取你的所有目标区间（我已经帮你整理好） =====================
region_list <- read.table(text = "
chr1 176693484 176784176
chr2 73021605 73175631
chr4 1246972788 1247060365
chr4 460799005 460856479
chr4 74233656 74258851
chr8 1173586992 1173651964
chr1 1242500093 1242672786
chr5 43301051 43316787
chr10 319014490 319015283
chr11 299648144 299800123
chr1 1306839664 1306963608
chr1 1767755699 1767914966
chr1 176784276 176827147
chr1 734046993 734121308
chr3 1466316772 1466392061
chr6 19682071 19859017
chr10 1737061712 1737148127
chr10 1959351501 1959461436
chr11 277045495 277122842
chr7 1085660486 1085683859
chr7 23595365 23778406
chr9 1241896772 1241908154
chr9 246672295 246754573
chr9 57029759 57154375
chr9 847142499 847241610
chr10 1750318896 1750357755
chr3 92906933 92936805
chr4 1117894753 1117945591
chr4 620705166 620739253
chr6 309482833 309484607
chr7 955316189 955463255
chr8 1430561950 1430658113
chr8 495344113 495466444
chr3 587553104 587668233
chr6 1292482231 1292498452
chr7 1215768527 1215873687
chr10 1503740373 1503761902
chr7 1247466817 1247521726
chr9 296597210 296647659
chr2 1144225112 1144228601
chr4 1356534103 1356655597
chr4 721692727 721834551
chr6 172650108 172704181
chr11 1667408672 1667409920
chr11 1789625382 1789751630
chr2 169533784 169658940
chr3 46898832 46936312
chr11 1698608243 1698683197
chr11 973083616 973147864
chr1 1675149734 1675179088
chr1 1728621262 1728743829
chr2 1018440755 1018484555
chr2 271574643 271672639
chr3 1466415486 1466598092
chr3 1758058989 1758080411
chr7 1354195202 1354304232
chr7 613694601 613711356
chr9 1260771716 1260920996
chr9 342136311 342178623
chr9 435573794 435665317
chr9 826644128 826757342
chr1 239857157 240029902
chr1 855815678 855976112
chr4 1215818014 1216115604
chr5 1855349829 1855425572
chr5 1969065026 1969069805
chr6 1321698759 1321778338
chr6 372559260 372616373
chr7 1349300017 1349369698
chr9 349492701 349625773
chr10 1304273635 1304336683
chr10 402878769 402909675
chr1 1620059080 1620170455
chr1 1675192234 1675216832
chr5 4920251 4984069
chr7 110166387 110192827
chr7 576556456 576592932
chr8 1431615534 1431618202
chr1 1489379487 1489411208
chr5 329881447 329922732
chr8 552878192 552965744
", header=F, stringsAsFactors=F)
colnames(region_list) <- c("chr", "start", "end")
# ===================== 2. 读取Fst / Pi / TajimaD 数据（和你原来一样） =====================
file_info <- data.frame(
type = c("142_F0_vs_CK_fst","142_CK_pi_50k","142_FO_pi_50k","tajima_142_CK","tajima_142_FO"),
filename = c("142_F0_vs_CK_fst_50k.windowed.weir.fst","142_CK_pi_50k.windowed.pi","142_FO_pi_50k.windowed.pi",
"tajima_142_CK.Tajima.D","tajima_142_FO.Tajima.D")
)
# 统一读取并格式化
alldata <- data.frame()
for(i in 1:nrow(file_info)){
type <- file_info$type[i]
data <- read.table(file_info$filename[i], header=T, stringsAsFactors=F)
data_f <- data[, c(1,2,ncol(data))]
colnames(data_f) <- c("chr","start","val")
data_f <- data_f %>% filter(grepl("^Chr[1-9]$|^Chr1[0-1]$", chr))
data_f$type <- type
alldata <- rbind(alldata, data_f)
}
# 格式化
alldata_f <- alldata %>%
separate(type, into=c("TP","pop"), sep="_fo", extra="merge", fill="right") %>%
mutate(pop = ifelse(TP=="Fst", "Fst", pop))
# Fst 负值修正
alldata_f$val[alldata_f$TP=="Fst" & alldata_f$val<0] <- 0
# ===================== 3. 批量循环：每个区间画一张图！ =====================
for(i in 1:nrow(region_list)){
# 当前区间
chr_i <- region_list$chr[i]
start_i <- region_list$start[i]
end_i <- region_list$end[i]
# 匹配数据（Chr大写）
plot_data <- alldata_f %>%
mutate(chr_upper = toupper(chr)) %>%
filter(chr_upper == toupper(chr_i),
start >= start_i,
start <= end_i)
if(nrow(plot_data) < 3) next # 跳过无数据区间
# 绘图（和你原图完全一样风格）
p <- ggplot(plot_data, aes(x=start, y=val, color=pop)) +
geom_point(size=0.6, alpha=0.6) +
geom_smooth(method="loess", span=0.2, linewidth=0.8, se=F) +
facet_wrap(~TP, ncol=1, scales="free_y", strip.position="left") +
scale_color_manual(values=c("black","#E41A1C","#377EB8")) +
theme_bw() +
labs(x=paste0(chr_i, " Position"), y="") +
theme(strip.background = element_blank(),
strip.text = element_text(size=12, face="bold"),
strip.placement = "outside",
panel.grid = element_blank(),
legend.title = element_blank(),
legend.position = "top")
# 保存
outname <- paste0("region_", chr_i, "_", start_i, "_", end_i, ".pdf")
ggsave(outname, p, width=10, height=7, device="pdf")
message("✅ 已生成: ", outname)
}
message("\n🎉 全部完成！所有区间的图都已保存！")