群体选择信号分析
群体进化与选择信号：
●生活在世界不同区域的生物群体在历史长河中经历千万年的自然选择、人工驯化、迁徙、遗传漂变等事件的洗礼，有的物种灭绝，有的物种分化出不同的亚种或亚群,甚至逐渐演变成新的物种。群体进化研究就是用来追溯和揭露这个进化过程的。
●现代群体进化研究，开始基于全基因组的变异信息，分析群体的遗传背景、群体结构、遗传多样性、基因交流情况、物种的形成机制以及群体的进化方向等生物学问题。
●群体进化主要的应用方向有:
1.人工驯化机制研究:通过对野生型和驯化型群体的研究，发现人工驯化作用下受选择的区域，进而鉴定到与重要经济性状相关的基因。
2.自然选择机制分析:研究生存于不同地理环境下的群体，挖掘适应性进化过程中受到选择的区域，进而为育种提供环境适应性基因资源。
3.种群历史研究:通过分析物种的起源地及各个分布区域中群体的遗传变异信息，探索物种的进化过程。
●选择信号分析是群体进化这个大的研究领域中的一个研究方向/方法，多用于前两个应用方向的研究。
●由于人工选择或自然选择作用造成的基因组结构特征的变化称为选择信号。选择信号分析是生物信息分析中的“考古学”。它可以帮助我们寻找生物群体在千万年的进化过程中，在基因组上留下的痕迹。
●这些信号与动物的选育方向以及驯化适应机制紧密相关，因此对选择信号进行检测有助于挖掘与动物经济性状相关的基因，了解性状形成的潜在的遗传学基础，解析动物的驯化适应机制，解析群体间的差异机制。对于畜禽遗传改良具有重要意义。
选择扫荡（遗传搭车效应）.png
选择信号分析的特性和应用场景
●对于物种的群体遗传研究，本身就涉及到多个领域，如考古学、生态学、社会学，然后才是基因组学和生物信息学分析。所以涉及群体遗传的研究项目,需要尽量将这几类因素综合在一起讨论，才能形成完整的故事。
1.有明确的研究方向和研究目的(如:研究猪的高原适应性)
2.样本材料具有相关特性且与研究目的相契合(如:藏猪具有高原适应性)
3.对材料背景非常了解，包括起源地/栖息地、(地理位置..环境和生态)、 迁徙历史、驯化历史、育成历史等(四面环山的四川盆地成为冰川时期动物群体的避难所)
4.当定位到受选择的区域之后，把受选择区域的位点做基因注释，需要了解与所研究性状相关的基因、蛋白、通路等(受选择区域中发现'vitamin B6 binding相关基因:维生素B6可帮助血红蛋白合成和促进氧的绑定)
5.建议每个亚群个体数量不小于10个，且来自不同家系
选择信号分析的一般流程.png
遗传背景分析
一、主成分分析(Principal Component Analysis)
●PCA是一种线性代数中的数据处理方法，它利用降维的思想，从高维度数据(如测序得到的百万级别SNP位点数据) 中提取关键的信息，以便我们使用更少的变量(指标)就可以对样本进行有效区分。这些被提取出的信息按照其效应从大到小排列，我们称之为主成分1(Principal Component1)、主成分2、主成分3...
●PCA分析的应用场景:
1.检测离群样本
2.推断群体分层和亚群间的遗传距离
二、进化树
●又称为系统发生树，它利用样本间的差异度将样本进行聚类，用一种类似树状分支的图形来概括各物种之间的亲缘关系，可用来描述物种之间的进化关系和遗传距离远近。
●不同的构树方法
1.基于距离的方法:首先通过各个物种之间的比较，根据一 定的假设(进化距离模型)推导得出分类群之间的进化距离，构建一个进化距离矩阵。 进化树的构建则是基于这个矩阵中的进化距离关系。(UPGMA, NJ)
2.基于特征的方法:不计算序列间的距离，而是将序列中有差异的位点作为单独的特征，并根据这些特征来建树(ML MP)
●进化树的解读
1.枝长:枝长累积距离越近的样本差异越小
2.自展值:进化树分支可信度（蓝圈，百分比75%以上比较可信）
3.标尺:代表序列的差异程度
进化树的解读.png
三、群体结构(structure) 分析
●先预设群体由若干亚群(k=n)构成,通过模拟算法找出在k=n的情况下，最合理的样本分类方法。最后再根据每次模拟的最大似然值，找出最适用这群体的K值。
●应用场景:
1.推算亚群划分情况
2.推算群体基因交流程度
3.推算个体的血统构成比例
●主流软件
1.STRUCTURE
2.fastSTRUCTURE
3.Admixture
四、LD衰减分析
●连锁不平衡
●当位于某一座位的特定等位基因与另一座位的某一等位基因同时出现的概率大于群体中随机分布的两个等位基因同时出现的概率时，就称这两个座位处于连锁不平衡状态。
●一般而言，两个位点在基因组上离得越近，相关性就越强，LD系数就越大。反之，LD系数越小。也就是说，随着位点间的距离不断增加，LD系数通常情况下会慢慢下降。这个规律，通常就会使用LD衰减图来呈现。
连锁不平衡.png
●LD衰减图就是利用曲线图来呈现基因组上分子标记间的平均LD系数随着标记间距离增加而降低的过程。
●大概的计算原理就是先统计基因组上两两标记间的LD系数大小，再按照标记间的距离对D系数进行分类,最终可以计算出一定距离的分子标记间的平均LD系数大小。
LD衰减图.png
图中横轴为标记间的物理距离，纵轴为平均r方值。
群体多样性越低，衰减速度越慢。
●LD衰减分析的应用
1.评估群体特性和选择强度:驯化选择会导致群体遗传多样性下降，位点间的连锁程度更高。所以，通常驯化程度越高选择强度越大的群体,LD衰减速度也越慢。例如商品化群体比自然群体通常更大的LD衰减距离。类似的自然选择、遗传漂变导致的群体遗传多样性下降，也会减慢LD衰减的速度。
2.检测受选择基因组区域:与有利突变紧密连锁的中性位点会由于选择作用在基因组上形成高频率的核心单倍型，以其为中心向基因组两侧扩展会形成长范围的扩展单倍型。然而随着与有利突变间距的增加，连锁不平衡程度会相应衰减，在一定范围内各扩展单倍型纯合的总和占核心单倍型纯合的比例可以被用来检测基因组范围内的选择作用。
3.GWAS分析中评估标记密度是否足够: GWAS分析本质就是利用标记和功能突变的相关性(LD关系)，来检测与性状相关的功能突变的位置。一般而言LD系数大于0.8就是强相关。如果LD系数小于0.1，则可以认为没有相关性。如果LD衰减到0.1这么大的区间内都没有标记覆盖的话即使这个区间有一一个效应很强的功能突变，也是检测不到关联信号的。所以通常可以通过比较LD衰减(到0.1)距离和标记间的平均距离来判断标记是否对全基因组有足够的覆盖度。(GWAS最低标记量≈基因组大小/LD衰减距离)
选择信号检测方法
常用群体内检测指标的计算方法大致分为三种：1.基于核苷酸多态性降低的π、θw；2.基于分离位点频率的Tajima’D；3.基于连锁不平衡增加的EHH、iHS。以上三类指标对应于基因组受选择特征的三个维度，而后才有了群体间的选择指标：1.由π衍生的π ratio、ROD、Fst；2.由EHH衍生的XPEHH。
https://zhuanlan.zhihu.com/p/52064863
对于单个物种，基于选择的效应，选择信号检测的方法可以被分为4大类:
1.基于等位基因频率谱的方法
2.基于连锁不平衡增加的方法
3.基于群体分化的方法
4.基于基因组杂合度的方法
一、基于等位基因频率谱
●基因型频率和基因频率的改变是选择作用在基因组上最直接的体现。基因频谱(site-frequency spectrum)就是指某种等位基因在基因组上某个目标区域内出现的频繁程度。
●符合中性模型的群体，其群体中存在广泛的遗传多态，当突变发生时总能够维持在一个较低的频率，只有当群体基因组上出现或存在有利突变时，选择才会发生作用，从而产生所谓的选择清除或搭车效应。
●代表性的检测方法: Tajima's D, Fu andLi'sD, Fay and Wu'sH, CLR, Hp
●Tajima's D检验的目的是区分随机演变的DNA序列(“中性”)和在非随机过程中演化的DNA序列，包括定向选择或平衡选择。
●Tajima's D的计算原理:多态位点数量和平均非匹配数量的差值。
●D=0时，符合中性假设，群体未受到选择; D<0时，受到定向向选择; D>0时，受到平衡选择。
二、基于连锁不平衡
●基于连锁不平衡理论，位点间的连锁不平衡程度会随标记间距离的增加而逐渐降低。因此，在基因组上可以观察到选择作用造成的不同长度的扩展单倍型纯合(Extended Haplotype Homozygousity)。
●该方法的基本原理是:在中性条件下，基因组很难形成长范围的连锁不平衡的单倍型，因为新突变需要经历漫长的遗传漂变才能达到较高频率，而在漫长的时间里会发生大量基因重组事件，使得这种连锁不断被打破。而当群体处于正向选择作用下时，致因突变及其连锁位点在正选择的作用下，在短时间内会达到较高频率，形成大片段的纯合单倍型。扩展单倍型纯合度检验正是基于这样的特征来筛选受选择基因。
●代表性的检测方法: EHH, XP-EHH, iHS, nSL, OmegaPlus
三、基于群体分化
●同一物种不同群体之间由于环境不同或选择目标不同，其基因组等位基因频率会表现出歧化选择的效应。这种现象在相同基因座位不同等位基因均受到选择时表现尤为明显，即选择加速群体分化。因此，基于群体分化的方法，不同群体同一等位基因频率存在的差异程度大于两个群体处于中性条件下的期望时，就推断该位点存在选择作用。
●代表性的检测方法: Weir and Cockerhan's Fst, LSBL, di
●Fst的取值范围为0-1，1表示群体间完全分化的位点，0表示在群体间完全没有分化的位点。
●基于Fst的的检测方法多采用基因组单位点扫描的策略，而这样的方式容易受到遗传漂变等因素的影响，产生假阳性的显著位点。为尽量减少假阳性的发生，通常采用滑动窗口的策略，降低这些干扰因素，增加选择信号检测的准确性。
四、基于基因组杂合度
●当基因组上特定区域受到选择时，由于“选择性清除”作用的存在，该区域及其连锁的区域表现为多态性降低，同时纯和度增加。因此对基因组的杂合度进行检测，可以推断出基因组中受到选择的区域。基因组上受选择程度越高，则杂合度程度越低。
●代表性的检测方法: θπRatio, ROH
●核苷酸多态性θπ比率越偏离1，受选择程度越高。θπ比率的检测公式如下:θπratio=θπA/θπB
其中，θπA和θπB分别代表A群体和B群体的θπ值。θπ比率大于1， 反映A群体的基因组杂合度大于B群体的杂合度，则B群体相应基因组区域受到选择。θπ 比率小于1,则A群体的基因组杂合度低于B群体，则选择发生在A群体对应的基因组区域。
方法汇总.png
选择信号的应用策略
●由于单纯一种选择信号检测方法容易造成假阳性选择信号的产生，目前选择信号检测普遍采用两种或多种检测方法进行组合策略，不同检测方法相互验证，用以避免假阳性，例如基于群体分化和基于基因组杂合度检测组合的策略、基于连锁不平衡的检测方法与基于基因组杂合度检测方法组合的策略等。也有文献使用3种以上的方法同时进行选择信号检测，筛选2种或以上方法检测重叠的位点。
●此外，在功能基因挖掘方面，选择信号检测可以与其他功能基因挖掘策略进行组合，提高基因定位的效率和准确性，如选择信号检测与全基因组关联分析的组合，选择信号与拷贝数变异检测的组合。
参考文献.png
选择信号的计算（脚本）
群体结构分析样本量差异基本都会有些影响。
选择信号分析的时候没过滤，是因为在目前的ibs分析中，fst的频率抽样，受影响小些。群体结构的时候比较担心连锁的影响。选择信号主要就是要找这些连锁的地方。
一、Fst（基于SNP位点的多态性的选择信号的检测）
Fst/pi都是基于SNP位点的多态性来检测潜在的选择信号区域。Fst/pi ratio基于pi值
按照点来计算：
# 使用vcftools对每一个SNP变异位点进行计算
nohup /home/software/vcftools/src/cpp/vcftools --vcf test.vcf --weir-fst-pop 1_population.txt --weir-fst-pop 2_population.txt --out p_1_2-single &
按照窗口区域来计算（逐个位点计算FST时，可能会出现FST值很高的假阳性信号（中性选择导致），所以考虑到搭载效应同时计算了滑窗FST）
# 计算Fst （此步骤生成的XXX.windowed.weir.fst 可用于画曼哈顿图）
nohup /home/software/vcftools/src/cpp/vcftools --vcf test.vcf --weir-fst-pop 1_population.txt --weir-fst-pop 2_population.txt --out p_1_2_bin --fst-window-size 500000 --fst-window-step 50000 &
# test.vcf是SNP calling 过滤后生成的vcf文件；
# p_1_2_3 生成结果的prefix
# 1_population.txt是一个文件包含同一个群体中所有个体，一般每行一个个体。个体名字要和vcf的名字对应。
# 2_population.txt 包含了群体二中所有个体。
# 计算的窗口是500kb，而步长是50kb （根据你的需求可以作出调整）。我们也可以只计算每个点的Fst，去掉参数（--fst-window-size 500000 --fst-window-step 50000）即可。
# 为Fst排序
nohup sort -k 5 -g XXX.windowed.weir.fst > XXX.windowed.weir.sorted.fst &
# 窗口计数
wc -l XXX.windowed.weir.sorted.fst # 假设为10000
# 取前百分之一（需过滤空集）
nohup tail -n 100 XXX.windowed.weir.sorted.fst > XXX.windowed.weir.sorted.1%.fst &
# 计算Fst后输出表格的格式: 1.CHROM（染色体号） 2.BIN_START （开头）3.BIN_END（结尾） 4.N_VARIANTS（SNP数） 5.WEIGHTED_FST（加权Fst） 6.MEAN_FST（平均Fst）
# SnpEff注释
Fst（群体间分化指数）的值在0-1之间，分化指数越大，两个群体的差异越大。Fst值为0表示两个群体是随机交配的，基因型完全相似。如果Fst值为1，表示两个群体完全隔离。
当Fst出现负值时，表示位点在群体里无分化，一般处理为0再做后续分析。
曼哈顿图可视化
https://www.jianshu.com/p/db932369b2e8
1.制作csv表格
计算Fst后输出表格的格式为: 1.CHROM（染色体号） 2.BIN_START （开头）3.BIN_END（结尾） 4.N_VARIANTS（SNP数） 5.WEIGHTED_FST（加权Fst） 6.MEAN_FST（平均Fst）
awk -F '\t' '{print $1"\t"$2"\t"$5}' XXX.windowed.weir.fst > XXX.windowed.weir-125.fst
# 制作表格共四列:加表头1.SNP 2.CHR 3.BP（BIN_START） 4.P（WEIGHTED_FST加权Fst）
# 第一列为SNP的名字，第二列CHR为所在染色体，第三列BP为染色体上所在位置。要注意如果CHR中存在X，Y，需要转化为数字如赋予23，24等，其中第一列SNP的名字是可选择的，后三列是必须提供的
# 以表格形式打开，加上第一列SNP：snp1、snp2、snp3...
2.R可视化
install.packages('qqman')
library('qqman')
mydata2<-read.table("XXX.windowed.weir-125.fst",head=T)
color_set <- c("#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#66C2A5","#FC8D62")
par(cex=0.9)
manhattan(mydata2,col=color_set,logp = F,ylab='Fst',genomewideline = 0.6) # 在表格中查找前1%的最小值，画出阈值线，本例0.6
3.图解析
y坐标：Fst值
x坐标：染色体号
二、Tajima's D
这个是选择相关的一个参数，大于0代表群体观测杂合度高于预期杂合度，稀有等位基因频率降低（群体收缩或者平衡选择），小于0说明群体观测杂合位点少于预期值，稀有等位基因频率增加（群体扩张或者低频选择）。 也就是说，只有0是正常的，其他都是选择发生。
nohup /home/software/vcftools/src/cpp/vcftools --vcf XXX.vcf --TajimaD 500000 --out TajimaD &
三、π，核苷酸多样性，越大说明核苷酸多样性越高，越低说明两个座位DNA序列差异越小。(在群体内一段一段进行比较，则杂合度越低越受选择)
https://blog.csdn.net/yangl7/article/details/109546077
π用来分析碱基多态性，多态性越低，受选择程度越高。取值时与Fst相反，需要取数据的后1%。
# A，B两个群，仅关注A群，使用 --keep
nohup /home/software/vcftools/src/cpp/vcftools --gzvcf /home/important_data/lzxhebing/cattle_328_chr1-30_resequencing_final_vcf/cattle-328-chr1-30-snp-indel.filter.pass.vcf.gz --keep bohai.txt --window-pi 500000 --window-pi-step 50000 --out bohai_pi &
# 生成两个文件:XXX.windowed.pi、XXX.log，将性染色体删除，PI列降序排序
nohup sort -n -k 5r bohai_pi.windowed.pi > bohai_pi.windowed.sorted.pi &
wc -l bohai_pi.windowed.sorted.pi #49770行
nohup tail -n 500 bohai_pi.windowed.sorted.pi > bohai_pi.windowed.sorted.1%.pi &
# 使用bohai_pi.windowed.sorted.1%.pi文件通过R进行绘图
pi（θπ） 选择消除分析图可视化
library(ggplot2)
data<-read.table("allsample_500K.windowed.pi",header=T)
for (i in 1:29){
result_name = paste("chr",i,".png",sep="")
png(result_name,width=1500,height=300)
chrom = subset(data,CHROM==i)
xlab = paste("Chromosome",i,"(MB)",sep="")
p <- ggplot(chrom,aes(x=BIN_START/1000000,y=PI,group=Phenotype,color=Phenotype,shape=Phenotype)) + geom_line(size=0.5) + xlab(xlab)+ ylab("Pi") + theme_bw()
print(p)
dev.off()
}
https://www.plob.org/article/21645.html
选择信号结合.png
四、iHS检测群体内的选择信号(基于单体型haplotype的选择信号的检测)
在selective sweeps选择过程中，有些强烈受到选择的位点variants由于LD的因素会连带着其附近的位点variants一起被保留，并且不会受到重组recombination的打断。一些低重组区域的haplotypes的长度会高于那些高重组区域的haplotypes的长度。因此，对比同一genomic区域在不同群体中的haplotype的长度可以用来判断是否受到选择。例如：在一个群体内部，如果某一个体强烈受到选择，其haplotype的长度会远长于其它个体；同理，对于两个群体之间的比较，某一群体受到选择，则其基因组中的受选择区域的haplotypes会比未受到选择群体中的haplotypes更长。
例如：使用selscan软件计算了澳洲野犬的iHS，并通过常染色体上20 kb的滑动窗口通过规范（在selscan的软件中）对分数进行归一化。如果其中30%的站点的iHS绝对值高于阈值(或iHS绝对值的前1%)，我们将窗口确定为候选区域。（参考Genomic regions under selection in the feralization of the dingoes）
https://www.jianshu.com/p/23ab344d66f7
1.数据文件准备
# 提取vcf 压缩
nohup /home/software/bcftools/bcftools view -S id.txt test.vcf.gz -Ov > sampleid.vcf &
nohup bgzip -c -f sampleid.vcf > sampleid.vcf.gz &
# 过滤，去多等位基因位点，最大缺失率（max-missing）<10%即检出率>90%的SNP位点，最小等位基因频率（MAF）>0.03（哈代温伯格平衡检验的P值大于10^-6 -hwe）；--min-alleles 2 --max-alleles 2 去除多等位基因和等位基因缺失
nohup /home/software/vcftools/src/cpp/vcftools --vcf sampleid.vcf --remove-indels --max-missing 0.9 --maf 0.03 --min-alleles 2 --max-alleles 2 --recode --recode-INFO-all --out sampleid.miss.snp.recode.vcf &
# 去其他染色体，由于参考基因组版本问题，这里存在很多位于随机或未知染色体上的变异位点，所以需要先过滤掉这些位点
nohup /home/software/vcftools/src/cpp/vcftools --vcf sampleid.miss.snp.recode.vcf --recode --recode-INFO-all --not-chr XXX --out test_noincludeXXX &
# 转格式
nohup /home/software/plink --vcf test_noincludeXXX.miss.snp.recode.vcf --recode --allow-extra-chr --out test_noincludeXXX.miss.snp.recode &
# 换map文件位置名称 ，空格分隔
awk '{print$1,$1":"$4,$3,$4}' test_noincludeXXX.misssnp.recode.map > test_noincludeXXX.misssnp.recode.pos.map
# 转换成vcf格式
nohup /home/software/plink --allow-extra-chr --file test_noincludeXXX.misssnp.recode.pos --make-bed --out test_noincludeXXX.misssnp.recode.pos &
# 更换染色体，加头，压缩
2.vcf文件的phasing，beagle转换格式
由于后面计算时selscan软件对vcf文件格式要求：前九列同标准的vcf格式一样，后面紧接着的个体的基因型尽量是phasing后的。文件格式如下：< quality ><individual 1 genotype > ... < individual N genotype >，因此首先使用beagle软件转换格式
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
nohup perl -pe "s/\s\.:/\t.\/.:/g" test.miss.snp.recode.vcf | bgzip -c > out.vcf.gz &（重测序数据会遇到这种问题）
nohup java -jar /home/software/beagle.25Nov19.28d.jar gt=out.vcf.gz out=test.beagle ne=44 > test.beagle.vcf.log &
#有报错：java.lang.OutOfMemoryError: Java heap space 解决办法：
nohup java -jar -Xmn12G -Xms24G -Xmx48G /home/software/beagle.25Nov19.28d.jar gt=test_noincludeXXX.misssnp.recode.pos-chr-36head.vcf out=test.beagle ne=44 > test.beagle.vcf.log &
3.准备分染色体的vcf文件，单倍型以及基于单倍型的分析都需要逐条染色体运行，所以需要把vcf文件以及map按照染色体ID进行拆分。
vi vcftools.sh
for i in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30;
do /home/software/vcftools/src/cpp/vcftools --vcf test.beagle.vcf --chr ${i} --recode --recode-INFO-all --out ${i}.test;
done
nohup bash vcftools.sh & # 得到${i}.test.recode.vcf文件
4.准备分染色体的map文件
vi map.bash
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30;
do /home/software/vcftools/src/cpp/vcftools --vcf ${k}.test.recode.vcf --plink --out ${k}.MT;
done
nohup bash map.bash & # 得到${k}.MT.map文件
5.map计算遗传距离
vi map2.sh
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 ;
do awk 'BEGIN{OFS=" "} {print 1,".",$4/1000000,$4}' ${k}.MT.map > ${k}.MT.map2;
done
nohup bash map2.sh & # 得到${k}.MT.map2文件
6. 计算iHS
vi iHS.sh
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30;
do /home/software/selscan/bin/linux/selscan --ihs --vcf ${k}.test.recode.vcf --map ${k}.MT.map2 --out ${k}.iHS;
done
nohup bash iHS.sh & # 得到 ${k}.iHS.ihs.out文件，运行时间久！
（1）# 重测序数据直接计算iHS（对每个100K窗口单倍型进行打分，标准化使其正态分布，并统计绝对值>2的SNP所占的比例）
# 提取结果
vi hu.sh
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 ;
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
nohup bash hu.sh &
# 标准化（将norm文件第7列数据即iHS数据进行统计，统计绝对值>2的SNP所占的比例）
vi hu2.sh
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 ;
do /home/software/selscan/bin/linux/norm --ihs --files ${k}.hu.ihs.out --bp-win --winsize 100000 ;
done
nohup bash hu2.sh &
# 提取结果
vi hu3.sh
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 ;
do awk '{print '${k}',$1,$2,$4}' ${k}.hu.ihs.out.100bins.norm.100kb.windows > ${k}.hu3.ihs.out.100kb.windows;
done
nohup bash hu3.sh &
# 得到两个文件 XXX.hu.ihs.out.100bins.norm.100kb.windows与XXX.hu.ihs.out.100bins.norm
cat ./*.hu3.ihs.out.100kb.windows > all.hu3.ihs.out.100kb.windows
sort -k 4n,4 all.hu3.ihs.out.100kb.windows > all.hu3.ihs.out.100kb.windows.sort
# 按照第四列取前5%的位点，画图（此处也可使用XXX.hu.ihs.out.100bins.norm文件画图，但数据量过大，推荐使用XXX.hu.ihs.out.100bins.norm.100kb.windows）
制作表格，格式如下，逗号分隔保存：
X SNP CHR BP P
1 snp1 1 2400001 0.265873
2 snp2 1 2500001 0.449057
3 snp3 1 2600001 0.596059
# R
install.packages('qqman')
library('qqman')
data1 <- read.table("all.iHS.csv",header = T,sep = ',')
color_set <- c("#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#66C2A5","#FC8D62")
par(cex=0.9)
manhattan(data1,logp = F,col=color_set)
# SNPEFF注释，参考下面
（2）以下芯片数据操作
# 提取结果 ( 可更改为 awk '{print '${k}',$2,$6}' )
vi awk2.sh
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30;
do awk '{print '${k}',$2,$6}' ${k}.iHS.ihs.out > ${k}._.ihs.out ;
done
nohup bash awk2.sh & # 得到{k}._.ihs.out
cat ./*._.ihs.out > all.test.ihs.out # 合并所有染色体结果
# 根据第四列（iHS值）从大到小排列：
sort -k 1n,1 -k 4n,4 all.test.ihs.out > all.test.ihs.sort.out
7.制作窗口文件（芯片数据）
设置窗口滑动区间，计算落到各个区间的iHS总和，计算平均值，提取iHS绝对值为top 1%的窗口区间，认定为选择信号强的区间，根据文献，设置500kb的窗口；计算各个窗口内的均值，计算期望和标准差，利用以下公式得到标准化的iHS：
iHS.png
# 使用excel表格中的函数制作50k窗口位置文件，格式如下:
1 0 50000
1 50001 100001
1 100002 200002
...
# 由于窗口文件是从表格中粘贴的，需要注意分隔符的设置，将空格键替换为tab键，并删除最末尾的tab键：
vi sed.sh
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30;
do sed -i 's/ /\t/g' ${k}.window.bed;
sed -i 's/[ \t]*$//g' ${k}.window.bed;
cat ${k}.window.bed | tr -s [:space:] > ${k}._window.bed；
done
得到{k}._window.bed
8.制作区间文件（芯片数据）
# 制作位置区间bed文件
vi ihs.bed.sh
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30;
do awk '{print $1,$2,$2+1,$3}' ${k}._.ihs.out > ${k}._.ihs.bed;
done
nohup bash ihs.bed.sh &
# 取交集,需要统计各个SNP落在哪个窗口，用bedtools统计两个区间（SNP位置，SNP+1为位置区间；窗口为第二个区间）的交集。
vi intersect.sh
for k in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27;
do /home/software/bedtools2/bin/bedtools intersect -a ${k}._window.bed -b ${k}._.ihs.bed -wao > ${k}.window.ihs.intersect.bed;
done
nohup bash intersect.sh &
9.用ihs_window_mean_xunhuan.py脚本计算各个窗口（即第2，3列区间）内iHS的均值（芯片数据）
vi ihs_window_mean_xunhuan.py
# -*- coding: utf-8 -*-
"""
Created on Sun Jul 11 22:16:50 2021
@author: Administrator
"""
import numpy as np
import glob
dirname ="/*.txt"
files = glob.glob(dirname)
for filename in files:
#删除不需要的数据
data = []
fp=open(filename)
lines = fp.readlines()
for line in lines:
if line[-2] != '0':
data.append(line)
fp.close()
fp=open(filename, 'w')
for row in data:
fp.write(row)
fp.close()
#计算平均值
data = np.loadtxt(filename, dtype=float, delimiter='\t')
result = np.empty((0,8))
i = 0
count = 0
temp = data[0]
for row in data:
if temp[1] == row[1]:
count += row[6]
i += 1
else:
temp[6] = count / i
result = np.concatenate((result, np.expand_dims(temp,0)), axis=0)
temp = row
i = 1
count = row[6]
temp[6] = count / i
result = np.concatenate((result, np.expand_dims(temp,0)), axis=0)
filename = filename[:-3] + 'csv'
np.savetxt(filename, result, delimiter=',')
10.用得到的窗口ihs均值（all.window.ihs.intersect.csv），在R中计算标准iHS值（芯片数据）
setwd("")
sink("window_biaozhun_ihs.csv")
data=read.table("all_window_ihs.csv")
options(max.print=30000)
scale(data,center=TRUE,scale=TRUE)
attr(1,"scaled:center")
attr(1,"scaled:scale")
sink()
# 得到window_biaozhun_ihs.csv文件，包含一列数据，即标准化的iHS值
11.画曼哈顿图
install.packages('qqman')
library('qqman')
data1 <- read.table("/window_biaozhun_ihs3.csv",header = T,sep = ',')
color_set <- c("#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#66C2A5","#FC8D62")
par(cex=0.9)
manhattan(data1,col=color_set,ylab='fst'，genomewideline = 1.96)
SnpEff注释
https://www.jianshu.com/p/b2b45d2523db
https://www.cnblogs.com/zhanmaomao/p/10964636.html
1.安装 ，构建数据库
#会产生一个snpEff目录 所有的程序都在这里面
下载 https://pcingola.github.io/SnpEff/
unzip snpEff_latest_core.zip
2.配置自己的基因组和注释文件
以绵羊（sheep）参考基因组为例：打开snpEFF文件夹下的snpEff.contig，在Third party databases下面增加新的物种信息：
#-------------------------------------------------------------------------------
# Third party databases （原来的）
#-------------------------------------------------------------------------------
#sheep genome, version sheepv1.0 (新加的)
sheepv1.0.genome : sheep （新加的）
# Databases for human GRCh38 (hg38) （原来的）
3.注释文件为 gff 格式
# 在SnpEff目录下，创建目录 data, data/sheepv1.0, data/genomes
mkdir data && cd data
mkdir sheepv1.0
mkdir genomes
# 将 gff3 (.gff) 文件放入sheepv1.0目录下，并改名为 genes.gff
# 将基因组序列文件（.fasta）放入 genomes 目录下，并改名为 sheepv1.0.fa
nohup java -jar snpEff.jar build -gff3 -v sheepv1.0 & # 在 snpEff 目录下，执行命令
4.注释文件为 gtf 格式
# 如果注释文件为.gtf，可参考 gff 中的步骤，要将注释文件重新命名为 genes.gtf
nohup java -jar snpEff.jar build -gtf22 -v sheepv1.0 &
5.开始注释 (提取的要注释的文件中染色体应为NC格式且文件以tab键分隔，提取三列即可：染色体号、起始位置、终止位置，若为SNP位点则第三列用1代替，且位于snpEff目录下操作命令）
（1）若文件为按照窗口计算Fst后输出的文件，则提取为bed文件格式，bed格式(取染色体号，起始位置和结束位置和Fst 的值)
# 计算Fst后输出表格的格式: 1.CHROM（染色体号） 2.BIN_START （开头）3.BIN_END（结尾） 4.N_VARIANTS（SNP数） 5.WEIGHTED_FST（加权Fst） 6.MEAN_FST（平均Fst）
awk -F '\t' '{print $1"\t"$2"\t"$3"\t"$5}' XXX.windowed.weir.sorted.1%.fst > XXX.windowed.weir.sorted.1%.bed
#如果bed文件中染色体为1.2.3格式，需替换染色体
nohup java -Xmx4g -jar ./snpEff.jar sheepv1.0 -i bed XXX.windowed.weir.sorted.1%.bed-chr.bed > XXX.windowed.weir.sorted.1%.bed-chr.anno.bed &
（2）若文件为按照位点计算（例如：fst按位点计算、重测序SNP数据、重测序INDEL数据）但重测序使用中，得到的注释文件不理想
# 用awk提取fst中的snp 位置信息，并用vcftools对应回vcf文件中
nohup awk '{print $1,$4,$4+1}' test.map > test-awk.bed &
# 注释
nohup java -Xmx4g -jar ./snpEff.jar sheepv1.0 -i bed test-awk.bed > test-awk.anno.bed &
（3）若文件为vcf文件 (重测序SNP数据、重测序INDEL数据)
比如只关注CDS中的注释信息，不考虑上游、下游、UTR、基因间区等信息
可以选择以下参数简化输出
-no-downstream
-no-upstream
-no-utr
-no-intergenic
-no-intron
nohup java -jar snpEff.jar ann -no-utr -no-downstream -no-upstream -no-intergenic sheepv1.0 test.vcf.gz > snpeff.vcf &
# 此处可将变异信息删除做一个小的vcf文件，只保留到vcf的info信息即可，则最后出现的文件小一些
6.结果说明
(1)SnpEff结果解读 - 简书 (jianshu.com)
第一列：CHROM 发生突变的染色体ID。
第二列：POS：发生突变的染色体上的具体位置。
第三列：ID 可以在后面的注释信息中找到geneID。
第四列：REF 参考基因组上的碱基或者序列。
第五列：ALT 发生突变后的碱基或者序列。
第六列：QUAL 得分，Phred格式的数值。代表着此为点是纯和的概率。此值越大，概率越低，代表着此为点是变异位点的可能性越大。
第七列：FILTER 过滤情况 ：一般分析后的结果都为PASS，则表示该位点是变异位点。
第八列：INFO 变异位点的相关信息。
第九列：FORMAT 变异位点的格式：比如 GT:PL:ADF:ADR:AD:GP:GQ
第十列：SAMPLEs 各个样本的值，这些值对应着第9列的各个部分，不同部分之间的值使用冒号分隔。
而SnpEff结果文件中，在INFO这一列中，增添了一个字段，ANN
ANN=A|upstream_gene_variant|MODIFIER|AT1G69210|AT1G69210|transcript|AT1G69210.1|protein_coding||c.-3686C>T|||||3686|,A|upstream_gene_variant|MODIFIER|AT1G69210|AT1G69210|transcript|AT1G69210.2|protein_coding||c.-3686C>T|||||3686|,A|downstream_gene_variant|MODIFIER|SP1L2|AT1G69230|transcript|AT1G69230.1|protein_coding||c.2799C>T|||||2509|,A|downstream_gene_variant|MODIFIER|MES15|AT1G69240|transcript|AT1G69240.1|protein_coding||c.4371C>T|||||4138|,A|downstream_gene_variant|MODIFIER|SP1L2|AT1G69230|transcript|AT1G69230.2|protein_coding||c.*2799C>T|||||2568|,A|intron_variant|MODIFIER|SIK1|AT1G69220|transcript|AT1G69220.1|protein_coding|8/17|c.1323+52C>T||||||,A|intron_variant|MODIFIER|SIK1|AT1G69220|transcript|AT1G69220.2|protein_coding|8/17|c.1242+52C>T||||||
新增的字段由|进行间隔，并且这个字段中包含了突变位点的注释信息，因此非常重要。
重点关注以下几点：
Allele ANN=A 说明：突变后的碱基是A。
Annotation upstream_gene_variant 造成的基因上游的突变。downstream_gene_variant 造成的基因下游的突变。
Annotation_Impact 突变位点造成的影响：一般可以划分为四类 HIGH，MODERATE，LOW，MODIFILER 一般突变位点造成HIGH最好，往后依次效果越低。
Gene_Name 基因名称。
Gene_ID 基因ID。
Feature_Type 想要分析的特征类型，transcript, motif, miRNA 等。
Feature_ID 根据Feature Type指定的特征，给出对应的ID。
Transcript_BioType 转录本类型, 通常采用Ensembl数据库的转录本类型。
(2)snpEff_genes.txt和snpEff_summary.html这两个文件记录总结性信息比较简单。
*.ann.vcf 是一个注释结果文件，其就在vcf的info信息新添加了anno一列信息，其具体每个值含义如下：
1.Allele：突变之后的碱基，第一个突变位点由T碱基突变成了C碱基，对应Allel的值为C；
2.Annotation：由sequence ontology定义的突变类型；
3.Annotation_Impact：对变异位点有害程度的简单评估，取值有HIGH, MODERATE, LOW, MODIFIER 4种；
4.Gene_Name：基因名字；
5.Gene_ID：基因ID；
6.Feature_Type:想要分析的特征类型，transcript, motif, miRNA 等;
7.Feature_ID:根据Feature Type指定的特征，给出对应的ID;
8.Transcript_BioType:转录本类型, 通常采用Ensembl数据库的转录本类型;
9.Rank:只有当变异位点位于基因区域时才有值，会给出变异位点所处的exon/intron的编号和该基因的exon/intron的总数，比如一个突变位点位于基因的第3个exon上，该基因一共有12个exon, 对应的Rank的值为3/12,当变异位点位于基因区域以外时，该字段的值为空;
10.HGVS.c:采用HGVS标准命名的基因水平的变异情况;
11.HGVS.p:采用HGVS标准命名的蛋白质水平的变异情况，只有当突变位点位于编码区是才会有值;
12.cDNA.pos/cDNA.length:突变位点在cDNA上的位置/cDNA的总长度;
13.CDS.pos/CDS.length:突变位点在CDS上的位置/CDS的总长度;
14.AA.pos/AA.length:突变位点在氨基酸序列上的位置/氨基酸序列的总长度;
15.Distance:变异位点与最近的特征的距离，当变异位点位于基因间区时，会给出与最近的基因之间的距离；当变异位点位于exon区域时，会给出与最近的内含子边界的距离，不同的情况，距离的定义不同;
16.ERRORS/WARNINGS/INFO:对注释结果的可靠程度进行评估。
其他注释方法
一、ensembl数据库BioMart注释
可观察到所需的目标物种位点注释信息
1.ensembl——BioMart
dateset.png
2.Attributes #按需选择
attributes1.png
attributes2.png
3.Results
results.png
二、ensembl vep注释
1.通过位点计算Fst，再将利用Fst文件对应出的那1%的vcf文件的前5列，用excel打开复制前5列，（从染色体那一行开始复制）使用NC染色体
1.png
