Hi all,
I have haploidized whole-genome fasta sequenses for both extant and extinct species. I want to apply D-statistic to assess possible introgression events, but since I deal with ancient sequences It seems better to restrict my analysis to transversion site only to avoid bias resulting from cytosin deamination. How I can perform it? Is there some scripts, allowed to do in after alignment?
I've attached a script that will do what you want. You'll need Python (>3.6) installed, with Numpy and Scipy libraries (see: https://docs.anaconda.com/anaconda/install/ for an easy to install python environment).
The script will read a PHYLIP formatted alignment, and calculate the D-statistic using transversions only. It will compute significance using a specified number of bootstraps, and calculating a Z-score.
Usage:
python3 calcD.py example.phy P1 P2 P3 P4 100
Where P1 (etc) are the sample names for the P1 though P4 lineages for the D-statistic, and "100" is the number of bootstraps. P4 is the outgroup, and we are testing introgression from P3 P1 and P3 P2
For example if your phylip file is:
4 15
A ATGTTACTGTAAAAA
B GTGGCTGAGATCCCT
C GCGTCTCTCAAAAAA
D ACGGTAGACTTCCCT
You would call it like so:
python3 calcD.py example.phy A B C D 100
Let me know if you have any issues. I wrote it just now, so it isn't extensively tested. If you want to change to using transitions as well, you can change "calcD_transversions_only" to "calcD" at lines 135 and 143.