K-mer子序列的生成

作者: 何物昂 | 来源:发表于2019-10-27 12:12 被阅读0次

K-mer子序列的生成
基因组组装----k-mer
Python应用——分析序列k-mer
Kraken使用手册
基因组分析 K-mer 第3回估计1000bp的全基因长度
基因组组装和注释理论
KE_TG
基因组分析 K-mer 第1回理解K-mer和Coverage
Kettle案例三:添加序列和选择/改名值
C++生成随机数据方法

简要介绍

在生信中，k-mer指生物序列中长度为k的子序列。 $k$ -mer包含着生物序列的两个基本特征：
1. 单体组分信息
2. 序列顺序信息

通过两个信息可以基本确定一个序列。在许多生物相关领域都有着广泛应用，如用于序列组装(构建De Bruijn图)，生物序列特征提取(机器学习的输入特征)，序列分析等领域。

同时在不属于生物相关的其他领域，如计算语言学中，它被称为n-gram。

一条长度为 $L$ 的序列，有 $L-k+1$ 个k-mer。比如序列: "ATGCA"有5个1-mer(A,T,G,C,A)，4个2-mer(AT,TG,GC,CA)，2个4-mer(ATGC,TGCA)。对于一个由n种单体组成的序列，理论上可能存在 $n^k$ 种k-mer, 比如对于DNA序列来说可能存在 $4^4$ 种不同的4-mer。对于蛋白质一级结构序列来说，可能存在 $20^3$ 种不同3-mer。

计算k-mer，当前已有优化过的软件来计算k-mer。这里代码k-mer，仅供学习使用。

代码

从这里看到一份有点技巧性的的代码

from itertools import islice

seq = "ATGCGGCAA"
n = 3
def n_grams(a, n):
    z = (islice(a, i, None) for i in range(n))
    return zip(*z)

for i in n_grams(seq, n):
    print(i)
#('A', 'T', 'G')
#('T', 'G', 'C')
#('G', 'C', 'G')
#('C', 'G', 'G')
#('G', 'G', 'C')
#('G', 'C', 'A')
#('C', 'A', 'A')

下文分解代码~

islice(iterable, start, stop[, step])

islice创建一个生成器，返回从 iterable 里选中的元素。start,stop为下标索引值，start表示从什么位置开始，stop 表示到什么位置结束，step为步伐，默认为1。

from itertools import islice
seq = "ATGCGGCAA"
k1 = islice(seq, 0, None)  # 从索引0开始，返回一个生成器
for i in k1: 
    print(i, end=" ")  # 默认是换行符，改成空格，一行输出
#  A T G C G G C A A

k2 = islice(seq, 1, None)  # 从索引1开始
for i in k2: 
    print(i, end=" ")  
#  T G C G G C A A

k3 = islice(seq, 2, None)  # 从索引2开始，
for i in k3: 
    print(i, end=" ")  
#  G C G G C A A

z = (islice(seq, i, None) for i in range(n))

z 也是一个生成器，可生成 islice(seq, 0, None) ，islice(seq, 1, None)， islice(seq, 2, None)

zip(*z)

z = (islice(seq, i, None) for i in range(n))
z1 = zip(*z) ## *拆包，会将z迭代器拆成三个元素，每个元素是一个生成
器
### (islice(seq, 0, None) ,islice(seq, 1, None), islice(seq, 2, None) = *z
"""
kmer1 = [i for i in z1]

z2 = zip(('A', 'T', 'G', 'C', 'G', 'G', 'C', 'A', 'A'),
         ('T', 'G', 'C', 'G', 'G', 'C', 'A', 'A'),
         ('G', 'C', 'G', 'G', 'C', 'A', 'A'), )
kmer2 = [i for i in z2]

kmer1 等于kmer2。

zip 的作用是： 同时迭代zip() 里的每个可迭代对象，
并将其组合成一个新元组, 一旦其中某个对象迭代完成就停止。
如：
A T G
T G C
...
C A A
"""

for i in z1: 
    print(i) 
#('A', 'T', 'G')
#('T', 'G', 'C')
#('G', 'C', 'G')
#('C', 'G', 'G')
#('G', 'G', 'C')
#('G', 'C', 'A')
#('C', 'A', 'A')

技巧

假如我要生成3-ker的字符串, 
开始生成这样的三个字符串，利用islice的特性
ATGCGGCAA  # 即从0开始
TGCGGCAA  ## 去掉第1个字符，从1开始
GCGGCAA   ## 去掉第2个字符， 从2开始
之后同时迭代形成新的元组，利用zip的特性
(A, T, G)
(T, G, C)
(G, C, G)
...
这样就可保持这样的取值形式（seq[i], seq[i+1], seq[i+2]）

其它

然而。。。我以为弄那花里胡哨的的，会提升代码运行速度。实际并没有。还不如下面这种形式直接点。

seq = "ATGCGGCAA"
for i in range(len(seq)-3+1):
    print(seq[i:i+3])
#('A', 'T', 'G')
#('T', 'G', 'C')
#('G', 'C', 'G')
#('C', 'G', 'G')
#('G', 'G', 'C')
#('G', 'C', 'A')
#('C', 'A', 'A')