我正在尝试根据 SNP 数据(~4000SNPs)计算 tajD。我有 fasta 和 vcf 文件,最初尝试使用我的 fasta 文件。在帮助下,我明白了我有一个可以从文件中采样的脚本,但是我需要将 fasta 文件拆分为不同的群体。我很害怕这样做(尽管我必须给出我现在遇到的错误),所以我分阶段使用我的 vcf 文件,并希望使用它来代替。我正在使用 R 包 pegas,并收到以下错误。
> rm(list=ls(all=T))
> library("ape")
> library("ade4")
> library("adegenet")
> library("pegas")
> b8c18FromVcf <- read.vcf("b8c18_2phased.vcf")
文件显然尚未访问:扫描文件 b8c18_2phased.vcf 3.194102 / 3.194102 Mb Done。读取 4074 / 4074 位点。完毕。
获得单倍型
> b8c18haplos <- haplotype(b8c18FromVcf)
分析个人编号 186 / 186
来自 haplos 的 tajD
> tajd <- tajima.test(b8c18haplos)
警告信息:在 tajima.test(b8c18haplos) 中:Tajima 测试需要至少 4 个序列
我将在此处附加指向我的分阶段和非分阶段文件的链接。 https://drive.google.com/open?id=0B6qb8IlaQGFZTmQ1YXRVbnFSRzA https://drive.google.com/open?id=0B6qb8IlaQGFZQm9HZjZSUkE3NEU
有什么想法吗?
最后,我想知道是否有办法在 tajD 命令中对种群进行子集化。我在同一个 vcf 文件中有 7 个人口,我应该分别计算每个人口的 tajD。如果没有,什么是子集 vcfs 的最佳工具。我已经对此进行了大量的谷歌搜索,但似乎没有一个是直截了当的。
带着感谢,