Variants calling

【群体遗传】选择信号分析-Fst、Pi

2022-10-11  本文已影响0人  Making

1.选择信号分析—Fst选择信号分析

概念
Fst= (Ht-Hs）/ Ht
Hs:亚群体中的平均杂合度
Ht:复合群体中的平均杂合度

Fst值

理论上计算Fst的步骤
理论上要估算FST，需要以下步骤：

链接：Fst的计算原理与实战 - 简书 (jianshu.com)

image.png

两篇少样本数量不影响Fst文献
两篇少样本数量不影响Fst文献

FST计算方法

按照计算方法可分为：按照SNP单点计算和滑动窗口模式计算，一般使用的是滑动窗模式

vcftools计算Fst实战

SNP单点计算
##对每一个SNP变异位点进行计算
vcftools --vcf sample.vcf --weir-fst-pop population_1.txt --weir-fst-pop population_2.txt --out sample_1_2

       

滑动窗口模式计算
#按照窗口模式计算（更准确）
vcftools --vcf sample.vcf --weir-fst-pop population_1.txt --weir-fst-pop population_2.txt --out P_1_2 --fst-window-size 500000 --fst-window-step 50000
# sample.vcf是SNP calling 或芯片数据过滤后生成的vcf 文件；
# p_1_2_3 生成结果的prefix
# population_1.txt是一个文件包含同一个群体中所有个体，一般每行一个个体。个体名字要和vcf的名字对应。
# population_2.txt 包含了群体二中所有个体。
##--fst-window-size  # 设置计算Fst的窗口大小，根据自己的数据进行设置，看看别人文章里怎么用的
##--fst-window-step  # 设置计算Fst的步长长度，根据自己的数据进行设置
#计算的窗口是500kb，而步长是50kb （根据你的需其可以作出调整）。我们也可以只计算每个点的Fst，去掉参数（--fst-window-size 500000 --fst-window-step 50000）即可。
#注：建议参考已发表文献，根据相关物种的在基因组单位区段面积的位点数目、LD平均值或在r2=0.1的时候位置来定窗口大小

#根据加权fst为Fst排序
sort -k 5 -g p1_p2_window.windowed.weir.fst > p1_p2.sorted.fst
#窗口计数
wc -l p1_p2.sorted.fst #假设为20000
#取前1%
tail -n 200 p1_p2.sorted.fst > p1_p2.sorted.1%.fst
#找出前1%中最小的加权fst值

参数说明：
--vcf 输入vcf文件
--weir-fst-pop 输入群体的群体ID名，该文件必须是txt格式，每个ID占一行。
--fst-window-size 500000 --fst-window-step 50000 ，这里窗口设置为500kb，步长设置为50kb。根据情况调整窗口大小。

image.png
(18条消息) 【群体遗传】Fst（群体间分化指数）_陈有朴的博客-CSDN博客_fst遗传分化指数 
例如：
Sample1
Sample2
Sample3

分别对不同结果进行图形绘制

##图1
library(ggplot2)
data<-read.table("test1.out.windowed.weir.fst",header=T)
sc3 = subset(data,CHROM=="Gm01")
p <- ggplot(sc3,aes(x=BIN_END/1000000,y=WEIGHTED_FST)) + geom_point(size=0.5, colour="blue") + xlab("Physical distance (Mb)")+ ylab("Fst") + ylim(-1,1)
p + theme_bw()

##图2
library(ggplot2)
data<-read.table("test.out.weir.fst",header=T)
sc3 = subset(data,CHROM=="Gm01")
p <- ggplot(sc3,aes(x=POS,y=WEIR_AND_COCKERHAM_FST)) + geom_point(size=0.5, colour="blue") + xlab("Physical distance (Mb)")+ ylab("Fst") + ylim(-1,1)
p + theme_bw()

计算haploPS、XP-EHH、 Fst，正向选择分析方法寻找性状相关的位点
<meta charset="utf-8">

π的计算

π，核苷酸多样性，越大说明核苷酸多样性越高，越低说明两个座位DNA序列差异越小。
vcftools

LD-blockde 计算

计算和可视化，参考

分析方法 参考文献：中国农业科学 scientific reports

找到重测序的数据，基因分型，找到单倍型，call SNP,过滤，注释snp,可视化，分别计算haploPS,XP-EHH,Fst.求交集，公共基因进行注释，富集分析。

2.计算总体PI值和TajimaD

核酸多样性PI，值越大说明核苷酸多样性越高
TajimaD用于检验DNA序列在演化国产过程中是否遵循中性演化模型。D>0:平衡选择，突然群体收缩；D<0:最近的选择性清除，最近瓶颈后的群体扩张，与消除基因连锁；D=0：群体根据突变-漂移平衡演变，没有选择。

vcftools --vcf Filter.snp.vcf --window-pi 500000 --out total
vcftools --vcf Filter.snp.vcf --TajimaD 500000 --out total

计算每个亚群的PI值和TajimaD值

#创建bash脚本
vi test.sh
#写入以下内容：
   for i in {1..3};do
   #根据ID，从vcf文件中提取每个亚群的信息
   vcftools --vcf Filter.snp.vcf --keep ${i}-population.txt --recode -- 
   recode-INFO-all --out p${i}
   #根据提取的信息计算每个亚群的PI值
   vcftools --vcf p${i}.recode.vcf --out q${i}_pi_500kb --window-pi 
   500000 --window-pi-step 50000
   #根据提取的信息计算每个亚群的TajimaD值
   vcftools --vcf p${i}.recode.vcf --out q${i}_500kb.TajimaD --TajimaD 
   500000
   done
#给bash脚本添加可执行权限
chmod +x test.sh
#运行脚本
./test.sh

根据pi值画箱线图和小提琴图

#添加分组信息G1、G2、G3
data=data.frame(q1_500kb_pi_windowed)
data
Group=c("G1")
data1=data.frame(data,Group)
data=data.frame(q2_500kb_pi_windowed)
data
Group=c("G2")
data2=data.frame(data,Group)
data=data.frame(q3_500kb_pi_windowed)
data
Group=c("G3")
data3=data.frame(data,Group)
#将三组数据合并
library(dplyr)
total_data<-dplyr::bind_rows(data1,data2,data3)
#画箱线图+小提琴图
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

根据PI值画折线图

#输入数据时，每列必须是数值型数据，尤其是染色体那一列
data=data.frame(q1_500kb_pi_windowed)
Group=c("G1")
data1=data.frame(data,Group)
data=data.frame(q2_500kb_pi_windowed)
Group=c("G2")
data2=data.frame(data,Group)
data=data.frame(q3_500kb_pi_windowed)
Group=c("G3")
data3=data.frame(data,Group)
#将三组数据合并
library(dplyr)
library(ggplot2)
total_data<-dplyr::bind_rows(data1,data2,data3)
for (i in 1:7){ 
  result_name = paste("chr",i,".png",sep="") 
  png(result_name,width=1500,height=300) 
  chrom = subset(total_data,CHROM==i) 
  xlab = paste("Chromosome",i,"(MB)",sep="") 
  p <- ggplot(chrom,aes(x=BIN_END/1000000,y=PI,color=Group,shape=Group)) + geom_line(size=0.5) + xlab(xlab)+ ylab("Pi") + theme_bw() 
  print(p) 
  dev.off() 
  }

链接：https://www.jianshu.com/p/4dbde8533607

参考：https://www.jianshu.com/p/c520c36fe340
https://www.jianshu.com/p/b73a8d6233be
https://blog.csdn.net/u014182497/article/details/52672308
http://www.omicshare.com/forum/thread-3688-1-1.html
群体中的Fst值-学习篇 - 简书 (jianshu.com)
https://www.jianshu.com/p/4dbde8533607

上一篇 下一篇

猜你喜欢

热点阅读