Genepop格式可用在許多保護(hù)遺傳學(xué)研究中。它是Genepop應(yīng)用程序的格式，是許多群體遺傳學(xué)分析的事實(shí)格式。如果讀者來自其他領(lǐng)域（例如，有大量的測序經(jīng)驗(yàn)），你可能沒有聽說過這種格式，但是這種格式確被廣泛使用（正如它的引用記錄證明的），值得一看。這里將把一些以前的秘笈中數(shù)據(jù)集轉(zhuǎn)換成這種格式，并介紹來自Biopython的Genepop的解析器。
準(zhǔn)備好
讀者需要運(yùn)行前一個(gè)秘笈，因?yàn)樗妮敵鰧@個(gè)秘笈是必需的。
如果讀者沒有用Docker，可能你無法使用前面產(chǎn)生的部分代碼（大多數(shù)是處理數(shù)據(jù)轉(zhuǎn)換的基本食糧），你可以在https://github.com/tiagoantao/pygenomics找到它們的代碼，并用下列代碼下載：
pip install pygenomics
注意在此階段，我們還沒有用到Genepop應(yīng)用程序（在下一個(gè)秘笈將會(huì)改變），因此無需現(xiàn)在安裝它。
像往常一樣，有一個(gè)可用的notebook在 03_PopGen/Genepop_Format.ipynb，但是它仍需要讀者運(yùn)行前一個(gè)notebook來產(chǎn)生必需的文件。
如何做
參看下列步驟：
1，這里先加載元數(shù)據(jù)如下（我們將使用前一個(gè)秘笈的一個(gè)精簡版本）：
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，現(xiàn)在檢查PLINK數(shù)據(jù)文件和元數(shù)據(jù)的一致性，因?yàn)槲覀冃枰謇砣后w的映射文件來產(chǎn)生一個(gè)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))
前面的代碼生成來自元數(shù)據(jù)文件的所有個(gè)體的一個(gè)列表。然后，它會(huì)打開文件hapmap1.ped，取樣1%的譜系信息，（這里選擇了1%，因?yàn)槿?%的過程比10%或100%個(gè)樣品快得多；我們只需要譜系信息，而沒有遺傳信息）然后比較信息對。最后將報(bào)告所有的在PED文件中，但不在元數(shù)據(jù)文件中的個(gè)體。我們進(jìn)行對每個(gè)PED的行進(jìn)行替換過程，因?yàn)樽x者可以發(fā)現(xiàn)一些PLINK文件是空白符分離，而另一些制表符分隔的。在一個(gè)完美的情況下，它什么也不會(huì)輸出。但這里有一個(gè)不正確的數(shù)據(jù)項(xiàng)。該數(shù)據(jù)項(xiàng)（其中有家庭標(biāo)識(shí)為2469和個(gè)人身份標(biāo)識(shí)為na20281）不具有與家庭身份報(bào)道一致的元數(shù)據(jù)。
小技巧：
對讀者自己的數(shù)據(jù)，一定要徹底對數(shù)據(jù)和元數(shù)據(jù)檢查一致性問題。所有數(shù)據(jù)源（如果你有一個(gè)以上的），要確保它們之間的一致。如果沒有，至少對所有有問題的情況，更好的是采取行動(dòng)去了解和糾正任何潛在的問題。默認(rèn)的假設(shè)應(yīng)該有問題（不是一切的都沒問題）。雖然讀者可能是自己產(chǎn)生數(shù)據(jù)，但是檢查一下更好。疏漏和錯(cuò)誤是肯定會(huì)發(fā)生的，犯錯(cuò)誤是正常的，基本的方式是檢查他們。過分自信是缺乏經(jīng)驗(yàn)的表現(xiàn)。
3，這里從PLINK轉(zhuǎn)換一些數(shù)據(jù)集到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)
注意這里對所有前綴相同的文件進(jìn)行處理，特別是第一個(gè)對輸入文件的傳遞：對該前綴添加“.ped”和“.map”來添加輸入文件，第二個(gè)前綴添加“.gp”來產(chǎn)生Genepop文件，而“.pops”將包含Genepop文件的群體的順序。參看一下這兩個(gè)產(chǎn)生的文件，可以用來熟悉它們的內(nèi)容，雖然我們將在下個(gè)秘笈中更詳盡地刨析它們。作為PLINK格式，沒有群體結(jié)構(gòu)的信息，這里需要把pop_ind的字典傳遞過去。該字典將被用于創(chuàng)建一個(gè)具有群體結(jié)構(gòu)的Genepop文件。這里用我的包中提供的一個(gè)函數(shù)來轉(zhuǎn)換PLINK到Genepop數(shù)據(jù)，它需要一些時(shí)間來運(yùn)行。注意這里只是轉(zhuǎn)換二次取樣的數(shù)據(jù)，這是為了下游的分析的許多計(jì)算過程更有效率，但是要注意在許多讀者自己的分析中需要完整的數(shù)據(jù)集。這個(gè)函數(shù)將忽略沒有群體的個(gè)體 ，而這意味著它將排除那個(gè)在一致性步驟中的錯(cuò)誤的標(biāo)記了家族的個(gè)體。A將轉(zhuǎn)換成1，C到2，T到3，而G到4。雖然一個(gè)pops文件將給輸出文件創(chuàng)建群體的順序，它將總是按照詞典順序排序。如果你對這個(gè)函數(shù)如何工作感到好奇，歡迎參看這里
https://github.com/tiagoantao/pygenomics/blob/master/genomics/popgen/plink/convert.py。提前告知它只是包含文字處理過程。
4，Biopython提供了一個(gè)對Genepop文件的在內(nèi)存中的解析器，這里嘗試打開一個(gè)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
輸出如下：
關(guān)于Genepop人群名稱默認(rèn)的假設(shè)就是，最后一個(gè)個(gè)體是用來確定一個(gè)群體。這是稍微有些非正規(guī)，我們也會(huì)產(chǎn)生一個(gè)包含人群名稱的".pop"文件（像以前的秘笈一樣）。由于以前的秘笈中標(biāo)記的采樣過程是隨機(jī)的，讀者可能會(huì)看到一個(gè)稍微不同的基因座的數(shù)目（loci）。由于整個(gè)數(shù)據(jù)集在內(nèi)存中，這里可以直接訪問任何人群的任何個(gè)體。這里進(jìn)行打印的是最后一行，即我們訪問第二人群的第一個(gè)人和打印他的名字，還有第一個(gè)SNP等位基因（分別為3和2，因此編碼為T和C）。第一個(gè)SNP被稱為1/rs2710888/949705。在這里，1個(gè)代表的染色體序號，中間的ID是SNP的rs標(biāo)識(shí)（在NCBI的dbSNP數(shù)據(jù)庫的標(biāo)識(shí)符），和最后一個(gè)數(shù)是對人類Build 36染色體的位置。最后，我們刪除記錄，因?yàn)樗加昧舜罅康膬?nèi)存。
小技巧
請注意，這些輸出取決于Genepop是如何編碼的（這在我的to_genepop函數(shù)中）。還取決于例如，編碼1234是任意的（ACTG只是便利）或事實(shí)上的人群的詞典順序或基因座的名稱包括rs標(biāo)識(shí)和染色體的位置。如果讀者從另一個(gè)地方得到你的文件，就必須檢查他們使用任何約定（對你這可能會(huì)是方便的或不方便的）。如果讀者生成你自己的文件，一定要使用一些約定，將對下游的分析有用（像這里一樣）。當(dāng)然，這里的說法是一般性的；你可以將它應(yīng)用到其他格式的文件中，只要他們有任何形式的靈活性。
5，更符合實(shí)際情況的是，我們將對現(xiàn)代的數(shù)據(jù)使用大文件解析器，因?yàn)樗粫?huì)在內(nèi)存中加載全部文件，反而提供一個(gè)迭代器代替列表，代碼如下：
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))
這個(gè)count_individuals函數(shù)顯示了如何對一個(gè)Genepop文件用大文件解析器進(jìn)行遍歷；當(dāng)?shù)臅r(shí)候，如果找到一個(gè)空的元組，它是一個(gè)新人群的標(biāo)志。其他任何東西都是一個(gè)個(gè)體，由元組（一對）構(gòu)成，包括個(gè)體的名稱和一系列基因座（這里我們沒有讀?。Ｈ缜八?，每個(gè)人群的個(gè)體將返回完全相同的值。然后這里打印輸出三個(gè)不同文件的基因座的數(shù)目：10%抽樣，10%抽樣，只包含常染色體，以及10%抽樣的LD清洗的常染色體。輸出結(jié)果（由于隨機(jī)產(chǎn)生文件，會(huì)稍有不同）顯示出第一個(gè)文件有更多的標(biāo)記，第二個(gè)稍少些（因?yàn)榈诙€(gè)是第一個(gè)文件的子集，去除了性染色體和線粒體）和最后一個(gè)更少，因?yàn)樗荓D清洗后的第二個(gè)文件的子集。
附加信息
事實(shí)上在網(wǎng)絡(luò)上有一個(gè)Genepop的接口，讀者可以使用手動(dòng)的例子（特別是對于小規(guī)模的文件）：http://genepop.curtin.edu.au/
