生信自学笔记（九）智慧的长者与多序列联配之clustal全局联配算法

要不，还是先讲个黑暗的小故事吧。

国王愈来愈烦躁了，他觉得这个国家满哪儿都是人，大街上走着人，池塘里泡着人，屋顶上晾着人，自己去四下巡游，什么风景都看不着。

“这可不行，这么多人，东西哪够分啊，好东西都分给这么些人了，我还玩啥呢？”国王愤愤嘟囔了半个月，终于下定决心，张榜公告天下：地不分南北，人不分贵贱，全国除皇室外所有家庭，一律只能有一个孩子；多余者，自己亲自送入皇宫充当奴隶；如有瞒报，满门抄斩。

全体国民在皇宫卫士的刀枪下都十分顺从地把自己多余的孩子送到了奴隶集中营里面，但是宰相就不乐意了，他位高权重，只可惜前段时间周游列国没来得及阻止国王，现在回来也只能听任生米做成熟饭，但他膝下有六个儿子，手心手背都是肉，哪个也舍不得送走。手足无措之际，他突然想起自己这次在野草国请回一位一见如故的智者，据说他还掌握着神秘的魔法，赶紧叫来智者商量对策。

智者听罢扶了扶灰框眼镜，抚掌大笑：“易事，易事！”

智者首先把宰相的六个儿子召来，他双目一扫，立即先叫两个人站了出来。

宰相大惑不解：“他们是一对双胞胎，长者缘何先将他们找出来。”

智者呵呵一笑：“哈哈哈，难怪他们这么相似，原来是双胞胎。你少安毋躁，一会我自会解释。”

智者衣袖一挥，原本就十分相似的两人登时变得一模一样。智者又找了一对很像的兄弟出来，也让他们变得一模一样。智者将第二对人中的一个和第一对人中的一个比较一番，看到第一对人脸上和第二对人相比多了一点黑痣，衣袖一挥，第二对的第一个人脸上便多了一个黑痣，再一挥，第二对的第二个人脸上也多了一个黑痣。

智者如法炮制，很快，宰相的六个儿子就看上去几乎一模一样了。

智者拍拍手：“好啦，您可以对外宣称您只有一个儿子了。只是要小心，千万不要让其他人看见他们同时出现。”

宰相带着六个儿子一齐伏倒在地：“感谢大师！”

站起来后，宰相掸了掸身上的土：“大师我还是不明白，为何您要以这样一种奇怪的方式改变我儿子的面貌，您既然有易容之术，就不能随便挑一个儿子作为模板，把其他孩子依次变成他的模样吗？”

智者又是一笑：“你错啦。我没有你说的易容术。我身上只有两个本领，一是给我两个人，我能用肉眼迅速判断他们外貌的相似程度；二是我能给某个人添上某个特征。所以我要先两个两个地找出最相似的人，施法让他们尽量看上去一样。如果下次有必要对这二人的面貌进行修改，那么怎么改第一个人的面貌，就对第二个人如法炮制。如果人更多，比如我一开始选了四个人，那么我只要知道其中一个要怎样修改，剩下的人我也就知道了。”

宰相深自叹服，亲备黄金百两送老人离去。

春去秋来，王国的人果然不像以前那样多了，只是路上蹒跚的老年人越来越多。宰相也很老了，他什么事也不想管，整日静静坐在家里，听人家议论：“宰相家的孩子真是贪玩，这么大了还天天都能在外面看见他。”宰相这时就摇着头叹气：“哪里是天天出去啊，明明是孩子们六天才能出去一趟。”

这一年，城里突然流行起伤寒来，宰相家也没能幸免，一个孩子不幸去世了。国王闻之亲自前来，慰问失去了"唯一孩子"的老宰相。白天的葬礼很快举行完了，宰相摆出酒席，与国王大醉一场。到了深夜，国王想起白天他国进贡的珍品还没拆看，赶紧起身告辞。

宰相府很大，没来几次的国王很快在重重亭台中迷路了，他身为国王，自然不能做出问路这种有失体统的事情，便只得跌跌撞撞向前走。走着走着，突然一间房里响起音乐声，国王好奇心顿起，看看四下无人便扒着床沿往里看去：在里面轻拢丝竹的，正是自己在白天亲眼目睹着下葬的可怜的年轻人！

国王险些叫出声，一阵冰冷彻骨之后，他软绵绵地倒了下来。

clustal算法的核心部分，也就是智者将六个儿子变得外貌上极为相似的整个过程。智者的两个法宝，对应着碱基相似度的计算和前面的Needleman Wusch两条序列全局联配算法。

具体步骤

首先，导入Needleman Wusch算法，这里的引入相对路径还是颇费周折的。嗯，sys.path里面存放的就是python寻找包的路径的集合。

import os
import sys
# 添加两条序列联配算法所在的路径
sys.path.append(os.path.abspath('../两条序列联配算法'))
import Needleman_Wusch as NW

然后记录所有序列对之间的相似度得分，分数越小越相似。

# 得分越小，相似度越高
S = ["TAC","TACG", "TAACG", "TTGG", "TTATGG", "TTCC", "AATCG"]
T = []
for i in range(0, len(S)):for j in range(i+1, len(S)):F = NW.mark(S[i], S[j], 0, 1, 4)T.append({'score': F[-1][-1], 'seq1': S[i], 'seq2': S[j]})
T.sort(key = lambda x: x['score'])

完成“如法炮制”函数

# change_all函数的功能是计算如何能从L[0]通过逐步加空格(用-表示)推导出s，并把从L[0]推导到s的方法迁移应用于其他所有L中的元素
def change_all(L,s):for i in range(len(s)-1):if L[0][i] != s[i]:for j in range(len(L)):L[j] = L[j][0:i]+'-'+L[j][i:]if L[0] == s:breakif L[0] != s:k = len(s)-len(L[0])for j in range(len(L)):L[j] = L[j]+'-'*k

