Dsuite據(jù)說是相較于其他相似軟件，占用的資源很少，這就給了我在自己服務(wù)器上操作的空間。
官網(wǎng)：https://github.com/millanek/Dsuite
官網(wǎng)上的教程：https://github.com/millanek/tutorials/tree/master/analysis_of_introgression_with_snp_data

1 安裝

要進(jìn)行編譯，必須擁有相當(dāng)新的 GCC (>=4.9.0) 和 zlib 壓縮庫 (https://www.zlib.net），這也是為什么我裝了新系統(tǒng)就去整了個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統(tǒng)計的所有可能的種群/物種的三元
DtriosCombine：整合Dtrios在基因組區(qū)域運行的結(jié)果（如每條染色體）
Dinvestigate：對D統(tǒng)計量顯著提高的三元進(jìn)行后續(xù)分析：計算f4統(tǒng)計信息，以及沿基因組窗口中的f_d和f_dM
Fbranch（需進(jìn)行上述可選安裝）：計算與種群/物種相關(guān)的樹的分支的D和f統(tǒng)計量

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

準(zhǔn)備文件：

vcf文件，接受bgzip和gzip壓縮文件，可以包含Indel和多等位基因，但不會被計算（本文使用beagle文件是因為其他計算f4的軟件都要求填充）
sets.txt：兩列文件，個體名稱和種群名稱
注意：必須要有一個Outgroup；對于不需要的種群，將種群名改為xxx即可屏蔽而不用提取vcf文件里的樣本。
tree文件（可選，建議加上）：Newick 格式的樹。這棵樹應(yīng)該有與物種/種群名稱相對應(yīng)的標(biāo)簽，與sets.txt文件中對應(yīng)；分支長度可以存在但不會被使用（建議別加長度，可能出錯）。需要注意的是，如果選擇加上樹文件，sets.txt文件里種群名必須是英文+數(shù)字的形式，可以有下劃線。關(guān)于樹文件的構(gòu)建，教程放在末尾。

命令

~/software/Dsuite/Dsuite/Build/Dsuite Dtrios test.beagle.vcf SETS.txt 
#可選命令：
 -k  #默認(rèn)=20，將數(shù)據(jù)集劃分為Jackknife塊的數(shù)量，整個數(shù)據(jù)集至少應(yīng)該是20
-j  #默認(rèn)=NA，JJackknife塊中snp的數(shù)量，如果使用了該命令，則替換-k
-r，--region=start,length #只處理VCF文件的子集，如 --region=20001,10000 將處理 20001 到 30000 的變異位點
-t，--tree= TREE_FILE #一個帶有newick格式的樹的文件，該樹指定種群/物種之間的關(guān)系。根據(jù)樹排列的三組的D和f4-ratio值將輸出在一個后綴為_tree.txt的文件中
-o，--out-prefix=OUT_FILE_PREFIX #指定輸出文件的前綴，默認(rèn)情況下，前綴取自SETS.txt文件的名稱
-n，--run-name #run-name將包含在輸出文件名的PREFIX后面
--no-f4-ratio #不計算f4-ratio
-l NUMLINES #VCF輸入的行數(shù)，如果通過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數(shù)據(jù);也就是說，每個“個體”都是一個種群，然后從AD(等位基因深度)字段估計等位基因頻率，例如MIN_DEPTH=5可能是合理的;當(dāng)讀取次數(shù)較少時，等位基因頻率被設(shè)置為缺失
-c, --no-combine #不輸出"_combine.txt"和"_combine_stderr.txt"文件，這些文件只有在DtriosCombine有用
--KS-test-for-homoplasy #檢測強(qiáng)abba信息位點是否沿基因組聚集

并行程序

此外，還提供了一個DtriosParallel命令，用法類似于Dtrios。允許同時輸入多個vcf文件，該命令會調(diào)用DtriosCombine自動合并vcf文件，所有的vcf文件應(yīng)包含相同的樣本（如不同染色體的vcf文件）

./utils/DtriosParallel SETS.txt INPUT_FILE.vcf [INPUT_FILE.vcf ...]
-j，-k，-t，-n，-l，-g，-p命令同上
-c，--no-combine #不運行DtriosCombine來獲取單個綜合結(jié)果文件
--cores CORES，#(default=CPU count)  設(shè)置進(jìn)程數(shù)
--keep-intermediate #保留區(qū)域Dsuite Dtrios結(jié)果。
 --logging_level {DEBUG,INFO,WARNING,ERROR,CRITICAL}, -v {DEBUG,INFO,WARNING,ERROR,CRITICAL} #最小級別的輸出日志
--dsuite-path DSUITE_PATH  #顯式地將路徑設(shè)置為Dsuite所在的目錄。默認(rèn)情況下，腳本將首先執(zhí)行檢查Dsuite是否可以從$PATH訪問。如果不是，它將嘗試在../Build/Dsuite中定位Dsuite。
--environment-setup ENVIRONMENT_SETUP, #應(yīng)該運行該命令來為Dsuite設(shè)置環(huán)境。例如，'模塊加載GCC'或'conda

相較于Dtrios最大的優(yōu)勢應(yīng)該是可以通過"--cores"設(shè)置進(jìn)程數(shù)，其他大部分有用參數(shù)與上文類似。

結(jié)果文件

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

繪熱圖

首先準(zhǔn)備一個plot_order.txt文件，單列，包含每個種群名稱?；蛘呖梢允褂迷撁睿?/p>

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

其次準(zhǔn)備兩個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

結(jié)果圖.png

該圖是不考慮p1物種，列出p2和p3物種之間最大的D值和f4-radio值，
選擇的種群有無基因流對結(jié)果有很大影響

4 DtriosCombine 整合上一步運行的結(jié)果（如不同染色體的運算結(jié)果）

Dsuite DtriosCombine [OPTIONS] DminFile1 DminFile2 DminFile3 ....
-o，--out-prefix=OUT_FILE_PREFIX #指定輸出文件的前綴，默認(rèn)情況下，前綴為“out”
-t，--tree= TREE_FILE #一個帶有newick格式的樹的文件，該樹指定種群/物種之間的關(guān)系。根據(jù)樹排列的三組的D和f4-ratio值將輸出在一個后綴為_tree.txt的文件中
-s , --subset=start,length #只處理VCF文件的子集，如 --subset=20001,10000 將處理 20001 到 30000 的變異位點

5 Dinvestigate 對D統(tǒng)計量顯著升高的組進(jìn)行后續(xù)分析

準(zhǔn)備文件

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是窗口和步長，對于動物來說，這個值似乎有點小了，建議自己摸索合適的參數(shù)

畫圖

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

結(jié)果圖.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 #輸出文件名，默認(rèn)Fbranch
--outgroup OUTGROUP #newick文件里Outgroup組的名字，默認(rèn)Outgroup
--use_distances #使用newick文件的實際節(jié)點距離繪制樹。(默認(rèn)值:False)
--ladderize #在繪圖之前對輸入樹進(jìn)行階梯化。(默認(rèn)值:False)

結(jié)果圖.png

樹文件的準(zhǔn)備

有兩種方法

第一種：利用treemix

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

詳細(xì)教程見：http://m.itdecent.cn/p/f33bbd8e0996

第二種：PAUP（官網(wǎng)提供的教程中提到的）

這種方法有個問題，小數(shù)據(jù)集還好，但是大數(shù)據(jù)集會分析不出來。比如我的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...”

設(shè)定詳情.png