【问题标题】:Plotting the score matrix from a Needleman-Wunsch pairwise sequence alignment in matplotlib从 matplotlib 中的 Needleman-Wunsch 成对序列比对中绘制得分矩阵
【发布时间】:2017-04-29 02:03:28
【问题描述】:

我正在尝试根据 Python 中的全局对齐算法(或 Needleman–Wunsch algorithm)绘制矩阵。

我不知道 matplotlib 是否是这种情况下的最佳工具。我尝试使用Bokeh,但结构很难适应我想要的矩阵。

我使用Bio.SeqIO(BioPython 的标准序列输入/输出接口)来存储两个序列。

我怎样才能得到与这张图片相似的结果:

这在 Matplotlib 中可能吗?我该怎么做?

更新

最后,我能够根据 ImportanceOfBeingErnest 给出的答案构建算法。结果如下:

这里是这个实现的要点:plot_needleman_wunsch.py

这是整个项目(正在进行中):bmc-sequence-alignment

【问题讨论】:

  • 如果你添加你的两个数据序列,我会试一试。
  • seq_alpha = Seq("GATCCA")seq_beta = Seq("GTGCCT")
  • 对不起,我不是生物统计学专家。你能重现数字和箭头(颜色和方向)的来源吗? matplotlib 绝对没有内置的东西来解决这个问题。但如果你能提供,我可以尝试绘制这个
  • 我也不是生物统计学家,但我会尝试用例子来解释。我对箭头的颜色不感兴趣,但对方向感兴趣。要了解箭头的来源,请阅读tutorial 中的 3.1.13.1.2 部分。简而言之,我有两个大小为 m 和 n 的数组。矩阵为(m + 1 x n + 1)。然后用“空白”填充第一行和第一列。为了获得方向,询问对角线、左侧和向上之间的最大值是多少
  • “要了解箭头的来源,请阅读本教程中的 3.1.1 和 3.1.2 部分。”我对阅读如何重现您的数据不感兴趣。如果您提供了您想要绘制的数据,那么有人可能会为您绘制它。

标签: python matplotlib bioinformatics biopython


【解决方案1】:

在问题中放置箭头的算法没有明确的描述;因此,这个答案的重点是在 matplotlib 中产生类似情节的方式。

这里的想法是将数字放置在图中的整数位置并在n+0.5 处绘制次要网格线以获得类似表格的外观。箭头绘制为 4 列数组中定义的位置之间的注释(前 2 列:箭头起点的 x 和 y,第三和第四列:箭头终点的 x,y)。

import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import numpy as np;np.random.seed(5)
plt.rcParams["figure.figsize"] = 4,5
param = {"grid.linewidth" : 1.6,
         "grid.color"     : "lightgray",
         "axes.linewidth" : 1.6,
         "axes.edgecolor" : "lightgray"}
plt.rcParams.update(param)

#Data
headh = list("GATCCA")
headv = list("GTGCCT")

v = np.zeros((7,7), dtype=int)
v[1:,1:] = np.random.randint(-2,7, size=(6,6))

arrows = np.random.randint(0,v.shape[1], size=(14,4))
opt = np.array([(0,1),(1,0),(1,1)])
arrows[:,2:] = arrows[:,:2] + opt[np.random.randint(0,3,size=14 )]

arrowsb = np.random.randint(0,v.shape[1], size=(7,4))
optb = np.array([(0,1),(1,0),(1,1)])
arrowsb[:,2:] = arrowsb[:,:2] + optb[np.random.randint(0,3,size=7 )]

#Plot
fig, ax=plt.subplots()
ax.set_xlim(-1.5, v.shape[1]-.5 )
ax.set_ylim(-1.5, v.shape[0]-.5 )
ax.invert_yaxis()
for i in range(v.shape[0]):
    for j in range(v.shape[1]):
        ax.text(j,i,v[i,j], ha="center", va="center")
for i, l in enumerate(headh):
    ax.text(i+1,-1,l, ha="center", va="center", fontweight="semibold")
for i, l in enumerate(headv):
    ax.text(-1,i+1,l, ha="center", va="center", fontweight="semibold")

ax.xaxis.set_minor_locator(ticker.FixedLocator(np.arange(-1.5, v.shape[1]-.5,1)))
ax.yaxis.set_minor_locator(ticker.FixedLocator(np.arange(-1.5, v.shape[1]-.5,1)))
plt.tick_params(axis='both', which='both', bottom='off', top='off', 
                left="off", right="off", labelbottom='off', labelleft='off')
ax.grid(True, which='minor')


arrowprops=dict(facecolor='crimson',alpha=0.5, lw=0, 
                shrink=0.2,width=2, headwidth=7,headlength=7)
for i in range(arrows.shape[0]):
    ax.annotate("", xy=arrows[i,2:], xytext=arrows[i,:2], arrowprops=arrowprops)
arrowprops.update(facecolor='blue')
for i in range(arrowsb.shape[0]):
    ax.annotate("", xy=arrowsb[i,2:], xytext=arrowsb[i,:2], arrowprops=arrowprops)
plt.show()

【讨论】:

  • 哇,非常感谢!这就是我想要的。我将使用适合算法的解决方案上传更新。再次感谢您。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2013-11-24
  • 2012-09-22
  • 1970-01-01
相关资源
最近更新 更多