基于單體型haplotype的選擇信號的檢測。在selective sweeps選擇過程中擂错，有些強烈受到選擇的位點variants由于LD的因素會連帶著其附近的位點variants一起被保留缕陕，并且不會受到重組recombination的打斷粱锐。一些低重組區(qū)域的haplotypes的長度會高于那些高重組區(qū)域的haplotypes的長度。因此扛邑，對比同一genomic區(qū)域在不同群體中的haplotype的長度可以用來判斷是否受到選擇怜浅。例如：在一個群體內(nèi)部，如果某一個體強烈受到選擇蔬崩，其haplotype的長度會遠長于其它個體恶座；同理，對于兩個群體之間的比較沥阳，某一群體受到選擇奥裸，則其基因組中的受選擇區(qū)域的haplotypes會比未受到選擇群體中的haplotypes更長。

選擇壓力分析基本原理

image

原始群體中沪袭，遺傳多樣性是十分高的湾宙，整個序列的核酸diversity都高。而在受到選擇之后冈绊，diversity會發(fā)生波動侠鳄。核酸多樣性下降 可能就是由于under selection導致的。

在演化/馴化過程中死宣，如果某一基因X占優(yōu)勢伟恶，即X的基因型占據(jù)主導地位，則基因X所在區(qū)域的雜合率/多樣性會顯著下降毅该。本質(zhì)就是 比較基因組不同區(qū)域多樣性（雜合率）的變化

• 群體遺傳關(guān)心的問題：
  • 遺傳結(jié)構(gòu)（phylogeny+structure）
  • 基因組上受選擇區(qū)域：群體水平基因組不同位置的區(qū)域遺傳多樣性變化的規(guī)律（例如：Pi博秫、Tajima's D, Fst）
• 變異類型：
  • 中性突變（同義、相同類型的氨基酸眶掌、不影響環(huán)境適應(yīng)性）：平衡選擇挡育，這種基因型頻率是大致恒定的
  • 有利突變（正選擇）：選擇掃蕩（Selective sweep），與有利突變的中性突變的頻率會顯著提升
    
    selective sweeps
  • 有害突變（負選擇）：背景選擇（negative selection/background selection/ purifying selection） 是潛在的噪音

負選擇會對正選擇有一定的干擾作用朴爬，都能產(chǎn)生大量的低頻突變即寒，但是正選擇會產(chǎn)生相對較多的高頻突變。

兩個亞群體之間的比較

多樣性水平在亞群間比較召噩，一般包括線性相關(guān)分析母赵、亞群體間的差異比較兩類。動植物重測序多是后者具滴。Fst/pi ratio基于pi值凹嘲。

1. 群體分化程度Fst (Fixation index): 比較兩個亞群體間的Pi值和亞群體內(nèi)的Pi值的差異。
   • 由PI值計算演變來（序列兩兩差異取均值）
   • 兩個亞群體在某一段seq區(qū)域的差異度构韵。0是無差異周蹭，數(shù)值越大溯革，則說明兩個亞群體之間已經(jīng)發(fā)生了明顯的分化（亞群內(nèi)個體相似，亞群間差異大）

Fst=(\pi(between) - \pi(within))/ \pi(between) 

2. 多樣性變化倍數(shù)Pi ratio：某區(qū)間在亞群間的多樣性差異的倍數(shù)谷醉，簡單粗暴，就關(guān)注多樣性值的高低變化冈闭。
   • 例如野生群體A/栽培群體B俱尼；野生群體A的多樣性較高，而栽培群體B的多樣性較低萎攒，所以多樣性降低最顯著的基因組區(qū)域遇八，就與馴化改良基因相關(guān)
3. 其它比較值：ROD值、XP-CLR值等耍休。而多個品種間的比較分化差異的di值

基于haplotype單體型的比較

前面pi/fst等都是基于SNP位點的多態(tài)性來檢測潛在的選擇信號區(qū)域刃永。另一種方法是基于單體型haplotype的選擇信號的檢測。在selective sweeps選擇過程中羊精，有些強烈受到選擇的位點variants由于LD的因素會連帶著其附近的位點variants一起被保留斯够，并且不會受到重組recombination的打斷。一些低重組區(qū)域的haplotypes的長度會高于那些高重組區(qū)域的haplotypes的長度喧锦。因此读规，對比同一genomic區(qū)域在不同群體中的haplotype的長度可以用來判斷是否受到選擇。例如：在一個群體內(nèi)部燃少，如果某一個體強烈受到選擇束亏，其haplotype的長度會遠長于其它個體；同理阵具，對于兩個群體之間的比較碍遍，某一群體受到選擇，則其基因組中的受選擇區(qū)域的haplotypes會比未受到選擇群體中的haplotypes更長阳液。

