生信學(xué)習(xí)筆記：使用SNP data做基因滲入分析 (4)

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

直接給出樹文件snapp.as_newick.tre現(xiàn)在應(yīng)該具有以下內(nèi)容八匠，其中分支長度由冒號后面的數(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絮爷，我們需要手動將此物種添加到樹中。要將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。由于這個數(shù)字在第四列中找到勉耀，這代表著BBAA位點(diǎn)的數(shù)量(C-BBAA)指煎，因此指定了該行中前兩個物種共享的派生位點(diǎn)數(shù)量，“altfas”和“neocan”便斥。這意味著“neocan”似乎與“altfas”的關(guān)系比與數(shù)據(jù)集中的其他物種關(guān)系更密切至壤。

要將“neocan”插入樹文件作為“altfas”的姐妹物種，將“altfas”替換為“（altfas：1.0椭住，neocan：1.0）”崇渗，包括括號。這可以使用文本編輯器完成京郑，或者在命令行上使用以下命令完成宅广，將修改后的樹字符串保存在名為的新文件中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)該與先前寫入的輸出相同跟狱，除了多了一個名為samples__tree.txt的文件。

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

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）小于其他兩個文件中的D -statistic （在兩種情況下均為0.778765）叼丑。這并不奇怪关翎，因?yàn)槿缟纤觯珼 min值是給定三個物種中漸滲的保守度量鸠信。

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

 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

你會看到最高D -statistic為0.778765纵寝，其中對應(yīng)的P1 =“altfas”，P2 =“neocan”星立，P3 =“telvit”爽茴。

使用以下命令從文件中查找此給定三個物種的共享站點(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ì)胞表示在這兩個物種中發(fā)現(xiàn)的最顯著的D-統(tǒng)計(jì)學(xué)，在P1中的所有可能的物種中碎紊。然而佑附，為了準(zhǔn)備這個圖，我們首先需要準(zhǔn)備一個文件仗考，列出沿?zé)釄D軸列出P2和P3物種的順序音同。根據(jù)我們使用的樹文件（snapp.complete.tre），對物種進(jìn)行排序是有意義的秃嗜。因此权均，將以下列表寫入名為的新文件species_order.txt：

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

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

要繪制熱圖必搞，請使用Ruby腳本plot_d.rb必指，指定Dsuite的輸出文件之一，文件species_order.txt恕洲，繪制的最大D值以及輸出文件的名稱作為命令行選項(xiàng)塔橡。以samples__Dmin.txt文件作為輸入啟動腳本，最大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等瀏覽器泌类。熱圖圖如下圖所示：

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

在進(jìn)一步解釋熱圖中顯示的漸滲模式之前让簿，我們還將為其他兩個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

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

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

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

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

 Dsuite Dinvestigate

需要的輸入文件有除了VCF輸入文件NC_031969.f5.sub1.vcf.gz和樣本表之外samples.txt鞠评，還需要兩條信息。一個是僅包含來自一個或多個三元組的物種名稱的文件壕鹉，另一個是應(yīng)該用于滑動窗口分析的窗口大刑昊稀（加上窗口遞增的步長）。

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

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

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

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

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

看看該文件：

  less altfas_neocan_telvit_localFstats__2500_500.txt

該文件中的行對應(yīng)于各個染色體窗口，文件中包含的六列包含有關(guān)染色體ID以及每個窗口的起始和結(jié)束位置的信息杈抢。接下來是每行上的三個數(shù)字数尿，它們代表第四列中窗口的D-統(tǒng)計(jì)量，以及另外兩個統(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個窗口歼捐，另外該窗口的數(shù)量可以通過指定更大或更小的窗口大小來調(diào)整信息密度-w

為了可視化5號染色體上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()

作圖會生成名為altfas_neocan_telvit_localFstats__2500_500.pdf的文件。

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

重復(fù)這個分析佳魔，但是現(xiàn)在使用不同的三個物種的組合,“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)在顯示更多的變化宁脊，甚至在幾個窗口中變?yōu)樨?fù)數(shù)，表明在這些窗口中贤姆，neooli與neobri的相似性高于neopul榆苞。總體而言霞捡，兩項(xiàng)統(tǒng)計(jì)數(shù)據(jù)都低于第一個三個物種的組合语稠。值得注意的是，這里在5Mbp的地方弄砍，再次可以看到明顯的下降，由于同一地區(qū)不太可能在兩個物種不同三組合中具有更多的漸滲输涕，這可能表明這實(shí)際上是由技術(shù)的缺陷導(dǎo)致的而非生物過程引起的音婶。