转载自https://mp.weixin.qq.com/s/1IM6KSCGIq0cUJhKKXVcIw

背景

选择清除分析是常用于鉴定群体基因组水平受选择区域的一种方法，其背后的理论是：在受选择群体中，会出现有利等位基因的频率快速上升，与其相邻位点的多态性快速下降的现象。基于此，我们可以通过计算群体内基因组区域的等位基因频率的变化情况来鉴定受选择基因组区域。常用群体内检测指标的计算方法大致分为三种：1.基于核苷酸多态性降低的π、θw；2.基于分离位点频率的Tajima’D；3.基于连锁不平衡增加的EHH、iHS。以上三类指标对应于基因组受选择特征的三个维度，而后才有了群体间的选择指标：1.由π衍生的π ratio、ROD、Fst；2.由EHH衍生的XPEHH。本文主要介绍XPEHH的实操。

XPEHH分析

XPEHH是基于群体单倍型的选择清除分析，有很多工具可以实现XPEHH方法，比如R包rehh、Linux平台下selscan软件。在这里我用selscan软件完成。

1.软件下载与安装

selscan下载网址：https://github.com/szpiech/selscan/releases选择最新版本，复制链接地址

Linux下输入

2.查看分析所用数据文件及格式

selscan软件帮助文档可以在此下载：https://github.com/szpiech/selscan/blob/master/manual/selscan-manual.pdf

通过软件文档我们可以了解到XPEHH分析有三种输入方式：

wget -c https://github.com/szpiech/selscan/releases/download/v1.3.0/selscan-linux-1.3.0.tar.gz

tar -zxvf selscan-linux-1.3.0.tar.gz #解压压缩包

cd ./selscan-linux-1.3.0

pwd

selscan=软件所在绝对路径

2.查看分析所用数据文件及格式

selscan软件帮助文档可以在此下载：https://github.com/szpiech/selscan/blob/master/manual/selscan-manual.pdf

通过软件文档我们可以了解到XPEHH分析有三种输入方式：

selscan --xpehh --hap <pop1 haplotypes> --ref <pop2 haplotypes> --map <mapfile> --out <outfile> 

#hap文件格式输入

selscan –xpehh –vcf <pop1 vcf> --vcf-ref <pop2 vcf> --map <mapfile> --out <outfile> #vcf文件格式输入

selscan –xpehh –tped <pop1 tped> --tped-ref <pop2 tped> --map <mapfile> --out <outfile> #tped文件格式输入

因为我的数据文件是所有个体总的vcf文件，所以选用vcf文件作为输入文件格式最为方便。软件对vcf文件格式要求：前九列同标准的vcf格式一样，后面紧接着的个体的基因型尽量是phasing后的。文件格式如下：< quality ><individual 1 genotype > ... < individual N genotype >

1 100 rs1 0 1 . . . . 0|1 0|0

1 200 rs2 0 1 . . . . 1|1 1|0

1 300 rs3 0 1 . . . . 0|0 1|0

1 400 rs4 0 1 . . . . 0|0 1|0

1 500 rs5 0 1 . . . . 1|0 0|1

所以我们需要对原始的vcf文件进行phasing

另外我们还需要准备一个map文件，这个文件主要用来存放每个标记的ID以及所在的染色体、物理位置、遗传距离（cM）。文件格式如下：

1 rs1 0.01 100

1 rs2 0.02 200

1 rs3 0.03 300

1 rs4 0.04 400

1 rs5 0.05 500

3.数据文件准备

3.1 map文件的准备

除了做模式植物的，大多数情况下我们都没有一个现成的能和自己当下vcf文件一一对应的map文件。我们可以根据物种总的遗传距离和基因组大小，算平均的遗传距离，而后根据每个标记的物理位置推断遗传距离。

这是原始的vcf文件。

由于参考基因组版本问题，这里存在很多位于随机或未知染色体上的变异位点，所以需要先过滤掉这些位点。

grep -v "random" file.vcf > file1.vcf 

grep -v "Un" file1.vcf > filtered.vcf

另外需要注意的是，selscan软件要求输入的基因型只能是0或1，而不能是其它。所以我们需要在这里把多等位基因位点过滤掉，避免后续分析出现问题（我当时就吃了这个亏，导致分析重来

Vcftools --vcf filtered.vcf --min-alleles 2 --max-alleles 2 --recode --out Vitis_all_sample

接下来我们就从过滤后的vcf文件提取chr/id/physical position这三列

cut -f 1,3,2 Vitis_all_sample.recode.vcf > dis.map

直接cut vcf文件会带#号的注释行一同输出到dis.map文件中，所以需要去掉以#号开头的注释行

sed -i '/^#/d' dis.map

接下来算遗传距离的均值，并把物理距离转换成遗传距离：vitis基因组大小 ~500M，查到16年一篇文献vitis总的遗传距离为2424cM，相当于每单位物理位置为0.00000485cM 生成最终的map文件

awk '{ print $1" "$3" "$2*0.00000485" "$2}' file.map > new.map

注意：map文件要求是空格分割

3.2 vcf文件的phasing

在phasing之前我需要根据研究目的从总的样本的vcf文件中提取两个亚群的子vcf文件。准备两个群体分类表，每行是一个样本名称（注意这个样本名称要和vcf文件记录的样本号保持一致） 我使用bcftools提取vcf文件的子集

#首先对vcf文件进行压缩

bcftools view Vitis_all_sample.recode.vcf -Oz -o Vitis_all_sample.recode.vcf.gz

#然后建立索引

bcftools index Vitis_all_sample.recode.vcf.gz

#最后根据群体分类表提取子vcf文件

