群体遗传进化与选择分析：Pi、Fst、TajimaD、XP-CLR计算方法

生信分析笔记

在遗传学中，群体指的是一组具有共同遗传特征的个体，而个体则是指单个生物体。群体中的个体之间可以存在遗传交流和基因流动，这会导致群体中的基因频率发生变化。今天分享的笔记是群体进化与选择分析，包括Pi、Fst、TajimaD、XP-CLR的介绍和计算方法。

首先，咱们都知道时间不会停止，也就意味着历史的车轮不会停止，自然界一直在不断地演化，不管是动物还是植物，都在不停的选择和分化。例如玉米小麦等植物，在很久以前可能发源于杂草，后来经过人的驯化，才改良为现在适宜栽培的品种。

什么是正选择？

正选择可以用自然选择来解释：假如一个基因或位点能够使个体有着更强的生存力，这样就会使个体的后代更多，如此一来，这个基因或位点在群体中就越来越多。

举个例子：以前非洲大草原有短颈鹿，后来偶然的突变导致长颈鹿产生，由于长颈鹿能吃到更多的食物，有着更高的存活率，所以导致这个突变受到正选择。

负选择

如果群体中的某个个体出现一个致命的突变，从而使自己或者是后代从群体中被淘汰，这也导致群体中该位点的多态性的降低。

举个例子：假如正常生长中的玉米群体中偶然发生了一个突变，导致水分吸收受阻，这种缺陷导致后代被淘汰，因此该突变位点的多态性会降低。

平衡选择

平衡选择指多个等位基因在一个群体的基因库中以高于遗传漂变预期的频率被保留，如杂合子优势。

核酸多样性Pi

Pi指的是核苷酸多样性，Pi值越大说明核苷酸多样性越高。通常用于衡量群体内的核苷酸多样性，也可以用来推演进化关系，可以理解成先在群体内两两求Pi，再计算群体的均值,常用软件是vcftools。

vcftools --vcf input.vcf --window-pi 200000 --window-pi-step 100000 --keep 1.sample.list --out pi_window_1.sample.list
# 检查文件的行数
wc -l pi_window_1.sample.list.windowed.pi

批量计算Pi的脚本：

#!/bin/bash
#定义所有以.txt结尾的sample.list文件
sample_files=(*.txt)
#循环执行命令
for sample in "${sample_files[@]}";do
    #生成输出文件名
    output="pi_window_${sample}"

    #输出调试信息
    echo "Running vcftools for $sample,output file: $output"
    #确认文件存在
    if [ ! -f "$sample" ];then
        echo "Error:File $sample does not exist"
        continue
    fi
    #运行vcftools命令
    vcftools --vcf input.vcf --window-pi 200000 --window-pi-step 100000 --keep "$sample" --out "$output"
done

利用ggplot2画图的脚本：

# 读取数据
pi_window_1 <- read.table("pi_window_1.sample.list.windowed.pi",header = TRUE,sep = "\t")
pi_window_2 <- read.table("pi_window_2.sample.list.windowed.pi",header = TRUE,sep = "\t")
pi_window_3 <- read.table("pi_window_3.sample.list.windowed.pi",header = TRUE,sep = "\t")
pi_window_1$Group <- 1
pi_window_2$Group <- 2
pi_window_3$Group <- 3
df_plot_all <- rbind(pi_window_1,pi_window_2,pi_window_3)

# 绘制结果
library(ggplot2)
p1 <- ggplot(df_plot_all)+
  geom_smooth(aes(BIN_START,PI,color=as.factor(Group),group=Group),method = "loess",span = 0.1,se = F,linewidth = 1)+
  # geom_smooth(aes(BIN_START,PI,color=as.factor(Group),group=Group),method = "gam",span = 0.1,se = F,linewidth = 1)+
  scale_color_manual(values = c("red","green","blue"))+
  xlab("Physical Position")+
  ylab("Pi")+
  theme_bw()+
  theme(legend.position = "right")

# 输出图片
ggsave(filename="pi_200kwindow_100kstep.pdf",plot = p1,width = 16,height = 9)
p1

群体内选择检验TajimaD

Tajima's D是日本学者Tajima Fumio 1989年提出的一种统计检验方法，用于检验DNA序列在演化过程中是否遵循中性演化模型。

D > 0: 平衡选择，突然收缩，稀有等位基因以低频率存在。 D < 0: 经历瓶颈效应，随后群体扩张，稀有等位基因以高频率存在。 D = 0: 平衡演变，没有太大频率差异。

计算方法：

vcftools --vcf input.vcf --TajimaD 100000 --keep 1.sample.list --out TajimaD_1.sample.list

批量计算TajimaD的脚本：

