绘制序列图及其反向补码

Question

我是一名生物信息学家，最近开始学习 python，我有兴趣绘制一个 Graph.I 具有一组节点和边的图。

节点

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

边

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

当我使用以上信息构造法线图时，我得到 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. As a bioinformatician, you should really be working with a library suited for bioinformatics, namely BioPython个。

然后进行反向补码如下所示：

from Bio.Seq import Seq

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

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

class Node(object):
    def __init__(self, seq):
        self.seq = seq
        self.revcompl = self.revcompl()

    def revcompl(self):
        complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
        return "".join(complement.get(base, base) for base in reversed(self.seq))

    def __eq__(self, other):
        return self.seq == other.seq or self.revcompl == other.seq

    def __hash__(self):
        return hash(self.seq) ^ hash(self.revcompl)

    def __str__(self):
        return '({}, {})'.format(self.seq, self.revcompl)

所以现在我们可以将我们的集合或原始节点转换为我们的新节点列表及其反向补充。

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 的 DOT 代码。

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 个边。

绘制序列图及其反向补码

Plot a graph of sequences and their reverse complement

python

graph

bioinformatics

graphviz