1

在 Python 中，我使用EggLib。我正在尝试计算在 VCF 文件中找到的每个 SNP 的 Jost 的 D 值。

数据

此处的数据采用 VCF 格式。数据集很小，有 2 个种群，每个种群 100 个人和 6 个 SNP（都在 1 号染色体上）。

每个个体被命名为Pp.Ii，其中p是它所属的人口指数，i是个体指数。

代码

我的困难在于人口结构的规范。这是我的审判

### Read the vcf file ###
vcf = egglib.io.VcfParser("MyData.vcf") 

### Create the `Structure` object ###
# Dictionary for a given cluster. There is only one cluster.
dcluster = {}            
# Loop through each population 
for popIndex in [0,1]:  
    # dictionnary for a given population. There are two populations
    dpop = {}            
    # Loop through each individual
    for IndIndex in range(popIndex * 100,(popIndex + 1) * 100):     
            # A single list to define an individual
        dpop[IndIndex] = [IndIndex*2, IndIndex*2 + 1]
    dcluster[popIndex] = dpop

struct = {0: dcluster}

### Define the population structure ###
Structure = egglib.stats.make_structure(struct, None) 

### Configurate the 'ComputeStats' object ###
cs = egglib.stats.ComputeStats()
cs.configure(only_diallelic=False)
cs.add_stats('Dj') # Jost's D

### Isolate a SNP ###
vcf.next()
site = egglib.stats.site_from_vcf(vcf)

### Calculate Jost's D ###
cs.process_site(site, struct=Structure)

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Library/Python/2.7/site-packages/egglib/stats/_cstats.py", line 431, in process_site
    self._frq.process_site(site, struct=struct)
  File "/Library/Python/2.7/site-packages/egglib/stats/_freq.py", line 159, in process_site
    if sum(struct) != site._obj.get_ning(): raise ValueError, 'invalid structure (sample size is required to match)'
ValueError: invalid structure (sample size is required to match)

该文档表明here

[结构对象] 是一个包含两个项目的元组，每个项目都是一个字典。第一个代表内群，第二个代表外群。

内组词典本身就是一本包含更多词典的词典，每个群体都有一个词典。每个集群字典都是人口字典，人口本身由字典表示。人口字典再次是个人字典。幸运的是，个人由列表表示。

个体列表包含属于该个体的所有样本的索引。对于单倍体数据，个体将是单项列表。在其他情况下，所有单独的列表都必须具有相同数量的项目（一致的倍性）。请注意，如果倍性多于一个，则不会强制将给定个体的样本分组在原始数据中。

ingroup 字典的键是标识每个集群的标签。在集群字典中，键是人口标签。最后，在人口字典中，键是单独的标签。

第二个字典代表外群。它的结构更简单：它有单独的标签作为键，对应的样本索引列表作为值。外群字典类似于任何内群人口字典。倍性需要匹配所有内群和外群个体。

但我无法理解它。提供的示例适用于 fasta 格式，我不明白将逻辑扩展到 VCF 格式。

4

1 回答 1

0

有两个错误

第一个错误

该函数make_structure返回 Structure 对象，但不将其保存在stats. 因此，您必须保存此输出并在函数中使用它process_site。

Structure = egglib.stats.make_structure(struct, None)

第二个错误

Structure 对象必须指定单倍体。因此，将字典创建为

dcluster = {}            
for popIndex in [0,1]:  
    dpop = {}            
    for IndIndex in range(popIndex * 100,(popIndex + 1) * 100):     
        dpop[IndIndex] = [IndIndex]
    dcluster[popIndex] = dpop

struct = {0: dcluster}
于 2017-03-15T22:25:25.650 回答