轉(zhuǎn)載自https://mp.weixin.qq.com/s/1IM6KSCGIq0cUJhKKXVcIw

背景

選擇清除分析是常用于鑒定群體基因組水平受選擇區(qū)域的一種方法晤柄，其背后的理論是：在受選擇群體中氏仗，會出現(xiàn)有利等位基因的頻率快速上升股冗，與其相鄰位點(diǎn)的多態(tài)性快速下降的現(xiàn)象葛闷。基于此憋槐，我們可以通過計算群體內(nèi)基因組區(qū)域的等位基因頻率的變化情況來鑒定受選擇基因組區(qū)域。常用群體內(nèi)檢測指標(biāo)的計算方法大致分為三種：1.基于核苷酸多態(tài)性降低的π淑趾、θw秦陋；2.基于分離位點(diǎn)頻率的Tajima’D；3.基于連鎖不平衡增加的EHH治笨、iHS驳概。以上三類指標(biāo)對應(yīng)于基因組受選擇特征的三個維度，而后才有了群體間的選擇指標(biāo)：1.由π衍生的π ratio旷赖、ROD顺又、Fst；2.由EHH衍生的XPEHH等孵。本文主要介紹XPEHH的實(shí)操稚照。

XPEHH分析

XPEHH是基于群體單倍型的選擇清除分析，有很多工具可以實(shí)現(xiàn)XPEHH方法俯萌，比如R包rehh果录、Linux平臺下selscan軟件。在這里我用selscan軟件完成咐熙。

1.軟件下載與安裝

selscan下載網(wǎng)址：https://github.com/szpiech/selscan/releases選擇最新版本弱恒，復(fù)制鏈接地址

Linux下輸入

2.查看分析所用數(shù)據(jù)文件及格式

selscan軟件幫助文檔可以在此下載：https://github.com/szpiech/selscan/blob/master/manual/selscan-manual.pdf

通過軟件文檔我們可以了解到XPEHH分析有三種輸入方式：

wget?-c?https://github.com/szpiech/selscan/releases/download/v1.3.0/selscan-linux-1.3.0.tar.gz

tar?-zxvf?selscan-linux-1.3.0.tar.gz?#解壓壓縮包

cd?./selscan-linux-1.3.0

pwd

selscan=軟件所在絕對路徑

2.查看分析所用數(shù)據(jù)文件及格式

selscan軟件幫助文檔可以在此下載：https://github.com/szpiech/selscan/blob/master/manual/selscan-manual.pdf

通過軟件文檔我們可以了解到XPEHH分析有三種輸入方式：

selscan?--xpehh?--hap?<pop1?haplotypes>?--ref?<pop2?haplotypes>?--map?<mapfile>?--out?<outfile>?

#hap文件格式輸入

selscan?–xpehh?–vcf?<pop1?vcf>?--vcf-ref?<pop2?vcf>?--map?<mapfile>?--out?<outfile>?#vcf文件格式輸入

selscan?–xpehh?–tped?<pop1?tped>?--tped-ref?<pop2?tped>?--map?<mapfile>?--out?<outfile>?#tped文件格式輸入

因?yàn)槲业臄?shù)據(jù)文件是所有個體總的vcf文件，所以選用vcf文件作為輸入文件格式最為方便棋恼。軟件對vcf文件格式要求：前九列同標(biāo)準(zhǔn)的vcf格式一樣返弹，后面緊接著的個體的基因型盡量是phasing后的。文件格式如下：< quality ><individual 1 genotype > ... < individual N genotype >

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

所以我們需要對原始的vcf文件進(jìn)行phasing

另外我們還需要準(zhǔn)備一個map文件爪飘，這個文件主要用來存放每個標(biāo)記的ID以及所在的染色體义起、物理位置、遺傳距離（cM）师崎。文件格式如下：

1 rs1 0.01 100

1 rs2 0.02 200

1 rs3 0.03 300

1 rs4 0.04 400

1 rs5 0.05 500

3.數(shù)據(jù)文件準(zhǔn)備

3.1 map文件的準(zhǔn)備

除了做模式植物的默终，大多數(shù)情況下我們都沒有一個現(xiàn)成的能和自己當(dāng)下vcf文件一一對應(yīng)的map文件。我們可以根據(jù)物種總的遺傳距離和基因組大小犁罩，算平均的遺傳距離齐蔽，而后根據(jù)每個標(biāo)記的物理位置推斷遺傳距離。

這是原始的vcf文件昼汗。

