Dsuite据说是相较于其他相似软件，占用的资源很少，这就给了我在自己服务器上操作的空间。
官网：https://github.com/millanek/Dsuite
官网上的教程：https://github.com/millanek/tutorials/tree/master/analysis_of_introgression_with_snp_data

1 安装

要进行编译，必须拥有相当新的 GCC (>=4.9.0) 和 zlib 压缩库 (https://www.zlib.net），这也是为什么我装了新系统就去整了个gcc。

git clone https://github.com/millanek/Dsuite.git #直接克隆编译就行
cd Dsuite
make
#可选安装
cd utils
python3 setup.py install --user --prefix= #=后面无任何符号，包括空格，用于画图

2 命令:

Dtrios：计算D (ABBA-BABA)和f4-ratio统计的所有可能的种群/物种的三元
DtriosCombine：整合Dtrios在基因组区域运行的结果（如每条染色体）
Dinvestigate：对D统计量显着提高的三元进行后续分析：计算f4统计信息，以及沿基因组窗口中的f_d和f_dM
Fbranch（需进行上述可选安装）：计算与种群/物种相关的树的分支的D和f统计量

a) Dsuite Dtrios [OPTIONS] INPUT_FILE.vcf SETS.txt
b) Dsuite Dquartets [OPTIONS] INPUT_FILE.vcf SETS.txt
c) Dsuite Dinvestigate [OPTIONS] INPUT_FILE.vcf.gz SETS.txt test_trios.txt
d) Dsuite Fbranch [OPTIONS] TREE_FILE.nwk FVALS_tree.txt

3 Dtrios

准备文件：

vcf文件，接受bgzip和gzip压缩文件，可以包含Indel和多等位基因，但不会被计算（本文使用beagle文件是因为其他计算f4的软件都要求填充）
sets.txt：两列文件，个体名称和种群名称
注意：必须要有一个Outgroup；对于不需要的种群，将种群名改为xxx即可屏蔽而不用提取vcf文件里的样本。
tree文件（可选，建议加上）：Newick 格式的树。这棵树应该有与物种/种群名称相对应的标签，与sets.txt文件中对应；分支长度可以存在但不会被使用（建议别加长度，可能出错）。需要注意的是，如果选择加上树文件，sets.txt文件里种群名必须是英文+数字的形式，可以有下划线。关于树文件的构建，教程放在末尾。

命令

~/software/Dsuite/Dsuite/Build/Dsuite Dtrios test.beagle.vcf SETS.txt 
#可选命令：
 -k  #默认=20，将数据集划分为Jackknife块的数量，整个数据集至少应该是20
-j  #默认=NA，JJackknife块中snp的数量，如果使用了该命令，则替换-k
-r，--region=start,length #只处理VCF文件的子集，如 --region=20001,10000 将处理 20001 到 30000 的变异位点
-t，--tree= TREE_FILE #一个带有newick格式的树的文件，该树指定种群/物种之间的关系。根据树排列的三组的D和f4-ratio值将输出在一个后缀为_tree.txt的文件中
-o，--out-prefix=OUT_FILE_PREFIX #指定输出文件的前缀，默认情况下，前缀取自SETS.txt文件的名称
-n，--run-name #run-name将包含在输出文件名的PREFIX后面
--no-f4-ratio #不计算f4-ratio
-l NUMLINES #VCF输入的行数，如果通过unix管道读取VCF，则必须输入
-g，--use-genotype-probabilities       # (optional) use probabilities (GP tag) or calculate them from likelihoods (GL or PL tags) using a Hardy-Weinberg prior，the probabilities are used to estimate allele frequencies in each population/species（原文贴在这，因为我看不懂）
-p, --pool-seq=MIN_DEPTHVCF #包含pool-seq数据;也就是说，每个“个体”都是一个种群，然后从AD(等位基因深度)字段估计等位基因频率，例如MIN_DEPTH=5可能是合理的;当读取次数较少时，等位基因频率被设置为缺失
-c, --no-combine #不输出"_combine.txt"和"_combine_stderr.txt"文件，这些文件只有在DtriosCombine有用
--KS-test-for-homoplasy #检测强abba信息位点是否沿基因组聚集

并行程序

此外，还提供了一个DtriosParallel命令，用法类似于Dtrios。允许同时输入多个vcf文件，该命令会调用DtriosCombine自动合并vcf文件，所有的vcf文件应包含相同的样本（如不同染色体的vcf文件）

