基因组索引

基因组

基因组序列是以染色体或者scafford形式排布的， read比对回基因组时开销很大。为了解决这个问题，需要对参考基因组建立索引，加速比对的过程。

建立索引，一般采用kmer对基因组的序列进行处理，然后针对特地的kmer记录其在染色体索引上的位置信息。

Target(T): CGTGCTGGCTT

# 首先针对Target序列，建立索引,kmer=5
# index T         sort index T
# CGTGC: 0        CGTGC: 0,4
# GTGCG: 1        GCGTG: 3
# TGCGT: 2   ==>  GTGCG: 1
# GCGTG: 3        GTGCT: 5
# CGTGC: 4        TGCGT: 2
# GTGCT: 5        TGCTT: 6
# TGCTT: 6

# Pattern(P): GCGTGC 比对回T

# 第一种情况
# P前五个字符 GCGTG比对到T索引为3的位置，同时确认P剩下的一个字符C，发现也能比对上T，
# 因此P能比对上T索引为3的位置。

# 第二种情况
# P的后五个字符CGTGC，能比对上T的0和4位置，但是对于0位置，P的前一个字符G比对不上；
# 对于4位置，G能比对上T，所以P的比对上T位置是4-1=3

# 一共2 index hits

1
2
3

# P: GCGTGA
# 第一种情况
# P的前5个字符 GCGTG比对上T索引为3的位置，但是剩下的字符A比对不上T的C字符, 所以P不在T中

1 2	# P: GCGTAC # P的前五个字符不在索引中，比对不上

数据结构

针对上面需要对基因组进行k-mer索引的情况，需要采用合适的数据结构来存储基因组的这些信息。 Python中的dict可以实现。只需要对Pattern的kmer在T的索引中查找就行，查找可以使用二分法查找。二分法查找的时间复杂度是 $O(log_2(n))$

python 实现

import sys
import bisect 

class Index(object):
    def __init__(self, t, k):
        """创建不同长度的索引"""
        self.k = k # kmer 长度
        self.index = []
        for i in range(len(t) - k + 1): # 循环每个k-mer
            self.index.append((t[i:i+k], i)) # 追加(k-mer, offset)
        self.index.sort() # 按字母顺序排列
    def query(self, p):
        """返回P第一个k-mer匹配的索引"""
        kmer = p[:self.k]
        i = bisect.bisect_left(self.index, (kmer, -1)) # 二分法查找
        hits = []
        while i < len(self.index):
            if self.index[i][0] != kmer:
                break
            hits.append(self.index[i][1])
            i += 1
        return hits

def queryIndex(p, t, index):
    k = index.k
    offsets = []
    for i in index.query(p):
        if p[k:] == t[i+k:i+len(p)]: # 验证P剩下的字符 len(p) - k
            offsets.append(i)
    return offsets

t = 'ACTTGGAGATCTTTGAGGCTAGGTATTCGGGATCGAAGCTCATTTCGGGGATCGATTACGATATGGTGGGTATTCGGGA'
p = 'GGTATTCGGGA'

index = Index(t, 4)
print(queryIndex(p, t, index))