#!/bin/bash
#定义所有以.txt结尾的sample.list文件
sample_files=(*.txt)
#循环执行命令
for sample in "${sample_files[@]}";do
    #生成输出文件名
    output="TajimaD_${sample}"
    #输出调试信息
    echo "Running vcftools for $sample,output file: $output"
    #确认文件存在
    if [ ! -f "$sample" ];then
        echo "Error: File $sample does not exist"
        continue
    fi
    #运行vcftools 命令
    vcftools --vcf input.vcf --TajimaD 100000 --keep "$sample" --out "$output"
done

群体间分歧度Fst

Fst叫固定分化指数，用于估计亚群间平均多态性大小与整个种群平均多态性大小的差异，反映的是群体结构的变化。Fst的取值范围是[0,1]

当Fst=1时，表明亚群间有着明显的种群分化，值越高表示分化程度越高。在中性进化条件下，Fst的大小主要取决于遗传漂变和迁移等因素的影响。

假设种群中的某个等位基因对特定环境的适应度较高而经历适应性选择，那该基因的频率在种群中会升高，种群的分化水平增大，群体Fst升高。

计算方法：

 vcftools --vcf all_2229.vcf --fst-window-size 200000 --fst-window-step 100000 --weir-fst-pop 1.sample.list --weir-fst-pop 2.sample.list --out fst.1.sample.list.2.sample.list

批量计算Fst的脚本：

#!/bin/bash
#定义所有的sample.list文件
sample_files=("W.sample.list.txt" "S.sample.list.txt")
#获取样本文件的数量
num_samples=${#sample_files[@]}
#确保有足够的文件进行成对比较
if ((num_samples < 2));then
    echo "Error: Not enough sample files.At least two are required."
    exit 1
fi
#循环执行命令，两两比较sample.list文件
for ((i=0;i<num_samples; i++ ));do
    for ((j=i+1;j<num_samples; j++ ));do
        sample1="${sample_files[i]}"
        sample2="${sample_files[j]}"
        #生成输出文件名
        output="fst.${sample1%.sample.list.txt}.${sample2%.sample.list.txt}"
        #输出调试信息
        echo "Running vcftools for $sample1 and $sample2,output file: $output"
        #确认文件存在
        if [ ! -f "$sample1"];then
            echo "Error: File $sample1 does not exist"
            continue
        fi
        if [ ! -f "$sample2" ];then
            echo "Error: File $sample2 does not exist"
            continue
        fi
        #运行vcftools命令
        vcftools --vcf input.vcf --fst-window-size 2000 --fst-window-step 1000 --weir-fst-pop "$sample1" --weir-fst-pop "$$
    done
done

跨群体复合似然比检验XP-CLR

XP-CLR 是陈华老师等在 2010 年开发的方法，全称叫 the cross-population composite likelihood ratio test（跨群体复合似然比检验），是一种是基于选择扫荡（selective sweeep）的似然方法。

选择扫荡可以增加群体之间的遗传分化，导致等位基因频率偏离中性条件下的预期值。XP-CLR利用两个群体之间的多基因座等位基因频率差异（multilocus allele frequency differentiation）建立模型，使用布朗运动来模拟中性下的遗传漂移，并使用确定性模型来近似地对附近的单核苷酸多态性（SNPs）进行选择性扫描，以下是计算方法：

#!/bin/bash
# 定义样本文件数组
sample_files=("1.sample.list.txt" "2.sample.list.txt" "3.sample.list.txt")
# 获取所有染色体编号
chromosomes=$(bcftools view -h input.vcf | grep '^##contig' | cut -d '=' -f3 | cut -d ',' -f1)
echo "Chromosomes: $chromosomes"
# 定义输出目录
output_dir="./out"
mkdir -p $output_dir
# 循环遍历样本文件两两组合
for i in 0 1 2; do
    for j in $(seq $((i+1)) 2); do
        samplesA=${sample_files[$i]}
        samplesB=${sample_files[$j]}
        # 获取无后缀样本名
        samplesA_name=$(basename "${samplesA%.sample.list.txt}")
        samplesB_name=$(basename "${samplesB%.sample.list.txt}")

        # 定义组合编号
        combo="${samplesA_name}${samplesB_name}"   
        # 定义组合输出文件
        output_file="${output_dir}/chr_all_${combo}.txt"
        # 循环遍历每个染色体，并运行 xpclr，输出到临时文件
        for chr in $chromosomes
        do
            echo "Running XPCLR for chromosome $chr with samples ${samplesA_name} and ${samplesB_name}..."
            temp_output="${output_dir}/temp_${chr}_${combo}.txt"
            xpclr --out "$temp_output" --format vcf --input input.vcf --samplesA $samplesA --samplesB $samplesB --size 200000 --step 100000 --chr $chr
            # 将临时文件的内容追加到组合输出文件
            cat "$temp_output" >> "$output_file"
            # 删除临时文件
            rm "$temp_output"
        done
    done
done
echo "XPCLR analysis completed for all combinations and chromosomes."

以上就是今天分享的全部内容，如果对您有所帮助，欢迎转发收藏。

微信扫一扫分享文章