Tajima D
这个是选择相关的一个参数，大于0代表群体观测杂合度高于预期杂合度，稀有等位基因频率降低（群体收缩或者平衡选择），小于0说明群体观测杂合位点少于预期值，稀有等位基因频率增加（群体扩张或者低频选择）。 也就是说，只有0是正常的，其他都是选择发生。
π
π，核苷酸多样性，越大说明核苷酸多样性越高，越低说明两个座位DNA序列差异越小。
Fst
Fst, 分化系数，从0到1说明亲缘关系越来越远。接近于0说明两个个体亲缘关系近，接近1说明亲缘关系远。
Hardy-Weinberg
Hardy-Weinberg平衡检测，这个主要是检测基因型频率是否等于基因频率乘积。比如A:0.3,a:0.7那么Aa的频率是否为0.42
####计算pi
vcftools --vcf merged.vcf --window-pi 500000 --out pi
####计算TajimaD
vcftools --vcf marged.vcf --TajimaD 500000 --out TajimaD
###计算HW
vcftools --vcf merged.vcf --hardy --out HW
###fst需要两个以上pop
vcftools --vcf SNP.vcf --weir-fst-pop 1.txt --weir-fst-pop 2.txt --out 1_VS_2 --fst-window-size 500000
### 将亚群拆分出来
bcftools view -S ARO merged.vcf -O v -o ARO.vcf
绘图：
library(ggplot2)
###读取相应文件即可
mydata<-read.table("pi.windowed.pi",header = T)
mydata1<-na.omit(mydata)
###画pi
pdf('pi.pdf')
ggplot(mydata1,aes(x=BIN_START/1000000,y=PI,group=factor(CHROM),colour=CHROM))+geom_line()+facet_wrap(mydata1$CHROM)+xlab("Physical distance(Mb)")+ylab("PI")+theme(legend.position = "none")+ theme_bw() +theme(panel.grid.major=element_line(colour=NA), panel.background = element_rect(fill = "transparent",colour = NA),
plot.background = element_rect(fill = "transparent",colour = NA),
panel.grid.minor = element_blank())
dev.off()
###画TJM D
pdf('TJM.pdf')
ggplot(mydata1,aes(x=BIN_START/1000000,y=TajimaD,group=factor(CHROM),colour=CHROM))+geom_line()+facet_wrap(mydata1$CHROM)+xlab("Physical distance(Mb)")+ylab("Tajima's D")+theme(legend.position = "none")+ theme_bw() +theme(panel.grid.major=element_line(colour=NA), panel.background = element_rect(fill = "transparent",colour = NA),
plot.background = element_rect(fill = "transparent",colour = NA),
panel.grid.minor = element_blank())
dev.off()
###画HW
ggplot(a,aes(x=POS/1000000,y=P_HWE,group=factor(CHR),color=CHR))+geom_point()+facet_wrap(a$CHR)+xlab("Physical distance(Mb)")+theme_bw() +theme(panel.grid.major=element_line(colour=NA), panel.background = element_rect(fill = "transparent",colour = NA),
plot.background = element_rect(fill = "transparent",colour = NA),
panel.grid.minor = element_blank())
结果：
TJM
pi
HW
#添加分组信息NIP、wild、aus、tro、tem、aro
data1=read.table('NIP.windowed.pi',header = T)
Group=c("NIP")
NIP=data.frame(data1,Group)
data2=read.table('wild.windowed.pi',header = T)
Group=c("wild")
wild=data.frame(data2,Group)
data3=read.table('aus.windowed.pi',header = T)
Group=c("aus")
aus=data.frame(data3,Group)
data4=read.table('tro.windowed.pi',header = T)
Group=c("tro")
tro=data.frame(data4,Group)
data5=read.table('tem.windowed.pi',header = T)
Group=c("tem")
tem=data.frame(data5,Group)
data6=read.table('aro.windowed.pi',header = T)
Group=c("aro")
aro=data.frame(data6,Group)
data7=read.table('XI.windowed.pi',header = T)
Group=c("XI")
XI=data.frame(data7,Group)
#将6组数据合并
library(ggplot2)
library(dplyr)
total_data<-dplyr::bind_rows(NIP,wild,aus,tro,tem,aro,XI)
#画箱线图+小提琴图
pdf("all_pi.pdf")
p <- ggplot(total_data,aes(Group,PI,fill=Group))+geom_boxplot(width=.2)+geom_violin()
mytheme <- theme(plot.title = element_text(face = "bold.italic", size = "14", colour = "brown"), axis.title = element_text(face = "bold.italic", size = "10",color = "blue"), axis.text.x = element_text(face = "bold", size = 8, angle = 45, hjust = 1, vjust = 1), axis.text.y = element_text(face = "bold",size = 8), panel.background = element_rect(fill = "white", color = "black"), panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank(), legend.text = element_text(size = 8), legend.title = element_text(size = 10, face = "bold"), panel.grid.minor.x = element_blank())
p + mytheme
dev.off()
Fst结果
参考链接：
Hardy-Weinberg平衡，Tajima's D ,π和Fst检测,还有XP-CLR - 知乎 (zhihu.com)
网友评论