群体适应性分析——选择信号
●选择信号：生物在自然选择或人工选择过程中，由于选择作用在基因组上留下的遗传多态性降低、连锁不平衡等选择印记。
●选择性清除（selective sweep）：受选择基因座位的遗传多样性随着选择作用的发生而降低，并且由于相邻基因间的高度连锁不平衡，多态性一并降低的过程。
●搭车效应：受到正向选择作用的基因座位不仅本身基因频率会升高，并且会借助连锁不平衡的作用使得附近位点的基因频率一并升高的现象。为去除这种效应，可比较相关基因上、下游各50Kb长的区域，以及邻近基因的FST，pi和XPEHH值。
为了排除这个效应，在对选择信号分析过程中，通常需要使用基于单倍体连锁不平衡的检验方法来校验受选择的位点。
●选择信号方法检测出的功能基因中，不同群体的某方面功能候选基因的非同义突变和同义突变的比值的比较可用于排除选择压力的放松。
搭车效应.png
正选择信号的检测方法
常用群体内检测指标的计算方法大致分为三种：1.基于核苷酸多态性降低的π、θw；2.基于分离位点频率的Tajima’D；3.基于连锁不平衡增加的EHH、iHS。以上三类指标对应于基因组受选择特征的三个维度，而后才有了群体间的选择指标：1.由π衍生的π ratio、ROD、Fst；2.由EHH衍生的XPEHH。
一、基于等位基因频率谱
●基因型频率和基因频率的改变是选择作用在基因组上最直接的体现。基因频谱(site-frequency spectrum)就是指某种等位基因在基因组上某个目标区域内出现的频繁程度。
●符合中性模型的群体，其群体中存在广泛的遗传多态，当突变发生时总能够维持在一个较低的频率，只有当群体基因组上出现或存在有利突变时，选择才会发生作用，从而产生所谓的选择清除或搭车效应。
●代表性的检测方法: Tajima's D, Fu andLi'sD, Fay and Wu'sH, CLR, Hp
●Tajima's D检验的目的是区分随机演变的DNA序列(“中性”)和在非随机过程中演化的DNA序列，包括定向选择或平衡选择。
●Tajima's D的计算原理:多态位点数量和平均非匹配数量的差值。
●D=0时，符合中性假设，群体未受到选择; D<0时，受到定向向选择; D>0时，受到平衡选择。
二、基于连锁不平衡
●基于连锁不平衡理论，位点间的连锁不平衡程度会随标记间距离的增加而逐渐降低。因此，在基因组上可以观察到选择作用造成的不同长度的扩展单倍型纯合(Extended Haplotype Homozygousity)。
●该方法的基本原理是:在中性条件下，基因组很难形成长范围的连锁不平衡的单倍型，因为新突变需要经历漫长的遗传漂变才能达到较高频率，而在漫长的时间里会发生大量基因重组事件，使得这种连锁不断被打破。而当群体处于正向选择作用下时，致因突变及其连锁位点在正选择的作用下，在短时间内会达到较高频率，形成大片段的纯合单倍型。扩展单倍型纯合度检验正是基于这样的特征来筛选受选择基因。
●代表性的检测方法: EHH, XP-EHH, iHS, nSL, OmegaPlus
三、基于群体分化
●同一物种不同群体之间由于环境不同或选择目标不同，其基因组等位基因频率会表现出歧化选择的效应。这种现象在相同基因座位不同等位基因均受到选择时表现尤为明显，即选择加速群体分化。因此，基于群体分化的方法，不同群体同一等位基因频率存在的差异程度大于两个群体处于中性条件下的期望时，就推断该位点存在选择作用。
●代表性的检测方法: Weir and Cockerhan's Fst, LSBL, di
●Fst的取值范围为0-1，1表示群体间完全分化的位点，0表示在群体间完全没有分化的位点。
●基于Fst的的检测方法多采用基因组单位点扫描的策略，而这样的方式容易受到遗传漂变等因素的影响，产生假阳性的显著位点。为尽量减少假阳性的发生，通常采用滑动窗口的策略，降低这些干扰因素，增加选择信号检测的准确性。
四、基于基因组杂合度
●当基因组上特定区域受到选择时，由于“选择性清除”作用的存在，该区域及其连锁的区域表现为多态性降低，同时纯和度增加。因此对基因组的杂合度进行检测，可以推断出基因组中受到选择的区域。基因组上受选择程度越高，则杂合度程度越低。
●代表性的检测方法: θπRatio, ROH
●核苷酸多态性θπ比率越偏离1，受选择程度越高。θπ比率的检测公式如下:θπratio=θπA/θπB
其中，θπA和θπB分别代表A群体和B群体的θπ值。θπ比率大于1， 反映A群体的基因组杂合度大于B群体的杂合度，则B群体相应基因组区域受到选择。θπ 比率小于1,则A群体的基因组杂合度低于B群体，则选择发生在A群体对应的基因组区域。
链接：https://www.jianshu.com/p/7f0e3897f818
选择信号的计算（脚本）
一、FST
●基于等位基因频率的常用检测方法，衡量群体间的分化程度进而推测可能存在的正选择。但是无法判断受选择方向。
●全基因组 FST检验可以通过比较单个位点FST值找到 局部受到选择的基因，但由于一些物种的基因组比 较大，在研究中会使用滑框(大动物常用 50 k 的滑框) 进行 FST 值的计算。在中性突变下，FST 的大小主要受迁移和遗传漂变等因素的影响。但当一个突变受 到人工选择或自然选择时，其频率的升高就会增大 种群的分化水平，FST 值也越大(0＜FST＜1，1 表示该区域种群完全分化)。通过揭示变异的“离群”模式、多样性的缺失以及连锁标记受到选择的影响， 就可以实现对选择痕迹的检测。FST 的计算具有多种方式，目前最常用的方式是通过均方误差
●计算公式：
计算公式
●其中 MSG 是群体内的均方误差，MSP 是群体间的均方误差，nc是校正后整个群体的平均样本大小，此方法主要适用于不同群体之间选择信号的检测。
1 Fst(50K窗口，50K步长)
vcftools --vcf test.vcf \
--weir-fst-pop JBC.txt \
--weir-fst-pop EU.txt \
--out Fst_JBC-EU \
--fst-window-size 50000 \
--fst-window-step 50000
# test.vcf是SNP calling 过滤后生成的vcf文件；
# Fst_JBC-EU生成结果的prefix
# JBC.txt是一个文件包含同一个群体中所有个体，一般每行一个个体。个体名字要和vcf的名字对应。
# EU.txt 包含了群体二中所有个体。
# 计算的窗口是50kb，而步长是50kb （根据你的需求可以作出调整）。
2 根据fst值排序
sort -k 5 -g Fst_JBC-EU.windowed.weir.fst > Fst_JBC-EU.windowed.weir.sorted.fst
3 提取前5%
tail -n 2487 Fst_JBC-EU.windowed.weir.sorted.fst > Fst_JBC-EU.sort.5%
二、XPEHH分析
●基于连锁不平衡和单倍型
●计算公式：
XPEHH.png
●XPEHH值具有方向性，正的较大XPEHH值表明A群体可能在该位点发生了选择事件，负的较大XPEHH值表明B群体可能在改位点发生了选择事件。
●只能筛选正选择信号, xpehh>2表示目标群体的受选择方向，xpehh<-2表示参考群体的受选择方向
前期准备：
1 beagle填充，还有所谓的phasing
java -jar -Xmn12G -Xms24G -Xmx48G /home/software/beagle.25Nov19.28d.jar gt=QC.JBC_EU.filter-geno005-maf005.vcf out=QC.JBC_EU.filter-geno005-maf005.beagle ne=82
2 提取群体样品
vcftools --gzvcf QC.JBC_EU.filter-geno005-maf005.beagle.vcf.gz \
--keep JBC.txt \
--recode \
--recode-INFO-all \
--out JBC.beagle
vcftools --gzvcf QC.JBC_EU.filter-geno005-maf005.beagle.vcf.gz \
--keep EU.txt \
--recode \
--recode-INFO-all \
--out EU.beagle
由于是群体间的，因此需要将两个群体的vcf文件分别提取出来，后面只用其中一个做展示，另一个只是修改输入输出文件名即可
！！！注意：第3-5步只用其中一个群体计算遗传距离即可，是通用的
开始分析：
1 拆分染色体(1_JBC.chr.sh )
for i in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do vcftools --vcf JBC.beagle.recode.vcf \
--recode \
--recode-INFO-all \
--chr ${i} \
--out JBC.chr${i};
done
2 准备分染色体的map文件(2_map.sh)
用的上一步的生成文件，也就是每个染色体的vcf的单倍型文件
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do vcftools --vcf JBC.chr${k}.recode.vcf \
--plink \
--out chr${k}.MT;
done
3 map计算遗传距离(3_map.distance.sh)
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do awk 'BEGIN{OFS=" "} {print 1,".",$4/1000000,$4}' chr${k}.MT.map > chr${k}.MT.map.distance;
done
4 计算XPEHH(4_xpehh.sh)
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do /home/software/selscan/bin/linux/selscan --xpehh \
--vcf JBC.chr${k}.recode.vcf \
--vcf-ref EU.chr${k}.recode.vcf \
--map chr${k}.MT.map.distance \
--threads 10 \
--out chr${k}.EU_JBC;
done
得到 chr${k}.EU_JBC文件29个
重测序数据直接计算XPEHH（对每个100K窗口单倍型进行打分，标准化使其正态分布，并统计值>2的SNP所占的比例） 相对芯片来说步骤少点
最终得到的标准化的xpehh值表示的是xpehh>2的SNP在此区间所占的比例
具有方向性，xpehh>2表示目标群体的受选择方向，xpehh<-2表示参考群体的受选择方向
5 第一列改为染色体号(5_add.chr.sh)
并将空格改为tab键分割
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do awk '{print '${k}',$2,$3,$4,$5,$6,$7,$8}' chr${k}.EU_JBC.xpehh.out > Chr${k}.EU_JBC.xpehh.out;
sed -i 's/ /\t/g' Chr${k}.EU_JBC.xpehh.out;
done
6 加50kb窗口且标准化(6_norm.xpehh.sh)
可按需求修改窗口大小
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do /home/software/selscan/bin/linux/norm --xpehh \
--files Chr${k}.EU_JBC.xpehh.out \
--bp-win --winsize 50000;
done
产生两种文件.100bins.norm文件和.100bins.norm.50kb.windows文件
.100bins.norm.50kb.windows文件的表头为
<win start> <win end> <scores in win> <gt the fraction of XP-EHH scores >2> <lt the fraction of XP-EHH scores < -2> <approx percentile for gt threshold wins> <approx percentile for lt threshold wins> <max score> <min score>
7 提取并合并结果文件(7_extract.merge.all.sh)
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do awk '{print '${k}',$1,$2,$4,$5,$8,$9}' Chr${k}.EU_JBC.xpehh.out.norm.50kb.windows > Chr${k}.EU_JBC.chart.xpehh.norm.50kb;
cat ./*.EU_JBC.chart.xpehh.norm.50kb > all.xpehh.norm.50kb;
done
三、XPCLR分析
●XP-CLR 是陈华老师、Nick Patterson 和 David Reich 在 2010 年发表的方法，全称叫 the cross-population composite likelihood ratio test（跨群体复合似然比检验），是一种是基于选择扫荡（selective sweeep）的似然方法。
●选择扫荡可以增加群体之间的遗传分化，导致等位基因频率偏离中性条件下的预期值。XP-CLR 利用了两个群体之间的多基因座等位基因频率差异建立模型，使用布朗运动来模拟中性下的遗传漂移，并使用确定性模型来近似地对附近的单核苷酸多态性（SNPs）进行选择性扫描。
●该方法利用已有的数据集通过极大似然法估计等位基因频率等群体参数，然后预测等位基因频率在中性模型下的“失真”程度进而判断是否有选择发生。极大似然法估计群体参数的一个明显缺点就是过度依赖现有的数据集，可能会造成假阳性的选择信号，作者利用成对SNP连锁不平衡的方法，降低了可能存在连锁的SNP的权重。
1 拆分染色体(1_JBC_EU.chr.sh)
for i in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do vcftools --vcf QC.filter-geno005-maf005.beagle.vcf.gz \
--recode \
--recode-INFO-all \
--chr ${i} \
--out QC.chr${i};
done
2准备分染色体的map文件(2_map.sh)
用的上一步的生成文件
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do vcftools --vcf QC.chr${k}.recode.vcf \
--plink --out chr${k}.MT;
done
3 map计算遗传距离(3_map.distance.sh)
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do awk 'BEGIN{OFS=" "} {print 1,".",$4/1000000,$4}' chr${k}.MT.map > chr${k}.MT.map.distance;
done
4 计算XPCLR(4_xpclr.sh)
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do xpclr --out chr${k} \
--format vcf \
--input QC.JBC_EU.chr${k}.recode.vcf \
--samplesA EU.txt \
--samplesB JBC.txt \
--map chr${k}.MT.map.distance \
--chr ${k} \
--gdistkey None \
--phased \
--size 50000 \
--step 50000;
done
单个染色体做：
xpclr --out chr9 \
--format vcf \
--input QC.chr9.recode.vcf \
--samplesA EU.txt \
--samplesB JBC.txt \
--map chr9.MT.map.distance \
--chr 9 \
--gdistkey None \
--size 50000 --step 25000
A为参考群体，B为目标群体
--out: 输出文件
--format： 输入格式，vcf，hdf5，zarr，txt（对应2种基因型，和一个snp位点文件）
--input：输入vcf，或者hdf5， 选择txt时，不选择，
--samplesA： 样本A名称文件； txt时，不选择
--samplesB：样本B名称文件； txt时，不选择
--map: 基因位点文件
--popA： 样本A基因型文件，txt时使用
--popB： 样本B基因型文件，txt时使用
--chr： 染色体，和输入文件一致即可，每次分析一个
--ld： LD值，进行权重分析
--maxsnps： 一个窗口最大的SNP
--minsnps： 一个窗口最小的SNP
--size： 窗口大小
--step： 步长
5 提取并合并结果文件(5_extract.merge.all.sh)
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do awk '{print $2,$3,$4,$12,$13}' Chr${k} > Chr${k}.EU_JBC.chart.xpclr.50kb.windows;
cat ./*.EU_JBC.chart.xpclr.50kb.windows > all.xpclr.50kb.windows;
done
6 结果处理
没有norm_xpclr值的按照0处理（个人认为）,放在excel中处理
四、iHS分析
●基于连锁不平衡和单倍型
●
iHS.png
●iHS=0时，表明该位点处于平衡选择状态。
●当 iHS 值为较大的正值时，表示单倍型可能与祖先一致，当 iHS 值为较大的负值时，表示可能具有新衍生的单倍型。
●重测序数据进行计算的时候，最后标准化后产生的iHS表示的是每个窗口内|iHS| > 2的比例，因此最高不超过1。
●优势与局限：
iHS 检验 是基于祖先等位基因的一种检验方法，其优势是通 过祖先的信息能更准确地观察到迚化中受选择的情 况。而其劣势是现实中很难获得祖先的群体信息， 因为很多古老的品种已经灭绝或野生品种很难获得。
beagle以及分染色体、计算遗传距离与XPEHH相同
1拆分染色体(1_chr.sh)
2准备分染色体的map文件(2_map.sh)
3 map计算遗传距离(3_map.distance.sh)
4 计算iHS(4_ihs.sh)
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do /home/software/selscan/bin/linux/selscan --ihs \
--vcf chr${k}.vcf \
--map chr${k}.MT.map.distance \
--out chr${k}.iHS;
done
得到 ${k}.iHS.ihs.out文件29个
重测序数据直接计算iHS（对每个100K窗口单倍型进行打分，标准化使其正态分布，并统计绝对值>2的SNP所占的比例） 相对芯片来说步骤少点
最终得到的ihs值表示的是|iHS|>2的SNP在此区间所占的比例
5 第一列改为染色体号(5_add.chr.sh)
并将空格改为tab键分割
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do awk '{print '${k}',$2,$3,$4,$5,$6}' chr${k}.iHS.ihs.out > Chr${k}.ihs.out ;
sed -i 's/ /\t/g' Chr${k}.ihs.out;
done
6 加50kb窗口且标准化(6_norm.ihs.sh)
可按需求修改窗口大小
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do /home/software/selscan/bin/linux/norm --ihs \
--files Chr${k}.ihs.out \
--bp-win --winsize 50000;
done
产生两种文件.100bins.norm文件和.100bins.norm.50kb.windows文件
7 提取并合并结果文件(7_extract.merge.all.sh)
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do awk '{print '${k}',$1,$2,$4}' Chr${k}.ihs.out.100bins.norm.50kb.windows > Chr${k}.chart.ihs.out.50kb.windows;
cat ./*.chart.ihs.out.500kb.windows > all.ihs.out.50kb.windows;
done
按照第四列排序
sort -k 4n,4 all.ihs.out.50kb.windows > all.ihs.out.50kb.windows.sort
筛选出绝对值大的前5%注释
五、CLR分析(SweeD软件)
1拆分染色体(1_chr.sh)
nohup bash chr.sh &
for i in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do vcftools --vcf QC.filter-geno005-maf005.vcf \
--recode \
--recode-INFO-all \
--chr ${i} \
--stdout > QC.JBC.filter-geno005-maf005.chr${i}.vcf;
done
2 CLR(2_CLR.sh)
需要chr_grid.txt文件，第一列为染色体名称，第二列为每条染色体长度除以要设置的窗口长度
1 1580
2 1360
3 1200
sed -i 's/\t/,/g' chr_grid.txt
##############
for line in `cat chr_grid.txt`
do
chr=${line%%,*}
grid=${line##*,}
/home/ywf/software/sweed/SweeD -name ${chr}.100kb \
-input QC.JBC.filter-geno005-maf005.chr${chr}.vcf \
-grid ${grid} \
-minsnps 200 \
-maf 0.05 \
-missing 0.1;
done
-name 指定程序运行名和输出结果文件的后缀
-input 输入需要分析的数据文件
-sampleList 指定参与分析的样本文件
-grid 每条染色体应分割成多少个网格，数值越大，分割的片段越多
-minsnps 参与计算的窗口中最少的SNP数量
-maf 小于该频率的次等位基因将被过滤掉
-missing 缺失率小于该值的位点将被过滤
3对结果第一列加上染色体(3_add.chr.delet.merge.sh)
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29;
do awk '{print '${k}',$1,$2,$3,$4,$5}' SweeD_Report.${k}.100kb > SweeD_Report.chr.${k}.100kb;
sed -i 's/ /\t/g' SweeD_Report.chr.${k}.100kb;
# 删除结果中前三行
sed -i '1,3d' SweeD_Report.chr.*.100kb;
# 合并结果文件
cat ./SweeD_Report.chr.*.100kb > all.CLR.100k.txt;
done
其中，第二列为位置，第三列为似然法计算出的值
4 按照第三列排序，取前5%，并对其位置进行取整（使用excel表格处理）
24710*0.05 = 1236
取整：=ROUND(B1,0)
注意：注释前请换染色体号
加50k窗口
