之前在做叶绿体基因组核酸多样的时候，先是用全基因组做，选择窗口和步长使用dnasp5来求pi值。后来发现文章里放的都是编码序列和非编码序列，或者每个基因的pi值。叶绿体的编码序列一般在80个多个，如果再加上非编码区，那么有一百多个序列需要计算，手动算的话挺费时费劲的，一分钟算三个的话得近一个小时。。。作为懒人的我肯定是不愿意每天花一个小时做这样的重复工作，于是……
  我找了文章里的计算方法，发现使用的竟然全部是dnasp（难道都是一个个算的么？dnasp是windows的一款GUI软件，似乎并不能批量计算），后来发现vcftools好像可以计算pi，而且是linux版本的，通过脚本很容易实现批量计算，于是抓紧试了试。但是vcftools只支持vcf格式的文件，没关系，把比对的序列转成vcf格式就行了，最终实现批量计算多个基因的pi值，小小开心下。但是有个问题，同样的序列用vcftools计算得到的pi和用dnasp计算得到的pi不一样。开始是觉得我转换得到vcf的时候不对。于是我找到了另一款计算pi的工具：perl中有个Bio::PopGen::Statistics模块，可以计算pi，而且输入的是比对后的序列。但是！结果仍然和dnasp计算的不一样！令人惊讶的是，popgene和vcftools的计算结果是相同的，我百思不得其解。难道是dnasp算的有问题。于是找到了计算pi的公式，手动算了下，结果和dnasp的值相同。难道。。。另外两个算错了？？？于是研究了Bio::PopGen::Statistics模块中计算pi值的代码，终于找到了原因。
   Bio::PopGen::Statistics在计算pi的时候和vcftools类似，先把比对结果转换成vcf格式，然后计算。但是不同的是，如果两个序列比如AAA,AAT，直接计算的话差异碱基数为1。当转成vcf的时候是用 3 A T 0/0 1/1，来表示(3号位点其中一个是A，另一个是T)。这时候的差异碱基数就变成了4（把0/0，1/1分别看成两个位点，如果还原成序列的话就是A，A，T，T，及从原来的2条序列变成了4条。。。。），实际计算的时候是把样品数当成4而不是2。也就是说，差异碱基数变成了原先的四倍，样品（序列）数变成了原先的两倍，所以最终算的值和dnasp的结果不同。
   考虑到vcftools是用于vcf文件的计算，如果有个位点的信息是0/1，那么确实要表示成两个碱基，所以样品数量要变成原先的两倍，计算的也没毛病。dnasp是直接计算序列的pi，没有考虑太多，和直接套公式计算的结果相同，结果没得说。Bio::PopGen::Statistics呢，我给的是比对的序列，然后你按照vcf的来算，似乎并不友好（没详细看它的文档），按理说我给它的序列和给dnasp的序列是相同的，结果应该也是相同的，但是实际情况并不是。
   综上所述，如果想计算的是比对后序列的pi，dnasp首选，在这里也是唯一的选择，因为即使你把序列给Bio::PopGen::Statistics算，它算的结果依旧不是你想要的。。。如果你的结果是群体的变异位点信息，那么就选择vcftools把~（当然，不怕麻烦的话也可以把变异位点还原成序列，然后用dnasp来做……）。
   到这里就结束了？？怎么可能，我还没找到批量计算对齐后序列的pi（和dnasp计算的结果相同）方法呢。于是自己根据公式造了小轮子。这里注意到有的对齐序列是有indel的（用“-”字符表示），计算的时候是不把含有indel的位点算进去的。这也就是为什么在设置步长和窗口的时候，dnasp给的结果中，窗口的数值有的不是我们设置的倍数，如下图，设置的窗口时200，最后竟然出现了706，这就是由于501到706中间有6个位点是“-”，这个计算的时候要把它去掉，窗口200再加上6个“-”，所以出现了706。此外，最后一个窗口往往要小于设置的窗口，程序计算的时候也需要注意下。计算脚本用perl完成（只会perl…），支持步长和窗口的设置，计算结果和dnasp一样，但是比其灵活：便于做批量的计算。

   最后，代码就不放到这了，有需要的话可以加群：636121790。