利用Python读取fasta文件并进行一系列操作（三）

利用Python读取fasta文件并进行一系列操作（三）

---

概述：

---

本节目标：计算智人与猩猩ABO基因的相对熵
语言：python3.8
模块：pysam scipy
整体思路：先计算出序列中AG CT AC AT GC GT六种组合序列所占比，再计算相对熵

步骤：

---

1. 利用pysam模块分别读取智人ABO基因所有序列和猩猩ABO基因所有序列

import pysam as sam
from scipy import stats
hfasta = sam.FastaFile('homo_ABO.fasta')
cfasta = sam.FastaFile('chimp_ABO.fasta')
c_all = cfasta.fetch('NC_036888.1:c105542142-105516107', 0, 26036)
h_all = hfasta.fetch('NG_006669.2', 0, 42144)

2. 利用split函数切序列并用len函数返回结果-1，分别得出六种组合序列数量，再求占比

cAG = len(c_all.split('AG')) - 1
cCT = len(c_all.split('CT')) - 1
cAC = len(c_all.split('AC')) - 1
cAT = len(c_all.split('AT')) - 1
cGC = len(c_all.split('GC')) - 1
cCT = len(c_all.split('CT')) - 1
p = cCT/sum([cAG,cCT,cAC,cAT,cGC,cCT])
h_all = hfasta.fetch('NG_006669.2', 0, 42144)
hAG = len(h_all.split('AG')) - 1
hCT = len(h_all.split('CT')) - 1
hAC = len(h_all.split('AC')) - 1
hAT = len(h_all.split('AT')) - 1
hGC = len(h_all.split('GC')) - 1
hCT = len(h_all.split('CT')) - 1
q = hCT/sum([hAG,hCT,hAC,hAT,hGC,hCT])

3. 利用scipy模块的entropy函数求相对熵

result = stats.entropy([p,q])
print(result)

结果：

---

0.6930868205144018