debruijn-tp/debruijn/debruijn.py

"""
Small assembly module based on de bruijn graphs
"""
import os
import statistics
import networkx as nx
from networkx import algorithms

def read_fastq(fichier):
    """
    Returns an iterator object that retrieves only the nucleic sequences of a
    fastq file.
    Arguments:
        fichier, str: path to fastq file

    Returns:
        a str generator, generator of sequences
    """
    with open(fichier, 'r') as filin:
        for line in filin:
            yield filin.readline().strip()
            filin.readline()
            filin.readline()

def cut_kmer(seq, k):
    """
    Returns an iterator that returns k-mers of k-size of a sequence
    Arguments:
        seq, str: a sequence
        k, int: k-mer size, must be shorter than len(seq)

    Returns:
        an iterator returning str
    """
    for i in range(len(seq)-(k-1)):
        yield seq[i:i+k]

def build_kmer_dict(fichier, k):
    """
    Arguments:
        fichier, str: path to fastq file
        k, int: k-mer size, must be shorter than len(seq)

    Return:
        hash_table, dict: dictionnary with key = k-mer as str
                          and value count of k-mer occurence
    """
    hash_table = {} # initialise empty hash table
    it_fastq = read_fastq(fichier)
    for seq in it_fastq: # for each sequence
        it_kmer = cut_kmer(seq, k) # count each occurence of k-mer
        for kmer in it_kmer:
            try:
                hash_table[kmer]
            except KeyError:
                hash_table[kmer] = 1
            else:
                hash_table[kmer] += 1

    return hash_table

def build_graph(hash_table):
    """
    Returns a graph from a hash table

    Arguments:
        hash_table, dict: dictionnary obtained with build_kmer_dict() function
    Return:
        graph, nx.DiGraph: the de Bruijn tree corresponding to hash_table
    """
    graph = nx.DiGraph()
    for key in hash_table:
        graph.add_edge(key[:-1], key[1:], weight=hash_table[key])

    return graph

def get_starting_nodes(graph):
    """
    Returns the list of starting nodes of a graph

    Arguments:
        graph, nx.DiGraph: de Bruijn tree

    Return:
        starting_nodes, list of strings: list of starting nodes 
    """
    starting_nodes = []
    for node in graph:
        if graph.in_degree(node) == 0: # if count of input edge == 0
            starting_nodes.append(node)

    return starting_nodes

def std(values):
    """
    Computes standard deviation from a list of value
    Arguments:
        values, list: list of values

    Returns :
        standard deviation of the 'values' data list
    """
    return statistics.stdev(float(values))



def get_sink_nodes(graph):
    """

    Arguments:
        graph, nx.DiGraph: de Bruijn tree

    Return:
        sink_nodes, list of strings: list of terminal nodes 
    """
    sink_nodes = []
    for node in graph:
        if graph.out_degree(node) == 0: # if count of output edge == 0
            sink_nodes.append(node)

    return sink_nodes


def path_average_weight(graph, path):
    """
    Arguments:
        graph, nx.DiGraph: a de bruijn graph
        path, list of str(nodes): list of nodes constituing a path

    Return:
        mean weight, float: the mean weight of the path
    """
    weight = 0
    for i in range(len(path)-1):
        weight += graph[path[i][i+i][weight]
    
    return weight/(len(path)-1)



def remove_paths(graph, paths, delete_entry_node=False, delete_sink_node=False):
    """
    Removes paths in 'paths' list from 'graph' object. Don't delete entry and
    sink nodes unless specified.

    Arguments:
        graph, nd.DiGraph(): de bruijn graph
        paths, list of string lists: list of paths composed of nodes
        delete_entry_node, boolean: Delete entry if True
        delete_sink_node, boolean: Delete sink node if True

    Return:
        graph, nd.DiGraph(): de bruijn graph with deletes paths.
    """
    entry = 1
    sink = -1
    if delete_entry_node == True:
        entry = 0
    if delete_sink_node == True:
        sink = None
    for path in paths:
        graph.remove_nodes_from(path[entry:sink]) 

    return graph


def select_best_path(graph, paths, path_len, mean_weights,
                     delete_entry_node=False, delete_sink_node=False):
    """
    
    """
    max_weight = max(mean_weights)
    heaviest = [i for i, j in enumerate(mean_weights) if j == mean_weights]
    if len(heaviest) > 1:
        max_len = max(path_lengths)
        longest = [i for i in heaviest if path_len[i] == max_len]
        if len(longest) > 1:
            Random.seed(9001)
            best = random.choice[longest]
        else:
            best = longest[0]
    else:
        best = heaviest[0]
    paths.pop(best)

    return remove_paths(graph, paths, delete_entry_node, delete_sink_node)


def fill(text, width=80):
    """Split text with a line return to respect fasta format"""
    return os.linesep.join(text[i:i+width] for i in range(0, len(text), width))

def save_contigs(tuples, outname):
    """
    Arguments:
        tuples, tuple: Obtained from get_contigs()
        outname, str: name of the file to be written
    """
    i = 0
    with open(outname, "w") as outfile:
        for duo in tuples:
            i += 1
            outfile.write(">contig_{} len={}\n".format(i, duo[1]))
            outfile.write("{}\n".format(fill(duo[0])))
    return


def get_contigs(graph, starting_nodes, sink_nodes):
    """
    Arguments:
        graph, nx.DiGraph: de Bruijn tree
        starting_nodes, list of strings: list of starting nodes 
        sink_nodes, list of strings: list of terminal nodes 

    Return:
        contigs, list of tupple: list of tupple (contigs, len(contigs))
    """
    contigs = []
    for starting_node in starting_nodes:
        for sink_node in sink_nodes:
            if algorithms.has_path(graph, starting_node, sink_node) == True:
                path = algorithms.shortest_path(graph, starting_node, sink_node)
                contig = path[0] # base of the contig is seq of the first node
                for i in range(len(path)-1):
                    contig += path[i+1][-1] # adds last char of node
                contigs.append((contig, len(contig)))
    
    return contigs

def solve_bubble():
    pass


def simplify_bubbles():
    pass


def solve_entry_tips():
    pass


def solve_out_tips():
    pass