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

Dsuite軟件文章:Malinsky, M., Matschiner, M. and Svardal, H. (2021) Dsuite ‐ fast D‐statistics and related admixture evidence from VCF files. Molecular Ecology Resources 21, 584–595. doi: https://doi.org/10.1111/1755-0998.13265

1. 軟件安裝

git clone https://github.com/millanek/Dsuite.git
cd Dsuite
make

Dsuite主要包含三個:“Dtrios”破讨,“DtriosCombine”和“Dinvestigate”不同命令尽狠。
Dtrios計算所有可能的種群/物種的D(ABBA-BABA)統(tǒng)計和f4統(tǒng)計曹货,主要的參數(shù)是-j谊娇,用于設(shè)定計算窗口的大小范圍:

Usage: Dsuite Dtrios [OPTIONS] INPUT_FILE.vcf SETS.txt
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 SETS.txt 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_FILE.nwk               (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 _tree.txt suffix
-n, --run-name                          run-name will be included in the output file name

2. 輸入文件

(1)VCF文件

可以用gzip或bgzip壓縮鼓拧,它可以包含多等位基因座和插入/缺失,但僅使用雙等位基因SNP。

(2)個體/群體SETS.txt

一個文本文件洼哎,每行一個樣本烫映,一個制表符將該樣本的名字與其所屬的物種/種群的名稱分開沼本,如下所示:

Ind1    Species1
Ind2    Species1
Ind3    Species2
Ind4    Species2
Ind5    Species3
Ind6    Outgroup
Ind7    Outgroup

Dtrios需要將至少一個樣本指定為outgroup。
Dquartets對所有物種/種群均一視同仁锭沟,不應(yīng)指定任何外群抽兆。

(3)Newick格式的樹文件(可選的)

樹文件應(yīng)具有與物種/種群名稱相對應(yīng)的葉子標(biāo)簽,不使用枝長族淮,樹必需有根辫红。

(Species2:6.0,(Species1:5.0,(Species3:3.0,Species4:4.0)));
(4)test_trios.txt文件(可選的)

每行三個種群/物種,以制表符按P1 P2 P3的順序分隔

Species1    Species2    Species3
Species1    Species4    Species2

3. 軟件運行

Dsuite Dtrios INPUT_FILE.vcf.gz SETS.txt
#或者添加過濾:最小深度至少為1000
NUMLINES=$(bcftools view -i 'INFO/DP>1000' INPUT_FILE.vcf.gz | wc -l)  # to get NUMLINES
bcftools view -i 'INFO/DP>1000' INPUT_FILE.vcf.gz | Dsuite Dtrios -l $NUMLINES stdin SETS.txt

4. 結(jié)果解讀

帶有后綴BBAA.txt祝辣,Dmin.txt和可選的tree.txt(如果使用了-t選項)的輸出文件包含以下結(jié)果:D統(tǒng)計量贴妻,Zscore,未矯正的p值蝙斜,f4-比率以及 BBAA名惩,BABA和ABBA模式的計數(shù)情況。
帶有后綴Combine.txt和Combine_stderr.txt的輸出文件用作DtriosCombine的輸入孕荠。 如果不需要使用DtriosCombine娩鹉,則可以刪除這些文件。
(1) samples_combine.txt文件

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

在這里稚伍,每一行顯示分析一個三個物種的分析結(jié)果弯予,例如在第一行中,altfas用作P1个曙,neobri被認為是P2锈嫩,而neocan被放置在P3。然后該行的第五和第六列中的數(shù)字分別代表著:ABBA位點在該三個物種中的數(shù)量(C-ABBA)(其中衍生的等位基因由“neobri”和“neocan”共享)和BABA位點的計數(shù)垦搬,C-BABA(衍生的等位基因由“altfas”和“neocan”共享)呼寸。除了第5和第6列中BABA和ABBA位點的數(shù)量之外,第4列列出了“BBAA”位點的數(shù)量(C-BBAA)悼沿,P1和P2共享衍生的等位基因(因此通過“altfas”和“ neobri“共享)等舔。

ABBA-BABA測試基于P1和P2是姊妹物種的假設(shè),當(dāng)計算給定三個物種的D-統(tǒng)計量時糟趾,Dsuite首先重新排列分配給P1慌植,P2和P3的物種(因此ABBA甚牲,BABA和BBAA位點的數(shù)量也重新排列),這是根據(jù)某些規(guī)則:

  • 在第一組計算中蝶柿,測試所有可能的重新排列丈钙,并且報告每給定三個物種的最低D-統(tǒng)計量,稱為D min交汤。因此雏赦,D min是給定三個物種中D-統(tǒng)計量的保守估計。此輸出將寫入具有__Dmin.txt結(jié)尾的文件芙扎。
  • 這樣的排列使得P1和P2是給定三個物種的兩種物種星岗,它們共享最大數(shù)量的衍生地點。換句話說戒洼,進行重排以使C-BBAA > C-ABBA和C-BBAA > C-BABA俏橘。另外,旋轉(zhuǎn)P1和P2圈浇,使得在P2和P3之間共享的派生位點的數(shù)量大于P1和P3之間的派生站點的數(shù)量寥掐。總之磷蜀,這意味著C-BBAA > C-ABBA >C-BABA召耘。這種類型的重排背后的想法是,共享最大數(shù)量的衍生地點的兩個物種很可能是給定三個物種中真正的兩個姊妹物種褐隆,因此正確地放置在位置P1和P2污它。預(yù)計這種假設(shè)在某些條件下會持續(xù)存在(例如時鐘式演化,缺乏同質(zhì)性妓灌,缺乏漸滲轨蛤,以及泛化的祖先種群),但有時候很難說真實數(shù)據(jù)集的可靠性虫埂∠樯剑基于這種類型的重新排列的D-統(tǒng)計被寫入具有結(jié)尾__BBAA.txt的文件。
  • 要直接告訴Dsuite哪個給定三個物種應(yīng)該被視為姊妹物種(因此掉伏,應(yīng)該將其分配給P1和P2)缝呕,我們可以使用參數(shù)-t或提供包含所有數(shù)據(jù)集種類的樹文件--tree。然后輸出將被寫入帶有結(jié)尾_tree.txt的附加文件斧散。

samples_Dmin.txt

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

第四列現(xiàn)在顯示每一個給定三個物種的的D-統(tǒng)計量供常,第五列顯示基于對D = 0 的零假設(shè)的歸一化的p值。

讓我們選擇一個給定三個物種的例子鸡捐,看看它是如何出現(xiàn)在三個不同的文件samples_combine.txt中samples_Dmin.txt栈暇,和samples_BBAA.txt。我們將選擇第一個給定三個物種的例子箍镜,包括“altfas”源祈,“neocan”和“neobri”的那一行煎源。要僅從三個文件中查看此給定三個物種的例子的行,可以使用此命令:

cat samples_combine.txt | grep altfas | grep neocan | grep neobri
cat samples_Dmin.txt | grep altfas | grep neocan | grep neobri
cat samples_BBAA.txt | grep altfas | grep neocan | grep neobri

這應(yīng)該會分別輸出以下三行:

###samples_combine.txt
altfas  neobri  neocan  1378.2  13479.6 4066.95
###samples_Dmin.txt
altfas  neocan  neobri  0.493787    0
###samples_BBAA.txt
altfas  neocan  neobri  0.493787    0

這里簡單解析一下結(jié)果香缺,首先samples_combine.txt品種名字母排序的方式與其它兩個文件有點不同手销,P1在三個文件中保持相同(“altfas”),但是P2和P3的順序被交換了(“neocan”和“neobri”)图张。此交換還暗示ABBA锋拖,BABA和BBAA模式的計數(shù)相應(yīng)地交換。因此在交換之后并且P1 =“altfas”祸轮,P2 =“neocan”兽埃,P3 =“neobri”,計數(shù)如下:C-ABBA = 4066.95倔撞,C-BABA = 1378.2讲仰,C-BBAA = 13479.6慕趴。因此痪蝇,“neocan”和“neobri”共享4066.95個衍生位點,“altfas”和“neobri”共享1378.2個衍生位站冕房,“altfas”和“neocan”共享13479.6個衍生位站躏啰。有了這些數(shù)量,D=(4066.95 - 1378.2)/(4066.95 + 1378.2)= 0.493787耙册。這個數(shù)字與Dsuite在這兩個文件(samples_Dmin.txt和samples_BBAA.txt给僵。)生成的報告一致。

重復(fù)與上一步相同的操作详拙,但這一次使用給定開頭為neo的三個物種的帝际,文件samples__Dmin.txt和重新排列不同samples__BBAA.txt。三個“neobri”饶辙,“neocra”和“neogra”就是這樣一個例子蹲诀。使用這些命令可以在所有三個文件中查看此給定開頭為neo的行:

cat samples_combine.txt | grep neobri | grep neocra | grep neogra
cat samples_Dmin.txt | grep neobri | grep neocra | grep neogra
cat samples_BBAA.txt | grep neobri | grep neocra | grep neogra

結(jié)果會分別輸出以下的行:

###samples_combine.txt
neobri  neocra  neogra  3788.23 3552.38 2992.93
###samples_Dmin.txt
neogra  neocra  neobri  0.0321294   0.145298
###samples_BBAA.txt
neocra  neobri  neogra  0.0854723   8.58201e-07

<article class="_2rhmJa">

文件中的結(jié)果samples_BBAA.txt表明,當(dāng)P1 =“neocra”弃揽,P2 =“neobri”脯爪,P3 =“neogra”時,則C-BBAA = 3788.23矿微,C-ABBA = 3552.38痕慢,C-BABA = 2992.93,因此D =(3552.38 - 2992.93)/(3552.38 + 2992.93)= 0.0854723涌矢。

然而掖举,文件中的結(jié)果samples_Dmin.txt顯示,這次發(fā)生了另一次重新排列(因此C-BBAA不大于其他兩個計數(shù)的重新排列)產(chǎn)生較低的D-統(tǒng)計:P1 =“neogra”娜庇,P2 =“neocra” 塔次,并且P3 =“neobri”滨巴,則C-BBAA = 2992.93,C-ABBA = 3788.23俺叭,并且C-BABA = 3552.38是鬼,因此D =(3788.23-3552.38)/(3788.23 + 3552.38)= 0.0321294。

這說明了D-min值沈撞,報告了可能的最低D值對于給定三物種的統(tǒng)計灌曙,有時選擇該三個物種的重新排列,其中P1和P2實際上彼此共享較少的派生位點裕照,而不是它們兩者與P3共享攒发。這與ABBA-BABA測試的原始假設(shè)相沖突,即P1和P2彼此之間的關(guān)系比P3更緊密晋南。在解釋Dsuite分析的結(jié)果時惠猿,應(yīng)該記住,在文件中報告的D min值_Dmin.txt實際上是D -statistic的保守估計负间。文件中報告的值以_BBAA.txt結(jié)尾為基礎(chǔ)偶妖,這些值基于確保C-BBAA > C-ABBA > C-BABA的重新排列,通痴#可以更好地測量D -statistic趾访,但是,最好的選擇可能是使用--tree選項運行分析董虱,提供一個輸入樹扼鞋,直接告訴Dsuite如何重新排列所有三個物種。
參考:
https://github.com/millanek/Dsuite
http://www.reibang.com/p/e97eb7b4b2ca

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末愤诱,一起剝皮案震驚了整個濱河市云头,隨后出現(xiàn)的幾起案子,更是在濱河造成了極大的恐慌淫半,老刑警劉巖溃槐,帶你破解...
    沈念sama閱讀 217,185評論 6 503
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件,死亡現(xiàn)場離奇詭異撮慨,居然都是意外死亡竿痰,警方通過查閱死者的電腦和手機,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 92,652評論 3 393
  • 文/潘曉璐 我一進店門砌溺,熙熙樓的掌柜王于貴愁眉苦臉地迎上來影涉,“玉大人,你說我怎么就攤上這事规伐⌒非悖” “怎么了?”我有些...
    開封第一講書人閱讀 163,524評論 0 353
  • 文/不壞的土叔 我叫張陵,是天一觀的道長鲜棠。 經(jīng)常有香客問我肌厨,道長,這世上最難降的妖魔是什么豁陆? 我笑而不...
    開封第一講書人閱讀 58,339評論 1 293
  • 正文 為了忘掉前任柑爸,我火速辦了婚禮,結(jié)果婚禮上盒音,老公的妹妹穿的比我還像新娘表鳍。我一直安慰自己,他們只是感情好祥诽,可當(dāng)我...
    茶點故事閱讀 67,387評論 6 391
  • 文/花漫 我一把揭開白布譬圣。 她就那樣靜靜地躺著,像睡著了一般雄坪。 火紅的嫁衣襯著肌膚如雪厘熟。 梳的紋絲不亂的頭發(fā)上,一...
    開封第一講書人閱讀 51,287評論 1 301
  • 那天维哈,我揣著相機與錄音绳姨,去河邊找鬼。 笑死笨农,一個胖子當(dāng)著我的面吹牛就缆,可吹牛的內(nèi)容都是我干的。 我是一名探鬼主播谒亦,決...
    沈念sama閱讀 40,130評論 3 418
  • 文/蒼蘭香墨 我猛地睜開眼,長吁一口氣:“原來是場噩夢啊……” “哼空郊!你這毒婦竟也來了份招?” 一聲冷哼從身側(cè)響起,我...
    開封第一講書人閱讀 38,985評論 0 275
  • 序言:老撾萬榮一對情侶失蹤狞甚,失蹤者是張志新(化名)和其女友劉穎锁摔,沒想到半個月后,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體哼审,經(jīng)...
    沈念sama閱讀 45,420評論 1 313
  • 正文 獨居荒郊野嶺守林人離奇死亡谐腰,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 37,617評論 3 334
  • 正文 我和宋清朗相戀三年,在試婚紗的時候發(fā)現(xiàn)自己被綠了涩盾。 大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片十气。...
    茶點故事閱讀 39,779評論 1 348
  • 序言:一個原本活蹦亂跳的男人離奇死亡,死狀恐怖春霍,靈堂內(nèi)的尸體忽然破棺而出砸西,到底是詐尸還是另有隱情,我是刑警寧澤,帶...
    沈念sama閱讀 35,477評論 5 345
  • 正文 年R本政府宣布芹枷,位于F島的核電站衅疙,受9級特大地震影響,放射性物質(zhì)發(fā)生泄漏鸳慈。R本人自食惡果不足惜饱溢,卻給世界環(huán)境...
    茶點故事閱讀 41,088評論 3 328
  • 文/蒙蒙 一、第九天 我趴在偏房一處隱蔽的房頂上張望走芋。 院中可真熱鬧理朋,春花似錦、人聲如沸绿聘。這莊子的主人今日做“春日...
    開封第一講書人閱讀 31,716評論 0 22
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽熄攘。三九已至兽愤,卻和暖如春,著一層夾襖步出監(jiān)牢的瞬間挪圾,已是汗流浹背浅萧。 一陣腳步聲響...
    開封第一講書人閱讀 32,857評論 1 269
  • 我被黑心中介騙來泰國打工, 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留哲思,地道東北人洼畅。 一個月前我還...
    沈念sama閱讀 47,876評論 2 370
  • 正文 我出身青樓,卻偏偏與公主長得像棚赔,于是被迫代替她去往敵國和親帝簇。 傳聞我的和親對象是個殘疾皇子,可洞房花燭夜當(dāng)晚...
    茶點故事閱讀 44,700評論 2 354

推薦閱讀更多精彩內(nèi)容