2015-01-09 38 views
3

我想用Python中的位編碼編碼核苷酸「A」,「G」,「C」和「T」。例如:Python:使用位。用零和1編碼核苷酸

'A' = 00 
'G' = 01 
'C' = 10 
'T' = 11 

爲了構建一個包含k聚體巨大的字典,如:

dic = { 'ATGACTGACT':231, 'AAATGACGGAC':500 ... } 

我認爲這可能減少,因爲「ATGC」由字典所需要的內存量將需要4字節,但同一個字將需要8位和一位編碼。

我不知道這是可以做到,如果是這樣,它使用Python

在此先感謝我怎麼會呢!

編輯:對不起,我沒有解釋自己的權利。

我想要的是步行穿過由ATGC組成的序列,其大小爲k的滑動窗口,並計算每個k聚體在該序列中出現的次數。例如:

'ATGAATGAA' # with a sliding window of 5 would be 
dic = { 'ATGAA':2, 'TGAAT':1, 'GAATG':1, 'AATGA':1, } 

因爲我想與「聯合運輸線協定」的所有可能組合來構建字典大小k開始前閱讀的順序,以訪問字典與每個k-mer的關鍵並總結1的價值,我想知道是否有可能k-mers在該字典存儲位編碼。或多或少:

dic = {1011001010: 3, 0000110011: 666, ... etc } 

目前我正在用itertools構建該詞典。

# k-mers of size 8 
{''.join(x):0 for x in itertools.product('ATGC', repeat=8)} 

我想另一個問題將是每個k-mer的需要,以訪問字典

+1

現在還不清楚:在你的dic例子中,你仍然有完整的k-mer版本,那麼在哪裏減少內存?另一個重要的問題是,在這之後你打算怎麼做?部分搜索?使用指標? – pmod 2015-01-09 22:03:48

+6

'00001100'會編碼什麼? AAT? TA? AATA? – 2015-01-09 22:08:43

+2

對於可變長度字符串,您可以使用像[bitstring](https://pythonhosted.org/bitstring/index.html)(特別是ConstBitStream用作字典鍵)的模塊。我不知道這將如何影響性能/內存使用。順便說一句,你可以[過早優化](http://ubiquity.acm.org/article.cfm?id=1513451)?我不是說你是;問題是,恆定的4倍內存改進是否值得額外的複雜性,並且現在應該做出改變? – Lack 2015-01-09 22:31:19

回答

1

被轉換爲位編碼這不正是你問什麼

In [11]: d={'A':'00','G':'01','C':'10','T':'11'} 

In [12]: int('0B'+''.join([d[c] for c in 'ATGACTGACT']),2) 
Out[12]: 215883 

In [13]: int('0B'+''.join([d[c] for c in 'ATGACTGACT'[::-1]]),2) 
Out[13]: 925212 

In [14]: 

但他們的意見中提出的異議通過pmodIgnacio Vazquez-Abrams是非常重要的,我認爲你應該認真考慮你的方法。

4

您可以將您的kmers轉換爲二進制文件,但是Ignacio指出您仍然需要知道它們的長度,因此您可能還需要將其存儲。所以,對於很長的序列,這仍然會節省存儲空間。

下面是一些示例代碼,這需要序列,它們編碼並再次對其進行解碼:

encoding_map = {'A': 0, 'G': 1, 'C': 2, 'T': 3} 
decoding_lst = ['A', 'G', 'C', 'T'] 


def encode(k): 
    code = 0 
    for ch in k: 
     code *= 4 
     code += encoding_map[ch] 
    return code, len(k) 


def decode(enc): 
    code, length = enc 
    ret = '' 
    for _ in range(length): 
     index = code & 3 
     code >>= 2 
     ret = decoding_lst[index] + ret 
    return ret 


kmers = ['ATGACTGACT', 'ATGC', 'AATGC'] 

kmerdict = {k: encode(k) for k in kmers} 

print(kmerdict) 

for key, enc in kmerdict.items(): 
    print(enc, decode(enc)) 

典型輸出:

{'AATGC': (54, 5), 'ATGC': (54, 4), 'ATGACTGACT': (215883, 10)} 
(54, 5) AATGC 
(54, 4) ATGC 
(215883, 10) ATGACTGACT 

順便說一句,它不應該不管多久了序列是,Python應該能夠處理編碼和解碼,因爲整數擴展到足以保存數字的位。

+0

嗯,適合我。 – quamrana 2015-01-09 22:35:02

+0

對不起,你是對的。我在我的終端上犯了一個錯字。 – Lack 2015-01-09 22:40:33

0

正如@ gbofi的回答所示,將一個k-mer轉換爲04**k - 1之間的整數非常簡單。另一種,做編碼的主要數學方法是:

def kmer_to_int(kmer): 
    return sum(4**i * "ATGC".index(x) for i, x in enumerate(kmer)) 

如果這是不是建立一個二進制字符串,然後將其轉換爲int更快我沒有測試。

此代碼給出第一個字符的輸入中的最低位的位置,所以"AT"變得0b0100,或4"TA"變得0b00011。如果要將編碼視爲最重要的第一個字母,請在生成器表達式中使用enumerate(reversed(kmer))而不是enumerate(kmer)

正如其他人所評論的,這些整數只對給定長度k是唯一的。如果它們僅以尾隨A(例如,"ATG","ATGA""ATGAA""ATGAAA"等全部編碼爲36)的數量相同,則將針對不同長度的字符串的編碼給出相同的整數。至於你更廣泛的計算特定k-mers在更大序列中出現的更廣泛的目標,我不確定你是否會看到以這種方式編碼k-mers的好處。好處可能取決於數據集的細節。

整數鍵的一個優點是它們允許您使用列表而不是字典來保存您的計數。你可以用lst = [0] * 4**k建立一個合適的列表,然後用lst[kmer_to_int(kmer)] += 1遞增你所看到的值。給定相同數量的條目,列表的開銷比字典的開銷要低,但我不確定這種差異是否會大到足以有所幫助。

如果您的數據稀疏分佈(也就是說,很多4 ** k個可能的k-mer序列從未出現在您的輸入中),使用列表可能仍會浪費大量內存,因爲列表總是4**k元素長。更好的方法可能是使用一些其他方法來簡化對稀疏數據的代碼dict

一個選項是dict類的一些方法,以避免將結果集中的所有值初始化爲0。如果您將您的增量代碼更改爲d[key] = d.get(key, 0) + 1,則無論key是否已在字典中,它都將起作用。

另一種選擇是使用collections.Counter而不是常規的dictCounter類專門用於對輸入序列中的項目實例進行計數,這似乎正是您正在做的。它認爲它沒有看到的任何密鑰的值爲0

相關問題