繼續(xù)上一周的內(nèi)容，上次講到了如何使用--tree Dsuite的功能。要使用--treeDsuite的功能，我們顯然需要一個(gè)樹文件。一般樹文件可以通過SNP數(shù)據(jù)畫取進(jìn)化樹來獲得，這里就不詳細(xì)講解這方面內(nèi)容，留著以后再回顧。

直接給出樹文件snapp.as_newick.tre現(xiàn)在應(yīng)該具有以下內(nèi)容，其中分支長(zhǎng)度由冒號(hào)后面的數(shù)字編碼：

 ((altfas:2.594278046475669,((((((neobri:0.4315861683285048,(neooli:0.3619356529967183,neopul:0.3619356529967182):0.06965051533178657):0.12487340944557368,(neogra:0.5002021127811682,neohel:0.5002021127811681):0.05625746499291029):0.046536377715870936,(neocra:0.3843764668581544,neomar:0.3843764668581544):0.21861948863179498):0.13696107245777023,neosav:0.7399570279477196):0.3048349276159825,telvit:1.044791955563702):0.06307730474269935,(neochi:0.1368491602607037,neowal:0.1368491602607037):0.9710201000456977):1.4864087861692676):3.570689300019641,astbur:6.16496734649531);

由于這里給出的樹不包括Neolamprologus cancellatus，我們需要手動(dòng)將此物種添加到樹中。要將Neolamprologus cancellatus（“neocan”）添加到樹中，我們需要決定放置它的位置。也許最好的方法是進(jìn)行單獨(dú)的系統(tǒng)發(fā)育分析，但正如我們后面將要看到的，無論如何，Neolamprologus cancellatus（“neocan”）在物種樹中沒有單一的真實(shí)位置。因此，現(xiàn)在使用來自文件的共享派生位點(diǎn)的計(jì)數(shù)samples__BBAA.txt作為最佳放置物種的指標(biāo)就足夠了。因此，我們需要找到neocan與任何其他物種共享的最大衍生位點(diǎn)的數(shù)量。為此，我們可以使用以下命令：

cat samples__combine.txt | grep neocan | sort -n -k 4 -r | head -n 1
  cat samples__combine.txt | grep neocan | sort -n -k 5 -r | head -n 1
  cat samples__combine.txt | grep neocan | sort -n -k 6 -r | head -n 1

這應(yīng)該產(chǎn)生以下幾行：

altfas  neocan  neooli  14834.4 1408.33 4274.5
  altfas    neobri  neocan  1378.2  13479.6 4066.95
  neocan    neochi  neowal  801.883 754.94  12863.7

所有這些行中最大的數(shù)字是14834.4。由于這個(gè)數(shù)字在第四列中找到，這代表著BBAA位點(diǎn)的數(shù)量(C-BBAA)，因此指定了該行中前兩個(gè)物種共享的派生位點(diǎn)數(shù)量，“altfas”和“neocan”。這意味著“neocan”似乎與“altfas”的關(guān)系比與數(shù)據(jù)集中的其他物種關(guān)系更密切。

要將“neocan”插入樹文件作為“altfas”的姐妹物種，將“altfas”替換為“（altfas：1.0，neocan：1.0）”，包括括號(hào)。這可以使用文本編輯器完成，或者在命令行上使用以下命令完成，將修改后的樹字符串保存在名為的新文件中snapp.complete.tre：

cat snapp.as_newick.tre | sed "s/altfas/(altfas:1.0,neocan:1.0)/g" > snapp.complete.tre

現(xiàn)在樹文件應(yīng)該包含以下信息:

(((altfas:1.0,neocan:1.0):2.594278046475669,((((((neobri:0.4315861683285048,(neooli:0.3619356529967183,neopul:0.3619356529967182):0.06965051533178657):0.12487340944557368,(neogra:0.5002021127811682,neohel:0.5002021127811681):0.05625746499291029):0.046536377715870936,(neocra:0.3843764668581544,neomar:0.3843764668581544):0.21861948863179498):0.13696107245777023,neosav:0.7399570279477196):0.3048349276159825,telvit:1.044791955563702):0.06307730474269935,(neochi:0.1368491602607037,neowal:0.1368491602607037):0.9710201000456977):1.4864087861692676):3.570689300019641,astbur:6.16496734649531);

再次運(yùn)行Dsuite，這次添加-t（或--tree）選項(xiàng)以指定新準(zhǔn)備的樹文件snapp.complete.tre：

