用於計算分子生物學的 Python 工具
用於計算分子生物學的 Python 工具
Genepop 模組允許使用 Python 介面存取 Genepop 功能。這表示現在可以從 Python 存取 Genepop 的絕大多數方法（哈迪-溫伯格平衡的精確檢定、族群分化、基因型不平衡、F 統計、無效等位基因頻率、微衛星的等位基因大小統計等等）。由於程式碼只是一個封裝器，因此需要安裝 Genepop 才能使用。
一個 Genepop 檔案的解析器也可用（並在教學文件中說明）。可以處理 Genepop 格式的檔案，而無需 Genepop。
提供了兩個介面：一個通用、更複雜且更有效率的介面（GenePopController），以及一個簡化、更易於使用、不完整且效率較低的介面（EasyController）。由於其介面特性，EasyController 可能無法處理非常大的檔案。另一方面，它提供了實用函式來計算一些非常簡單的統計資料，例如等位基因計數，這些在通用介面中無法直接使用。
較複雜的介面假設使用者是更熟練的 Python 開發人員（例如，透過使用迭代器），而且目前尚未有相關文件。但即使對於經驗豐富的 Python 開發人員來說，只要所需的功能在 EasyController 中公開且其效能被認為可以接受，EasyController 仍然可以很方便。
有關計算方法詳細資訊，請參閱 Genepop 文件，其中提供了所有推導計算的論文連結。
為了使用控制器，必須在系統中安裝 Genepop，可以從 此處 下載。
在開始之前，讓我們測試安裝（為此您需要一個 Genepop 格式的檔案）
from Bio.PopGen.GenePop.EasyController import EasyController
ctrl = EasyController(your_file_here)
print(ctrl.get_basic_info())
將 your_file_here 替換為您檔案的名稱和路徑。如果您收到 IOError: Genepop not found，則表示 Biopython 找不到您的 Genepop 執行檔。如果 Genepop 不在 PATH 中，您可以將其新增到建構函式行中，即：
ctrl = EasyController(your_file_here, path_to_genepop_here)
如果一切正常，現在我們可以繼續使用 Genepop。對於以下範例，我們將使用 Genepop 檔案 big.gen，該檔案隨單元測試一起提供。我們也假設有一個 ctrl 物件已使用所選的相關檔案初始化。
我們先取得一些基本資訊
pop_names, loci_names = ctrl.get_basic_info()
返回檔案中可用的族群名稱和基因座名稱清單。
注意：大多數現有的 Genepop 檔案都提供關於族群名稱的錯誤資料。在許多情況下，該資訊可能不可信。評估族群資訊大多數時候是根據族群在檔案中的相對位置，而不是名稱。因此，檔案中的第一個族群索引為 0，第二個索引為 1，依此類推…
讓我們取得特定族群和特定等位基因的異質合子性資訊
(exp_homo, obs_homo, exp_hetero, obs_hetero) = ctrl.get_heterozygosity_info(0, "Locus2")
將取得族群 0 和基因座 2 的預期和觀察到的同質合子性和異質合子性（如果您使用另一個檔案，請相應調整族群位置和基因座名稱）。
可以取得特定族群中特定基因座的所有等位基因清單
allele_list = ctrl.get_alleles(0, "Locus2")
allele_list 將為 [3, 20]（即等位基因 3 和 20 在該族群中）。
等位基因的數量只需使用 len(allele_list) 即可取得。
也可以取得所有族群中特定基因座的所有等位基因清單
all_allele_list = ctrl.get_alleles_all_pops("Locus2")
all_allele_list 將為 [3, 20]。
可以取得特定族群中等位基因的頻率
allele_data = ctrl.get_allele_frequency(0, "Locus2")
allele_data 將為 (62, {3: 0.88700000000000001, 20: 0.113})。也就是說，有 62 個基因。88.7% 是 3，11.3% 是 20。
我們可以取得類似的基因型資訊（雙倍體資料）。也會回報預期頻率
genotype_list = ctrl.get_genotype_frequency(0, "Locus2")
genotype_list 將為：[(3, 3, 24, 24.3443), (20, 3, 7, 6.3114999999999997), (20, 20, 0, 0.34429999999999999)]
讓我們解釋第一個元素：有 24 個個體的基因型為 (3, 3)，而預期具有該基因型的個體數為 24.2443。
讓我們從一般多基因座 F 統計開始
Fis, Fst, Fit = ctrl.get_multilocus_f_stats()
這將取得多基因座 Fis、Fst 和 Fit。
讓我們取得每個基因座的統計資料（以及更多）
Fis, Fst, Fit, Qintra, Qinter = ctrl.get_f_stats("Locus2")
這將取得單基因座 Fis、Fst 和 Fit、Qintra 和 Qinter。
下方有 Fst 和 Fis 的特定章節（其中引入了成對變體和族群特定變體）。在 Fis 章節中會解釋 Qintra 和 Qinter。
讓我們取得特定基因座的成對 Fst
pair_fst = ctrl.get_avg_fst_pair_locus("Locus4")
將返回一個對應表，其中鍵是由 population1、population2（族群 ID）組成的對。population2 始終小於 population1。例如：族群 0 和族群 3 之間 Locus4 的成對 Fst 由 pair_fst[(3,0)] 給出。
您也可以取得多基因座成對 Fst 估計值
multilocus_fst = ctrl.get_avg_fst_pair()
這將返回與上述相同的資料結構，但具有多基因座成對 Fst。
我們現在將取得特定基因座/族群的 Fis 以及其他一些統計資料
allele_dict, summary_fis = ctrl.get_fis(0, "Locus2")
讓我們詳細查看 get_fis 的輸出
summary_fis = (62, -0.1111, -0.11269999999999999)
allele_dict = {3: (55, 0.8871, -0.1111), 20: (7, 0.1129, -0.1111)}
summary_fis 保存一個三元組，包含：等位基因總數、Cockerham 和 Weir Fis、Robertson 和 Hill Fis。
allele_dict 保存每個等位基因（每個等位基因都是鍵）的等位基因重複次數、頻率和 Cockerham 和 Weir Fis。
因此，從上述結果可以讀出以下資訊：有 62 個具有 2 個不同等位基因的基因（55 個是 3 型，7 個是 20 型）。3 型的頻率為 0.89，20 型的頻率為 0.11。所有 CW Fis 都是 -0.111，RH Fis 是 -0.112。
現在讓我們取得多基因座 Fis
pop_list = ctrl.get_avg_fis()
pop_list 將為每個族群返回一個元素。每個元素都是一個包含以下內容的四元組
我們可以取得遷徙數量的估計值
samp_size, priv_allele_freq, mig10, mig25, mig50, migcorr = ctrl.estimate_nm()
samp_size 是平均樣本大小，priv_allele_freq 是私有等位基因的平均頻率，mig10 是 Ne=10 的遷徙數量，mig25 是 Ne=25 的遷徙數量，mig50 是 Ne=50 的遷徙數量，migcorr 是在校正預期大小後的遷徙數量。
檢定通常是計算密集型的，因為它們通常基於馬可夫鏈演算法。在某些情況下，可以使用完整列舉方法，但這些方法只能應用於等位基因數量非常少的基因座。這表示大多數檢定將需要相當長的時間才能完成。
有關以下馬可夫鏈參數（預熱、分批和迭代）的更多詳細資訊，請參閱 Genepop 手冊。也請參閱手冊以了解何時適用完整列舉。
讓我們開始測試每個族群中每個基因座的哈迪-溫伯格平衡
loci_map = ctrl.test_hw_pop(1, "excess")
第二個參數可以是 probability、excess 或 deficiency。probability 是標準的 Haldane HW 檢定。如果您對異質合子缺乏感興趣，請使用 deficiency；如果您對過剩感興趣，請使用 excess。
輸出是一個對應表，其中鍵是基因座名稱。內容是一個元組，包含 P 值、標準誤差、Fis (Weir 和 Cockerham)、Fis (Robertson 和 Hill) 以及步驟。
pop_test, loc_test, all_test = ctrl.test_hw_global("deficiency")
如果您對異質合子缺乏感興趣，請使用 deficiency；如果您對過剩感興趣，請使用 excess。probability 在這裡不適用，就像在 test_hw_pop 中一樣。
輸出是一個三元組
我們可以利用對數似然比統計量 (G 檢定) 來檢定 2 個基因座是否處於連鎖不平衡狀態。
chi2, df, pval = ctrl.test_ld_all_pair(
"Locus1", "Locus2", dememorization=1000, batches=10, iterations=100
)
返回 G 統計的卡方值、自由度和 P 值。
距離隔離 (IBD) 分析需要 Genepop 檔案的特殊形式
範例
...
Pop
0 15, 0201 0303 0102 0302 1011
Pop
0 30, 0202 0301 0102 0303 1111
Pop
0 45, 0102 0401 0202 0102 1010
Pop
0 60, 0103 0202 0101 0202 1011
Pop
0 75, 0203 0204 0101 0102 1010
POP
15 15, 0102 0202 0201 0405 0807
...
請注意，我們正在使用的範例檔案不能用於此情況。
只有單次呼叫 IBD 分析
estimate, distance, (a, b), (bb, bblow, bbhigh) = ctrl.calc_ibd(
self, is_diplo=True, stat="a", scale="Log", min_dist=0.00001
)
is_diplo 指定資料是雙倍體 (True) 還是單倍體 (False)。stat 是 'a' 或 'e'（詳情請參閱 Genepop 手冊）。scale 是 'Log' 或 'Linear'。'Log' 用於 2D 座標，而 'Linear' 用於 1D。只有高於 min_dist 的成對比較才會用於計算回歸係數。
該方法返回
三角矩陣的解釋方式應該如下：Python 方式，矩陣是使用數字列表的列表實現的，如下所示
[[0.1], [0.2, 0.3], [0.4, 0.5, 0.6]]
上述資料結構對應於以下三角矩陣
1 2 3
2 0.1
3 0.2 0.3
4 0.4 0.5 0.6