Постройте график последовательностей и их обратное дополнение

Question

Я биоинформатик и недавно начал собирать питон, и меня интересует построение графика. У меня есть набор узлов и ребер.

узлы

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' 

В биологии мы имеем эту концепцию комплементарного спаривания оснований, где А = Т, Т = А, С = С и С = С. Таким образом, бесплатный «ACG» - это «ТГК» и обратная поддержка «ACG» = «CGT», то есть мы отменим дополнение.

В списке узлов мы видим 12 узлов, в которых 6 узлов являются обратными дополнением друг к другу т.е.

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

Теперь я хотел бы построить график, где есть шесть узлов, узел должен иметь свой собственный значение и его значение обратного дополнения и общее 10 ребер, т. е. график не должен быть отключен. Как визуализировать этот график, используя graphviz в python.? Если есть что-то еще, кроме графика, который может помочь мне визуализировать этот тип графика, пожалуйста, дайте мне знать.?

Answer 1

Я не совсем уверен, что вы пытаетесь достичь здесь (так и быть в курсе, что вы можете иметь в XY problem), но давайте рассмотрим ваш вопрос и посмотреть, где он получает нас.

Узла должен иметь свое собственное значение и его обратное значение комплемента

Так что нам нужно какой-нибудь объект, чтобы сохранить последовательность и обратный комплемент последовательности.

есть different ways of making the reverse complement of a sequence. Как bioinformatician, вы действительно должны работать с библиотекой, подходящей для биоинформатики, а именно 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, я бы попарно сравнить все последовательности, есть ли у них какие-либо совпадения между ними, создать список смежности и от генерации кода DOT для 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 ребер.