利用python为探索进化路径的实验构建‘dna’库
最近在阅读了
Wirkshark网络分析的艺术以及Wireshark网络分析就这么简单以后，提起了继续写技术博客的想法。之前总是觉得很多的小的细节可能不用被记录，所以大概把自己弄的很没有话说的样子，但是后来深受启发。毕竟每一次项目中的每一次疑问加解答都时可以成为一篇技术博客的文章的。于是重新开始整理自己所作的工作。
问题
前阵子一个组由于课题的目的，希望可以设计一连串的类似于primer的dna序列，不过由于出发点的问题，所以是想要基于一段原始的DNA序列进行所有突变可能的枚举的。原始的长度有17bp那么长，那么所有突变的可能大概就有4^17-1这么多种可能，这么大的内容我用python是很难进行枚举的。每一次枚举到了一半的时候就占满了服务器500G的内存。所以我只能想别的方法。
涉及解决方案
- 枚举方面。
作为一个程序员，当然不可能说写17个循环从而来枚举所有的可能。但是不这样又很难实现这样的枚举所有序列问题。后来想到了一个解决方案，就是类似于树的方法。
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)的一个树的遍历的方法。因为这个穷举过程归根结底其实就是对一个树的表现的过程。所以应该也会有两种思路，我用的就是类似于BFS的构建过程。这样可以保证遍历所有的节点。
*但是正如我之前说的，这个穷举实在太大，所以没有什么意义，所以我后来也深入了解了一下对方的课题目的。从而修改了自己的目标，这也是一个 生信人 非常重要的一个要求。不要盲目的跑程序和写程序，一定要跟对方沟通，让对方知道你能做什么，你发现的问题。如果对方本身也目的模糊，你一定要从一个逻辑思考的角度给他梳理清楚，否则就是在坑自己了。当然也在 坑别人就是。 *
- 回到最开始的目的。他们想做应该是设计这么一个library，从而去寻找进化的路径。当然这个其实是非常困难的，但是如果实验设计的好的话，应该也是可以做的。那么一个生信的能做些什么呢。
因为我们现在有一个野生型的序列，我们能做的应该是利用计算机把所有进化的可能性尽可能均匀分布，从而构建这么一个库，剩下的可能是一个碰运气的过程了。
我的思路大概是这样的。
- 按照与野生型的差异(我们可以简单的用变了多少个碱基来描述)作为每一列的特征。这样我们最多应该是有17列的，也就是从突变1个碱基一直到每一个碱基都不一样了。
- 由于他们需要探索进化的途径，那么我们的每一列，应该都是前一列的一个衍生。(当然这个模型十分简单，因为也许突变会在几次突变后重新往野生型的序列去变异)
- 明确下来后，这也就非常的简单了。因为每一次其实也就是只突变了一个。
- 接下来就是明确代码所需要实现的模块。
a. 给定序列，突变位置，突变的数量，即随机的在突变位置的候选中选择突变的数量的数字，然后在这些数字的基础上进行突变。
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. 给定序列，突变位置，结合a函数递归生成结果。重点在于我们的递归每一次只需要基于前一次变化，往后再变化1个新的位置的碱基(自然界中当然可能会发生在已发生过的位置，但是这里的模型简单不考虑)，这样就形成了一个递归的链条，其实也就是我刚刚说的思路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. 给定需要变化的区域，突变数量的范围进行构建一个符合要求的dna库。由于上一个函数的原因，其实我们仅仅需要给定一个初始点，就可以得出一行。那么有多少个初始点，就可以有多少行。
从实验目的来说，如果突变数量从1开始，初始点的可能性只有
17*3 == 51种，其中可能还有一部分需要被硬性筛选条件过滤的。所以只剩下有限的24种初始点。但是如果突变数量从2开始，其实就十分的多了，所以也就不加赘述和计算。
从以上的内容中，我们实现了为一个探索进化路径实验设计DNA库的编程脚本。其中还有一部分我之前没有涉及到的内容我也简单进行一下描述。
设计一段DNA序列，其实类似于设计一段primer，有一些硬性的条件需要我们去满足。例如：序列中不能有5bp以上的寡聚核苷酸链。序列中不能有自身互补的结构。序列的GC含量需要满足40%-60%间。
这些过程都可以很轻易的用Python进行实现，这里我直接贴上代码，具体的使用和为什么这么写由大家自己去思考。
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