檢測haplotype的選擇信號最好利用定相phased后的數(shù)據(jù)集怕敬。方法有EHH和CLR法。這里利用R包中的rehh包進行分析帘皿。rehh有強大的說明和教程文檔赖捌，后續(xù)深入了解其原理時值得進一步學習研究。rehh tutorial

rehh的實踐

1. 讀取數(shù)據(jù)矮烹。分別讀取兩個群體的phased的vcf數(shù)據(jù)
   • polarize_vcf設(shè)為FALSE. because we have not used an outgroup genome to set our alleles as derived or ancestral
   • 根據(jù)maf進行過濾位點

# read in data for each species# house
house_hh <- data2haplohh(hap_file = "./house_chr8.vcf.gz",polarize_vcf = FALSE)
# bactrianus
bac_hh <- data2haplohh(hap_file = "./bac_chr8.vcf.gz",polarize_vcf = FALSE)

## filter on MAF - 0.05
house_scan <- scan_hh(house_hh_f, polarized = FALSE)
bac_scan <- scan_hh(bac_hh_f, polarized = FALSE)

image

2. 計算計算單體型的iHS值越庇。
   • polarized = FALSE
   • freqbin =1 if we know the ancestral allels or derived allels. rehh can apply weights to different bins of allele frequencies in order to test whether there is a significant deviation in the iHS statistic.
   • log Pvalue用于檢測outliers值點，

## perform haplotype genome scan- iHS
house_scan <- scan_hh(house_hh_f, polarized = FALSE)
bac_scan <- scan_hh(bac_hh_f, polarized = FALSE)

## perform iHS on house
house_ihs <- ihh2ihs(house_scan, freqbin = 1)
### plot the iHS statistics
ggplot(house_ihs$ihs, aes(POSITION, IHS)) + geom_point()
### plot the log P-value 
ggplot(house_ihs$ihs, aes(POSITION, LOGPVALUE)) + geom_point()

iHS statistics plot

iHS Pvalue for outliers

3. 計算群體之間的EHH值 xpEHH
   • 計算 cross-population 的 EHH test奉狈。利用之前iHS檢測出的iES值卤唉。
   • include_freq=T,we get the frequencies of alleles in our output, which might be useful if we want to see how selection is acting on a particular position

## perform xp-ehh
house_bac <- ies2xpehh(bac_scan, house_scan,
                       popname1 = "bactrianus", popname2 = "house",
                       include_freq = T)
#PLOT the xpEHH values
ggplot(house_bac, aes(POSITION, XPEHH_bactrianus_house)) + geom_point()

image

負數(shù)值代表在pop2(house in this case)中的強烈的選擇信號。

image

4. 檢測受選擇區(qū)域中的haplotype structure

# find the highest hit 以最顯著的位點做示例
hit <- house_bac %>% arrange(desc(LOGPVALUE)) %>% top_n(1)
# get SNP position
x <- hit$position # POSITION 19191935

# detect the position which the selection occured in the haplotype objects
marker_id_h <- which(house_hh_f@positions == x)
marker_id_b <- which(bac_hh_f@positions == x)

# calculate the bifurcation of haplotypes around our site of selection
house_furcation <- calc_furcation(house_hh_f, mrk = marker_id_h)
bac_furcation <- calc_furcation(bac_hh_f, mrk = marker_id_b)
# plot the bifurcation plot
plot(house_furcation, xlim = c(19.18E+6, 19.22E+6))
plot(bac_furcation, xlim = c(19.18E+6, 19.22E+6))

image

house_furcation

image

bac_furcation

# calculate the haplotype length around our signature of selection.
house_haplen <- calc_haplen(house_furcation)
bac_haplen <- calc_haplen(bac_furcation)

# see how haplotype structure differs between our two populations.

plot(house_haplen)
plot(bac_haplen)

image

house_furcation

image

bac_furcation

the blue haplotype is much larger around this target and is also more numerous in the European house sparrow.

5. 輸出數(shù)據(jù)

# write out house bactrianus xpEHH
house_bac <- as_tibble(house_bac)
colnames(house_bac) <- tolower(colnames(house_bac))
write_tsv(house_bac, "./house_bac_xpEHH.tsv")

https://speciationgenomics.github.io/