有许多不同的统计量可以衡量差异选择的区域或基因渗入障碍的区域。这里介绍常用的四个统计量的计算：
- pi, 衡量遗传多样性
- Fst, 衡量群体分化程度
- dxy, 对核苷酸多样性绝对差异的量度
- fd, 对基因流/基因渗入的量度
这里使用 Simon Martin写的脚本，下载here，于Python2的环境下运行。
1. 计算每个物种的pi以及每对物种的Fst和dxy
#转换vcf文件为Simons格式geno.gz
parseVCF.py -i subset.vcf.gz | bgzip > subset.geno.gz
popgenWindows.py -g subset.geno.gz -o subset.Fst.Dxy.pi.csv.gz \
-f phased -w 20000 -m 10000 -s 20000 \
-p PundPyt -p NyerPyt -p NyerMak -p PundMak -p kivu 64253 \
--popsFile pop.info
参数解读：
- -w 20000指定的窗口大小为20 kb
- -s 20000可滑动20 kb
- -m 10000请求这些窗口至少覆盖10 kb的位点
- -f 指定基因型为phased(A|T)而不是unphased(A/T)
- -p 指定群体名
- --popsFile 指定一个文件，该文件一行表示一个个体。
2. 计算fd：
我们测试从群体Pundamilia nyererei（NyerMak）到群体P.sp.的渗入。以基伍湖丽鱼科鱼（NyerPyt）为外群。
ABBABABAwindows.py -w 20000 -m 10 -s 20000 -g subset.geno.gz \
-o subset.fd.PundP.NyerP.NyerM.Kivu.csv.gz \
-f phased --minData 0.5 --writeFailedWindow \
-P1 PundPyt -P2 NyerPyt -P3 NyerMak -O kivu \
--popsFile pop.info
参数：-T可以指定多线程运行。
3. 画图
rm(list = ls())
#读取通过滑动窗口计算得到的 FST, pi 和 dxy
windowStats<-read.csv("./genome_scan/Pundamilia_kivu_div_stats.csv",header=T)
#大体看一眼数据
names(windowStats)
length(windowStats$scaffold)
#添加染色体的长度信息
chrom<-read.table("./genome_scan/chrEnds.txt",header=T)
chrom$add<-c(0,cumsum(chrom$end)[1:21])
windowStats$mid<-windowStats$mid+chrom[match(windowStats$scaffold,chrom$chr),3]
#读取fd输出文件
fd<-read.delim(file = "./genome_scan/Pundamilia_ABBABABA.csv",sep=",",header=T,na.strings = "NaN")
#去除没有位点的窗口
fd<-fd[!is.na(fd$mid),]
#添加位置信息
fd$mid<-fd$mid+chrom[match(fd$scaffold,chrom$chr),3]
#画基因组分布图
par(mfrow=c(2,1),oma=c(3,0,0,0),mar=c(1,5,1,1))
#Fst: Makobe Island
plot(windowStats$mid,windowStats$Fst_NyerMak_PundMak,cex=0.5,pch=19,xaxt="n",
ylab="Fst Makobe",ylim=c(0,1),col=windowStats$scaffold)
abline(h=mean(windowStats$Fst_NyerMak_PundMak,na.rm=T),col="grey")
#Fst: Python Island
plot(windowStats$mid,windowStats$Fst_PundPyt_NyerPyt,cex=0.5,pch=19,xaxt="n",
ylab="Fst Python",ylim=c(0,1),col=windowStats$scaffold)
abline(h=mean(windowStats$Fst_PundPyt_NyerPyt,na.rm=T),col="grey")
axis(1,at=chrom$add+chrom$end/2,tick = F,labels = 1:22)
#Dxy: Makobe Island
plot(windowStats$mid,windowStats$dxy_NyerMak_PundMak,cex=0.5,pch=19,xaxt="n",
ylab="dxy Makobe",ylim=c(0,0.03),col=windowStats$scaffold)
abline(h=mean(windowStats$dxy_NyerMak_PundMak,na.rm=T),col="grey")
#Dxy: Python Island
plot(windowStats$mid,windowStats$dxy_PundPyt_NyerPyt,cex=0.5,pch=19,xaxt="n",ylab="dxy Python",ylim=c(0,0.03),col=windowStats$scaffold)
abline(h=mean(windowStats$dxy_PundPyt_NyerPyt,na.rm=T),col="grey")
#fd
plot(fd$mid,fd$fdM,cex=0.5,pch=19,xaxt="n",
ylab="fdM",ylim=c(0,1),col=windowStats$scaffold)
abline(h=mean(fd$fdM,na.rm=T),col="grey")
axis(1,at=chrom$add+chrom$end/2,tick = F,labels = 1:22)
###############################--------探索性统计---------############################
#dxy和fst相关性
par(mfrow=c(1,1),mar=c(5,5,1,1),oma=c(1,1,1,1))
plot(windowStats$dxy_NyerMak_PundMak,windowStats$Fst_NyerMak_PundMak,
xlab="dxy",ylab="Fst",pch=19,cex=0.3)
#忽略异常值
plot(windowStats$dxy_NyerMak_PundMak,windowStats$Fst_NyerMak_PundMak,
xlab="dxy",ylab="Fst",pch=19,cex=0.3,xlim=c(0,0.01))
#dxy是否与pi相关
plot(rowMeans(cbind(windowStats$pi_PundMak,windowStats$pi_NyerMak),na.rm=T),
windowStats$dxy_NyerMak_PundMak,cex=0.3,
xlab="pi",ylab="dxy")
abline(a=0,b=1,col="grey")
