sequencer.py

import fileinput
from StringIO import StringIO
import networkx as nx
import subprocess



# https://neil.fraser.name/news/2010/11/04/
def commonOverlapIndexOf(text1, text2):
  text1_length = len(text1)
  text2_length = len(text2)
  if text1_length == 0 or text2_length == 0:
    return 0
  if text1_length > text2_length:
    text1 = text1[-text2_length:]
  elif text1_length < text2_length:
    text2 = text2[:text1_length]
  if text1 == text2:
    return min(text1_length, text2_length)
  best = 0
  length = 1
  while True:
    pattern = text1[-length:]
    found = text2.find(pattern)
    if found == -1:
      return best
    length += found
    if text1[-length:] == text2[:length]:
      best = length
      length += 1


kmers = []
for line in fileinput.input():
  kmer = line.strip()
  kmers.append(kmer)

minimum_overlap_length = 10
g = nx.DiGraph()
for kmer1 in kmers:
  for kmer2 in kmers:
    if kmer1 == kmer2 or kmer1 in kmer2 or kmer2 in kmer1:
      continue
    best = commonOverlapIndexOf(kmer1, kmer2)
    if best > minimum_overlap_length:
      g.add_edge(kmer1, kmer2, overlap_idx=best)
    best = commonOverlapIndexOf(kmer2, kmer1)
    if best > minimum_overlap_length:
      g.add_edge(kmer2, kmer1, overlap_idx=best)

g = nx.line_graph(g)
dot_file = "./graph.dot"
fname = "./graph.png"
png_file = open(fname, 'w')
nx.write_dot(g, dot_file)
proc = subprocess.Popen(["dot", "-Tpng", dot_file], stdout=png_file)
proc.wait()