Dsuite Dtrios -t snapp.complete.tre NC_031969.f5.sub1.vcf.gz samples.txt

Dsuite應(yīng)該在幾分鐘內(nèi)再次完成分析。輸出應(yīng)該與先前寫入的輸出相同，除了多了一個(gè)名為samples__tree.txt的文件。

為了進(jìn)一步探索基于給定三個(gè)物種不同排列的規(guī)則的輸出文件之間的差異，找到每個(gè)文件中報(bào)告的最高D-統(tǒng)計(jì)數(shù)據(jù)samples__Dmin.txt），samples__BBAA.txt）和samples__tree.txt ）。一種方法是執(zhí)行以下三個(gè)命令：

cat samples__Dmin.txt | tail -n +2 | sort -n -k 4 -r | cut -f 4 | head -n 1
  cat samples__BBAA.txt | tail -n +2 | sort -n -k 4 -r | cut -f 4 | head -n 1
  cat samples__tree.txt | tail -n +2 | sort -n -k 4 -r | cut -f 4 | head -n 1

文件中報(bào)告的最高D -statistic samples__Dmin.txt（0.504355）小于其他兩個(gè)文件中的D -statistic （在兩種情況下均為0.778765）。這并不奇怪，因?yàn)槿缟纤?，D min值是給定三個(gè)物種中漸滲的保守度量。

使用以下命令找到這些極端D-statistics 的所對(duì)應(yīng)給出的三個(gè)物種：

 cat samples__Dmin.txt | tail -n +2 | sort -n -k 4 -r | head -n 1
  cat samples__BBAA.txt | tail -n +2 | sort -n -k 4 -r | head -n 1
  cat samples__tree.txt | tail -n +2 | sort -n -k 4 -r | head -n 1

你會(huì)看到最高D -statistic為0.778765，其中對(duì)應(yīng)的P1 =“altfas”，P2 =“neocan”，P3 =“telvit”。

使用以下命令從文件中查找此給定三個(gè)物種的共享站點(diǎn)的計(jì)數(shù)：

cat samples__combine.txt | grep altfas | grep neocan | grep telvit

你應(yīng)該看到C-BBAA = 12909.6（在第四列中），C-BABA = 893.302（在第五列中），并且C-ABBA = 7182.28（在第六列中）。因此，“altfas”和“neocan”共享12909.6派生位點(diǎn)，“neocan”和“telvit”共享7182.28派生位點(diǎn)，但“altfas”和“telvit”僅分享893.302派生位點(diǎn)。

因此，D -statistic =（7182.28 - 893.302）/（7182.28 + 893.302）= 0.778765 ，和上面計(jì)算的一樣。

為了更好地概述D -statistics 支持的漸滲模式，我們將以熱圖的形式繪制這些模型，其中位置P2和P3中的物種在水平軸和垂直軸上排序，以及相應(yīng)熱圖的顏色細(xì)胞表示在這兩個(gè)物種中發(fā)現(xiàn)的最顯著的D-統(tǒng)計(jì)學(xué)，在P1中的所有可能的物種中。然而，為了準(zhǔn)備這個(gè)圖，我們首先需要準(zhǔn)備一個(gè)文件，列出沿?zé)釄D軸列出P2和P3物種的順序。根據(jù)我們使用的樹文件（snapp.complete.tre），對(duì)物種進(jìn)行排序是有意義的。因此，將以下列表寫入名為的新文件species_order.txt：

  altfas
  neocan
  neochi
  neowal
  telvit
  neosav
  neocra
  neomar
  neogra
  neohel
  neobri
  neooli
  neopul

請(qǐng)注意，即使outgroup“astbur”包含在樹文件中，我們也會(huì)將其排除在外，因?yàn)樗贒suite分析中從未放置在P2或P3的位置。

要繪制熱圖，請(qǐng)使用Ruby腳本plot_d.rb，指定Dsuite的輸出文件之一，文件species_order.txt，繪制的最大D值以及輸出文件的名稱作為命令行選項(xiàng)。以samples__Dmin.txt文件作為輸入啟動(dòng)腳本，最大D為0.7，并將輸出命名為samples__Dmin.svg：

ruby plot_d.rb samples__Dmin.txt species_order.txt 0.7 samples__Dmin.svg

