Genepop格式可用在许多保护遗传学研究中。它是Genepop应用程序的格式,是许多群体遗传学分析的事实格式。如果读者来自其他领域(例如,有大量的测序经验),你可能没有听说过这种格式,但是这种格式确被广泛使用(正如它的引用记录证明的),值得一看。这里将把一些以前的秘笈中数据集转换成这种格式,并介绍来自Biopython的Genepop的解析器。
准备好
读者需要运行前一个秘笈,因为它的输出对这个秘笈是必需的。
如果读者没有用Docker,可能你无法使用前面产生的部分代码(大多数是处理数据转换的基本食粮),你可以在https://github.com/tiagoantao/pygenomics找到它们的代码,并用下列代码下载:
pip install pygenomics
注意在此阶段,我们还没有用到Genepop应用程序(在下一个秘笈将会改变),因此无需现在安装它。
像往常一样,有一个可用的notebook在 03_PopGen/Genepop_Format.ipynb,但是它仍需要读者运行前一个notebook来产生必需的文件。
如何做
参看下列步骤:
1,这里先加载元数据如下(我们将使用前一个秘笈的一个精简版本):
from collections import defaultdict
f = open(relationships_w_pops_121708.txt)
pop_ind = defaultdict(list)
f.readline() # header
for line in f:
toks = line.rstrip().split(\t)
fam_id = toks[0]
ind_id = toks[1]
pop = toks[-1]
pop_ind[pop].append((fam_id, ind_id))
f.close()
2,现在检查PLINK数据文件和元数据的一致性,因为我们需要清理群体的映射文件来产生一个Genepop文件,代码如下:
all_inds = []
for inds in pop_ind.values():
all_inds.extend(inds)
for line in open(hapmap1.ped):
toks = line.rstrip().replace( , \t).split(\t)
fam = toks[0]
ind = toks[1]
if (fam, ind) not in all_inds:
print(Problems with %s/%s % (fam, ind))
前面的代码生成来自元数据文件的所有个体的一个列表。然后,它会打开文件hapmap1.ped,取样1%的谱系信息,(这里选择了1%,因为取1%的过程比10%或100%个样品快得多;我们只需要谱系信息,而没有遗传信息)然后比较信息对。最后将报告所有的在PED文件中,但不在元数据文件中的个体。我们进行对每个PED的行进行替换过程,因为读者可以发现一些PLINK文件是空白符分离,而另一些制表符分隔的。在一个完美的情况下,它什么也不会输出。但这里有一个不正确的数据项。该数据项(其中有家庭标识为2469和个人身份标识为na20281)不具有与家庭身份报道一致的元数据。
小技巧:
对读者自己的数据,一定要彻底对数据和元数据检查一致性问题。所有数据源(如果你有一个以上的),要确保它们之间的一致。如果没有,至少对所有有问题的情况,更好的是采取行动去了解和纠正任何潜在的问题。默认的假设应该有问题(不是一切的都没问题)。虽然读者可能是自己产生数据,但是检查一下更好。疏漏和错误是肯定会发生的,犯错误是正常的,基本的方式是检查他们。过分自信是缺乏经验的表现。
3,这里从PLINK转换一些数据集到Genepop格式:
from genomics.popgen.plink.convert import to_genepop
to_genepop(hapmap1_auto, hapmap1_auto, pop_ind)
to_genepop(hapmap10, hapmap10, pop_ind)
to_genepop(hapmap10_auto, hapmap10_auto, pop_ind)
to_genepop(hapmap10_auto_noofs_ld, hapmap10_auto_noofs_ld, pop_ind)
to_genepop(hapmap10_auto_noofs_2, hapmap10_auto_noofs_2, pop_ind)
注意这里对所有前缀相同的文件进行处理,特别是第一个对输入文件的传递:对该前缀添加“.ped”和“.map”来添加输入文件,第二个前缀添加“.gp”来产生Genepop文件,而“.pops”将包含Genepop文件的群体的顺序。参看一下这两个产生的文件,可以用来熟悉它们的内容,虽然我们将在下个秘笈中更详尽地刨析它们。作为PLINK格式,没有群体结构的信息,这里需要把pop_ind的字典传递过去。该字典将被用于创建一个具有群体结构的Genepop文件。这里用我的包中提供的一个函数来转换PLINK到Genepop数据,它需要一些时间来运行。注意这里只是转换二次取样的数据,这是为了下游的分析的许多计算过程更有效率,但是要注意在许多读者自己的分析中需要完整的数据集。这个函数将忽略没有群体的个体 ,而这意味着它将排除那个在一致性步骤中的错误的标记了家族的个体。A将转换成1,C到2,T到3,而G到4。虽然一个pops文件将给输出文件创建群体的顺序,它将总是按照词典顺序排序。如果你对这个函数如何工作感到好奇,欢迎参看这里
https://github.com/tiagoantao/pygenomics/blob/master/genomics/popgen/plink/convert.py。提前告知它只是包含文字处理过程。
4,Biopython提供了一个对Genepop文件的在内存中的解析器,这里尝试打开一个1%取样的常染色体的文件:
from Bio.PopGen.GenePop import read
rec = read(open(hapmap1_auto.gp))
print(Number of loci %d % len(rec.loci_list))
print(Number of populations %d % len(rec.pop_list))
print(Population names: %s % , .join(rec.pop_list))
print(Individuals per population %s % , .join([str(len(inds)) for inds in rec.populations]))
ind = rec.populations[1][0]
print(Individual %s, SNP %s, alleles: %d %d % (ind[0], rec.loci_list[0], ind[1][0][0], ind[1][0][1]))
del rec
输出如下:
关于Genepop人群名称默认的假设就是,最后一个个体是用来确定一个群体。这是稍微有些非正规,我们也会产生一个包含人群名称的".pop"文件(像以前的秘笈一样)。由于以前的秘笈中标记的采样过程是随机的,读者可能会看到一个稍微不同的基因座的数目(loci)。由于整个数据集在内存中,这里可以直接访问任何人群的任何个体。这里进行打印的是最后一行,即我们访问第二人群的第一个人和打印他的名字,还有第一个SNP等位基因(分别为3和2,因此编码为T和C)。第一个SNP被称为1/rs2710888/949705。在这里,1个代表的染色体序号,中间的ID是SNP的rs标识(在NCBI的dbSNP数据库的标识符),和最后一个数是对人类Build 36染色体的位置。最后,我们删除记录,因为它占用了大量的内存。
小技巧
请注意,这些输出取决于Genepop是如何编码的(这在我的to_genepop函数中)。还取决于例如,编码1234是任意的(ACTG只是便利)或事实上的人群的词典顺序或基因座的名称包括rs标识和染色体的位置。如果读者从另一个地方得到你的文件,就必须检查他们使用任何约定(对你这可能会是方便的或不方便的)。如果读者生成你自己的文件,一定要使用一些约定,将对下游的分析有用(像这里一样)。当然,这里的说法是一般性的;你可以将它应用到其他格式的文件中,只要他们有任何形式的灵活性。
5,更符合实际情况的是,我们将对现代的数据使用大文件解析器,因为它不会在内存中加载全部文件,反而提供一个迭代器代替列表,代码如下:
from Bio.PopGen.GenePop.LargeFileParser import read as read_large
def count_individuals(fname):
rec = read_large(open(fname))
pop_sizes = []
for line in rec.data_generator():
if line == ():
pop_sizes.append(0)
else:
pop_sizes[-1] += 1
return pop_sizes
print(Individuals per population %s % , .join([str(ninds) for ninds in count_individuals(hapmap1_auto.gp)]))
print(len(read_large(open(hapmap10.gp)).loci_list))
print(len(read_large(open(hapmap10_auto.gp)).loci_list))
print(len(read_large(open(hapmap10_auto_noofs_ld.gp)).loci_list))
这个count_individuals函数显示了如何对一个Genepop文件用大文件解析器进行遍历;当迭代它的时候,如果找到一个空的元组,它是一个新人群的标志。其他任何东西都是一个个体,由元组(一对)构成,包括个体的名称和一系列基因座(这里我们没有读取)。如前所述,每个人群的个体将返回完全相同的值。然后这里打印输出三个不同文件的基因座的数目:10%抽样,10%抽样,只包含常染色体,以及10%抽样的LD清洗的常染色体。输出结果(由于随机产生文件,会稍有不同)显示出第一个文件有更多的标记,第二个稍少些(因为第二个是第一个文件的子集,去除了性染色体和线粒体)和最后一个更少,因为它是LD清洗后的第二个文件的子集。
附加信息
事实上在网络上有一个Genepop的接口,读者可以使用手动的例子(特别是对于小规模的文件):http://genepop.curtin.edu.au/