画一张序列图和它们的反向完全数

Question

我是一个生物信息学家，最近开始学习python，我对绘制一个图很感兴趣。在

节点

set(['ACG', 'ATC', 'GAT', 'CGG', 'CGT', 'AAT', 'ATT', 'GGA', 'TTC', 'CCG', 'TCC', 'GAA'])

边缘

^{pr2}$

当我用上述信息构造正规图时，我得到了12个节点和10个边，即使用下面的函数得到两个断开的图。在

def visualize_de_bruijn():
    """ Visualize a directed multigraph using graphviz """
    nodes = set(['ACG', 'ATC', 'GAT', 'CGG', 'CGT', 'AAT', 'ATT', 'GGA',      'TTC', 'CCG', 'TCC', 'GAA'])
    edges = [('ACG', 'CGG'), ('CGG', 'GGA'), ('GGA', 'GAA'), ('GAA', 'AAT'), ('AAT', 'ATC'), ('GAT', 'ATT'), ('ATT', 'TTC'), ('TTC', 'TCC'), ('TCC', 'CCG'), ('CCG', 'CGT')]
    dot_str = 'digraph "DeBruijn graph" {\n'
    for node in nodes:
        dot_str += '  %s [label="%s"] ;\n' % (node, node)
    for src, dst in edges:
        dot_str += '  %s -> %s ;\n' % (src, dst)
    return dot_str + '}\n'

在生物学中，我们有一个互补碱基配对的概念，其中A=T，T=A，G=C，C=G。所以，ACG的互补词是TGC，ACG的反互补词是CGT，也就是说，我们颠倒补语。在

在节点列表中，我们看到12个节点，其中6个节点是彼此的反向补码，即

ReverseComplement('ACG') = CGT 
ReverseComplement('CGG') = CCG 
ReverseComplement('GGA') = TCC   
ReverseComplement('GAA') = TTC 
ReverseComplement('AAT') = ATT 
ReverseComplement('ATC') = GAT

现在我想构造一个有六个节点的图，一个节点应该有它自己的值和它的逆补值，以及总共10条边，即图不应该断开。如何在python中使用graphviz可视化这个图形。？如果除了graphviz之外还有什么可以帮助我可视化这种图形的，请告诉我。？在

Answer 1

我真的不确定你想在这里完成什么（所以要知道你可能有一个XY problem)），但让我们来看看你的问题，看看它会给我们带来什么。在

a node should have its own value and its reverse complement value

所以我们需要一些对象来存储序列和该序列的反向补码。在

有different ways of making the reverse complement of a sequence。作为一个生物信息学家，你应该真正使用一个适合生物信息学的库，即BioPython。在

那么做一个反向补码如下所示：

from Bio.Seq import Seq

seq = 'ATCG'
print str(Seq(seq).reverse_complement()) # CGAT

但是对于这个问题，生成Seq对象可能有点过头了，所以我只使用下面的标准字典。我们还需要比较Node对象，因此我们需要重写__eq__。因为我们要生成一个唯一的set对象，所以我们还需要实现__hash__。对于漂亮的打印，我们还实现了__str__。在

^{pr2}$

所以现在我们可以把我们的集合或原始节点转换成我们的新节点列表，并用它们的反向补码。在

nodes = set(['ACG', 'ATC', 'GAT', 'CGG', 'CGT', 'AAT', 'ATT', 'GGA', 'TTC', 'CCG', 'TCC', 'GAA'])
newnodes = set(Node(seq) for seq in nodes)
assert len(newnodes) == 6

现在我们需要连接节点。在你的问题中，你没有真正说明你是如何生成带有边的列表的。您发布的内容的可视化看起来确实与您描述的一样：两个断开连接的图。在

然而，当我创建一个DeBruijn图时，我会成对比较所有的序列，看看它们之间是否有重叠，创建一个邻接列表，然后从中生成graphviz的点代码。在

from itertools import product

def suffix(seq, overlap):
    return seq[-overlap:]

def prefix(seq, overlap):
    return seq[:overlap]

def has_overlap_seq(seq1, seq2, overlap=2):
    if seq1 == seq2:
        return False
    return suffix(seq1, overlap) == prefix(seq2, overlap)

def get_adjacency_list_seqs(sequences, overlap=2):
    for seq1, seq2 in product(sequences, repeat=2):
        if has_overlap_seq(seq1, seq2, overlap):
            yield seq1, seq2

def make_dot_plot(adjacency_list):
    """Creates a DOT file for a directed graph."""
    template = """digraph "DeBruijn graph"{{
{}
}}""".format
    edges = '\n'.join('"{}" -> "{}"'.format(*edge) for edge in adjacency_list)
    return template(edges)

如果我为你的原稿做这个nodes

seq_adjacency_list = get_adjacency_list_seqs(nodes)
print make_dot_plot(seq_adjacency_list)

我得到一个单连通图：

所以我不确定在生成edges列表的原始实现中是否存在错误，或者您是否正在尝试完全其他的操作。在

现在向前看，我们可以调整前面的代码来处理序列字符串，从而也可以使用我们之前创建的Node对象。在

def has_overlap_node(node1, node2, overlap=2):
    if node1 == node2:
        return False
    return suffix(node1.seq, overlap) == prefix(node2.seq, overlap) or \
           suffix(node1.seq, overlap) == prefix(node2.revcompl, overlap) or \
           suffix(node1.revcompl, overlap) == prefix(node2.seq, overlap) or \
           suffix(node1.revcompl, overlap) == prefix(node2.revcompl, overlap)

def get_adjacency_list_nodes(nodes, overlap=2):
    for node1, node2 in product(nodes, repeat=2):
        if has_overlap_node(node1, node2, overlap):
            yield node1, node2

应用这个

nodes_adjacency_list = get_adjacency_list_nodes(newnodes)
print make_dot_plot(nodes_adjacency_list)

产生

它确实有6个节点，但有12个而不是请求的10个边。在