文件中的熱圖繪制samples__Dmin.svg是可縮放矢量圖形（SVG）格式。使用能夠以SVG格式讀取文件的程序打開此文件，例如使用Firefox或Adobe Illustrator等瀏覽器。熱圖圖如下圖所示：

正如你可以從右下角的顏色圖例中看到的那樣，這個(gè)熱圖的顏色同時(shí)顯示了兩個(gè)東西，D -statistic以及它的p值。因此，紅色表示較高的D統(tǒng)計(jì)，通常更飽和的顏色表示更重要。因此，用于漸滲的最強(qiáng)信號(hào)用飽和紅色顯示，如顏色圖例的右下角。雖然該圖中沒有一個(gè)細(xì)胞實(shí)際上具有該顏色，但“neocan”的行或列中的幾乎所有細(xì)胞都具有紅色 - 紫色，對(duì)應(yīng)于0.3左右的高度顯著的D-統(tǒng)計(jì)。

在進(jìn)一步解釋熱圖中顯示的漸滲模式之前，我們還將為其他兩個(gè)Dsuite的輸出文件生成熱圖，samples__BBAA.txt并且samples__tree.txt：

ruby plot_d.rb samples__BBAA.txt species_order.txt 0.7 samples__BBAA.svg
  ruby plot_d.rb samples__tree.txt species_order.txt 0.7 samples__tree.svg

兩個(gè)熱圖應(yīng)如下所示（samples__BBAA.svg頂部，samples__tree.svg下方）：

正如您所看到的，上述兩個(gè)熱圖總體上比基于D min的第一個(gè)熱圖更飽和，與D min作為漸滲的保守估計(jì)的解釋一致。然而，最顯著的差異在于Telmatochromis vittatus（“telvit”）所示的模式。后兩個(gè)圖中最飽和的紅色是連接“neocan”和“telvit”的細(xì)胞的顏色，相比之下，基于D min的第一個(gè)圖中只有淺藍(lán)色。因此，這些圖支持了Neolamprologus cancellatus（“neocan”）和Telmatochromis vittatus之間發(fā)生漸滲的假設(shè)。（ “telvit”）。然而，“neocan”和“telvit”的行和列中的其他紅紫色單元可能只是這種漸滲的副作用。samples__tree.svg例如，在文件的最后一個(gè)熱圖中，紅紫色細(xì)胞可能是由Neolamprologus cancellatus（“neocan”）與其他物種共享的位點(diǎn)引起的，因?yàn)檫@些其他物種與Telmatochromis vittatus（“telvit”）密切相關(guān)，是推斷的基因滲入捐贈(zèng)物種。

到目前為止，我們只計(jì)算d -statistics整個(gè)5號(hào)染色體（包括在VCF唯一的染色體），但我們還沒有測(cè)試過是否d t-統(tǒng)計(jì)是均勻或可變的整個(gè)染色體。這些信息對(duì)于確定最近的漸滲是如何發(fā)生有價(jià)值的，因?yàn)轭A(yù)期年齡漸滲事件會(huì)產(chǎn)生比舊的漸滲事件更大規(guī)模的D-統(tǒng)計(jì)變異。該信息還可以幫助識(shí)別個(gè)體獨(dú)特的滲入的基因座。

為了量化基因組中滑動(dòng)窗口的D-統(tǒng)計(jì)量，可以使用Dsuite的Dinvestigate命令，我們將它應(yīng)用于具有最強(qiáng)信號(hào)漸滲信號(hào)的三個(gè)組合，由三種物種Altolamprologus fasciatus（“altfas”），Neolamprologus cancellatus組成。（“neocan”）和Telmatochromis vittatus（“telvit”）。只需輸入以下內(nèi)容，即可查看此命令的可用選項(xiàng)：

 Dsuite Dinvestigate

需要的輸入文件有除了VCF輸入文件NC_031969.f5.sub1.vcf.gz和樣本表之外samples.txt，還需要兩條信息。一個(gè)是僅包含來自一個(gè)或多個(gè)三元組的物種名稱的文件，另一個(gè)是應(yīng)該用于滑動(dòng)窗口分析的窗口大?。由洗翱谶f增的步長(zhǎng)）。

準(zhǔn)備一個(gè)指定三物種“altfas”，“neocan”和“telvit”的文件，并命名該文件test_trios.txt。這可以使用以下命令在命令行上完成：

 echo -e "altfas\tneocan\ttelvit" > test_trios.txt

