最近在閱讀了
Wirkshark網(wǎng)絡(luò)分析的藝術(shù)以及Wireshark網(wǎng)絡(luò)分析就這么簡單以后琼锋，提起了繼續(xù)寫技術(shù)博客的想法挣跋。之前總是覺得很多的小的細(xì)節(jié)可能不用被記錄，所以大概把自己弄的很沒有話說的樣子庆寺，但是后來深受啟發(fā)走净。畢竟每一次項目中的每一次疑問加解答都時可以成為一篇技術(shù)博客的文章的。于是重新開始整理自己所作的工作邮旷。
問題
前陣子一個組由于課題的目的黄选，希望可以設(shè)計一連串的類似于primer的dna序列，不過由于出發(fā)點的問題婶肩，所以是想要基于一段原始的DNA序列進(jìn)行所有突變可能的枚舉的办陷。原始的長度有17bp那么長，那么所有突變的可能大概就有4^17-1這么多種可能律歼，這么大的內(nèi)容我用python是很難進(jìn)行枚舉的民镜。每一次枚舉到了一半的時候就占滿了服務(wù)器500G的內(nèi)存。所以我只能想別的方法险毁。
涉及解決方案
- 枚舉方面制圈。
作為一個程序員们童，當(dāng)然不可能說寫17個循環(huán)從而來枚舉所有的可能。但是不這樣又很難實現(xiàn)這樣的枚舉所有序列問題鲸鹦。后來想到了一個解決方案慧库，就是類似于樹的方法。
init_start = ['A','T','C','G']
while len(init_start) <= 17:
current = init_start.pop(0)
for _i in ['A','T','C','G']:
init_start += [current + _i]
用一個很簡單的類似于BFS(Breadth-First-Search)的一個樹的遍歷的方法馋嗜。因為這個窮舉過程歸根結(jié)底其實就是對一個樹的表現(xiàn)的過程齐板。所以應(yīng)該也會有兩種思路，我用的就是類似于BFS的構(gòu)建過程葛菇。這樣可以保證遍歷所有的節(jié)點覆积。
*但是正如我之前說的，這個窮舉實在太大熟呛，所以沒有什么意義宽档，所以我后來也深入了解了一下對方的課題目的。從而修改了自己的目標(biāo)庵朝，這也是一個 生信人 非常重要的一個要求吗冤。不要盲目的跑程序和寫程序，一定要跟對方溝通九府，讓對方知道你能做什么椎瘟，你發(fā)現(xiàn)的問題。如果對方本身也目的模糊侄旬，你一定要從一個邏輯思考的角度給他梳理清楚肺蔚，否則就是在坑自己了。當(dāng)然也在 坑別人就是儡羔。 *
- 回到最開始的目的宣羊。他們想做應(yīng)該是設(shè)計這么一個library，從而去尋找進(jìn)化的路徑汰蜘。當(dāng)然這個其實是非常困難的仇冯，但是如果實驗設(shè)計的好的話，應(yīng)該也是可以做的族操。那么一個生信的能做些什么呢苛坚。
因為我們現(xiàn)在有一個野生型的序列，我們能做的應(yīng)該是利用計算機(jī)把所有進(jìn)化的可能性盡可能均勻分布色难，從而構(gòu)建這么一個庫泼舱，剩下的可能是一個碰運(yùn)氣的過程了。
我的思路大概是這樣的枷莉。
- 按照與野生型的差異(我們可以簡單的用變了多少個堿基來描述)作為每一列的特征娇昙。這樣我們最多應(yīng)該是有17列的，也就是從突變1個堿基一直到每一個堿基都不一樣了依沮。
- 由于他們需要探索進(jìn)化的途徑涯贞，那么我們的每一列，應(yīng)該都是前一列的一個衍生危喉。(當(dāng)然這個模型十分簡單宋渔，因為也許突變會在幾次突變后重新往野生型的序列去變異)
- 明確下來后，這也就非常的簡單了辜限。因為每一次其實也就是只突變了一個皇拣。
- 接下來就是明確代碼所需要實現(xiàn)的模塊。
a. 給定序列薄嫡，突變位置氧急，突變的數(shù)量，即隨機(jī)的在突變位置的候選中選擇突變的數(shù)量的數(shù)字毫深，然后在這些數(shù)字的基礎(chǔ)上進(jìn)行突變吩坝。
def randome_shuffle(dna_seq, change_area, num_alt):
'''
:param dna_seq:
:param change_area: 0-coordinate. consist of all pos to change.Not continuous.
:param num_alt:
:return:
'''
if len(change_area) < num_alt:
raise SyntaxError, 'alt allele num is bigger than the area range.'
count_alt = 0
changed_pos_bucket = []
while count_alt < num_alt:
pos = random.choice(change_area)
if pos in changed_pos_bucket:
continue
else:
changed_pos_bucket += [pos]
count_alt += 1
_cache = base[::]
_cache.remove(dna_seq[pos].upper())
alt = _cache[random.randint(0, 2)]
if pos != len(dna_seq) - 1:
dna_seq = dna_seq[:pos] + alt + dna_seq[pos + 1:]
else:
dna_seq = dna_seq[:pos] + alt
return dna_seq, tuple(sorted(changed_pos_bucket))
b. 給定序列，突變位置哑蔫，結(jié)合a函數(shù)遞歸生成結(jié)果钉寝。重點在于我們的遞歸每一次只需要基于前一次變化，往后再變化1個新的位置的堿基(自然界中當(dāng)然可能會發(fā)生在已發(fā)生過的位置闸迷，但是這里的模型簡單不考慮)嵌纲，這樣就形成了一個遞歸的鏈條，其實也就是我剛剛說的思路1中的每一行腥沽。逮走。
def recursive_generate_lib(start_seq, area):
bucket = []
if not area:
return bucket
new_seq, changed_pos = randome_shuffle(start_seq, area, 1)
bucket.append(new_seq)
area.remove(changed_pos[0])
if area:
bucket += recursive_generate_lib(new_seq, area)
return bucket
c. 給定需要變化的區(qū)域，突變數(shù)量的范圍進(jìn)行構(gòu)建一個符合要求的dna庫今阳。由于上一個函數(shù)的原因师溅，其實我們僅僅需要給定一個初始點，就可以得出一行盾舌。那么有多少個初始點险胰，就可以有多少行。
從實驗?zāi)康膩碚f矿筝，如果突變數(shù)量從1開始起便，初始點的可能性只有
17*3 == 51種，其中可能還有一部分需要被硬性篩選條件過濾的窖维。所以只剩下有限的24種初始點榆综。但是如果突變數(shù)量從2開始，其實就十分的多了铸史，所以也就不加贅述和計算鼻疮。
從以上的內(nèi)容中，我們實現(xiàn)了為一個探索進(jìn)化路徑實驗設(shè)計DNA庫的編程腳本琳轿。其中還有一部分我之前沒有涉及到的內(nèi)容我也簡單進(jìn)行一下描述判沟。
設(shè)計一段DNA序列耿芹，其實類似于設(shè)計一段primer，有一些硬性的條件需要我們?nèi)M足挪哄。例如：序列中不能有5bp以上的寡聚核苷酸鏈吧秕。序列中不能有自身互補(bǔ)的結(jié)構(gòu)。序列的GC含量需要滿足40%-60%間迹炼。
這些過程都可以很輕易的用Python進(jìn)行實現(xiàn)砸彬，這里我直接貼上代碼，具體的使用和為什么這么寫由大家自己去思考斯入。
base = ['A', 'C', 'T', 'G']
complementary = {'A': 'T',
'T': 'A',
'C': 'G',
'G': 'C'}
from collections import Counter
def cal_GC_content(dna_seq):
bucket = []
for _i in dna_seq:
bucket.append(_i)
caled_dict = Counter(bucket)
for _base in base:
if _base not in caled_dict.keys():
caled_dict[_base] = 0
gc_sum = sum([caled_dict.get('G'), caled_dict.get('C')])
return float(gc_sum) / sum(caled_dict.values())
def filter_poly(dna_seq, poly_num):
for _i in base:
if _i * poly_num in dna_seq:
return False
return True
def filter_self_complementary(dna_seq, length):
for _i in base:
if _i * length in dna_seq:
if complementary[_i] * length in dna_seq:
return False
return True
