距離上次更新這個專題后披诗，該教程的原作者更新了一下原教程甲献，使用了一個新的宰缤，非常好用的工具，Dsuite進(jìn)行基因滲入分析，這期的推文就和大家一起回顧并學(xué)習(xí)一下該工具的使用撵溃，加深對基因滲入分析的了解疚鲤。
首先如果你是第一次或者沒有看到之前的文章锥累，先去看看滲入分析的背景介紹和該教程所用的數(shù)據(jù)的相關(guān)背景知識：
然后還有就是舊版的教程缘挑，大家有興趣可以看看，比較一下使用的方法的異同：
準(zhǔn)備工作
測試數(shù)據(jù)的下載
wget https://raw.githubusercontent.com/mmatschiner/tutorials/master/analysis_of_introgression_with_snp_data/data/NC_031969.f5.sub1.vcf.gz
保存好我們的樣本名文件桶略，這里除了outgroup物種外有13個物種：
vi samples.txt
###將下面的樣本名復(fù)制進(jìn)去
IZA1 Outgroup
IZC5 Outgroup
AUE7 altfas
AXD5 altfas
JBD5 telvit
JBD6 telvit
JUH9 neobri
JUI1 neobri
LJC9 neocan
LJD1 neocan
KHA7 neochi
KHA9 neochi
IVE8 neocra
IVF1 neocra
JWH1 neogra
JWH2 neogra
JWG8 neohel
JWG9 neohel
JWH3 neomar
JWH4 neomar
JWH5 neooli
JWH6 neooli
ISA6 neopul
ISB3 neopul
ISA8 neosav
IYA4 neosav
KFD2 neowal
KFD4 neowal
對Dsuite進(jìn)行安裝：
git clone https://github.com/millanek/Dsuite.git
cd Dsuite
make
Dsuite主要包含三個：“Dtrios”语淘，“DtriosCombine”和“Dinvestigate”不同命令，這里我們計算D統(tǒng)計量主要用到 Dtrios际歼，下面是該工具的幫助文檔惶翻，大家可以熟悉一下，主要的參數(shù)就是-j,主要用于設(shè)定計算窗口的大小范圍：
Usage: Dsuite Dtrios [OPTIONS] INPUT_FILE.vcf SETS.txt
Calculate the Dmin-statistic - the ABBA/BABA stat for all trios of species in the dataset (the outgroup being fixed)
the calculation is as definded in Durand et al. 2011
The SETS.txt should have two columns: SAMPLE_ID SPECIES_ID
The outgroup (can be multiple samples) should be specified by using the keywork Outgroup in place of the SPECIES_ID
-h, --help display this help and exit
-j, --JKwindow (default=20000) Jackknife block size in SNPs
-r , --region=start,length (optional) only process a subset of the VCF file
-t , --tree=TREE_FILE.nwk (optional) a file with a tree in the newick format specifying the relationships between populations/species
D values for trios arranged according to these relationships will be output in a file with _tree.txt suffix
-n, --run-name run-name will be included in the output file name
實戰(zhàn)演示
如Dsuite幫助文本上面所示鹅心，該程序可以運(yùn)行Dsuite Dtrios [OPTIONS] INPUT_FILE.vcf SETS.txt吕粗。在我們的例子中，輸入文件是NC_031969.f5.sub1.vcf.gz旭愧，而sets文件是包含之前寫入的樣本和物種ID的表samples.txt颅筋。因此，啟動我們的測試數(shù)據(jù)集的Dsuite分析：
Dsuite Dtrios NC_031969.f5.sub1.vcf.gz samples.txt
該運(yùn)行應(yīng)該在幾分鐘內(nèi)就能完成输枯。
現(xiàn)在使用該ls命令查看當(dāng)前目錄中的文件议泵。應(yīng)該能看到Dsuite已經(jīng)編寫了一些以sample-table文件命名的文件：
samples__BBAA.txt
samples__Dmin.txt
samples__combine.txt
samples__combine_stderr.txt
接下來，查看文件的內(nèi)容samples__combine.txt桃熄，例如使用less命令先口。您將看到此文件的第一部分具有以下內(nèi)容
altfas neobri neocan 1378.2 13479.6 4066.95
altfas neobri neochi 2281.67 2322.45 8154.16
altfas neobri neocra 1610.98 1495.17 14174.6
altfas neobri neogra 1618.4 1559.8 14312.1
altfas neobri neohel 1549.77 1416.71 14741.8
altfas neobri neomar 1980.94 1728.7 13168.8
altfas neobri neooli 1668.03 1460.84 14734.4
altfas neobri neopul 1279.66 1164.73 14283.7
altfas neobri neosav 1729.49 1616.15 12626.2
altfas neobri neowal 2525.99 2440.28 8572.34
altfas neobri telvit 2581.12 2688.36 8258.12
altfas neocan neochi 12306.4 1310.06 3745.6
altfas neocan neocra 14200.7 1395.33 4097.46
altfas neocan neogra 14825.3 1436.9 4286.4
altfas neocan neohel 14063.7 1361.83 4102.08
altfas neocan neomar 14794.6 1422.08 4181.25
altfas neocan neooli 14834.4 1408.33 4274.5
...
在這里，每一行顯示分析一個三個物種的分析結(jié)果瞳收，例如在第一行中碉京，Altolamprologus fasciatus（“altfas”）用作P1，Neolamprologus brichardi被認(rèn)為是P2螟深，而Neolamprologus cancellatus被放置在P3谐宙。然后該行的第五和第六列中的數(shù)字分別代表著:ABBA位點(diǎn)在該三個物種中的數(shù)量(C-ABBA)（其中衍生的等位基因由“neobri”和“neocan”共享）和BABA位點(diǎn)的計數(shù)，C-BABA（衍生的等位基因由“altfas”和“neocan”共享）血崭。除了第5和第6列中BABA和ABBA位點(diǎn)的數(shù)量之外卧惜，第4列列出了“BBAA”位點(diǎn)的數(shù)量(C-BBAA)，P1和P2共享衍生的等位基因（因此通過“altfas”和“ neobri“共享）夹纫。
你可能會注意到咽瓷，在此文件中，所有三元組都按字母順序排序; 因此舰讹，P1在P2之前按字母順序排列茅姜，P2在P3之前。然而，如上所述钻洒，ABBA-BABA測試基于P1和P2是姊妹物種的假設(shè)奋姿。當(dāng)計算給定三個物種的D-統(tǒng)計量時，Dsuite首先重新排列分配給P1素标，P2和P3的物種（因此ABBA称诗，BABA和BBAA位點(diǎn)的數(shù)量也重新排列），這是根據(jù)某些規(guī)則：
- 在第一組計算中头遭，測試所有可能的重新排列寓免，并且報告每給定三個物種的最低D-統(tǒng)計量，稱為D min计维。因此袜香，D min是給定三個物種中D-統(tǒng)計量的保守估計。此輸出將寫入具有__Dmin.txt結(jié)尾的文件鲫惶。
- 這樣的排列使得P1和P2是給定三個物種的兩種物種蜈首，它們共享最大數(shù)量的衍生地點(diǎn)。換句話說欠母，進(jìn)行重排以使C-BBAA > C-ABBA和C-BBAA > C-BABA欢策。另外，旋轉(zhuǎn)P1和P2艺蝴，使得在P2和P3之間共享的派生位點(diǎn)的數(shù)量大于P1和P3之間的派生站點(diǎn)的數(shù)量猬腰。總之猜敢，這意味著C-BBAA > C-ABBA >C-BABA姑荷。這種類型的重排背后的想法是，共享最大數(shù)量的衍生地點(diǎn)的兩個物種很可能是給定三個物種中真正的兩個姊妹物種缩擂，因此正確地放置在位置P1和P2鼠冕。預(yù)計這種假設(shè)在某些條件下會持續(xù)存在（例如時鐘式演化，缺乏同質(zhì)性胯盯，缺乏漸滲懈费，以及泛化的祖先種群），但有時候很難說真實數(shù)據(jù)集的可靠性博脑≡饕遥基于這種類型的重新排列的D-統(tǒng)計被寫入具有結(jié)尾__BBAA.txt的文件。
- 要直接告訴Dsuite哪個給定三個物種應(yīng)該被視為姊妹物種（因此叉趣，應(yīng)該將其分配給P1和P2）泞边，我們可以使用參數(shù)
-t或提供包含所有數(shù)據(jù)集種類的樹文件--tree。然后輸出將被寫入帶有結(jié)尾__tree.txt的附加文件疗杉。這里后面會講到阵谚。
那么，讓我們看看如何基于給定三個物種的前兩條規(guī)則計算D統(tǒng)計量∩沂玻看看文件samples__Dmin.txt奠蹬，例如使用less命令：
less samples__Dmin.txt
你應(yīng)該看到文件的第一部分，如下所示：
P1 P2 P3 Dstatistic p-value
altfas neocan neobri 0.493787 0
neobri neochi altfas 0.00885755 0.3836
neocra neobri altfas 0.0372849 0.013765
neogra neobri altfas 0.0184394 0.258714
neohel neobri altfas 0.0448571 0.113967
neomar neobri altfas 0.0679957 0.0195241
neooli neobri altfas 0.0662179 0.0133793
neopul neobri altfas 0.0470166 0.10957
neosav neobri altfas 0.0338796 0.103889
neowal neobri altfas 0.0172581 0.266909
neobri telvit altfas 0.0203511 0.199163
altfas neocan neochi 0.481745 0
altfas neocan neocra 0.491942 0
altfas neocan neogra 0.497877 0
altfas neocan neohel 0.501518 0
altfas neocan neomar 0.492416 0
altfas neocan neooli 0.504355 0
...
如標(biāo)題行所示嗡午，第四列現(xiàn)在顯示每一個給定三個物種的的D-統(tǒng)計量囤躁，第五列顯示基于對D = 0 的零假設(shè)的歸一化的p值。
讓我們選擇一個給定三個物種的例子翼馆，看看它是如何出現(xiàn)在三個不同的文件samples__combine.txt中samples__Dmin.txt割以，和samples__BBAA.txt金度。我們將選擇第一個給定三個物種的例子应媚，包括“altfas”，“neocan”和“neobri”的那一行猜极。要僅從三個文件中查看此給定三個物種的例子的行中姜，可以使用此命令：
cat samples__combine.txt | grep altfas | grep neocan | grep neobri
cat samples__Dmin.txt | grep altfas | grep neocan | grep neobri
cat samples__BBAA.txt | grep altfas | grep neocan | grep neobri
這應(yīng)該會分別輸出以下三行：
###samples__combine.txt
altfas neobri neocan 1378.2 13479.6 4066.95
###samples__Dmin.txt
altfas neocan neobri 0.493787 0
###samples__BBAA.txt
altfas neocan neobri 0.493787 0
這里簡單解析一下結(jié)果，首先samples__combine.txt品種名字母排序的方式與其它兩個文件有點(diǎn)不同跟伏，P1在三個文件中保持相同（“altfas”）丢胚，但是P2和P3的順序被交換了（“neocan”和“neobri”）。此交換還暗示ABBA受扳，BABA和BBAA模式的計數(shù)相應(yīng)地交換携龟。因此在交換之后并且P1 =“altfas”，P2 =“neocan”勘高，P3 =“neobri”峡蟋，計數(shù)如下：C-ABBA = 4066.95，C-BABA = 1378.2华望，C-BBAA = 13479.6蕊蝗。因此，“neocan”和“neobri”共享4066.95個衍生位點(diǎn)赖舟，“altfas”和“neobri”共享1378.2個衍生位站蓬戚，“altfas”和“neocan”共享13479.6個衍生位站。有了這些數(shù)量宾抓，D=（4066.95 - 1378.2）/（4066.95 + 1378.2）= 0.493787子漩。這個數(shù)字與Dsuite在這兩個文件（samples__Dmin.txt和samples__BBAA.txt。）生成的報告一致石洗。
重復(fù)與上一步相同的操作幢泼，但這一次使用給定開頭為neo的三個物種的，文件samples__Dmin.txt和重新排列不同samples__BBAA.txt劲腿。三個“neobri”旭绒，“neocra”和“neogra”就是這樣一個例子。使用這些命令可以在所有三個文件中查看此給定開頭為neo的行：
cat samples__combine.txt | grep neobri | grep neocra | grep neogra
cat samples__Dmin.txt | grep neobri | grep neocra | grep neogra
cat samples__BBAA.txt | grep neobri | grep neocra | grep neogra
結(jié)果會分別輸出以下的行：
###samples__combine.txt
neobri neocra neogra 3788.23 3552.38 2992.93
###samples__Dmin.txt
neogra neocra neobri 0.0321294 0.145298
###samples__BBAA.txt
neocra neobri neogra 0.0854723 8.58201e-07
文件中的結(jié)果samples__BBAA.txt表明，當(dāng)P1 =“neocra”挥吵，P2 =“neobri”重父，P3 =“neogra”時，則C-BBAA = 3788.23忽匈，C-ABBA = 3552.38房午，C-BABA = 2992.93，因此D =（3552.38 - 2992.93）/（3552.38 + 2992.93）= 0.0854723丹允。
然而郭厌，文件中的結(jié)果samples__Dmin.txt顯示，這次發(fā)生了另一次重新排列（因此C-BBAA不大于其他兩個計數(shù)的重新排列）產(chǎn)生較低的D-統(tǒng)計：P1 =“neogra”雕蔽，P2 =“neocra” 折柠，并且P3 =“neobri”，則C-BBAA = 2992.93批狐，C-ABBA = 3788.23扇售，并且C-BABA = 3552.38，因此D =（3788.23-3552.38）/（3788.23 + 3552.38）= 0.0321294嚣艇。
這說明了D-min值承冰，報告了可能的最低D值對于給定三物種的統(tǒng)計，有時選擇該三個物種的重新排列食零，其中P1和P2實際上彼此共享較少的派生位點(diǎn)困乒，而不是它們兩者與P3共享。這與ABBA-BABA測試的原始假設(shè)相沖突贰谣，即P1和P2彼此之間的關(guān)系比P3更緊密娜搂。在解釋Dsuite分析的結(jié)果時，應(yīng)該記住冈爹，在文件中報告的D min值__Dmin.txt實際上是D -statistic的保守估計涌攻。文件中報告的值以__BBAA.txt結(jié)尾為基礎(chǔ)，這些值基于確保C-BBAA > C-ABBA > C-BABA的重新排列频伤，通晨一眩可以更好地測量D -statistic，但是憋肖，最好的選擇可能是使用--tree選項運(yùn)行分析因痛，提供一個輸入樹，直接告訴Dsuite如何重新排列所有三個物種岸更，這部分就留下次再講吧鸵膏，要不然信息量太大大家也不好吸收。
