Dfoil – 推断祖先基因流

-2.

在一篇讲泥炭藓的MBE看到了dfoil的使用

《Extensive Genome-Wide Phylogenetic Discordance Is Due to Incomplete Lineage Sorting and Not Ongoing Introgression in a Rapidly Radiated Bryophyte Genus》

文章的逻辑是这样的：

基因树长得和物种树不一样，核基因树长得和细胞器基因树不一样，为了知道为啥长得不一样就做了基因流

D统计的值都不是很高，但是很显著。QuIBL推断的结果是这种系统发育的不一致主要是因为不完全谱系分选导致，而不是物种形成后的渐渗(关于QuIBL的这句话是文章原话，翻译过来的)。

假设基因流动是正在进行的或最近的，它可能会被限制在密切相关的物种，因为自分歧以来，生殖障碍随着时间的推移而积累-CoyneJA, Orr HA. 2004. Speciation. Sunderland (MA): Sinauer Associates

D统计值与物种分化呈负相关，而作者做出来D与Fst呈正相关，再加上之前做的admixture看不出来有近期基因流

综上，推断 可能是祖先基因流。于是要验证这个推断

所以做了Dfoil (https://github.com/jbpease/dfoil)，这个软件的优点是能推断P1/P2的祖先和P3、P4之间有没有基因流(但不能判断祖先基因流的方向)，以及非祖先基因流的方向

James B Pease, Matthew W. Hahn. 2015."Detection and Polarization of Introgression in a Five-taxon Phylogeny" Systematic Biology. 64 (4): 651-662.http://www.dx.doi.org/10.1093/sysbio/syv023 doi: 10.1093/sysbio/syv023

dfoil要求物种满足下图所示系统发生关系

((P1,P2),(P3,P4),O)且P3和P4分化的时间不能晚于P1P2

结果如下：

不同颜色背景的，是不同属的物种，蓝色这个属正好也是前面核基因树和细胞器基因树对不上的物种

数字表示窗口中支持祖先基因流的窗口的比例，总之就是推出了蓝色物种的祖先和红色物种有基因流

当然除了这个证据，作者还列举了第二个证据：

“渐渗区域的长度(length of consecutive introgressed region)可以用来推断渐渗发生的相对时间[1,2,3]，假设在近期发生了渐渗，预期在基因组中会发现更长的渐渗片段，随着时间推移由于重组，这些片段会分裂成更小的片段[4]”

[1]Barlow A, Cahill JA, Hartmann S, Theunert C, Xenikoudakis G, Fortes GG, Paijmans JLA, Rabeder G, Frischauf C, Grandal-d’Anglade A, et al. 2018. Partial genomic survival of cave bears in living brown bears.Nat Ecol Evol. 2(10):1563–1570.

[2]Moodley Y, Westbury MV, Russo I-RM,Gopalakrishnan S, RakotoariveloA, Olsen R-A, Prost S, Tunstall T, RyderOA, Dal? en L, et al. 2020.Interspecific gene flow and the evolution of specialisation in black and white rhinoceros. Mol Biol Evol. 37(11):3105–3117

[3]Westbury MV , Hartmann S .Barlow A ,Preick M , Ridush B , Nagel D ,Rathgeber T, Ziegler R, Baryshnikov G, Sheng G,et al. 2020. Hyena paleogenomes reveal a complex evolutionary history of cross-continental gene flow between spotted and cave hyena. Sci Adv.6(11):eaay0456

[4]We Rcek K, Hartmann S, Paijmans JLA, Taron U, Xenikoudakis G, Cahill JA, Heintzman PD, Shapiro B, Baryshnikov G, Bunevich AN, et al. 2017.Complex admixture preceded and followed the extinction of wisent in the wild. Mol Biol Evol. 34:598–612.

作者是咋做的呢，他看那些表示有渐渗的窗口是不是连在一块的

这个图是拿表的数据画的

惊人的发现绝大部分窗口都是单个单个的，更映照了是早期基因流，发生在祖先谱系中

综上所述, these results strongly support the hypothesis of ancient introgression among the ancestral Sphagnum species.

泥潭藓广泛的系统发育不一致是由于ILS和不持续基因流导致的，点题。

-1  Dfoil的下载与安装

Dfoil下载下来就能直接用了，主体是几个脚本，详细参数参见user manual

git clone https://www.github.com/jbpease/dfoil

0. 输入文件的准备

Dfoil要求的输入文件是这样的，每个物种一行fasta序列，并且要对齐，上面那篇文章是每个物种选了一个sequcencing coverage最高的个体来做。

用ANGSD生成每个物种的共有序列，然后进行过滤，过滤小于1M - 输出结果每个物种都会对齐，也用不着mafft()，然后按100kb窗口切分，接着把窗口里面N超过50%的过滤掉，最后再把一个窗口的放到一块。(这一步我也不清楚，输入文件是很久之前师姐就准备好的)

-3.具体做法是，每个物种挑选一个个体，用bwa mem比到同一个参考基因组上，每个物种就获得了一个bam文件

-2. angsd -doFasta2 -i 1.bam -o 1.fasta 把每个bam转成(共有)序列，实际上也是比对齐的序列

-1. 把比对齐的序列划窗口

seqkit可以实现. 这里表示窗口为4.步长也为4

最后数据有600个窗口，每个窗口(每个fasta文件)有9个物种(准确的说9个个体)

用dfoil的时候fasta文件里面只允许有5个个体，所以

0-1  samtools faidx input.fasta P1 P2 P3 P4 O > tmp.fasta  #提取5个物种

0-2  awk '!/^>/ { printf "%s", $0; n = "/n" } /^>/ { print n $0; n = "" } END { printf "%s", n }' tmp.fasta > clean.fasta  #把多行fasta变成一行fasta，并且保留原来的头

1.

python2  /PATH/fasta2dfoil.py --out xx.counts --names P1,P2,P3,P4,O 1.fa 2.fa 3.fa ......

/PATH/dfoil.py --infile xx.counts --out xx.out --plot xx.pdf  #plot这个参数不太重要，它画的图我感觉也用不上

这个xx.counts就是我们要的结果了，总共600行，因为我们有600个窗口

introgression这一列如果不是none，就表明这个窗口支持有基因流，后面几列就是基因流发生的方向与位置，introg23就是指P2到P3的基因流。introg123就是指P1P2祖先和P3之间是基因流。

2.可视化

这一块可以参照上文提到的MBE的泥炭藓

B图是这么多种物种组合，基因流都是啥款式的，是明显看到P12的祖先到P3的基因流很突出，C图就是这么突出的基因流，具体是那些P1 P2 P3。

C图说过了

我们可以统计每种基因流窗口的比例，然后参照着这两个图画出来，略、

附原文对B图 C图的注解