Python統(tǒng)計(jì)如何進(jìn)行DNA序列整理
Python統(tǒng)計(jì)在我們的使用中有很多的障礙，其中在DNA序列上的相關(guān)問題就需要我們不斷的去學(xué)習(xí)。下面我們就向大家介紹有關(guān)的問題，希望在以后Python統(tǒng)計(jì)的使用過程中有所收獲。
給定一堆DNA序列，即由字符A, C, G, T組成的字符串，統(tǒng)計(jì)所有長度為n的子序列出現(xiàn)的頻率。比如 ACGTACGT，子序列長度為2，于是 AC=2, CG=2, GT=2, TA=1，其余長度為2的子序列頻率為0.
***想到的就是建一個(gè)字典，key是所有可能的子序列，value是這個(gè)子序列出現(xiàn)的頻率。但是當(dāng)子序列比較長的時(shí)候，比如 n=8，需要一個(gè)有65536 (4的8次方) 個(gè)key-value pair的字典，且每個(gè)key的長度是8字符。這樣ms有點(diǎn)浪費(fèi)內(nèi)存。。
于是想到，所有的長度為n的子序列是有序且連續(xù)的，所以可以映射到一個(gè)長度為4的n次方的的list里。令 A=0, C=1, G=2, T=3，則把子序列 ACGT 轉(zhuǎn)換成 0*4^3 + 1*4^2 + 2*4 + 3 = 27, 映射到list的第27位。如此，list的index對應(yīng)子序列，而list這個(gè)index位置則儲存這個(gè)子序列出現(xiàn)的頻率。
于是我們先要建立2個(gè)字典，Python統(tǒng)計(jì)表示ACGT和0123一一對應(yīng)的關(guān)系：
- i2mD = {0:'A', 1:'C', 2:'G', 3:'T'}
- m2iD = dict(A=0,C=1,G=2,T=3)
- # This is just another way to initialize a dictionary
以及下面的子序列映射成整數(shù)函數(shù)：
- def motif2int(motif):
- '''convert a sub-sequence/motif to a non-negative integer'''
- total = 0
- for i, letter in enumerate(motif):
- total += m2iD[letter]*4**(len(motif)-i-1)
- return total
- Test:
- >>> motif2int('ACGT')
- 27
以上就是對Python統(tǒng)計(jì)的相關(guān)介紹。雖然我們內(nèi)部把子序列當(dāng)成正整數(shù)來存儲（確切地說，其實(shí)這個(gè)整數(shù)是沒有存在內(nèi)存里的，而是由其在list的index表示的），為了方便生物學(xué)家們看，輸出時(shí)還是轉(zhuǎn)換回子序列比較好。
【編輯推薦】