./utils/DtriosParallel SETS.txt INPUT_FILE.vcf [INPUT_FILE.vcf ...]
-j，-k，-t，-n，-l，-g，-p命令同上
-c，--no-combine #不运行DtriosCombine来获取单个综合结果文件
--cores CORES，#(default=CPU count)  设置进程数
--keep-intermediate #保留区域Dsuite Dtrios结果。
 --logging_level {DEBUG,INFO,WARNING,ERROR,CRITICAL}, -v {DEBUG,INFO,WARNING,ERROR,CRITICAL} #最小级别的输出日志
--dsuite-path DSUITE_PATH  #显式地将路径设置为Dsuite所在的目录。默认情况下，脚本将首先执行检查Dsuite是否可以从$PATH访问。如果不是，它将尝试在../Build/Dsuite中定位Dsuite。
--environment-setup ENVIRONMENT_SETUP, #应该运行该命令来为Dsuite设置环境。例如，'模块加载GCC'或'conda

相较于Dtrios最大的优势应该是可以通过"--cores"设置进程数，其他大部分有用参数与上文类似。

结果文件

head species_sets_no_geneflow_BBAA.txt
P1 P2 P3 Dstatistic Z-score p-value f4-ratio BBAA ABBA BABA
S01 S02 S00 0.00645161 0.228337 0.409692 7.6296e-05 40635 936 924
S01 S00 S03 0.0299321 1.08142 0.139754 0.000543936 28841 1724.75 1624.5
S01 S00 S04 0.0072971 0.308069 0.379015 0.00013279 28889.2 1691 1666.5
S01 S00 S05 0.00312175 0.133303 0.446977 5.6877e-05 28861.2 1687 1676.5

参考图.jpg

绘热图

首先准备一个plot_order.txt文件，单列，包含每个种群名称。或者可以使用该命令：

cut -f 2 sets.txt | uniq | head -n 20 > plot_order.txt

其次准备两个ruby脚本文件
链接：https://github.com/millanek/tutorials/blob/master/analysis_of_introgression_with_snp_data/src/plot_d.rb
https://github.com/millanek/tutorials/blob/master/analysis_of_introgression_with_snp_data/src/plot_f4ratio.rb

ruby plot_d.rb SETS_BBAA.txt plot_order.txt 0.7 species_sets_no_geneflow_BBAA_D.svg
ruby plot_f4ratio.rb SETS_BBAA.txt plot_order.txt 0.2 species_sets_no_geneflow_BBAA_f4ratio.svg

结果图.png

该图是不考虑p1物种，列出p2和p3物种之间最大的D值和f4-radio值，
选择的种群有无基因流对结果有很大影响

4 DtriosCombine 整合上一步运行的结果（如不同染色体的运算结果）

Dsuite DtriosCombine [OPTIONS] DminFile1 DminFile2 DminFile3 ....
-o，--out-prefix=OUT_FILE_PREFIX #指定输出文件的前缀，默认情况下，前缀为“out”
-t，--tree= TREE_FILE #一个带有newick格式的树的文件，该树指定种群/物种之间的关系。根据树排列的三组的D和f4-ratio值将输出在一个后缀为_tree.txt的文件中
-s , --subset=start,length #只处理VCF文件的子集，如 --subset=20001,10000 将处理 20001 到 30000 的变异位点

5 Dinvestigate 对D统计量显著升高的组进行后续分析

准备文件

vcf文件和sets.txt：同上
test_trios：列出想要研究的组，可以多行，每行三列（意味着可以直接做想要的组而不必f4）

#P1          #P2          #P3
mbuna   deep    Diplotaxodon

对于这组，描述为：研究deep和Diplotaxodon组的混合信号，并和mbuna比较

命令

Dsuite Dinvestigate -w 50,25 test.beagle.vcf SETS.txt test_trios.txt
其中50 25是窗口和步长，对于动物来说，这个值似乎有点小了，建议自己摸索合适的参数

画图

R脚本
plotInvestigateResults.R

# read in the results with 50 SNP windows and a step of 25 SNPs
bigStep <- read.table("mbuna_deep_Diplotaxodon_localFstats__50_25.txt",as.is=T,header=T) #修改文件名

