Questions regarding D-statistic calculation with qpDstat in ADMIXTURE

Dear Professor,

Recently, I encountered some issues while calculating D-statistics using qpDstat from ADMIXTURE. Following the methods described in the paper "Genomic architecture and introgression shape a butterfly radiation", I used one individual per species for the analysis. My pipeline was as follows:

msa_view ./maf/scaffold_9.maf --in-format MAF --out-format FASTA > ./fasta/scaffold_9.fa

snp-sites -cv -o VCF/scaffold_9.vcf ./fasta/scaffold_9.fa

perl -I /data/00/software/vcftools/vcftools-vcftools-954e607/src/perl /data/00/software/vcftools/src/perl/vcf-shuffle-cols -t template.vcf ./VCF/scaffold_9.vcf > VCF_2/scaffold_9.shuffle.vcf (the template.vcf was copied from a randomly selected VCF from the previous step as a reference)

perl -I /data/00/software/vcftools/vcftools-vcftools-954e607/src/perl /data/00/software/vcftools/vcftools-vcftools-954e607/src/perl/vcf-concat VCF_2/*.vcf > allGenomes.vcf

python2 /data/00/software/gdc/vcfToEigenstrat.py -v allGenomes.vcf -o allGenomes

python3 scafNamesToNums_new.py allGenomes.snp allGenomes.geno allGenomes_chr chr_ptg_rename_chr.map (to map to chromosome-level genome)

Finally, use these files for the D-statistic analysis: qpDstat -p qpDstat/file.par > allGenomes_chr.D.out

However, two of my species are actually different geographic subspecies of the same species. In the VCF file, the genetic distance between these two subspecies is extremely small, with a Pearson correlation coefficient as high as 97.12%. Nevertheless, as shown in the figure (I converted the D-statistic results into average D values between different species pairs), the consistency of D-statistics between the first and second columns (i.e., these two subspecies) is very low. Upon further inspection, I found that for the same set of four species, the significance of D-statistics differs greatly between the two subspecies, and even when both are significant, the inferred direction of gene flow often points to different species pairs.

Why does this happen? How can I improve the analysis? Should I use only one of the two subspecies? Or should I merge them into a single group? If merging is recommended, how should I do that?

I would greatly appreciate your guidance. I look forward to your reply.

Best regards,
Liuyj

Start by reviewing qpDstat/file.par and the generated allGenomes_chr.* files, then reproduce the comparison between the two subspecies using the stated MAF-to-VCF and Eigenstrat pipeline. Check whether input preparation and group definitions account for the differing results. Done means documenting the cause and providing a justified recommendation about separate, single, or merged treatment.