基于单倍型的选择信号检测方法

Posted on | In

选择的主要效应是改变种内或者种间的变异性水平。 当一个群体中出现有利变异时，正选择作用会使有利等位基因频率在短时间内提高，由于连锁不平衡作用，有利等位基因邻近的中性位点也会倾向于一起被选择， 这种现象称为“搭乘效应”(Hitchhiking Effect)。有利等位基因在群体中频率升高会消除或者减弱与其连锁的中性位点的变异，使其多态性也降低，这一过程被称为选择性清除(Selectivesweep) 。

选择信号的检测基于中性检验的思想， 其检测方法主要分为三类： 基于频率谱的方法(Frequency Spectrum–Based)， 如 Tajima’ D (Tajima, 1989)， Fay and Wu’s H(Fay and Wu, 2000) 等； 基于连锁不平衡的方法(Linkage Disequilibrium–Based)， 如EHH (Sabeti et al., 2002)， iHS (Voight et al., 2006)， XPEHH (Sabeti et al., 2007)等；基于群体分化的方法(Population Differentiation–Based)，如 FST(Weir and Cockerham,1984)， XP-CLR (Chen et al., 2010)， hapFLK(Fariello et al., 2013)等。

本文主要介绍基于连锁不平衡的选择信号检测方法。群体中不同基因座等位基因之间的非随机关联被称为连锁不平衡(Linkage disequilibrium, LD)。 正向选择作用下， 有利突变的等位基因频率会迅速升高， 并与其邻近的“搭乘”位点保持强连锁不平衡关系， 有利位点与其邻近位点会形成一个单倍型，直至重组将它们之间的 LD 打破。
基于 LD 的选择信号检测方法寻求强 LD区域的扩展（即长范围的单倍型），该方法适用于近期固定或者接近固定的选择位点，否则重组事件会打破 LD，单倍型范围也会缩短(Vitti et al., 2013)。 基于 LD 检测选择信号的方法中应用最广泛的是扩展单倍型纯合（ extended haplotypehomozygosity , EHH） (Sabeti et al., 2002)及其一系列衍生方法。 EHH 的定义如下：从一个核心区域（如推断的受选择位点）向两侧扩展距离 x， 两条随机选取的包含核心单倍型的染色体，从核心区域到距离 x 的整个区间为同源相同 (identical by descent， IBD)的概率。 在 EHH 基础上， Voight 等(2006)提出了单倍型纯合性积分(integrated haplotype score, iHS)的检验方法, iHS 通过计算衍生等位基因(derived)和祖先等位基因 (ancestral) EHH 积分的比值进行检验。在 EHH 和 iHS 的基础上， Sabeti 等(2007)又提出了与 iHS 类似的检验方法，即群体间扩展单倍型纯合(Cross Population Extended Haplotype Homozogysity,XPEHH)的方法， XPEHH 引入了群体的比较策略。

selscan是一款用于基于LD计算选择信号的简单易用的软件包，可以计算EHH、iHS、XPEHH、nSL等统计量，下面给出该软件进行选择信号计算的一般流程。

  • 准备input files
1
2
3
4
5
6
7
8
9
#plink转fastphase，注意selscan软件需要0/1编码的格式#
plink --bfile gwas_v2 --chr-set 26 --recode fastphase --out gwas_v2

#phasing#
for ((k=1;k<27;++k));do time fastPHASE -T10 -ogwas_v2.chr-${k} gwas_v2.chr-${k}.recode.phase.inp;done
#phased后的output文件处理*.recode.phase.inp#
for ((k=1;k<27;++k));do sed -e '1,21d' -e '/#/d' -e '$d' gwas_v2_chr${k}_hapguess_switch.out > v2_chr${k}_hapguess_switch.out;done
#按照染色体分割map文件，可在R中准备

  • iHS/XPEHH/nsl计算
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#基本用法
selscan --ihs --hap <haplotypes> --map <mapfile> --out <outfile> 
selscan --xpehh --hap <pop1 haplotypes> --ref <pop2 haplotypes> --map <mapfile> --out
selscan --nsl --hap <haplotypes> --map <mapfile> --out <outfile>

#XPEHH#
for ((k=1;k<27;++k));do ./selscan --xpehh --hap sel.case_${k}_hapguess_switch.out --ref sel.control_${k}_hapguess_switch.out --map chr${k}.map --out case_c_xpehh.chr${k};done

#iHS#
for ((k=1;k<27;++k));do ./selscan --ihs --hap sel.case_${k}_hapguess_switch.out --map chr${k}.map --out case_c_ihs.chr${k};done
for ((k=1;k<27;++k));do ./selscan --ihs --hap sel.control_${k}_hapguess_switch.out --map chr${k}.map --out control_c_ihs.chr${k};done

#nsl#
for ((k=1;k<27;++k));do ./selscan --nsl --hap sel.case_${k}_hapguess_switch.out --map chr${k}.map --out case_c_nsl.chr${k};done
for ((k=1;k<27;++k));do ./selscan --nsl --hap sel.control_${k}_hapguess_switch.out --map chr${k}.map --out sweep_out/control_c_nsl.chr${k};done

#合并文件
for(i in 1:26){
xpehh<-read.table(paste0("case_c_xpehh.chr",i,".xpehh.out"),header=T)
xpehh$chr<-i
write.table(xpehh,file=paste0("case_c_xpehh.chr",i,".xpehh.txt"),quote=F,col.names=F,row.names=T,sep=" ")
}
awk 1 case_c_xpehh.chr*.xpehh.txt > case-control.xpehh.txt

  • 作图
1
2
3
4
5
6
7
8
9
10
11
12
xp<-read.table("sweepout/case-control.xpehh.txt",header=T)
library(gap)
xp1<-data.frame(xp[,c("chr","pos","xpehh")])
xp2<-xp1[order(xp1$chr),]
mhtplot(xp2,control=mht.control(logscale=FALSE),pch=19)
write.table(xp2,file="xpehh.out",quote=F,col.names = T,row.names = F,sep=" ")

par(las=2, xpd=TRUE, cex.axis=1.8, cex=0.45)
ops<-mht.control(logscale=FALSE,col=rep(c("lightblue","hotpink2"),14),yline=2,xline=3,
                 cutoffs=0,labels=c(1:26),srt=0)
mhtplot(xp2,ops,pch=19,ylab="XP-EHH Score")
axis(2,pos=2,at=-4:6)