由于參考基因組版本問題肴熏，這里存在很多位于隨機(jī)或未知染色體上的變異位點(diǎn)，所以需要先過濾掉這些位點(diǎn)顷窒。

grep?-v?"random"?file.vcf?>?file1.vcf?

grep?-v?"Un"?file1.vcf?>?filtered.vcf

另外需要注意的是蛙吏，selscan軟件要求輸入的基因型只能是0或1源哩，而不能是其它。所以我們需要在這里把多等位基因位點(diǎn)過濾掉鸦做，避免后續(xù)分析出現(xiàn)問題（我當(dāng)時就吃了這個虧励烦，導(dǎo)致分析重來

Vcftools?--vcf?filtered.vcf?--min-alleles?2?--max-alleles?2?--recode?--out?Vitis_all_sample

接下來我們就從過濾后的vcf文件提取chr/id/physical position這三列

cut?-f?1,3,2?Vitis_all_sample.recode.vcf?>?dis.map

直接cut vcf文件會帶#號的注釋行一同輸出到dis.map文件中，所以需要去掉以#號開頭的注釋行

sed?-i?'/^#/d'?dis.map

接下來算遺傳距離的均值泼诱，并把物理距離轉(zhuǎn)換成遺傳距離：vitis基因組大小 ~500M坛掠，查到16年一篇文獻(xiàn)vitis總的遺傳距離為2424cM，相當(dāng)于每單位物理位置為0.00000485cM 生成最終的map文件

awk?'{?print?$1"?"$3"?"$2*0.00000485"?"$2}'?file.map?>?new.map

注意：map文件要求是空格分割

3.2 vcf文件的phasing

在phasing之前我需要根據(jù)研究目的從總的樣本的vcf文件中提取兩個亞群的子vcf文件治筒。準(zhǔn)備兩個群體分類表屉栓，每行是一個樣本名稱（注意這個樣本名稱要和vcf文件記錄的樣本號保持一致） 我使用bcftools提取vcf文件的子集

#首先對vcf文件進(jìn)行壓縮

bcftools?view?Vitis_all_sample.recode.vcf?-Oz?-o?Vitis_all_sample.recode.vcf.gz

#然后建立索引

bcftools?index?Vitis_all_sample.recode.vcf.gz

#最后根據(jù)群體分類表提取子vcf文件

bcftools?view?-S?pop.list?Vitis_all_sample.recode.vcf.gz?>?new.vcf

這里壓縮vcf文件以及建立索引這一步很重要，這樣才可以保證后面提取vcf文件不會出錯耸袜。如果你不壓縮以及建立索引友多，直接提取子集很容易出現(xiàn)no header..一類的報錯，這是因?yàn)関cf文件中的header部分并沒有指明全部的chr堤框，可以通過修改vcf文件header信息來解決這個問題域滥。但是我感覺這樣做有些麻煩，所以還是推薦先壓縮一下蜈抓，然后建立索引启绰，這樣能避免很多錯誤。

phasing軟件有很多沟使，常用的主要有SHAPEIT委可、Beagle。關(guān)于用SHAPEIT進(jìn)行phasing公眾號先前文章Haplotype phasing 常用軟件實(shí)操介紹和黃樹嘉如何使用Shapeit2對人類基因組數(shù)據(jù)進(jìn)行Phasing在這里我用Beagle軟件來做phasing格带。

Beagle軟件下載 官網(wǎng)：http://faculty.washington.edu/browning/beagle/beagle.html

右鍵復(fù)制.jar文件鏈接地址撤缴，然后轉(zhuǎn)到Linux界面下輸入

wget?-c?http://faculty.washington.edu/browning/beagle/beagle.18May20.d20.jar

pwd

beagle=path/beagle.18May20.d20.jar

因?yàn)檫@是用java語言寫的軟件，所以無須安裝叽唱，有java環(huán)境即可運(yùn)行。

常規(guī)的vcf文件就可以作為Beagle軟件的輸入微宝，此外還需要輸入一個map文件棺亭，其格式與selscan要求的map文件格式是一致的，前面已經(jīng)準(zhǔn)備過了蟋软，接下來就可以phasing了镶摘！等等，其實(shí)還有一件事沒有做.....單倍型以及基于單倍型的分析都需要逐條染色體運(yùn)行岳守，所以我們還需要把vcf文件以及map按照染色體ID進(jìn)行拆分凄敢。

拆分vcf文件

for?i?in?{1..19}

