基于VCF文件做基因渗入分析(Dsuite)

2023年12月7日发(作者：)

基于VCF文件做基因渗入分析（Dsuite）

1. 软件安装

git clone /millanek/

cd Dsuite

make

Dsuite主要包含三个：“Dtrios”，“DtriosCombine”和“Dinvestigate”不同命令。

Dtrios计算所有可能的种群/物种的D（ABBA-BABA）统计和f4统计，主要的参数是-j，用于设定计算窗口的大小范围：

Usage: Dsuite Dtrios [OPTIONS] INPUT_

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 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_ (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 _ suffix

-n, --run-name run-name will be included in the output file name

2. 输入文件

（1）VCF文件

可以用gzip或bgzip压缩，它可以包含多等位基因座和插入/缺失，但仅使用双等位基因SNP。

（2）个体/群体

一个文本文件，每行一个样本，一个制表符将该样本的名字与其所属的物种/种群的名称分开，如下所示：

Ind1 Species1

Ind2 Species1

Ind3 Species2

Ind4 Species2

Ind5 Species3

Ind6 Outgroup

Ind7 Outgroup

Dtrios需要将至少一个样本指定为outgroup。

Dquartets对所有物种/种群均一视同仁，不应指定任何外群。

（3）Newick格式的树文件(可选的)

树文件应具有与物种/种群名称相对应的叶子标签，不使用枝长，树必需有根。

(Species2:6.0,(Species1:5.0,(Species3:3.0,Species4:4.0)));

（4）test_文件(可选的)

每行三个种群/物种，以制表符按P1 P2 P3的顺序分隔

Species1 Species2 Species3

Species1 Species4 Species2

3. 软件运行

Dsuite Dtrios INPUT_

#或者添加过滤：最小深度至少为1000

NUMLINES=$(bcftools view -i 'INFO/DP>1000' INPUT_ | wc -l) # to get NUMLINES

bcftools view -i 'INFO/DP>1000' INPUT_ | Dsuite Dtrios -l $NUMLINES stdin

4. 结果解读

带有后缀，和可选的（如果使用了-t选项）的输出文件包含以下结果：D统计量，Zscore，未矫正的p值，f4-比率以

及 BBAA，BABA和ABBA模式的计数情况。

带有后缀和Combine_的输出文件用作DtriosCombine的输入。 如果不需要使用DtriosCombine，则可以删除这些文

件。

(1) samples_文件

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用作P1，neobri被认为是P2，而neocan被放置在P3。然后该行的

第五和第六列中的数字分别代表着：ABBA位点在该三个物种中的数量(C-ABBA)（其中衍生的等位基因由“neobri”和“neocan”共享）和

BABA位点的计数，C-BABA（衍生的等位基因由“altfas”和“neocan”共享）。除了第5和第6列中BABA和ABBA位点的数量之外，第4列

列出了“BBAA”位点的数量(C-BBAA)，P1和P2共享衍生的等位基因（因此通过“altfas”和“ neobri“共享）。

ABBA-BABA测试基于P1和P2是姊妹物种的假设，当计算给定三个物种的D-统计量时，Dsuite首先重新排列分配给P1，P2和P3的物种（因此

ABBA，BABA和BBAA位点的数量也重新排列），这是根据某些规则：

在第一组计算中，测试所有可能的重新排列，并且报告每给定三个物种的最低D-统计量，称为D min。因此，D min是给定三个物种中D-统计

量的保守估计。此输出将写入具有__结尾的文件。

这样的排列使得P1和P2是给定三个物种的两种物种，它们共享最大数量的衍生地点。换句话说，进行重排以使C-BBAA > C-ABBA和C-

BBAA > C-BABA。另外，旋转P1和P2，使得在P2和P3之间共享的派生位点的数量大于P1和P3之间的派生站点的数量。总之，这意味着C-

BBAA > C-ABBA >C-BABA。这种类型的重排背后的想法是，共享最大数量的衍生地点的两个物种很可能是给定三个物种中真正的两个姊妹

物种，因此正确地放置在位置P1和P2。预计这种假设在某些条件下会持续存在（例如时钟式演化，缺乏同质性，缺乏渐渗，以及泛化的祖先种

群），但有时候很难说真实数据集的可靠性。基于这种类型的重新排列的D-统计被写入具有结尾__的文件。

要直接告诉Dsuite哪个给定三个物种应该被视为姊妹物种（因此，应该将其分配给P1和P2），我们可以使用参数-t或提供包含所有数据集种类

的树文件--tree。然后输出将被写入带有结尾_的附加文件。

samples_

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

第四列现在显示每一个给定三个物种的的D-统计量，第五列显示基于对D = 0 的零假设的归一化的p值。

让我们选择一个给定三个物种的例子，看看它是如何出现在三个不同的文件samples_中samples_，和

samples_。我们将选择第一个给定三个物种的例子，包括“altfas”，“neocan”和“neobri”的那一行。要仅从三个文件中查看此

给定三个物种的例子的行，可以使用此命令：

cat samples_ | grep altfas | grep neocan | grep neobri

cat samples_ | grep altfas | grep neocan | grep neobri

cat samples_ | grep altfas | grep neocan | grep neobri

这应该会分别输出以下三行：

###samples_

altfas neobri neocan 1378.2 13479.6 4066.95

###samples_

altfas neocan neobri 0.493787 0

###samples_

altfas neocan neobri 0.493787 0

这里简单解析一下结果，首先samples_品种名字母排序的方式与其它两个文件有点不同，P1在三个文件中保持相同

（“altfas”），但是P2和P3的顺序被交换了（“neocan”和“neobri”）。此交换还暗示ABBA，BABA和BBAA模式的计数相应地交换。

因此在交换之后并且P1 =“altfas”，P2 =“neocan”，P3 =“neobri”，计数如下：C-ABBA = 4066.95，C-BABA = 1378.2，C-

BBAA = 13479.6。因此，“neocan”和“neobri”共享4066.95个衍生位点，“altfas”和“neobri”共享1378.2个衍生位

站，“altfas”和“neocan”共享13479.6个衍生位站。有了这些数量，D=（4066.95 - 1378.2）/（4066.95 + 1378.2）= 0.493787。这

个数字与Dsuite在这两个文件（samples_和samples_。）生成的报告一致。

重复与上一步相同的操作，但这一次使用给定开头为neo的三个物种的，文件samples__和重新排列不同samples__。三

个“neobri”，“neocra”和“neogra”就是这样一个例子。使用这些命令可以在所有三个文件中查看此给定开头为neo的行：

cat samples_ | grep neobri | grep neocra | grep neogra

cat samples_ | grep neobri | grep neocra | grep neogra

cat samples_ | grep neobri | grep neocra | grep neogra

结果会分别输出以下的行：

###samples_

neobri neocra neogra 3788.23 3552.38 2992.93

###samples_

neogra neocra neobri 0.0321294 0.145298

###samples_

neocra neobri neogra 0.0854723 8.58201e-07