접미사 트라이(Suffix trie)의 활용 #

NGS sequencing의 결과 적게는 수 천만 read 에서 수 억개의 read 가 생산되게 된다.

이 수많은 read를 genome에 align하기 위해 툴들이 가장 많이 사용하는 방법이 genome을 suffix-trie 자료구조로 변환하는 것이다. 변환된 genome의 자료구조를 바탕으로 read에 matching 시켜 위치를 찾는다.

아래 코드는 "" AATCGGGTTCAATCGGGGT ATCG GGGT "" 를 input으로 받을 경우 첫번째 string을 genome으로 인식하여 suffix-trie 자료구조를 만들고 두번째 행 이후로의 string을 pattern 혹은 read라고 생각하여 genome 상의 위치를 출력하는 코드이다.

class Node:
    def __init__(self, name):
        self.name = name
        self.posterior_node = {}
        self.pos = None

    def set_posterior_node(self, node, edge):
        self.posterior_node[edge] = node

    def get_posterior_node(self, edge):
        try:
            return self.posterior_node[edge]
        except KeyError:
            return -1
    def __str__(self):
        return str(self.name)

    def set_pos(self, pos):
        self.pos = pos

    def get_pos(self,pos):
        return self.pos


def depth_first(node):
    for edge in node.posterior_node:
        posterior_node = node.posterior_node[edge]
        print('{}->{}:{}'.format(node.name, posterior_node.name, edge))
        depth_first(posterior_node)


def check(node, text):
    next_node = node.get_posterior_node(text[0])
    if next_node == -1:
        return False
    elif next_node.get_posterior_node('$') != -1:
        return True
    else:
        return check(next_node, text[1:])


def main(file_):
    global number
    number = 0
    with open(file_) as f:
        lines = f.readlines()
        text = lines[0]
        patterns = [line.strip() for line in lines[1:]]

    first_node = Node(number)
    number += 1

    ## make trie
    suffix_pos = 0
    for pattern in patterns:
        for i in range(len(pattern)):
            edge = pattern[i]
            if i == 0:
                prior_node = first_node
            else:
                prior_node = next_node

            next_node = prior_node.get_posterior_node(edge)
            if next_node == -1:
                next_node = Node(number)
                number += 1
                prior_node.set_posterior_node(next_node, edge)
        last_node = Node(number)
        next_node.set_posterior_node(last_node, '$')
        number += 1

    #depth_first(first_node)

    ## find pattern match position
    result_match_pos = []

    for i in range(len(text)):
        if check(first_node, text[i:]):
            result_match_pos.append(i)


    print ' '.join([str(pos) for pos in result_match_pos])



if __name__ == '__main__':
    main(sys.argv[1])

Incoming Links #

Related Bioinformaticses (Bioinformatics 0) #

• FM-Index