bcftools view -S pop.list Vitis_all_sample.recode.vcf.gz > new.vcf

这里压缩vcf文件以及建立索引这一步很重要，这样才可以保证后面提取vcf文件不会出错。如果你不压缩以及建立索引，直接提取子集很容易出现no header..一类的报错，这是因为vcf文件中的header部分并没有指明全部的chr，可以通过修改vcf文件header信息来解决这个问题。但是我感觉这样做有些麻烦，所以还是推荐先压缩一下，然后建立索引，这样能避免很多错误。

phasing软件有很多，常用的主要有SHAPEIT、Beagle。关于用SHAPEIT进行phasing公众号先前文章Haplotype phasing 常用软件实操介绍和黄树嘉如何使用Shapeit2对人类基因组数据进行Phasing在这里我用Beagle软件来做phasing。

Beagle软件下载 官网：http://faculty.washington.edu/browning/beagle/beagle.html

右键复制.jar文件链接地址，然后转到Linux界面下输入

wget -c http://faculty.washington.edu/browning/beagle/beagle.18May20.d20.jar

pwd

beagle=path/beagle.18May20.d20.jar

因为这是用java语言写的软件，所以无须安装，有java环境即可运行。

常规的vcf文件就可以作为Beagle软件的输入，此外还需要输入一个map文件，其格式与selscan要求的map文件格式是一致的，前面已经准备过了，接下来就可以phasing了！等等，其实还有一件事没有做.....单倍型以及基于单倍型的分析都需要逐条染色体运行，所以我们还需要把vcf文件以及map按照染色体ID进行拆分。

拆分vcf文件

for i in {1..19}

do

vcftools --vcf file.vcf --chr $i --recode --out $i

done

注：葡萄有19条染色体，所以循环次数为1-19

拆分map文件

mkdir ./map#创建一个map文件夹

mv new.map ./map#把整理好的map文件移动到map文件夹

awk'{print > $1}'new.map#这行代码可以直接把map文件根据chr列拆分成多个文件

好了，现在可以开始phasing了！

Beagle软件的运行方式：java -Xmx[GB]g -jar beagle.jar [arguments]，需要注意的是输入的参数和值由=连接，中间不能有任何其它字符。

for i in{1..19}

do

Java -Xmx10g -jar$beaglegt=${i}.vcf map=../map/${i}ne=10 out=$inthreads=6

done

gt指输入的vcf文件，ne是指样品数量，nthreads就是线程数，out指定输出文件前缀。最后会生成.vcf.gz格式的文件。

4.selscan做XPEHH分析

等待两个亚群phasing结束后，我们就可以进行XPEHH分析了，运行代码为：

for i in{1..19}

do

$selscan--vcf"${i}.vcf.gz"--vcf-ref"../w/${i}.vcf.gz"--map"../map/${i}"--threads 6 --out$i

done

--vcf指想研究亚群的所有染色体的vcf文件，--vcf-ref是用作参考的亚群，map就是前面已经准备好的单条染色体的map文件。根据自己数据存放的位置，要灵活调整指定的文件路径。

最后生成结果

XPEHH分析结果可视化

展示选择清除分析结果的最常用的就是曼哈顿图。如下图，横轴代表染色体，纵轴是不同选择清除分析指标算的值，图中的每一个点代表一个sweep或marker；如果是GWAS每一个点就代表一个SNP。

所以要画这种图，我们先要把XPEHH分析所得染色体的文件合并为一个文件

cat *. > black.xpehh#在存放xpehh结果文件中直接运行这段命令

注意：因为每条染色体生成的xpehh文件都会有id\pos等表头，直接cat会造成生成的总文件中存在多个表头行，会影响整成绘图，所以我在R代码中增加了一个对文件的过滤步骤。

绘制曼哈顿图可以利用R包CMplot、qqman。我用CMplot展示下XPEHH的分析结果

绘图代码：

library(tidyverse)

library(CMplot)

input <- “path+file”

prefix <- “black”

#读入文件
xpehh <-

  read.table(input, header = T) %>%

  filter(!str_detect(id, 'id')) %>%

  separate(col = id, into = c("chr", 'position'), sep = "__", remove = F) %>%

  select(id, chr, pos, xpehh)

#作图
png(paste(prefix, "_XPEHH.png", sep = ""), width=960, height=480)

CMplot(xpehh,

      type = "p",

      plot.type = "m",

      LOG10 = F,

      cex = 0.2,

      band = 0.5,

      ylab.pos = 2,

      cex.axis = 1,

      ylab="XPEHH",

      file.output = F )

dev.off()

在R代码中需要注意以下几点：

1、filter(!str_detect(id, 'id')) 用来过滤掉文件中多余header行；

2、separate(col = id, into = c("chr", 'position'), sep = '__', remove = F) 用来分割id行，目的是为了生成chr列；

3、select(id, chr, pos, xpehh) 因为CMplot对输入的作图文件的要求：前三列必须为id\chr\pos。所以这里用select函数重新整理了以下作图文件。

参考资料：

基于LD衰减和等位基因频率检测自然选择-XPEHHhttps://mp.weixin.qq.com/s/-YydvSjsxTbRxjIPtUGLsA

Haplotype phasing 常用软件实操介绍https://mp.weixin.qq.com/s/YwXq_MA4SEXJhNzm-lF5-g

Beagle 5.1 http://faculty.washington.edu/browning/beagle/beagle.html

selscan v1.3.0 https://github.com/szpiech/selscan/releases/tag/v1.3.0

vcftools https://vcftools.github.io/man_latest.html

bcftools http://samtools.github.io/bcftools/