悄悄的题外话：上次说到的使用gatk提取变异，由于个体太多慢的不行（gatk4）,中途换成了gatk3，但问题不大，不过还是建议使用gatk4。

书归正题，我们这次一起来学习一下受选择基因的筛选和富集。做这一方面有很多方法，比如结合Fst,pi,tajimaD值筛选；HKA,PBS筛选。XPEHH啥啥的。当然啦，我最开始的时候使用的方法是Fst和tajimD来筛选的，具体可以参考这篇简书：（呃 找不到了 找到了我再发出来）。我们这次使用的方法是Fst与pi值结合筛选。话不多说，让我们愉快的开始吧！愉快的开始吧！快来开始吧！

筛选

既然使用的是Fst与pi，当然要计算它们啦，用vcftools

#!/bin/bash
v=sample.vcf
pop1=1.pop
pop2=2.pop
vcftools --vcf $v --weir-fst-pop $pop1 --weir-fst-pop $pop2 --fst-window-size 10000 --fst-window-step 2000 --out out
vcftools --vcf $v --window-pi 10000 --window-pi-step 2000 --keep $pop1 --out 1
vcftools --vcf $v --window-pi 10000 --window-pi-step 2000 --keep $pop2 --out 2

生成的文件后缀是 (为了方便后文不加后缀了)
.windowed.weir.fst 
.windowed.pi
窗口及步长分别是10000和2000

第二步，计算pi_ratio(即群体之前pi的比值)

perl merge2pi_ratio.pl 1 2 > 1_2.ratio
perl merge2pi_ratio.pl 2 1 > 2_1.ratio

pi_ratio

结果文件长这样：

ratio

第三步，将Fst和ratio结果整合，取Fst top10%, ratio10%

perl top_Fst_ratio_S.pl out 1_2.ratio 0.1 0.1 1.table 1.stat
perl top_Fst_ratio_S.pl out 2_1.ratio 0.1 0.1 2.table 2.stat
注：第一行得到的是2受选择，第二行得到的是1受选择

第四步 提取top对应的窗口

awk '$NF=="top"' 1.table > 2.top
bedtools merge -i 2.top > 2.top.merge
awk '$NF=="top"' 2.table > 1.top
bedtools merge -i 1.top > 1.top.merge

提取基因位置信息

awk '$3=="gene"' sample.gff > gene.gff

最后，生成我们需要的基因信息

bedtools  intersect -wo -F 0.1 -a 1.top.merge -b gene.gff | awk -F "\t"  '{print $4"\t"$7"\t"$8"\t"$10"\t"$12}' | sort -u  > 1.top.gene
bedtools  intersect -wo -F 0.1 -a 2.top.merge -b gene.gff | awk -F "\t"  '{print $4"\t"$7"\t"$8"\t"$10"\t"$12}' | sort -u  > 2.top.gene

可以看一下我们得到的结果：

output

因为gff文件的问题，可以看到基因的ID奇奇怪怪的，不能用这个ID去做富集，这时候我们就要把ID转换为我们通常的ID，我这里的做法呢是把基因对应的蛋白序列去和拟南芥进行比对，得到我们这个物种的ID所对应的拟南芥中基因的ID。

blast

这是我下载的示例物种的protein文件,同时也要下载拟南芥的蛋白信息(tair.fa)

fa

>后面的信息太繁杂了，只需要ID就好，因此做下处理

zcat protein.fa.gz | sed 's/.\+locus=\(.\+\)ID.\+/>\1/g' > pro.fa

拟南芥蛋白信息同理（根据实际情况处理哦）
提取受选择基因ID

cat 1.top.gene | awk '{print $5}' | awk -F ";" '{print $2}'| awk -F "=" '{print$2}' > 1.id

根据ID提取蛋白序列

perl stkfa_id.perl 1.id pro.fa > 1.fa

接下来就要运行blast啦
构建数据库

makeblastdb -in tair.fa -dbtype prot -out ta_db -parse_seqids -hash_index

比对

blastp -query 1.fa -db ta_db -out 1.output -evalue 1e-5 -outfmt 6 -num_alignments 1

提取最终ID信息

cat 1.output | awk '{print $2}' | uniq > finalid.txt

id.list

到这里提取受选择基因的前半部分就大功告成了，有很多细节还是要根据自己的情况来进行处理的，不可一概而论。正所谓师傅领进门，修行在个人，用到的perl脚本还得靠各位读者自行努力，当然你要是plmm我就给你了哦,嘿嘿嘿。