选择信号检测--基于单体型haplotype的选择信号的检测（iHS）

这篇播客主要借鉴了大佬的信号检测文章原文链接：

1.数据文件制备

一、首先手里拿到一份plink的ped、map文件：

(1)转vcf文件格式
plink --file test --recode vcf --out test
(2)利用vcf文件过滤掉一些质量较低的位点
 (最大缺失率（max-missing）<10%即检出率>90%的SNP位点，
 最小等位基因频率（MAF）>0.03；
(去除多等位基因和等位基因缺失: --min-alleles 2 --max-alleles 2）
 vcftools --vcf test --remove-indels --max-missing 0.9 --maf 0.03 --min-alleles 2 --max-alleles 2 --recode --recode-INFO-all --out test.zk.vcf
(3)质控结束后直接进行定相填充
java -Xmx1g -jar beagle... gt=test.zk.vcf out=test.gt.zk.vcf nthread=10
一般结束以后就有货了！！！
(3.1)顺手转个ped、map
plink --vcf test.gt.zk.vcf --recode --out test.gt

二、统一格式： map文件位置名称 ，空格分隔

awk '{print$1,$1":"$4,$3,$4}' test.gt.map > test.gt1.map
转换成vcf格式
plink --file test.gt1 --recode vcf --out test.gt1

三、使用perl代码处理vcf特异格式

perl -pe "s/\s\.:/\t.\/.:/g" test.miss.gt1.vcf | bgzip -c > test.gt1.vcf.gz
意思是：用perl执行"s/\s\.:/\t.\/.:/g"

四、再次定相填充

java -Xmx1g -jar beagle... gt=test.gt1.vcf.gz out=test.2gt 

五、分染色体：这个是重点

#vcf文件分割
for i in {1..27}
do
  vcftools --vcf 12.perl.gt.vcf.gz --chr ${i} --recode --recode-INFO-all --out ${i}.test
done
#map文件分割
for k in {1..27}
do
  vcftools --vcf ${k}.test.recode.vcf --plink --out ${k}.mp
done

六、遗传距离（Bp）的转化：这个是最重的重点！因为涉及到该物种的基因组大小。这个必须整明白。

for k in {1..27}
do
  awk 'BEGIN{OFS=" "} {print 1,".",$4/1000000,$4}' ${k}.mp.map > ${k}.MT.map2
done

看明白代码：只取map文件后两列！

七、计算iHS

for k in {1..27}
do
  selscan --ihs --vcf ${k}.test.recode.vcf --map ${k}.MT.map2 --out ${k}.iHS
done
# 得到 ${k}.iHS.ihs.out文件，运行时间久！




# 提取结果
for k in {1..27}
do
  awk '{print '${k}',$2,$3,$4,$5,$6}' ${k}.iHS.ihs.out > ${k}.hu.ihs.out
  sed -i 's/ /\t/g' ${k}.hu.ihs.out
done

七、搞大事
到了这一步就可以基本搞定了。一行代码，合并文件

cat *hu.ihs.out > all.ihs

接下来在all.ihs文件中取一列ihs绘制曼哈顿图就搞定了可别像我一样画一张难看的图能被老师骂死，话说人长的丑就不能图丑点嘛~~~