生信学习笔记：使用SNP data做基因渗入分析 (4)

继续上一周的内容，上次讲到了如何使用--tree Dsuite的功能。要使用--treeDsuite的功能，我们显然需要一个树文件。一般树文件可以通过SNP数据画取进化树来获得，这里就不详细讲解这方面内容，留着以后再回顾。

直接给出树文件snapp.as_newick.tre现在应该具有以下内容，其中分支长度由冒号后面的数字编码：

 ((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”）添加到树中，我们需要决定放置它的位置。也许最好的方法是进行单独的系统发育分析，但正如我们后面将要看到的，无论如何，Neolamprologus cancellatus（“neocan”）在物种树中没有单一的真实位置。因此，现在使用来自文件的共享派生位点的计数samples__BBAA.txt作为最佳放置物种的指标就足够了。因此，我们需要找到neocan与任何其他物种共享的最大衍生位点的数量。为此，我们可以使用以下命令：

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

这应该产生以下几行：

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

所有这些行中最大的数字是14834.4。由于这个数字在第四列中找到，这代表着BBAA位点的数量(C-BBAA)，因此指定了该行中前两个物种共享的派生位点数量，“altfas”和“neocan”。这意味着“neocan”似乎与“altfas”的关系比与数据集中的其他物种关系更密切。

要将“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

现在树文件应该包含以下信息:

(((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);

再次运行Dsuite，这次添加-t（或--tree）选项以指定新准备的树文件snapp.complete.tre：

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

Dsuite应该在几分钟内再次完成分析。输出应该与先前写入的输出相同，除了多了一个名为samples__tree.txt的文件。

为了进一步探索基于给定三个物种不同排列的规则的输出文件之间的差异，找到每个文件中报告的最高D-统计数据samples__Dmin.txt），samples__BBAA.txt）和samples__tree.txt ）。一种方法是执行以下三个命令：

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

文件中报告的最高D -statistic samples__Dmin.txt（0.504355）小于其他两个文件中的D -statistic （在两种情况下均为0.778765）。这并不奇怪，因为如上所述，D min值是给定三个物种中渐渗的保守度量。

使用以下命令找到这些极端D-statistics 的所对应给出的三个物种：

 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，其中对应的P1 =“altfas”，P2 =“neocan”，P3 =“telvit”。

使用以下命令从文件中查找此给定三个物种的共享站点的计数：

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

你应该看到C-BBAA = 12909.6（在第四列中），C-BABA = 893.302（在第五列中），并且C-ABBA = 7182.28（在第六列中）。因此，“altfas”和“neocan”共享12909.6派生位点，“neocan”和“telvit”共享7182.28派生位点，但“altfas”和“telvit”仅分享893.302派生位点。

因此，D -statistic =（7182.28 - 893.302）/（7182.28 + 893.302）= 0.778765 ，和上面计算的一样。

为了更好地概述D -statistics 支持的渐渗模式，我们将以热图的形式绘制这些模型，其中位置P2和P3中的物种在水平轴和垂直轴上排序，以及相应热图的颜色细胞表示在这两个物种中发现的最显着的D-统计学，在P1中的所有可能的物种中。然而，为了准备这个图，我们首先需要准备一个文件，列出沿热图轴列出P2和P3物种的顺序。根据我们使用的树文件（snapp.complete.tre），对物种进行排序是有意义的。因此，将以下列表写入名为的新文件species_order.txt：

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

请注意，即使outgroup“astbur”包含在树文件中，我们也会将其排除在外，因为它在Dsuite分析中从未放置在P2或P3的位置。

要绘制热图，请使用Ruby脚本plot_d.rb，指定Dsuite的输出文件之一，文件species_order.txt，绘制的最大D值以及输出文件的名称作为命令行选项。以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统计，通常更饱和的颜色表示更重要。因此，用于渐渗的最强信号用饱和红色显示，如颜色图例的右下角。虽然该图中没有一个细胞实际上具有该颜色，但“neocan”的行或列中的几乎所有细胞都具有红色 - 紫色，对应于0.3左右的高度显着的D-统计。

在进一步解释热图中显示的渐渗模式之前，我们还将为其他两个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

两个热图应如下所示（samples__BBAA.svg顶部，samples__tree.svg下方）：

正如您所看到的，上述两个热图总体上比基于D min的第一个热图更饱和，与D min作为渐渗的保守估计的解释一致。然而，最显着的差异在于Telmatochromis vittatus（“telvit”）所示的模式。后两个图中最饱和的红色是连接“neocan”和“telvit”的细胞的颜色，相比之下，基于D min的第一个图中只有浅蓝色。因此，这些图支持了Neolamprologus cancellatus（“neocan”）和Telmatochromis vittatus之间发生渐渗的假设。（ “telvit”）。然而，“neocan”和“telvit”的行和列中的其他红紫色单元可能只是这种渐渗的副作用。samples__tree.svg例如，在文件的最后一个热图中，红紫色细胞可能是由Neolamprologus cancellatus（“neocan”）与其他物种共享的位点引起的，因为这些其他物种与Telmatochromis vittatus（“telvit”）密切相关，是推断的基因渗入捐赠物种。

到目前为止，我们只计算d -statistics整个5号染色体（包括在VCF唯一的染色体），但我们还没有测试过是否d t-统计是均匀或可变的整个染色体。这些信息对于确定最近的渐渗是如何发生有价值的，因为预期年龄渐渗事件会产生比旧的渐渗事件更大规模的D-统计变异。该信息还可以帮助识别个体独特的渗入的基因座。

为了量化基因组中滑动窗口的D-统计量，可以使用Dsuite的Dinvestigate命令，我们将它应用于具有最强信号渐渗信号的三个组合，由三种物种Altolamprologus fasciatus（“altfas”），Neolamprologus cancellatus组成。（“neocan”）和Telmatochromis vittatus（“telvit”）。只需输入以下内容，即可查看此命令的可用选项：

 Dsuite Dinvestigate

需要的输入文件有除了VCF输入文件NC_031969.f5.sub1.vcf.gz和样本表之外samples.txt，还需要两条信息。一个是仅包含来自一个或多个三元组的物种名称的文件，另一个是应该用于滑动窗口分析的窗口大小（加上窗口递增的步长）。

准备一个指定三物种“altfas”，“neocan”和“telvit”的文件，并命名该文件test_trios.txt。这可以使用以下命令在命令行上完成：

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

我们将使用2,500个SNP的窗口大小，每次迭代增加500个SNP。要使用这些设置启动滑动窗口分析，请执行以下命令：

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

这分析应该只需要几分钟。然后Dsuite应该输出一个名为的文件altfas_neocan_telvit_localFstats__2500_500.txt。

看看该文件：

  less altfas_neocan_telvit_localFstats__2500_500.txt

该文件中的行对应于各个染色体窗口，文件中包含的六列包含有关染色体ID以及每个窗口的起始和结束位置的信息。接下来是每行上的三个数字，它们代表第四列中窗口的D-统计量，以及另外两个统计数据（Fd，Fdm），旨在量化受渐渗影响的窗口比例。

使用以下命令查找染色体5被划分的窗口数：

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

这应该表明在Dsuite的Dinvestigate分析中使用了509个窗口，另外该窗口的数量可以通过指定更大或更小的窗口大小来调整信息密度-w

为了可视化5号染色体上D和f d统计数据的变化，使用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，其估计混合比例大多接近0.5：并在大约%Mbp的地方有下降的信号，这是否是由参考基因组中缺失的数据或错误组装导致的假阳性结果，或者它是否显示出减少渐渗的生物信号，如果没有进一步的分析，很难说清楚。

重复这个分析，但是现在使用不同的三个物种的组合,“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

再次进行画图：

  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()

产生该图：

正如你所看到的，D -statistic（黑色）现在显示更多的变化，甚至在几个窗口中变为负数，表明在这些窗口中，neooli与neobri的相似性高于neopul。总体而言，两项统计数据都低于第一个三个物种的组合。值得注意的是，这里在5Mbp的地方，再次可以看到明显的下降，由于同一地区不太可能在两个物种不同三组合中具有更多的渐渗，这可能表明这实际上是由技术的缺陷导致的而非生物过程引起的。