又是一个目的，需要画图+理解
一，EHH的基本原理
参考文章
群体进化选择消除分析 - 组学大讲堂问答社区 (omicsclass.com)
基于LD衰减和等位基因频率检测自然选择-XPEHH (qq.com)
1.1基本概念介绍和理解
LD通常是指两个等位基因之间的关联程度，在这里，我们用LD来衡量一个位点的单个等位基因与不同距离处的多个位点之间的关联程度。由于重组，突变附近的LD将大幅衰减，大多数等位基因通常是比较古老的突变并且具有较短范围的LD。而一个新突变正在经历正选择（这里主要考虑Hard sweep）的关键特征就是等位基因频率在群体中异常快速上升，发生的时间足够短以至于重组基本上来不及破坏所选突变的单倍型，简而言之，正选择印记就是一个等位基因在给定群体频率下具有异常长范围的LD。LD范围的长短是相对的，主要取决于局部的重组速率，因此，基于LD衰减的检测方法必须控制重组率的局部差异，EHH是使用LD的衰减作为估计等位基因年龄的时钟来检测兴趣位点处的正选择。
Sabeti在2007年发表了跨群体延展单倍型纯合度（Cross-population Extended Haplotype Homozygosity, XPEHH）方法，适用于所选择的等位基因接近或已经在一个群体中固定，但是在种内所有群体中还保持多态性。
XPEHH 值是一种度量两个群体之间不平衡遗传连锁区域的指标。XPEHH 值通常在不同基因组位置上计算，正值表示一个群在某一基因座位上的连锁单倍型频率高于另一个人群，负值则表示相反情况。较高的 XPEHH 值可能暗示着此基因座位上存在自然选择的迹象。
1.2群体进化选择消除分析
还看到了一个有意思的现象群体选择分析|selscan 计算 XPEHH (qq.com)
1.3什么是Phasing
人类基因组的Phasing原理是什么？ - 知乎 (zhihu.com)
人类基因组的Phasing原理是什么？ (qq.com)
1.4查考文献
二，操作流程
在学习任何一个包之前，我们得先查清楚这个官网就跟snpeff一样
我觉得官网也说得很清楚，原理官网也有。selscan/manual/selscan-manual.pdf at master · szpiech/selscan (github.com)
去官网下载
szpiech/selscan: Haplotype based scans for selection (github.com)
tar -zxvf selscan-2.0.0.tar.gz
找到selscan所在的地方pwd
/mnt/e/绘图/xp_ehh/selscan-2.0.0/bin/linux
vim ~/.bashrc
export PATH=/mnt/e/绘图/xp_ehh/selscan-2.0.0/bin/linux:$PATH
source ~/.bashrc
根据文章说的，我们开始进行下一步
2021-01-04根据SNP验证自然选择—XPEHH分析 - 简书 (jianshu.com)
1.Beagle软件进行phasing
注意：前面都是我处理的过程，可以直接跳到后面去看代码。
selscan软件要求输入的基因型只能是0或1，而不能是其它。所以我们需要在这里把多等位基因位点过滤掉
从过滤后的vcf文件提取chr/id/physical position这三列，
需要去掉以#号开头的注释行sed -i '/^#/d' dis.map
需要对原始的vcf文件进行phasing
需要准备一个map文件，这个文件主要用来存放每个标记的ID以及所在的染色体、物理位置、遗传距离（cM）先对vcf进行原始处理，需要以上的东西，我大概整理一下，应该就是需要这些内容。
#VCF转ped/map
plink --vcf 文件.vcf --recode --out 文件--chr-set 24
#ped/map转bim/bed/fam
plink --file wenjain --make-bed --out wenajian --chr-set 24
#bim/bed/fam转VCF
plink --bfile wenjain --export vcf --out wenajin --chr-set 24
#bim/bed/fam转ped/map
plink -bfile xian1 -recode -out xian1 --chr-set 24
#ped/map格式的文件(mydata.ped/mydata.map)，就用**–file**
#是bed/fam格式的文件(mydata.bed/mydata.fam)，就用**–bfile**
我先把vcf格式的ped/map格式还有bim/bed/fam格式都转出来。
这里我们可以先查看一下map文件格式是什么
VCF 生成 plink 输入文件（ped和map） - 简书 (jianshu.com)
我已经看明白这个过程是干什么的了，根据文章说明selscan软件要求输入的基因型只能是0或1，需要在这里把多等位基因位点过滤掉，然后把vcf文件phasing
conda install -c bioconda vcftools
vcftools --vcf co --min-alleles 2 --max-alleles 2 --recode --out Vitis_all_sample
根据教程链接继续下载，这个教程的是BEAGLE2021-01-04根据SNP验证自然选择—XPEHH分析 - 简书 (jianshu.com)
Haplotype phasing 常用软件实操介绍 (qq.com)
wget -c https://faculty.washington.edu/browning/beagle/beagle.01Mar24.d36.jar先下载软件，把vcf分成24个，map也是（我寻思map应该也是需要分群的好了的map）
bcftools view Vitis_all_sample.recode.vcf -Oz -o Vitis_all_sample.recode.vcf.gz
bcftools view -S pop1.txt Vitis_all_sample.recode.vcf.gz > pop1.vcf
（这个地方的分群文件我是用的是1群体和2群体），pop1.txt文件格式可以参考使用bcftools提取指定样本的vcf文件（extract specified samples in vcf format） - 橙子牛奶糖 - 博客园 (cnblogs.com)
for i in {1..24}; do vcftools --vcf pop1.vcf --chr $i --recode --out $i.pop1; done
for i in {1..24}; do gzip $i.pop1.recode.vcf; done
plink --vcf Vitis_all_sample.recode.vcf --recode --out Vitis_all_sample_recode_pop1 --chr-set 24
mkdir ./map
mv Vitis_all_sample_recode_pop1.map ./map/
cd map
awk '{print > $1}' new.map
for i in {1..24}; do gzip $i; done
使用Beagle软件进行phasing，
for i in {1..24}; do gzip -d ${i}.gz; done（也可以不需要压缩文件）
for i in {1..24}; do java -Xmx30g -jar beagle.01Mar24.d36.jar gt=//mnt/e/绘图/xp_ehh/pop1/${i}.pop1.recode.vcf.gz map=/mnt/e/绘图/xp_ehh/map/${i} out=${i}.output; done
运行这个代码的时候我的map报错，说不能都是零
1 df-1-1791-1336867 0 1791
1 df-1-1857-1709605 0 1857
1 df-1-2569-710694 0 2569那我还是得安装这个步骤要求生成map文件
这里给两个计算的例子
接下来算遗传距离的均值，并把物理距离转换成遗传距离：vitis基因组大小 ~500M，查到16年一篇文献vitis总的遗传距离为2424cM，相当于每单位物理位置为0.00000485cM
[ \text{每单位物理位置（cM/bp）} = \frac{\text{总的遗传距离（cM）}}{\text{基因组大小（bp）}} ]
将给定的数值代入这个公式中：
[ \text{每单位物理位置（cM/bp）} = \frac{2424}{500,000,000} ]
[ \text{每单位物理位置（cM/bp）} = 0.000004848 ]
因此，每单位物理位置（每个碱基对）所对应的遗传距离约为0.00000485 cM/bp。
此例中基因组大小1G ，遗传距离为3950~3200cM，相当于每单位物理位置为0.0000032cM（按照3200cM算吧） 生成最终的map文件
首先，将基因组大小转换为以字节为单位：
[ 1G = 1 \times 10^9 \text{ 字节} ]
现在，我们将总的遗传距离范围和基因组大小代入公式中：
[ \text{每单位物理位置（cM/bp）} = \frac{3950 \text{~} 3200}{1 \times 10^9} ]
[ \text{每单位物理位置（cM/bp）} \approx 3.2 \times 10^{-6} \text{~} 3.95 \times 10^{-6} ]
所以我算出我的是0.00000166 cM/bp
由于我们的map文件是这样的格式，
第一列是染色体；第二列是位点ID；第三列是遗传位置；第四列是物理位置。cut -f 1,3,2 Vitis_all_sample.recode.vcf > dis.map
awk '{ print $1" "$3" "$2*0.00000166" "$2}' dis.map > dis.4.map
1 . 0.00431434 2599
1 . 0.00772066 4651
1 . 0.00772564 4654然后重新对我们map进行分割
mkdir ./map
mv dis.4.map ./map/
cd map
awk '{print > $1}' dis.4.map
使用Beagle软件进行phasing，
for i in {1..24}; do java -Xmx30g -jar beagle.01Mar24.d36.jar gt=//mnt/e/绘图/xp_ehh/pop1/${i}.pop1.recode.vcf.gz map=/mnt/e/绘图/xp_ehh/map/${i} out=${i}.output; done
这个是跑其中一个群体，接下来我们同样的跑第二个群，这里把前面探索的代码献上
需要的文件
总的vcf，分群文件1052_1052
1，selscan软件要求输入的基因型只能是0或1，需要在这里把多等位基因位点过滤掉
conda install -c bioconda vcftools
vcftools --vcf c.vcf --min-alleles 2 --max-alleles 2 --recode --out Vitis_all_sample
2，输出map格式的文件，分成24个
plink --vcf Vitis_all_sample --recode --out Vitis_all_sample --chr-set 24
cut -f 1,2,4 Vitis_all_sample.map > Vitis_all_sample.cut.map
awk '{ print $1" "$2" "$3*0.00000166" "$3}' Vitis_all_sample.cut.map > Vitis_all_sample.cut.awk.map
mkdir ./map
mv Vitis_all_sample.cut.awk.map ./map/
cd map
awk '{print > $1}' Vitis_all_sample.cut.awk.map
cd ..
3，先压缩vcf文件，分成24个
bcftools view Vitis_all_sample.recode.vcf -Oz -o Vitis_all_sample.recode.vcf.gz
bcftools view -S 3.txt Vitis_all_sample.recode.vcf.gz > pop3.vcf
bcftools view -S 1.txt Vitis_all_sample.recode.vcf.gz > pop1.vcf
for i in {1..24}; do vcftools --vcf pop1.vcf --chr $i --recode --out $i.pop1; done
for i in {1..24}; do vcftools --vcf pop3.vcf --chr $i --recode --out $i.pop3; done
for i in {1..24}; do gzip $i.pop1.recode.vcf; done
【for i in {1..24}; do gzip -d ${i}.gz; done（也可以不需要压缩文件）】
使用Beagle软件进行phasing，
for i in {1..24}; do java -Xmx30g -jar beagle.01Mar24.d36.jar gt=/mnt/e/绘图/xp_ehh/pop1/${i}.pop1.recode.vcf.gz map=/mnt/e/绘图/xp_ehh/map/${i} out=${i}.output1; done nthreads=62，使用selscan做XPEHH分析
终于到这一步了，普天同庆，前面真的好麻烦。根据文章的提示
使用 Vcftools 、Selscan、Rehh包进行全基因组选择信号（选择压力）分析 (qq.com)
基于LD衰减和等位基因频率检测自然选择-XPEHH (qq.com)
2021-01-04根据SNP验证自然选择—XPEHH分析 - 简书 (jianshu.com)
等待两个亚群phasing结束后，我们就可以进行XPEHH分析了，考虑到我之前对文件的存放进行了处理需要切换成自己的路径就行
taskset -c 0-3 nice -n 10 bash -c 'for i in {1..24}; do selscan --xpehh --vcf "/mnt/e/绘图/xp_ehh/output1/${i}.output1.vcf.gz" --vcf-ref "/mnt/e/绘图/xp_ehh/output3/${i}.output1.vcf.gz" --map "/mnt/e/绘图/xp_ehh/map/${i}" --threads 6 --out $i ; done'
cut -f 1 *.xpehh.out > 1_3.xpehh
可以得到和文章一样的文件了，接下来进行可视化。
3，XPEHH分析结果可视化
install.packages("tidyverse")
install.packages("CMplot")
library(tidyverse)
library(CMplot)
##########################
input <- "E:/绘图/xp_ehh/1_3xpehh/chuli.xpehh"
prefix <- "black"
#读入文件
xpehh <-
read.table(input, header = T) %>%
select(id, chr, pos, xpehh)
getwd()
setwd("E:/绘图/xp_ehh/1_3xpehh")
#作图
png(paste(prefix, "_1XPEHH.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()
三，报错分析
1，如果没有跑出output，就单独一个一个的跑
2，直接删掉，记得map里面也需要删掉
3，出现killed，直接上服务器，或者是
taskset -c 0-3 nice -n 10 java -Xmx10g -jar beagle.01Mar24.d36.jar gt=/mnt/e/绘图/xp_ehh/pop4/18.pop4.recode.vcf.gz map=/mnt/e/绘图/xp_ehh/map/18 out=18.output1
4，出现这个不匹配，因为selscan需要的是snp位点一一配对，需要重新跑一次看看是不是output文件有问题
2563
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