完成整个过程

# 存储原始序列
F = []
# 存储联配之后序列
FF = []
for i in T:if len(FF) == len(S):breakif i['seq1'] not in F and i['seq2'] not in F:F.append(i['seq1'])F.append(i['seq2'])S1 = []S2 = []t1 = ""t2 = ""H = NW.mark(i['seq1'], i['seq2'], 0, 1, 4)NW.traceback(" "+i['seq1'], " "+i['seq2'], 0, 1, 4, H, H.shape[0]-1, H.shape[1]-1, S1, S2, t1, t2)if len(S1) == 1:S1.append(S2[0])else:S1[1] = S2[0]if len(FF) != 0:S3 = []S4 = []t3 = ""t4 = ""K = NW.mark(S1[0], FF[0], 0, 1, 4)NW.traceback(" " + S1[0], " " + FF[0], 0, 1, 4, K, K.shape[0] - 1, K.shape[1] - 1, S3, S4, t3, t4)if FF[0] != S4[0]:change_all(FF,S4[0])if S3[0] != S1[0]:change_all(S1,S3[0])FF.append(S1[0])FF.append(S1[1])elif i['seq1'] not in F:F.append(i['seq1'])S3 = []S4 = []t3 = ""t4 = ""K = NW.mark(i['seq1'], FF[0], 0, 1, 4)NW.traceback(" " + i['seq1'], " " + FF[0], 0, 1, 4, K, K.shape[0] - 1, K.shape[1] - 1, S3, S4, t3, t4)if FF[0] != S4[0]:change_all(FF, S4[0])FF.append(S3[0])elif i['seq2'] not in F:F.append(i['seq2'])S3 = []S4 = []t3 = ""t4 = ""K = NW.mark(i['seq2'], FF[0], 0, 1, 4)NW.traceback(" " + i['seq2'], " " + FF[0], 0, 1, 4, K, K.shape[0] - 1, K.shape[1] - 1, S3, S4, t3, t4)if FF[0] != S4[0]:change_all(FF, S4[0])FF.append(S3[0])
print('原始序列\t 联配后序列')
for i, j in zip(F, FF):print(i,'\t',j)

最终结果

原始序列  联配后序列
TACG     T-A-CG
TTGG     T-T-GG
TTCC     T-T-CC
TAACG    T-AACG
AATCG    A-ATCG
TAC      T-A-C-
TTATGG   TTATGG

当然，还是有着一些问题的，主要是直接调用原来的函数就会很冗余，但又懒得改成新的函数。另外我觉得可以改进的一点就是，一对序列和已经比较好的一组序列进行联配时，一开始的两序列联配的对象并不是简单粗暴的两组序列的第一个序列，而是分别位于两组序列之中的最为相似的一对序列，不过这样又会很复杂，反正各有利弊吧hhh。

生信自学笔记（九）智慧的长者与多序列联配之clustal全局联配算法相关推荐

生信自学笔记（十二）：基因组序列与基因预测
基因组在生物学中,一个生物体的基因组是指包含在该生物的DNA(部分病毒是RNA)中的全部遗传信息,或者说是一套染色体中完整的DNA序列. 对于单倍体细胞,基因组是指编码序列和非编码序列在内的全部DN ...
生信自学笔记（五）计分矩阵的实例
氨基酸替换矩阵 PAM 替换矩阵 PAM(Point Accepted Multation) 是基于进化的点突变模型产生的,如果两种氨基酸替换频繁,说明自然界接受这种替换,那么这对氨基酸替换得分就高. ...
生信自学笔记（二）生物信息
基本类型 1. 核苷酸序列数据 DNA 或 RNA 当中四种碱基的排列顺序. DNA : A T C G RNA : A G C U 2. 蛋白质序列和结构数据蛋白质序列是指 20 种氨基酸的排列顺 ...
生信分析学习笔记：（2）GO KEGG分析
生信分析学习笔记:(2)GO KEGG分析介绍教程 1.富集分析 (Over-Representation Analysis ) 2.GSEA(Gene Set Enrichment Analys ...
计算机通路的基本概念,【生信学习笔记】KEGG分子通路数据库
原标题:[生信学习笔记]KEGG分子通路数据库首先什么是一个通路? 通路可以定义为a series of actions among molecules in a cell,细胞中的分子的一系列的行 ...
**生信自学记录1——获取Fastq格式的反向互补序列**
` 生信自学记录1--获取Fastq格式的反向互补序列总共分为三步 1.读取基因序列的str格式,返回反向互补序列str 2.打开fastq格式的文本提取基因序列,返回互补序列list 3.读取互补 ...
生信漫谈送你超级好用的多序列拼接软件
作为科研人员(生物专业从业者,医学生,临床医生等),接触的各种测序非常多,那么怎么快速的分析你测序的结果呢,测了正向,测了反向序列,是不是测通了,只有拼接在一起才知道,科研真是费脑 https://m ...
生信学习笔记：fastp质控处理生成的report结果解读
文章目录前言 raw data 和 fastq文件 reads Q20和Q30 N值 Adapters Duplication Insert fastp report summary Adapter ...
Ubuntu 20.0.4 linux生信服务器笔记
系统硬盘挂载情况 $ sudo root # df -h查看硬件raid信息 # lspci |grep -i raid 17:00.0 RAID bus controller: Broadcom / ...

生信自学笔记（九）智慧的长者与多序列联配之clustal全局联配算法

具体步骤

生信自学笔记（九）智慧的长者与多序列联配之clustal全局联配算法相关推荐

最新文章

热门文章