建站
咨詢
售后
建站咨詢
新聞中心
Python DNA序列在使用的時(shí)候有很多需要我們注意的東西，其實(shí)在不斷的學(xué)習(xí)中有很多問題存在，下面我們就詳細(xì)的看看如何進(jìn)行相關(guān)的技術(shù)學(xué)校。ms是我?guī)煹艿膔otation project：給定一堆Python DNA序列，即由字符A, C, G, T組成的字符串，統(tǒng)計(jì)所有長度為n的子序列出現(xiàn)的頻率。
站在用戶的角度思考問題，與客戶深入溝通，找到蕉嶺網(wǎng)站設(shè)計(jì)與蕉嶺網(wǎng)站推廣的解決方案，憑借多年的經(jīng)驗(yàn)，讓設(shè)計(jì)與互聯(lián)網(wǎng)技術(shù)結(jié)合，創(chuàng)造個(gè)性化、用戶體驗(yàn)好的作品，建站類型包括：網(wǎng)站建設(shè)、網(wǎng)站設(shè)計(jì)、企業(yè)官網(wǎng)、英文網(wǎng)站、手機(jī)端網(wǎng)站、網(wǎng)站推廣、主機(jī)域名、虛擬主機(jī)、企業(yè)郵箱。業(yè)務(wù)覆蓋蕉嶺地區(qū)。
比如 ACGTACGT，子序列長度為2，于是 AC=2, CG=2, GT=2, TA=1，其余長度為2的子序列頻率為0.
最先想到的就是建一個(gè)字典，key是所有可能的子序列，value是這個(gè)子序列出現(xiàn)的頻率。Python DNA序列但是當(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è)字典，表示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')
雖然我們內(nèi)部把子序列當(dāng)成正整數(shù)來存儲（確切地說，其實(shí)這個(gè)整數(shù)是沒有存在內(nèi)存里的，而是由其在list的index表示的），為了方便生物學(xué)家們看，輸出時(shí)還是轉(zhuǎn)換回子序列比較好。
于是有了下面的整數(shù)映射成子序列函數(shù)，其中調(diào)用了另外一個(gè)函數(shù)baseN()，來源在此，感謝作者～
- def baseN(n,b):
- '''convert non-negative decimal integer n to
- equivalent in another base b (2-36)'''
- return ((n == 0) and '0' ) or ( baseN(n // b, b).lstrip('0') + \
- "0123456789abcdefghijklmnopqrstuvwxyz"[n % b])
- def int2motif(n, motifLen):
- '''convert non-negative integer n to a sub-sequence/motif with length motifLen'''
- intBase4 = baseN(n,4)
- return ''.join(map(lambda x: i2mD[int(x)],'0'*(motifLen-len(intBase4))+intBase4))
- Test:
- >>> int2motif(27,4)
- 'ACGT'
以下代碼從命令行讀入一個(gè)存有DNA序列的fasta文件，以及子序列長度，并輸出子序列和頻率。注意以下代碼需要Biopython module。
- if __name__ == '__main__':
- import sys
- from Bio import SeqIO
- # read in the fasta file name and motif length
- # from command line parameters
- fastafile = sys.argv[1]
- motifLen = int(sys.argv[2])
- # list to store subsequence frequency
- frequencyL = [0]*4**motifLen
- # go over each DNA sequence in the fasta file
- # and count the frequency of subsequences
- it = SeqIO.parse(open(fastafile),'fasta')
- for rec in it:
- chrom = rec.seq.tostring()
- for i in range(len(chrom)-motifLen+1):
- motif = chrom[i:i+motifLen]
- frequencyL[motif2int(motif)] += 1
- # print frequency result to screen
- for i, frequency in enumerate(frequencyL):
- print int2motif(i, motifLen), frequency
以上就是創(chuàng)新互聯(lián)對Python DNA序列的相關(guān)介紹。
分享名稱：PythonDNA序列中子序列出現(xiàn)頻率
轉(zhuǎn)載注明：http://www.5511xx.com/article/djosgcc.html
