距離上次更新這個(gè)專題后，該教程的原作者更新了一下原教程，使用了一個(gè)新的，非常好用的工具，Dsuite進(jìn)行基因滲入分析，這期的推文就和大家一起回顧并學(xué)習(xí)一下該工具的使用，加深對(duì)基因滲入分析的了解。

首先如果你是第一次或者沒有看到之前的文章，先去看看滲入分析的背景介紹和該教程所用的數(shù)據(jù)的相關(guān)背景知識(shí)：

然后還有就是舊版的教程，大家有興趣可以看看，比較一下使用的方法的異同：

準(zhǔn)備工作

測(cè)試數(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個(gè)物種：

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

對(duì)Dsuite進(jìn)行安裝：

git clone https://github.com/millanek/Dsuite.git
cd Dsuite
make

Dsuite主要包含三個(gè)：“Dtrios”，“DtriosCombine”和“Dinvestigate”不同命令，這里我們計(jì)算D統(tǒng)計(jì)量主要用到 Dtrios，下面是該工具的幫助文檔，大家可以熟悉一下，主要的參數(shù)就是-j,主要用于設(shè)定計(jì)算窗口的大小范圍：

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

實(shí)戰(zhàn)演示

如Dsuite幫助文本上面所示，該程序可以運(yùn)行Dsuite Dtrios [OPTIONS] INPUT_FILE.vcf SETS.txt。在我們的例子中，輸入文件是NC_031969.f5.sub1.vcf.gz，而sets文件是包含之前寫入的樣本和物種ID的表samples.txt。因此，啟動(dòng)我們的測(cè)試數(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
...

在這里，每一行顯示分析一個(gè)三個(gè)物種的分析結(jié)果，例如在第一行中，Altolamprologus fasciatus（“altfas”）用作P1，Neolamprologus brichardi被認(rèn)為是P2，而Neolamprologus cancellatus被放置在P3。然后該行的第五和第六列中的數(shù)字分別代表著:ABBA位點(diǎn)在該三個(gè)物種中的數(shù)量(C-ABBA)（其中衍生的等位基因由“neobri”和“neocan”共享）和BABA位點(diǎn)的計(jì)數(shù)，C-BABA（衍生的等位基因由“altfas”和“neocan”共享）。除了第5和第6列中BABA和ABBA位點(diǎn)的數(shù)量之外，第4列列出了“BBAA”位點(diǎn)的數(shù)量(C-BBAA)，P1和P2共享衍生的等位基因（因此通過“altfas”和“ neobri“共享）。

你可能會(huì)注意到，在此文件中，所有三元組都按字母順序排序; 因此，P1在P2之前按字母順序排列，P2在P3之前。然而，如上所述，ABBA-BABA測(cè)試基于P1和P2是姊妹物種的假設(shè)。當(dāng)計(jì)算給定三個(gè)物種的D-統(tǒng)計(jì)量時(shí)，Dsuite首先重新排列分配給P1，P2和P3的物種（因此ABBA，BABA和BBAA位點(diǎn)的數(shù)量也重新排列），這是根據(jù)某些規(guī)則：

• 在第一組計(jì)算中，測(cè)試所有可能的重新排列，并且報(bào)告每給定三個(gè)物種的最低D-統(tǒng)計(jì)量，稱為D min。因此，D min是給定三個(gè)物種中D-統(tǒng)計(jì)量的保守估計(jì)。此輸出將寫入具有__Dmin.txt結(jié)尾的文件。
• 這樣的排列使得P1和P2是給定三個(gè)物種的兩種物種，它們共享最大數(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)的兩個(gè)物種很可能是給定三個(gè)物種中真正的兩個(gè)姊妹物種，因此正確地放置在位置P1和P2。預(yù)計(jì)這種假設(shè)在某些條件下會(huì)持續(xù)存在（例如時(shí)鐘式演化，缺乏同質(zhì)性，缺乏漸滲，以及泛化的祖先種群），但有時(shí)候很難說真實(shí)數(shù)據(jù)集的可靠性?；谶@種類型的重新排列的D-統(tǒng)計(jì)被寫入具有結(jié)尾__BBAA.txt的文件。
• 要直接告訴Dsuite哪個(gè)給定三個(gè)物種應(yīng)該被視為姊妹物種（因此，應(yīng)該將其分配給P1和P2），我們可以使用參數(shù)-t或提供包含所有數(shù)據(jù)集種類的樹文件--tree。然后輸出將被寫入帶有結(jié)尾__tree.txt的附加文件。這里后面會(huì)講到。

那么，讓我們看看如何基于給定三個(gè)物種的前兩條規(guī)則計(jì)算D統(tǒng)計(jì)量?？纯次募amples__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)在顯示每一個(gè)給定三個(gè)物種的的D-統(tǒng)計(jì)量，第五列顯示基于對(duì)D = 0 的零假設(shè)的歸一化的p值。

讓我們選擇一個(gè)給定三個(gè)物種的例子，看看它是如何出現(xiàn)在三個(gè)不同的文件samples__combine.txt中samples__Dmin.txt，和samples__BBAA.txt。我們將選擇第一個(gè)給定三個(gè)物種的例子，包括“altfas”，“neocan”和“neobri”的那一行。要僅從三個(gè)文件中查看此給定三個(gè)物種的例子的行，可以使用此命令：

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)該會(huì)分別輸出以下三行：