我們將使用2,500個(gè)SNP的窗口大小，每次迭代增加500個(gè)SNP。要使用這些設(shè)置啟動(dòng)滑動(dòng)窗口分析，請(qǐng)執(zhí)行以下命令：

Dsuite Dinvestigate -w 2500,500 NC_031969.f5.sub1.vcf.gz samples.txt test_trios.txt 

這分析應(yīng)該只需要幾分鐘。然后Dsuite應(yīng)該輸出一個(gè)名為的文件altfas_neocan_telvit_localFstats__2500_500.txt。

看看該文件：

  less altfas_neocan_telvit_localFstats__2500_500.txt

該文件中的行對(duì)應(yīng)于各個(gè)染色體窗口，文件中包含的六列包含有關(guān)染色體ID以及每個(gè)窗口的起始和結(jié)束位置的信息。接下來是每行上的三個(gè)數(shù)字，它們代表第四列中窗口的D-統(tǒng)計(jì)量，以及另外兩個(gè)統(tǒng)計(jì)數(shù)據(jù)（Fd，F(xiàn)dm），旨在量化受漸滲影響的窗口比例。

使用以下命令查找染色體5被劃分的窗口數(shù)：

cat altfas_neocan_telvit_localFstats__2500_500.txt | tail -n +2 | wc -l

這應(yīng)該表明在Dsuite的Dinvestigate分析中使用了509個(gè)窗口，另外該窗口的數(shù)量可以通過指定更大或更小的窗口大小來調(diào)整信息密度-w

為了可視化5號(hào)染色體上D和f d統(tǒng)計(jì)數(shù)據(jù)的變化，使用R作圖：

  table <- read.table("altfas_neocan_telvit_localFstats__2500_500.txt", header=T)
  windowCenter=(table$windowStart+table$windowEnd)/2
  pdf("altfas_neocan_telvit_localFstats__2500_500.pdf", height=7, width=7)
  plot(windowCenter, table$D, type="l", xlab="Position", ylab="D (gray) / fD (black)", ylim=c(0,1), main="Chromosome 5 (altfas, neocan, telvit)", col="gray")
  lines(windowCenter, table$f_d)
  dev.off()

作圖會(huì)生成名為altfas_neocan_telvit_localFstats__2500_500.pdf的文件。

D -statistic（灰色）幾乎在所有窗口中顯著高于f d -statistic（黑色），并且f d -statistic，其估計(jì)混合比例大多接近0.5：并在大約%Mbp的地方有下降的信號(hào)，這是否是由參考基因組中缺失的數(shù)據(jù)或錯(cuò)誤組裝導(dǎo)致的假陽性結(jié)果，或者它是否顯示出減少漸滲的生物信號(hào)，如果沒有進(jìn)一步的分析，很難說清楚。

重復(fù)這個(gè)分析，但是現(xiàn)在使用不同的三個(gè)物種的組合,“neooli”），N. pulcher（“neopul），和N. brichardi（” neobri“）：

 echo -e "neooli\tneopul\tneobri" > test_trios.txt
  Dsuite Dinvestigate -w 2500,500 NC_031969.f5.sub1.vcf.gz samples.txt test_trios.txt

再次進(jìn)行畫圖：

  table <- read.table("neooli_neopul_neobri_localFstats__2500_500.txt", header=T)
  windowCenter=(table$windowStart+table$windowEnd)/2
  pdf("neooli_neopul_neobri_localFstats__2500_500.pdf", height=7, width=7)
  plot(windowCenter, table$D, type="l", xlab="Position", ylab="D (gray) / fD (black)", ylim=c(0,1), main="Chromosome 5 (neooli, neopul, neobri)", col="gray")
  lines(windowCenter, table$f_d)
  dev.off()

產(chǎn)生該圖：

正如你所看到的，D -statistic（黑色）現(xiàn)在顯示更多的變化，甚至在幾個(gè)窗口中變?yōu)樨?fù)數(shù)，表明在這些窗口中，neooli與neobri的相似性高于neopul?？傮w而言，兩項(xiàng)統(tǒng)計(jì)數(shù)據(jù)都低于第一個(gè)三個(gè)物種的組合。值得注意的是，這里在5Mbp的地方，再次可以看到明顯的下降，由于同一地區(qū)不太可能在兩個(gè)物種不同三組合中具有更多的漸滲，這可能表明這實(shí)際上是由技術(shù)的缺陷導(dǎo)致的而非生物過程引起的。