Tajima D
這個(gè)是選擇相關(guān)的一個(gè)參數(shù)，大于0代表群體觀測(cè)雜合度高于預(yù)期雜合度秦陋，稀有等位基因頻率降低（群體收縮或者平衡選擇）蔓彩，小于0說(shuō)明群體觀測(cè)雜合位點(diǎn)少于預(yù)期值，稀有等位基因頻率增加（群體擴(kuò)張或者低頻選擇）驳概。 也就是說(shuō)赤嚼，只有0是正常的，其他都是選擇發(fā)生顺又。
π
π更卒，核苷酸多樣性，越大說(shuō)明核苷酸多樣性越高待榔，越低說(shuō)明兩個(gè)座位DNA序列差異越小逞壁。
Fst
Fst, 分化系數(shù)流济，從0到1說(shuō)明親緣關(guān)系越來(lái)越遠(yuǎn)。接近于0說(shuō)明兩個(gè)個(gè)體親緣關(guān)系近腌闯，接近1說(shuō)明親緣關(guān)系遠(yuǎn)绳瘟。
Hardy-Weinberg
Hardy-Weinberg平衡檢測(cè)，這個(gè)主要是檢測(cè)基因型頻率是否等于基因頻率乘積姿骏。比如A:0.3,a:0.7那么Aa的頻率是否為0.42
####計(jì)算pi
vcftools --vcf merged.vcf --window-pi 500000 --out pi
####計(jì)算TajimaD
vcftools --vcf marged.vcf --TajimaD 500000 --out TajimaD
###計(jì)算HW
vcftools --vcf merged.vcf --hardy --out HW
###fst需要兩個(gè)以上pop
vcftools --vcf SNP.vcf --weir-fst-pop 1.txt --weir-fst-pop 2.txt --out 1_VS_2 --fst-window-size 500000
### 將亞群拆分出來(lái)
bcftools view -S ARO merged.vcf -O v -o ARO.vcf
繪圖：
library(ggplot2)
###讀取相應(yīng)文件即可
mydata<-read.table("pi.windowed.pi",header = T)
mydata1<-na.omit(mydata)
a<-mydata1[sample(nrow(mydata1), 1000), ] #隨機(jī)選取其中的1000行數(shù)據(jù)糖声，這個(gè)Hardy-Weinberg平衡文件需要的，因?yàn)榻Y(jié)果太多分瘦，無(wú)法利用ggplot繪出全部結(jié)果
###畫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())
結(jié)果：
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組數(shù)據(jù)合并
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()
pi結(jié)果1
pi結(jié)果2
Fst結(jié)果
參考鏈接：
Hardy-Weinberg平衡去团，Tajima's D ,π和Fst檢測(cè),還有XP-CLR - 知乎 (zhihu.com)