do

vcftools?--vcf?file.vcf?--chr?$i?--recode?--out?$i

done

注：葡萄有19條染色體，所以循環(huán)次數(shù)為1-19

拆分map文件

mkdir?./map#創(chuàng)建一個map文件夾

mv?new.map?./map#把整理好的map文件移動到map文件夾

awk'{print?>?$1}'new.map#這行代碼可以直接把map文件根據(jù)chr列拆分成多個文件

好了湿痢，現(xiàn)在可以開始phasing了涝缝！

Beagle軟件的運(yùn)行方式：java -Xmx[GB]g -jar beagle.jar [arguments]扑庞，需要注意的是輸入的參數(shù)和值由=連接，中間不能有任何其它字符拒逮。

for i in{1..19}

do

Java?-Xmx10g?-jar$beaglegt=${i}.vcf?map=../map/${i}ne=10?out=$inthreads=6

done

gt指輸入的vcf文件罐氨，ne是指樣品數(shù)量，nthreads就是線程數(shù)滩援，out指定輸出文件前綴栅隐。最后會生成.vcf.gz格式的文件。

4.selscan做XPEHH分析

等待兩個亞群phasing結(jié)束后玩徊，我們就可以進(jìn)行XPEHH分析了租悄，運(yùn)行代碼為：

for i in{1..19}

do

$selscan--vcf"${i}.vcf.gz"--vcf-ref"../w/${i}.vcf.gz"--map"../map/${i}"--threads?6?--out$i

done

--vcf指想研究亞群的所有染色體的vcf文件，--vcf-ref是用作參考的亞群恩袱，map就是前面已經(jīng)準(zhǔn)備好的單條染色體的map文件恰矩。根據(jù)自己數(shù)據(jù)存放的位置，要靈活調(diào)整指定的文件路徑憎蛤。

最后生成結(jié)果

XPEHH分析結(jié)果可視化

展示選擇清除分析結(jié)果的最常用的就是曼哈頓圖外傅。如下圖，橫軸代表染色體俩檬，縱軸是不同選擇清除分析指標(biāo)算的值萎胰，圖中的每一個點(diǎn)代表一個sweep或marker；如果是GWAS每一個點(diǎn)就代表一個SNP棚辽。

所以要畫這種圖技竟，我們先要把XPEHH分析所得染色體的文件合并為一個文件

cat?*.?>?black.xpehh#在存放xpehh結(jié)果文件中直接運(yùn)行這段命令

注意：因?yàn)槊織l染色體生成的xpehh文件都會有id\pos等表頭，直接cat會造成生成的總文件中存在多個表頭行屈藐，會影響整成繪圖榔组，所以我在R代碼中增加了一個對文件的過濾步驟。

繪制曼哈頓圖可以利用R包CMplot联逻、qqman搓扯。我用CMplot展示下XPEHH的分析結(jié)果

繪圖代碼：

library(tidyverse)

library(CMplot)

input?<-?“path+file”

prefix?<-?“black”

#讀入文件
xpehh <-

? read.table(input, header = T) %>%

? filter(!str_detect(id, 'id')) %>%

? separate(col = id, into = c("chr", 'position'), sep = "__", remove = F) %>%

? select(id, chr, pos, xpehh)

#作圖
png(paste(prefix, "_XPEHH.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()

在R代碼中需要注意以下幾點(diǎn)：

1、filter(!str_detect(id, 'id')) 用來過濾掉文件中多余header行包归；

2锨推、separate(col = id, into = c("chr", 'position'), sep = '__', remove = F) 用來分割id行，目的是為了生成chr列公壤；

3换可、select(id, chr, pos, xpehh) 因?yàn)?b>CMplot對輸入的作圖文件的要求：前三列必須為id\chr\pos。所以這里用select函數(shù)重新整理了以下作圖文件厦幅。

參考資料：

基于LD衰減和等位基因頻率檢測自然選擇-XPEHHhttps://mp.weixin.qq.com/s/-YydvSjsxTbRxjIPtUGLsA

Haplotype phasing 常用軟件實(shí)操介紹https://mp.weixin.qq.com/s/YwXq_MA4SEXJhNzm-lF5-g

Beagle 5.1 http://faculty.washington.edu/browning/beagle/beagle.html

selscan v1.3.0 https://github.com/szpiech/selscan/releases/tag/v1.3.0

vcftools https://vcftools.github.io/man_latest.html

bcftools http://samtools.github.io/bcftools/