自从第一次了解选择信号之后就对其产生了浓厚的兴趣，通过读文章知道有许多的方法来对基因组上的选择信号进行检测，比如Fst，EHH，IHS，LD等等方法都可以检测选择信号，但是我通过阅读文献，大部分人都选择了xpehh方法来检测选择信号，我想应该是这个方法检测的比较准确和权威吧！下面就简单介绍一下xpehh的使用。
xpehh有许多方法来得以实现，比如基于linux就可以很好地实现，但是对于linux来讲，除了专业人士以外，大部分人都对其比较陌生。如果是初次接触生信，就更是雪上加霜。所以相对来讲，基于windows的r语言，就容易很多，也容易上手。而且r语言相对其他计算机语言来讲简单的不是一点半点，而且r语言还有Rstudio这个比较方便的窗口式编辑器。好了，废话不多说，接下来就简单的分享一下我自己 的学习心得。
#首先在打开Rstudio之后，确定固定路径,
> setwd("E:\\R语言\\test") #比如我这个路径就是在E盘里面R语言这个文件夹下的test这个文件夹。
#然后调用“LD”这个r语言的package，注意这个包如果没有的话需要提前安装。
> library(LD) #载入LD的包
#然后在上面的路径下创建几个文件夹
>ldstart() #构建目录，将会在上述目录中创建'haplotypes', 'populations', 'snps', 'rehh_in', and 'rehh_out'.五个文件夹，其中haplotypes文件夹中手动放入单倍型文件(.haps和.sample)这两个文件是由shapeit软件获得。生成的单倍型文件名要以chr开头，后面跟染色体ID，如：chr30.haps，chr30.sample。shapeit的下载和使用方法可以从https://mathgen.stats.ox.ac.uk/genetics_software/shapeit/shapeit.html网址获得；population文件夹中需要手动放入两个种群的个体ID，每个ID一行，文件名必须与填充中的名称相对应。'snp s'文件夹中放入snp 的ID（例如rs1042522）的文本文件（每行一个id）。
#上述准备文件都OK以后，往下进行
> make.binary.haplotypes(chr=30) #此步会生成一个hap_sample.Rdata文件，在haplotyprs的文件夹下。需要注意的是，chr=30这是我给出的例子，结合自己的数据想要计算那些染色体就写那些比如第2号染色体就可以写chr=2，或者想要计算1到10号染色体就写chr=（1：10）。
#然后下一步生成xpehh的准备文件。
> make.input.files(populations = c('high.txt', 'low.txt'), chrs = 30, haps.subset = T, output = c('high','low')) #该过程会在rehh_in文件夹中生成9个文件。注意我这里的high和low两个文件是我的例子文件，需结合自己的文件进行编译。
#然后继续往下
> make.scanhh(chrs = 30) #这一步会在路径下生成一个scanhh.RData文件，这个文件很重要哦。
#然后在运行xpehh之前需要先加载scanhh.RData文件
> load("./scanhh.RData") #载入scanhh.RData文件
#好了前面的所有准备工作做完了，马上进入正题。
> xpehh(pop1 = scanhh.list$high, pop2 = scanhh.list$low, snp.list = 'snp.txt', filter = 2,method = 'unilateral', annot = T, write.xls = 'ss.snps', plot = T) #会在rehh_out文件夹中生成xpehh结果文件。
pop1---具有要与第二个总体进行比较的第一个总体的名称的标量字符。
pop2---具有要与第一个总体进行比较的第二个总体的名称的标量字符。
snp.list---就是放进snps文件夹中的snp 的ID 的文件。
filter ---用于筛选重要xp-ehh分数的数值（默认值=2）
method ----就是想计算单端还是双端，我个人认为设置为“both”比较好，这样双端单端都会出来结果。
后面的值默认就可以。
这样做完之后，就可以去我们的路径下的rehh_out文件夹中寻找我们的结果啦。是不是很简单呢。
完整的代码
> setwd("E:\\R语言\\test") #给定固定的路径
> library(LD) #载入LD的包
> ldstart() #构建目录，将会在上述目录中创建'haplotypes', 'populations', 'snps', 'rehh_in', and 'rehh_out'.五个文件夹，
> make.binary.haplotypes(chr=30) #此步会生成一个hap_sample.Rdata文件，在haplotyprs的文件夹下。
> make.input.files(populations = c('high.txt', 'low.txt'), chrs = 30, haps.subset = T, output = c('high','low')) #该过程会在rehh_in文件夹中生成9个文件。
> make.scanhh(chrs = 30)
> load("./scanhh.RData") #载入scanhh.RData文件
> xpehh(pop1 = scanhh.list$high, pop2 = scanhh.list$low, snp.list = 'snp.txt', filter = 2,method = 'both', annot = T, write.xls = 'ss.snps', plot = T)#会在rehh_out文件夹中生成xpehh结果文件。
之后更多的细节可以参考https://rdrr.io/github/cmcouto-silva/LD/src/R/xpehh.R。上述的流程与代码都是我自己学的，有可能不是完全正确，但是这个流程确实可以做出结果。在这里分享给大家，如果有想做xpehh的，可以给你带来一点帮助，还有如果我有哪步写的不是很正确，请大神们指教，我再更正。
网友评论