精华内容
下载资源
问答
  • 对齐方式 已实现的DNA序列比对算法的集合,包括最佳全局比对,带状全局比对和用于多个序列比对的近似算法。
  • 生物DNA序列比对算法研究!!!生物DNA序列比对算法研究!!!生物DNA序列比对算法研究!!!
  • 一款DNA序列比对的软件,可用于DNA序列的分析、比对,酶切位点的查找等。
  • cljam:Clojure的DNA序列比对图(SAM)库
  • Pontos是一种易于使用的图形Java程序,用于根据PHYLIP格式的DNA序列比对计算未校正的距离(或相似性)矩阵。 它还会与常规对齐方式产生“差异”对齐(反之亦然)。 它可以以不同的方式处理差距和歧义。 缺口可以是:...
  • 用于序列比对器的多平台体系结构(MASA)是一种灵活的软件框架,可简化在多个硬件/软件平台中创建DNA序列比对应用程序的过程。 该框架支持将巨大的DNA序列与超过2亿个碱基对(MBP)进行比对。 MASA-Core库基于...
  • Swift是一个DNA序列比对程序,使用Smith-Waterman算法产生缺口比对。 它以查询文件(FASTA格式)和参考文件(FASTA格式)为输入。 它输出参考名称,读取名称,空位比对,比对得分,比对开始​​和结束位置以及比对...
  • 提出了一种新的基于统计的序列特征定义,并以此对进化树进行...同时对这一方法进行了分类测试和进化树构建测试,发现所得测试结果与同样是用进化树进行构建的基于两序列比对的传统算法比较,其性能和适用性都有明显改善.
  • 它使用高性能的比对模块完善了MashMap,该模块能够计算非常大的序列的基本水平比对。 过程 每个查询序列分为-s[N], --segment-length=[N]定义的非重叠部分。 然后使用MashMap的滑动minhash映射算法和后续的过滤步骤...
  • 第一个并行的Smith-Waterman算法利用Intel Xeon Phi群集来加速长DNA序列比对。 该算法是用C ++(具有一组SIMD内部扩展),OpenMP和MPI编写的。 性能评估表明,我们的算法实现了非常稳定的性能,在单个至强融核上...
  • 关于生物信息学的资料 上面是进化算法在基因序列中的应用
  • cout请依次输入要比较的DNA序列(‘#’结束):"; cout<<"DNA1:"; cin>>dna1[i]; while(dna1[i]!='#') { i++; cin>>dna1[i]; } cout; cin>>dna2[j]; while(dna2[j]!='#') { j++; cin>>...

    #include<iostream>
    #include<stack>
    #include<stdlib.h>
    using namespace std;
    
    stack<char> s;//当前搜索路径的LCS
    stack<char> lcs[100];//所有的LCS
    
    int count=0;//记录LCS的数量
    
    bool lcs_have_exist(stack<char> lcs[],int count,stack<char> s)
    {//判断当前搜索的LCS是否已存在于 stack_lcs[MAX]
    	int i;
    	bool exist=false;
    	for (i=0;i<count;i++)
    	{
    		if(lcs[i]==s)
    		{
             exist=true;
    		 break;
    		}
    	}
    	return exist;
    }
    
    
    
    void clear_stack(stack<char> &s)
    {
    	while(!s.empty())
    	{
    		s.pop();
    	}
    }
    
    void out_stack()//输出栈s内容
    {
    	stack<char> copy=s;
    	cout<<count+1<<":";
    	while(!copy.empty())
    	{   
    		cout<<copy.top()<<" ";
    		copy.pop();
    	}
    }
    
    
    void ALL_LCS(char X[],int b[][100],int i,int j)
    {
    	if(i==0||j==0)//当前搜索路径走到尽头
    	{   
    		if(!lcs_have_exist(lcs,count,s))//若当前搜索的LCS是新的LCS,则保存且输出
    		{
    			out_stack();
    			cout<<endl;
    			lcs[count]=s;
    			count++;//增加一个LCS
    			
    		}
    		
    	}
    	else if(i>0&&j>0)
    	{
    		if(b[i][j]==3)//沿着对角线向上走
    		{
    			s.push(X[i]);
    			ALL_LCS(X,b,i-1,j-1);
    			s.pop();
    		}
    		else if(b[i][j]==1)//水平向左走
    		{
    		    ALL_LCS(X,b,i,j-1);
    		}
    	    else if(b[i][j]==2)//垂直向上走
    		{
                ALL_LCS(X,b,i-1,j);
    		}
    	    else if(b[i][j]==4)//先向左走,再向右
    		{
    	    	ALL_LCS(X,b,i,j-1);
    			ALL_LCS(X,b,i-1,j);
    		}
    	}
    }
    
    
    
    void LCS(char dna1[],char dna2[],int m,int n,int b[][100],int tag[][100])
    {//生成矩阵tag[m][n],b[m][n]作为逆推的记录
    	int i,j,k;
    	for(k=0;k<=m;k++)
    		 tag[k][0]=0;//初始化矩阵
    	 for( k=0;k<=n;k++)
    		 tag[0][k]=0;
    	for(i=1;i<=m;i++)
    		for(j=1;j<=n;j++)
    		{
    			if(dna1[i]==dna2[j])
    			{
    				tag[i][j]=tag[i-1][j-1]+1;
    				b[i][j]=3;//"↖"
    			}
    			else if(tag[i][j-1]<tag[i-1][j])
    			{
                    tag[i][j]=tag[i-1][j];
    				b[i][j]=2;//"↑" 
    			}
    			else if(tag[i][j-1]>tag[i-1][j])
    			{  
    				tag[i][j]=tag[i][j-1];
    				b[i][j]=1;//"←"
    			}
    			else
    			{
    				tag[i][j]=tag[i][j-1];
    				b[i][j]=4;//"←↑"
    			}
    		}
    		for(i=0;i<=m;i++)
    		{
    	    	cout<<endl;
    		    for(j=0;j<=n;j++)
    			{
    		    	cout<<tag[i][j]<<" ";
    			}
    		}
    }
    
    int main()
    {
         const int max=100;//dna序列最大个数
         char dna1[max];
    	 char dna2[max];//定义dna
    	 int tag[max][max];//定义矩阵
    	 int b[max][max];//记录轨迹
    	 int i=1;
    	 int j=1;
    	 int k=0;
         cout<<"请依次输入要比较的DNA序列(‘#’结束):"<<endl;
    	 cout<<"DNA1:";
    	 cin>>dna1[i];
    	 while(dna1[i]!='#')
    	 {
    		  i++;
    	      cin>>dna1[i];
    	 }
    	 cout<<"DNA2:";
    	 cin>>dna2[j];
    	  while(dna2[j]!='#')
    	 {
    		  j++;
    	      cin>>dna2[j];
    	 }
    	 cout<<endl;
         cout<<"DNA1的序列:";
    	 for(k=1;k<i;k++)
    	 {
    	     cout<<dna1[k]<<"  ";//输入需要比较的DNA序列
    	 }
    	 cout<<endl;
    	 cout<<"DNA2的序列:";
    	 for( k=1;k<j;k++)
    	 {
    	     cout<<dna2[k]<<"  ";
    	 }
    	 cout<<endl;
        i--;
    	j--;//剔除#号
    
    	LCS(dna1,dna2,i,j,b,tag);
    
        cout<<endl;
    	cout<<"最长子序列:"<<endl;	
    	ALL_LCS(dna1,b,i,j);
    	cout<<endl;	
    	return 0;
    }
    
    


    展开全文
  • 序列比对是生物信息学中基本的信息处理方法之一.DNA序列比对分为多序列比对和双序列比对两种.本文首先分析了经典蚁群算法在双序列比对中的应用,然后针对经典蚁群算法收敛速度慢和陷入局部最优值等缺点进行改进,...
  • 动态序列比对 [Ctrl+A 全选 注:如需引入外部Js需刷新才能执行]
  • DNA序列局部比对算法

    2018-03-09 18:22:08
    Smith-Waterman 算法的实现 配对 错位 缺失后得到的分值
  • DNA序列比对是生物信息学中的最重要的任务之一。本文针对多序列比对的特点,提出一种渐进蚁群算法,即将渐进比对算法和蚁群算法相结合。在渐进蚁群算法中,既能克服蚁群算法易于陷入局部最优解、收敛速度慢的特点...
  • 生物信息原理作业第三弹:DNA序列局部比对,利用Smith–Waterman算法,python3.6代码实现。 实例以及原理均来自https://en.wikipedia.org/wiki/Smith%E2%80%93Waterman_algorithm。 DNA序列局部比对 转载请保留...

    生物信息原理作业第三弹:DNA序列局部比对,利用Smith–Waterman算法,python3.6代码实现。

    实例以及原理均来自https://en.wikipedia.org/wiki/Smith%E2%80%93Waterman_algorithm

    DNA序列局部比对

    转载请保留出处!

     1 import numpy as np
     2 import pandas as pd
     3 sequence1 = 'TGTTACGG'
     4 sequence2 = 'GGTTGACTA'
     5 s1 = ''
     6 s2 = ''
     7 gap = -2
     8 score_matrix = pd.read_excel('score.xlsx')      #匹配得分
     9 print(score_matrix)
    10 best_matrix = np.empty(shape= (len(sequence2)+1,len(sequence1)+1),dtype = int)
    11 def get_match_score(s1,s2):
    12     score = score_matrix[s1][s2]
    13     return score
    14 
    15 def get_matrix_max(matrix):                    #得到最大分数下标
    16     Max = matrix.max()
    17     for i in range(len(sequence2)+1):
    18         for j in range(len(sequence1)+1):
    19             if matrix[i][j] == Max:
    20                 return (i,j)
    21             
    22 for i in range(len(sequence2)+1):
    23     for j in range(len(sequence1)+1):
    24         if i == 0 or j == 0:
    25             best_matrix[i][j] = 0
    26         else:
    27             match = get_match_score(sequence2[i-1],sequence1[j-1])
    28             gap1_score = best_matrix[i-1][j] + gap
    29             gap2_score = best_matrix[i][j-1] + gap
    30             match_score = best_matrix[i-1][j-1]+match
    31             score = max(gap1_score,gap2_score,match_score)
    32             if score>0:
    33                 best_matrix[i][j] = score
    34             else:
    35                 best_matrix[i][j] = 0
    36 print(best_matrix)
    37 
    38 #traceback
    39 i,j = get_matrix_max(best_matrix)
    40 while(best_matrix[i][j]!= 0):
    41     match = get_match_score(sequence2[i-1],sequence1[j-1])
    42     if i>0 and j>0 and best_matrix[i][j] == best_matrix[i-1][j-1]+match:
    43         s1 += sequence1[j-1]
    44         s2 += sequence2[i-1]
    45         i-=1;j-=1
    46     elif i>0 and best_matrix[i,j] == best_matrix[i-1,j]+gap:
    47         s1+='-'
    48         s2+=sequence2[i-1]
    49         i-=1
    50     else:
    51         s1+=sequence1[j-1]
    52         s2+='-'
    53         j-=1
    54 print(s1[::-1]+'\n'+s2[::-1])

    感觉我的得分矩阵写成Excel不必要,等我熟悉一下Numpy和Python命令行之后会修改的。

    转载于:https://www.cnblogs.com/zxzhu/p/7930376.html

    展开全文
  • 生物信息学原理作业第二弹:利用Needleman–Wunsch算法进行DNA序列全局比对。 具体原理:https://en.wikipedia.org/wiki/Needleman%E2%80%93Wunsch_algorithm。 利用Needleman–Wunsch算法进行DNA序列全局比对 ...

    生物信息学原理作业第二弹:利用Needleman–Wunsch算法进行DNA序列全局比对。

    具体原理:https://en.wikipedia.org/wiki/Needleman%E2%80%93Wunsch_algorithm

    利用Needleman–Wunsch算法进行DNA序列全局比对

    转载请保留出处!

    贴上python代码:

     1 # -*- coding: utf-8 -*-
     2 """
     3 Created on Sat Nov 25 18:20:01 2017
     4 
     5 @author: zxzhu
     6 后需修改:
     7 1.加命令行参数
     8 2.给出多种比对结果
     9 """
    10 
    11 import numpy as np
    12 import pandas as pd
    13 sequence1 = 'AACGTACTCA'
    14 sequence2 = 'TCGTACTCA'
    15 s1 = ''
    16 s2 = ''
    17 gap = -4
    18 score_matrix = pd.read_excel('score.xlsx')        #score matrix
    19 best_matrix = np.empty(shape= (len(sequence2)+1,len(sequence1)+1),dtype = int)
    20 
    21 def get_match_score(s1,s2):
    22     score = score_matrix[s1][s2]
    23     return score
    24 
    25 for i in range(len(sequence2)+1):
    26     for j in range(len(sequence1)+1):
    27         if i == 0:
    28             best_matrix[i][j] = gap * j
    29         
    30         elif j == 0:
    31             best_matrix[i][j] = gap *i
    32         else:
    33             match = get_match_score(sequence2[i-1],sequence1[j-1])
    34             gap1_score = best_matrix[i-1][j]+gap
    35             gap2_score = best_matrix[i][j-1]+gap
    36             match_score = best_matrix[i-1][j-1]+match
    37             best_matrix[i][j] = max(gap1_score,gap2_score,match_score)
    38 print(best_matrix)
    39 i,j = len(sequence2),len(sequence1)
    40 while(i>0 or j>0):
    41     match = get_match_score(sequence2[i-1],sequence1[j-1])
    42     if i>0 and j>0 and best_matrix[i][j] == best_matrix[i-1][j-1]+match:
    43         s1 += sequence1[j-1]
    44         s2 += sequence2[i-1]
    45         i-=1;j-=1
    46     elif i>0 and best_matrix[i,j] == best_matrix[i-1,j]+gap:
    47         s1+='-'
    48         s2+=sequence2[i-1]
    49         i-=1
    50     else:
    51         s1+=sequence1[j-1]
    52         s2+='-'
    53         j-=1
    54 print(s1[::-1]+'\n'+s2[::-1])

    后面会加入命令行。

    多种结果这里只取了一种,这个问题有待解决。

    如果有其他的方法我会及时添加。

    转载于:https://www.cnblogs.com/zxzhu/p/7903706.html

    展开全文
  • 一维成对CNN用于两个DNA序列的整体比对
  • DNA序列翻译为氨基酸

    2017-10-14 13:50:22
    perl 代码写的 DNA序列比对,用的是smith比对的算法,只要提交输入比较比对即可
  • DNA序列动态规划法比对的C代码,生物信息学的内容,动态规划法
  • Ferox是一种DNA序列比对应用程序,它使用模糊k-mers快速,准确地比对参考基因组的序列读段进行比对。 Ferox还可以用于比对整个基因组。 Ferox使用的种子机制是高度可配置的,从而允许在XML配置文件中声明性地创建...
  • MashMap实现了一种快速且近似的算法,用于计算长DNA序列之间的局部比对边界。 对于将基因组装配或长读(PacBio / ONT)映射到参考基因组可能很有用。 给定最小对齐长度和所需局部对齐的同一性阈值,Mashmap使用k -...
  • 全局与局部序列比对的实现(DNA)程序能实现什么一、实现步骤1.用户输入步骤2.代码实现步骤二、实验结果截图1.全局比对2.局部比对总结 程序能实现什么 a.完成gap值的自定义输入以及两条需比对序列的输入 b.完成得分...


    程序能实现什么

    a.完成gap值的自定义输入以及两条需比对序列的输入
    b.完成得分矩阵的计算及输出
    c.输出序列比对结果
    d.使用matplotlib对得分矩阵路径的绘制


    一、实现步骤

    1.用户输入步骤

    a.输入自定义的gap值
    b.输入需要比对的碱基序列1(A,T,C,G)换行表示输入完成
    b.输入需要比对的碱基序列2(A,T,C,G)换行表示输入完成

    输入(示例):
    在这里插入图片描述

    2.代码实现步骤

    1.获取到用户输入的gap,s以及t
    2.调用构建得分矩阵函数,得到得分矩阵以及方向矩阵
    3.将得到的得分矩阵及方向矩阵作为参数传到回溯函数中开始回溯得到路径,路径存储使用的是全局变量,存的仍然是方向而不是坐标位置减少存储开销,根据全局变量中存储的方向将比对结果输出。
    4.根据全局变量中存储的方向使用matplotlib画出路径

    全局比对代码如下:

    import matplotlib.pyplot as plt
    import numpy as np
    
    #定义全局变量列表finalList存储最后回溯的路径  finalOrder1,finalOrder2存储最后的序列 finalRoad用于存储方向路径用于最后画图
    def createList():
        global finalList
        global finalOrder1
        global finalOrder2
        global finalRoad
        finalList = []
        finalOrder1 = []
        finalOrder2 = []
        finalRoad = []
    
    
    #创建A G C T 对应的键值对,方便查找计分矩阵中对应的得分
    def createDic():
        dic = {'A':0,'G':1,'C':2,'T':3}
        return dic
    #构建计分矩阵
    # A G C T
    def createGrade():
        grade = np.matrix([[10,-1,-3,-4],
                           [-1,7,-5,-3],
                           [-3,-5,9,0],
                           [-4,-3,0,8]])
        return grade
    
    #计算两个字符的相似度得分函数
    def getGrade(a,b):
        dic = createDic()  # 碱基字典 方便查找计分矩阵
        grade = createGrade()  # 打分矩阵grade
        return grade[dic[a],dic[b]]
    
    #构建得分矩阵函数 参数为要比较序列、自定义的gap值
    def createMark(s,t,gap):
        a = len(s)                         #获取序列长度a,b
        b = len(t)
        mark = np.zeros((a+1,b+1))         #初始化全零得分矩阵
        direction = np.zeros((a+1,b+1,3))  #direction矩阵用来存储得分矩阵中得分来自的方向   第一个表示左方 第二个表示左上 第三个表示上方 1表示能往哪个方向去
                                           #由于得分可能会来自多个方向,所以使用三维矩阵存储
    
        direction[0][0] = -1               #确定回溯时的结束条件 即能够走到方向矩阵的值为-1
        mark[0,:] = np.fromfunction(lambda x, y: gap * (x + y), (1, b + 1), dtype=int)     #根据gap值将得分矩阵第一行计算出
        mark[:,0] = np.fromfunction(lambda x, y: gap * (x + y), (1, a + 1), dtype=int)     #根据gap值将得分矩阵第一列计算出
        for i in range(1,b+1):
            direction[0,i,0] = 1
        for i in range(1, a + 1):
            direction[i, 0, 2] = 1
    
        for i in range(1,a+1):
            for j in range(1,b+1):
                threeMark = [mark[i][j-1],mark[i-1][j-1],mark[i-1][j]]         #threeMark表示现在所要计算得分的位置的左边 左上 上边的得分
                threeGrade = [gap,getGrade(s[i-1],t[j-1]),gap]                 #threeGrade表示经过需要计算得左边 左上 上边的空位以及相似度得分
                finalGrade = np.add(threeMark,threeGrade)                      #finalGrade表示最终来自三个方向上的得分
                mark[i][j] = max(finalGrade)                                   #选取三个方向上的最大得分存入得分矩阵
                #可能该位置的得分可以由多个方向得来,所以进行判断并循环赋值
                for k in range(0,len([y for y,x in enumerate(finalGrade) if x == max(finalGrade)])):
                    directionList = [y for y,x in enumerate(finalGrade) if x == max(finalGrade)]
                    direction[i][j][directionList[k]]  = 1
        return mark,direction
    
    #回溯函数 参数分别为 得分矩阵 方向矩阵 现在所处得分矩阵的位置 以及两个序列
    def remount(mark,direction,i,j,s,t):
        if direction[i][j][0] == 1 :
            if direction[i][j-1][0] == -1:            #如果该位置指向左边 先判断其左边是否是零点
                finalList.append(0)                   #如果是 将该路径存入路径列表
                finalList.reverse()                   #将列表反过来得到从零点开始的路径
                index1 = 0                            #记录现在所匹配序列s的位置 因为两个字符串可能是不一样长的
                index2 = 0                            #记录现在所匹配序列t的位置
                for k in finalList:
                    if k == 0 :
                        finalOrder1.append("-")
                        finalOrder2.append(t[index2])
                        index2 += 1
                    if k == 1 :
                        finalOrder1.append(s[index1])
                        finalOrder2.append(t[index2])
                        index1 += 1
                        index2 += 1
                    if k == 2 :
                        finalOrder1.append(s[index1])
                        finalOrder2.append("-")
                        index1 += 1
                finalList.reverse()  # 将原来反转的路径再返回来
                finalRoad.append(np.array(finalList))  # 将此次的路径添加到最终路径记录用于最后画图
                finalList.pop()                       #输出后将当前方向弹出 并回溯
                return
            else :
                finalList.append(0)                   #如果不是零点 则将该路径加入路径矩阵,继续往下走
                remount(mark,direction,i,j-1,s,t)
                finalList.pop()                       #该方向走完后将这个方向弹出  继续下一轮判断 下面两个大的判断同理
        if direction[i][j][1] == 1 :
            if direction[i-1][j-1][0] == -1:
    
                finalList.append(1)
                finalList.reverse()                   # 将列表反过来得到从零点开始的路径
                index1 = 0                            # 记录现在所匹配序列s的位置 因为两个字符串可能是不一样长的
                index2 = 0                            # 记录现在所匹配序列t的位置
                for k in finalList:
                    if k == 0 :
                        finalOrder1.append("-")
                        finalOrder2.append(t[index2])
                        index2 += 1
                    if k == 1 :
                        finalOrder1.append(s[index1])
                        finalOrder2.append(t[index2])
                        index1 += 1
                        index2 += 1
                    if k == 2 :
                        finalOrder1.append(s[index1])
                        finalOrder2.append("-")
                        index1 += 1
                finalList.reverse()  # 将原来反转的路径再返回来
                finalRoad.append(np.array(finalList))  # 将此次的路径添加到最终路径记录用于最后画图
                finalList.pop()
                return
            else :
                finalList.append(1)
                remount(mark,direction,i-1,j-1,s,t)
                finalList.pop()
        if direction[i][j][2] == 1 :
            if direction[i-1][j][0] == -1:
                finalList.append(2)
                finalList.reverse()                      # 将列表反过来得到从零点开始的路径
                index1 = 0                               # 记录现在所匹配序列s的位置 因为两个字符串可能是不一样长的
                index2 = 0                               # 记录现在所匹配序列t的位置
                for k in finalList:
                    if k == 0 :
                        finalOrder1.append("-")
                        finalOrder2.append(t[index2])
                        index2 += 1
                    if k == 1 :
                        finalOrder1.append(s[index1])
                        finalOrder2.append(t[index2])
                        index1 += 1
                        index2 += 1
                    if k == 2 :
                        finalOrder1.append(s[index1])
                        finalOrder2.append("-")
                        index1 += 1
                finalList.reverse()                    # 将原来反转的路径再返回来
                finalRoad.append(np.array(finalList))  # 将此次的路径添加到最终路径记录用于最后画图
                finalList.pop()
                return
            else :
                finalList.append(2)
                remount(mark,direction,i-1,j,s,t)
                finalList.pop()
    
    #画箭头函数
    def arrow(ax,sX,sY,aX,aY):
        ax.arrow(sX,sY,aX,aY,length_includes_head=True, head_width=0.15, head_length=0.25, fc='w', ec='b')
    
    #画图函数
    def drawArrow(mark, direction, a, b, s, t):
        #a是s的长度为4    b是t的长度为6
        fig = plt.figure()
        ax = fig.add_subplot(111)
        val_ls = range(a+2)
        scale_ls = range(b+2)
        index_ls = []
        index_lsy = []
        for i in range(a):
            if i == 0:
                index_lsy.append('#')
            index_lsy.append(s[a-i-1])
        index_lsy.append('0')
        for i in range(b):
            if i == 0:
                index_ls.append('#')
                index_ls.append('0')
            index_ls.append(t[i])
        plt.xticks(scale_ls, index_ls)           #设置坐标字
        plt.yticks(val_ls, index_lsy)
        for k in range(1,a+2):
            y = [k for i in range(0,b+1)]
            x = [x for x in range(1,b+2)]
            ax.scatter(x, y, c='y')
        for i in range(1,a+2):
            for j in range(1,b+2):
                ax.text(j,a+2-i,int(mark[i-1][j-1]))
        lX = b+1
        lY = 1
        for n in range(0,len(finalRoad)):
            for m in (finalRoad[n]):
                if m == 0:
                    arrow(ax,lX,lY,-1,0)
                    lX = lX - 1
                elif m == 1:
                    arrow(ax,lX,lY,-1,1)
                    lX = lX - 1
                    lY = lY + 1
                elif m == 2:
                    arrow(ax, lX, lY, 0, 1)
                    lY = lY + 1
            lX = b + 1
            lY = 1
        ax.set_xlim(0, b + 2)  # 设置图形的范围,默认为[0,1]
        ax.set_ylim(0, a + 2)  # 设置图形的范围,默认为[0,1]
        ax.set_aspect('equal')  # x轴和y轴等比例
        plt.show()
        plt.tight_layout()
    if __name__ == '__main__':
        createList()
        print("Please enter gap:")
        gap = int(input())              #获取gap值 转换为整型   tip:刚开始就是因为这里没有进行类型导致后面的计算部分报错
        print("Please enter sequence 1:")
        s = input()                     #获取用户输入的第一条序列
        print("Please enter sequence 2:")
        t = input()                     #获取用户输入的第二条序列
        a = len(s)                      #获取s的长度
        b = len(t)                      #获取t的长度
        mark,direction = createMark(s,t,gap)
        print("The scoring matrix is as follows:")         #输出得分矩阵
        print(mark)
        remount(mark,direction,a,b,s,t) #调用回溯函数
        c = a if a > b else b          #判断有多少种比对结果得到最终比对序列的长度
        total = int(len(finalOrder1)/c)
        for i in range(1,total+1):     #循环输出比对结果
            k = str(i)
            print("Sequence alignment results "+k+" is:")
            print(finalOrder1[(i-1)*c:i*c])
            print(finalOrder2[(i-1)*c:i*c])
        drawArrow(mark, direction, a, b, s, t)
    
    
    

    局部比对代码如下

    import matplotlib.pyplot as plt
    import numpy as np
    import operator
    #在局部比对中 回溯结束的条件是方向矩阵中该位置的值全为0
    #定义全局变量列表finalList存储最后回溯的路径  finalOrder1,finalOrder2存储最后的序列
    def createList():
        global finalList
        global finalOrder1
        global finalOrder2
        global finalRoad
        finalList = []
        finalOrder1 = []
        finalOrder2 = []
        finalRoad = []
    
    
    #创建A G C T 对应的键值对,方便查找计分矩阵中对应的得分
    def createDic():
        dic = {'A':0,'G':1,'C':2,'T':3}
        return dic
    #构建计分矩阵
    # A G C T
    def createGrade():
        grade = np.matrix([[10,-1,-3,-4],
                           [-1,7,-5,-3],
                           [-3,-5,9,0],
                           [-4,-3,0,8]])
        return grade
    
    #计算两个字符的相似度得分函数
    def getGrade(a,b):
        dic = createDic()  # 碱基字典 方便查找计分矩阵
        grade = createGrade()  # 打分矩阵grade
        return grade[dic[a],dic[b]]
    
    #构建得分矩阵函数 参数为要比较序列、自定义的gap值
    def createMark(s,t,gap):
        a = len(s)                         #获取序列长度a,b
        b = len(t)
        mark = np.zeros((a+1,b+1))         #初始化全零得分矩阵
        direction = np.zeros((a+1,b+1,3))  #direction矩阵用来存储得分矩阵中得分来自的方向   第一个表示左方 第二个表示左上 第三个表示上方 1表示能往哪个方向去
                                           #由于得分可能会来自多个方向,所以使用三维矩阵存
        for i in range(1,a+1):
            for j in range(1,b+1):
                threeMark = [mark[i][j-1],mark[i-1][j-1],mark[i-1][j]]         #threeMark表示现在所要计算得分的位置的左边 左上 上边的得分
                threeGrade = [gap,getGrade(s[i-1],t[j-1]),gap]                 #threeGrade表示经过需要计算得左边 左上 上边的空位以及相似度得分
                finalGrade = np.add(threeMark,threeGrade)                      #finalGrade表示最终来自三个方向上的得分
                if max(finalGrade) >= 0:                                   #如果该最大值是大于0的则 选取三个方向上的最大得分存入得分矩阵 否则不对矩阵进行修改
                    mark[i][j] = max(finalGrade)
                    for k in range(0,len([y for y,x in enumerate(finalGrade) if x == max(finalGrade)])): #可能该位置的得分可以由多个方向得来,所以进行判断并循环赋值
                        directionList = [y for y,x in enumerate(finalGrade) if x == max(finalGrade)]
                        direction[i][j][directionList[k]]  = 1
        return mark,direction
    
    #回溯函数 参数分别为 得分矩阵 方向矩阵 现在所处得分矩阵的位置 以及两个序列
    def remount(mark,direction,i,j,s,t):
        if direction[i][j][0] == 1 :
            if all(direction[i][j-1] == [0,0,0]):            #如果该位置指向左边 先判断其左边是否是零点
                finalList.append(0)                   #如果是 将该路径存入路径列表
                finalList.reverse()                   #将列表反过来得到从零点开始的路径
                index1 = i                            #记录现在所匹配序列s的位置 因为两个字符串可能是不一样长的
                index2 = j-1                            #记录现在所匹配序列t的位置
                for k in finalList:
                    if k == 0 :
                        finalOrder1.append("-")
                        finalOrder2.append(t[index2])
                        index2 += 1
                    if k == 1 :
                        finalOrder1.append(s[index1])
                        finalOrder2.append(t[index2])
                        index1 += 1
                        index2 += 1
                    if k == 2 :
                        finalOrder1.append(s[index1])
                        finalOrder2.append("-")
                        index1 += 1
                finalList.reverse()
                finalRoad.append(np.array(finalList))  # 将此次的路径添加到最终路径记录用于最后画图
                finalList.pop()                       #输出后将当前方向弹出 并回溯
                return
            else :
                finalList.append(0)                   #如果不是零点 则将该路径加入路径矩阵,继续往下走
                remount(mark,direction,i,j-1,s,t)
                finalList.pop()                       #该方向走完后将这个方向弹出  继续下一轮判断 下面两个大的判断同理
        if direction[i][j][1] == 1 :
            if all(direction[i-1][j-1] == [0,0,0]):
                finalList.append(1)
                finalList.reverse()                     # 将列表反过来得到从零点开始的路径
                index1 = i-1                            # 记录现在所匹配序列s的位置 因为两个字符串可能是不一样长的
                index2 = j-1                            # 记录现在所匹配序列t的位置
                for k in finalList:
                    if k == 0 :
                        finalOrder1.append("-")
                        finalOrder2.append(t[index2])
                        index2 += 1
                    if k == 1 :
                        finalOrder1.append(s[index1])
                        finalOrder2.append(t[index2])
                        index1 += 1
                        index2 += 1
                    if k == 2 :
                        finalOrder1.append(s[index1])
                        finalOrder2.append("-")
                        index1 += 1
                finalList.reverse()
                finalRoad.append(np.array(finalList))  # 将此次的路径添加到最终路径记录用于最后画图
                finalList.pop()
                return
            else :
                finalList.append(1)
                remount(mark,direction,i-1,j-1,s,t)
                finalList.pop()
        if direction[i][j][2] == 1 :
            if all(direction[i-1][j] == [0,0,0]):
                finalList.append(2)
                finalList.reverse()                      # 将列表反过来得到从零点开始的路径
                index1 = i-1                               # 记录现在所匹配序列s的位置 因为两个字符串可能是不一样长的
                index2 = j                               # 记录现在所匹配序列t的位置
                for k in finalList:
                    if k == 0 :
                        finalOrder1.append("-")
                        finalOrder2.append(t[index2])
                        index2 += 1
                    if k == 1 :
                        finalOrder1.append(s[index1])
                        finalOrder2.append(t[index2])
                        index1 += 1
                        index2 += 1
                    if k == 2 :
                        finalOrder1.append(s[index1])
                        finalOrder2.append("-")
                        index1 += 1
                finalList.reverse()
                finalRoad.append(np.array(finalList))  # 将此次的路径添加到最终路径记录用于最后画图
                finalList.pop()
                return
            else :
                finalList.append(2)
                remount(mark,direction,i-1,j,s,t)
                finalList.pop()
    #画箭头函数
    def arrow(ax,sX,sY,aX,aY):
        ax.arrow(sX,sY,aX,aY,length_includes_head=True, head_width=0.15, head_length=0.25, fc='w', ec='b')
    
    #画图函数
    def drawArrow(mark, direction, a, b, s, t,mx,my):
        #a是s的长度为4    b是t的长度为6
        fig = plt.figure()
        ax = fig.add_subplot(111)
        val_ls = range(a+2)
        scale_ls = range(b+2)
        index_ls = []
        index_lsy = []
        for i in range(a):
            if i == 0:
                index_lsy.append('#')
            index_lsy.append(s[a-i-1])
        index_lsy.append('0')
        for i in range(b):
            if i == 0:
                index_ls.append('#')
                index_ls.append('0')
            index_ls.append(t[i])
        plt.xticks(scale_ls, index_ls)           #设置坐标字
        plt.yticks(val_ls, index_lsy)
        for k in range(1,a+2):
            y = [k for i in range(0,b+1)]
            x = [x for x in range(1,b+2)]
            ax.scatter(x, y, c='y')
        for i in range(1,a+2):
            for j in range(1,b+2):
                ax.text(j,a+2-i,int(mark[i-1][j-1]))
        lX = my + 1
        lY = a - mx + 1
        for n in range(0,len(finalRoad)):
            for m in (finalRoad[n]):
                if m == 0:
                    arrow(ax,lX,lY,-1,0)
                    lX = lX - 1
                elif m == 1:
                    arrow(ax,lX,lY,-1,1)
                    lX = lX - 1
                    lY = lY + 1
                elif m == 2:
                    arrow(ax, lX, lY, 0, 1)
                    lY = lY + 1
            lX = b + 1
            lY = 1
        ax.set_xlim(0, b + 2)  # 设置图形的范围,默认为[0,1]
        ax.set_ylim(0, a + 2)  # 设置图形的范围,默认为[0,1]
        ax.set_aspect('equal')  # x轴和y轴等比例
        plt.show()
        plt.tight_layout()
    
    if __name__ == '__main__':
        createList()
        print("Please enter gap:")
        gap = int(input())              #获取gap值 转换为整型   tip:刚开始就是因为这里没有进行类型导致后面的计算部分报错
        print("Please enter sequence 1:")
        s = input()                     #获取用户输入的第一条序列
        print("Please enter sequence 2:")
        t = input()                     #获取用户输入的第二条序列
        a = len(s)                      #获取s的长度
        b = len(t)                      #获取t的长度
        mark,direction = createMark(s,t,gap)
        print("The scoring matrix is as follows:")         #输出得分矩阵
        print(mark)
        maxDirection = np.argmax(mark)  #获取最大值的位置
        i = int(maxDirection/(b+1))
        j = int(maxDirection - i*(b+1))
        remount(mark,direction,i,j,s,t) #调用回溯函数
        print(finalOrder1)
        print(finalOrder2)
        drawArrow(mark, direction, a, b, s, t, i, j)
    

    二、实验结果截图

    1.全局比对

    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    2.局部比对

    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述


    总结

    本次实验使用动态规划对全局序列比对进行了实现,自己卡的最久的地方是回溯以及画图的时候。刚开始在实现回溯的过程中,老是找不准回溯的条件以及将所有的路径都记录下来的方法,最后是使用的方向矩阵,也就是重新定义一个与得分矩阵等大的矩阵(但是这个矩阵是三维),存放的是每个位置能够回溯的方向,第一个数值表示左边,第二个表示左上,第三个表示上方,为0时表示当前方向不能回溯,没有路径,为1时表示能回溯,当该位置的所有能走的方向都走完时即可返回。将所有路径记录下来的方法是定义全局变量,当有路径能够走到终点时便将这条路径存放入该全局变量中。
    绘图的时候使用的是matplotlib中的散点图,然后将每个点的得分以注释的形式标记在该点的右上角,并用箭头将路径绘出。不得不说的是,这个图确实太丑了,我学识浅薄,也没想到能画出这个图的更好的方法,还希望老师指点。
    总的来说这次实验经历的时间还比较长,主要是因为python也没有很熟悉,很多函数也是查了才知道,然后可视化更是了解的少,所以画出来的图出奇的丑,还有回溯的时候也是脑子转不过弯来,所以要学习的东西还有很多,需要更加努力。
    本次实验还能够有所改进的地方是:
    1.把两个比对算法结合,让用户能够选择使用哪种比对方式。
    2.作出一个更好看的界面,增加用户体验感。
    3.把图画的更美观。
    (老丁已阅,USC的同学们谨慎借鉴)

    展开全文
  • 输出序列比对结果 d.使用matplotlib对得分矩阵路径的绘制 一、实现步骤 1.用户输入步骤 a.输入自定义的gap值 b.输入需要比对的碱基序列1(A,T,C,G)换行表示输入完成 b.输入需要比对的碱基序列2(A,T,C,G)换行表示...
  • alignTools是用于逐对DNA序列比对的C实现工具的集合。 它包括以下对齐功能: 整体一致性 局部语言 适合对齐(跳跃状态) 重叠对齐 编辑距离 不久的将来将添加更多功能。 版本 0.7.23-r15 开始使用 安装 $ git ...
  • 序列比对与Clustal的使用以及各类常见的序列分析工具介绍 内容提要 第一部分多序列比对 意义方法算法 Clustal的使用 1.Clustalx 2.Clustalw 第二部分常见的序列分析软件分类简介 第一部分 多序列比对及Clustal的...
  • JFinisher是用于比对,编辑和操纵生物序列的软件。 它旨在帮助完成基因组组装。 从参考序列开始,程序使用Smith-Waterman局部比对算法和辅助方法对重叠群进行比对,从而可以管理所产生的比对。 它具有图形界面,可...

空空如也

空空如也

1 2 3 4 5 ... 15
收藏数 290
精华内容 116
关键字:

dna序列比对