# plot D in windows along scaffold 18:
plot(bigStep$windowStart/1000000, bigStep$D,type="l",xlab="scaffold 18 coordinate (Mb)",ylab="D (ABBA-BABA)")
# plot f_dM in windows along scaffold 18:
plot(bigStep$windowStart/1000000, bigStep$f_dM,type="l",xlab="scaffold 18 coordinate (Mb)",ylab="f_dM")
# plot f_d in windows along scaffold 18:
plot(bigStep$windowStart/1000000, bigStep$f_d,type="l",xlab="scaffold 18 coordinate (Mb)",ylab="f_d")
# use f_dM and zoom-in on the region of interest 
plot(bigStep$windowStart/1000000, bigStep$f_dM,type="l",xlim=c(4,4.5),xlab="scaffold 18 coordinate (Mb)",ylab="f_dM",ylim=c(-0.2,0.8))

# Same window: 50SNPs, and the smallest 1 SNP step
smallestStep <- read.table("mbuna_deep_Diplotaxodon_localFstats__50_1.txt",as.is=T,header=T)
# plot f_dM again:
plot(smallestStep$windowStart/1000000, smallestStep$f_dM,type="l",xlim=c(4,4.5),xlab="scaffold 18 coordinate (Mb)",ylab="f_dM",ylim=c(-0.2,0.8))

# Smaller 10SNP window, and the smallest 1 SNP step
smallerWindow <- read.table("mbuna_deep_Diplotaxodon_localFstats__10_1.txt",as.is=T,header=T)
# plot f_dM again:
plot(smallerWindow$windowStart/1000000, smallerWindow$f_dM,type="l",xlim=c(4.0,4.5),xlab="scaffold 18 coordinate (Mb)",ylab="f_dM",ylim=c(-0.2,0.9),pch=16)

# Finally, a tiny window with 2 SNPs and a step of 1 SNP
smallWindow <- read.table("mbuna_deep_Diplotaxodon_localFstats__2_1.txt",as.is=T,header=T)
plot(smallWindow$windowStart, smallWindow$f_dM,type="l",xlim=c(4000000,4500000),xlab="scaffold 18 coordinate (Mb)",ylab="f_dM",ylim=c(-0.2,0.8))

结果图.png

6 Fbranch

命令

Dsuite Fbranch simulated_tree_with_geneflow.nwk species_sets_with_geneflow_tree.txt > species_sets_with_geneflow_Fbranch.txt #可以在这一步修改种群名为想要的
python3 /Users/milanmalinsky/Sanger_work/Dsuite/Dsuite/utils/dtools.py species_sets_with_geneflow_Fbranch.txt simulated_tree_with_geneflow.nwk
-n RUN_NAME, --run-name RUN_NAME #输出文件名，默认Fbranch
--outgroup OUTGROUP #newick文件里Outgroup组的名字，默认Outgroup
--use_distances #使用newick文件的实际节点距离绘制树。(默认值:False)
--ladderize #在绘图之前对输入树进行阶梯化。(默认值:False)

结果图.png

树文件的准备

有两种方法

第一种：利用treemix

treemix -i sample.treemix.in.gz -o sample.ML.tree -root wBGU 

详细教程见：https://www.jianshu.com/p/f33bbd8e0996

第二种：PAUP（官网提供的教程中提到的）

这种方法有个问题，小数据集还好，但是大数据集会分析不出来。比如我的180个个体。而且也找不到Linux下安装的教程。
首先下载PAUP：http://phylosolutions.com/paup-test/
convert_vcf_to_nexus.rb脚本文件：https://github.com/millanek/tutorials/blob/master/species_tree_inference_with_snp_data/src/convert_vcf_to_nexus.rb
taxpartitions.txt种群信息文件：如下，意味着1，2序列是astbur，以此类推

BEGIN SETS;
    TAXPARTITION SPECIES =
        astbur: 1-2,
        altfas: 3-4,
        telvit: 5-6,
        neobri: 7-8,
        neochi: 9-10,
  END;

然后运行，挺慢的

ruby convert_vcf_to_nexus.rb nc.test.beagle.vcf test.beagle.nex #不接受压缩格式

完事之后检查一下nex文件中的个体编号，不能有"-"，否则会报错（一个快速的修改方法：使用paup的编辑功能）。
然后合并种群和nex文件

  cat test.beagle.nex taxpartitions.txt > pop.test.beagle.nex

得到的文件用paup打开：File→open→pop.test.beagle.nex→Execute
指定外群：Data→Define Outgroup
分析：Analysis→“SVDQuartets...”

设定详情.png