美文网首页Variants calling
【群体遗传】选择信号分析-Fst、Pi

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

作者: Making | 来源:发表于2022-10-11 12:48 被阅读0次

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

    概念
    • Fst:群体间遗传分化指数,是种群分化和遗传距离的一种衡量方法,分化指数越大,差异越大。适用于亚群体间多样性的比较。

    • 用于衡量种群分化程度,取值从0到1,为0则认为两个种群间是随机交配的,基因型完全相似;为1则表示是完全隔离的,完全不相似。它往往从基因的多样性来估计,比如SNP或者microsatellites(串联重复序列一种,长度小于等于10bp)。是一种以哈温平衡为前提的种群遗传学统计方法。

    • 两个种群之间遗传差异的基本测量是统计量FST。在遗传学中,F一词通常代表“近亲繁殖”,它倾向于减少群体中的遗传变异。遗传变异可以用杂合度来衡量,所以F一般表示群体中杂合性的减少。 FST是与它们所属的总群体相比,亚群体中杂合性的减少量。
      具体可以下面的公式表示:

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

    Fst值

    理论上计算Fst的步骤
    理论上要估算FST,需要以下步骤:

    • 找出每个亚群的等位基因频率。
    • 查找复合群体的平均等位基因频率
    • 计算每个亚群的杂合度(2pq)
    • 计算这些亚群杂合度的平均值,这是HS。
    • 根据总体等位基因频率计算杂合度,这是HT。
    • 最后,计算FST =(HT-HS)/ HT

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

    image.png

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

    FST计算方法

    • Fst值的取值范围是[0,1],最大值为1表明两个群体完全分化,最小值为0表明群体间无分化。
    • 实际使用FST<0--0.05,表示群体分化很小;0.05--0.15,中等程度的分化;0.15--0.25,较大的分化;0.25以上,分化很大。
    • 其实Fst分析就是看两个群体之间分化程度的一种方法,Fst值越大(越接近1)表明两个群体间分化程度越高,亲缘关系越远;Fst值越小(越接近0)表明群体间分化程度越低,亲缘关系越近。
    • 由于负值Fst无意义,所以通常计算出的负值Fst需要将其替换为0

    按照计算方法可分为:按照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遗传分化指数
    • population_2.txt和population_2.txt格式一样,只有一列样品信息,个体名字要和vcf的名字对应**
    例如:
    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

    相关文章

      网友评论

        本文标题:【群体遗传】选择信号分析-Fst、Pi

        本文链接:https://www.haomeiwen.com/subject/ghctzrtx.html