###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ǎn)單解析一下結(jié)果，首先samples__combine.txt品種名字母排序的方式與其它兩個(gè)文件有點(diǎn)不同，P1在三個(gè)文件中保持相同（“altfas”），但是P2和P3的順序被交換了（“neocan”和“neobri”）。此交換還暗示ABBA，BABA和BBAA模式的計(jì)數(shù)相應(yīng)地交換。因此在交換之后并且P1 =“altfas”，P2 =“neocan”，P3 =“neobri”，計(jì)數(shù)如下：C-ABBA = 4066.95，C-BABA = 1378.2，C-BBAA = 13479.6。因此，“neocan”和“neobri”共享4066.95個(gè)衍生位點(diǎn)，“altfas”和“neobri”共享1378.2個(gè)衍生位站，“altfas”和“neocan”共享13479.6個(gè)衍生位站。有了這些數(shù)量，D=（4066.95 - 1378.2）/（4066.95 + 1378.2）= 0.493787。這個(gè)數(shù)字與Dsuite在這兩個(gè)文件（samples__Dmin.txt和samples__BBAA.txt。）生成的報(bào)告一致。

重復(fù)與上一步相同的操作，但這一次使用給定開頭為neo的三個(gè)物種的，文件samples__Dmin.txt和重新排列不同samples__BBAA.txt。三個(gè)“neobri”，“neocra”和“neogra”就是這樣一個(gè)例子。使用這些命令可以在所有三個(gè)文件中查看此給定開頭為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é)果會(huì)分別輸出以下的行：

###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”時(shí)，則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不大于其他兩個(gè)計(jì)數(shù)的重新排列）產(chǎn)生較低的D-統(tǒng)計(jì)：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值，報(bào)告了可能的最低D值對(duì)于給定三物種的統(tǒng)計(jì)，有時(shí)選擇該三個(gè)物種的重新排列，其中P1和P2實(shí)際上彼此共享較少的派生位點(diǎn)，而不是它們兩者與P3共享。這與ABBA-BABA測(cè)試的原始假設(shè)相沖突，即P1和P2彼此之間的關(guān)系比P3更緊密。在解釋Dsuite分析的結(jié)果時(shí)，應(yīng)該記住，在文件中報(bào)告的D min值__Dmin.txt實(shí)際上是D -statistic的保守估計(jì)。文件中報(bào)告的值以__BBAA.txt結(jié)尾為基礎(chǔ)，這些值基于確保C-BBAA > C-ABBA > C-BABA的重新排列，通?？梢愿玫販y(cè)量D -statistic，但是，最好的選擇可能是使用--tree選項(xiàng)運(yùn)行分析，提供一個(gè)輸入樹，直接告訴Dsuite如何重新排列所有三個(gè)物種，這部分就留下次再講吧，要不然信息量太大大家也不好吸收。