计算每列的 kmer 核苷酸频率

Question

我有一个序列列表：

CAGGTAGCCA
CCGGTCAGAT
AGGGTTTGAT
TTGGTGAGGC
CAAGTATGAG
ACTGTATGCA
CTGGTAACCA

每个序列长 10 个核苷酸。我想计算每个二核苷酸位置的 AA、CC、CG 等的二核苷酸频率。

例如，第一个两列将包含这些碱基 CA、CC、AG、TT 等，第二个两个应该是 AG、CG、GG、TG 等。虚拟输出应该是这样的：

AA   0.25    0.34    0.56    0.43    0.00    0.90    0.45    0.34    0.31    0.01
CC   0.45    0.40    0.90    0.00    0.40    0.90    0.30    0.25    0.2     0.10
GC   0.00    0.00    0.34    1.00    0.30    0.30    0.35    0.90    0.1     0.30
TC   0.24    0.56    0.00    0.00    1.00    0.34    0.45    0.35    0.36    0.45

AA、CC、GC 和 TC 可以是任何二核苷酸组合。

我的尝试：

import itertools
hl = []

#making all possible dinucleotide combinations
bases=['A','T','G','C']
k = 2

kmers = [''.join(p) for p in itertools.product(bases, repeat=k)]
kmers = dict.fromkeys(kmers, 0)
print(kmers)
for i in range(9):
    hl.append(kmers)

# my seq file
f = open("H3K27ac.seq")
nLines = 0

for line in f:
    id = 0
# trying to read dinucleotide at each loop
    for (fst, sec) in zip(line[0::2], line[1::2]):
        id = id + 1
        hl[id][fst+sec] += 1
    nLines += 1
f.close()

nLines = float(nLines)
for char in ['AA', 'TT', 'CC', 'GG', 'CG', 'GC']:   
    print("{}\t{}".format(char, "\t".join(["{:0.2f}".format(x[char]/nLines) for x in hl])))

通过某种方式估计整个文件中所有二核苷酸的总数。我试图让它更通用，以便它可以用于三核苷酸和任何更多数量的核苷酸。

Answer 1

您能否详细说明“每个二核苷酸位置的频率”是什么意思？以下代码不计算任何类型的百分比或频率，但它可能有助于遍历列：

seqs = """CAGGTAGCCA
CCGGTCAGAT
AGGGTTTGAT
TTGGTGAGGC
CAAGTATGAG
ACTGTATGCA
CTGGTAACCA""".splitlines()

def chunkwise(it, n=2):
    from itertools import tee
    its = tee(it, n)
    for index, it in enumerate(its):
        for _ in range(index):
            _ = next(it, None)
    yield from zip(*its)

for col in zip(*map(chunkwise, seqs)):
    print(col)

输出：

(('C', 'A'), ('C', 'C'), ('A', 'G'), ('T', 'T'), ('C', 'A'), ('A', 'C'), ('C', 'T'))
(('A', 'G'), ('C', 'G'), ('G', 'G'), ('T', 'G'), ('A', 'A'), ('C', 'T'), ('T', 'G'))
(('G', 'G'), ('G', 'G'), ('G', 'G'), ('G', 'G'), ('A', 'G'), ('T', 'G'), ('G', 'G'))
(('G', 'T'), ('G', 'T'), ('G', 'T'), ('G', 'T'), ('G', 'T'), ('G', 'T'), ('G', 'T'))
(('T', 'A'), ('T', 'C'), ('T', 'T'), ('T', 'G'), ('T', 'A'), ('T', 'A'), ('T', 'A'))
(('A', 'G'), ('C', 'A'), ('T', 'T'), ('G', 'A'), ('A', 'T'), ('A', 'T'), ('A', 'A'))
(('G', 'C'), ('A', 'G'), ('T', 'G'), ('A', 'G'), ('T', 'G'), ('T', 'G'), ('A', 'C'))
(('C', 'C'), ('G', 'A'), ('G', 'A'), ('G', 'G'), ('G', 'A'), ('G', 'C'), ('C', 'C'))
(('C', 'A'), ('A', 'T'), ('A', 'T'), ('G', 'C'), ('A', 'G'), ('C', 'A'), ('C', 'A'))
>>>

编辑 - 我会使用 collections.Counter 来计算列中的“kmers”，因为您不需要首先构建一个全为零的字典。如果您尝试访问计数器中不存在的键值对，默认情况下它将 return 归零：

from collections import Counter
from itertools import product

seqs = """CAGGTAGCCA
CCGGTCAGAT
AGGGTTTGAT
TTGGTGAGGC
CAAGTATGAG
ACTGTATGCA
CTGGTAACCA""".splitlines()

def chunkwise(it, n=2):
    from itertools import tee
    its = tee(it, n)
    for index, it in enumerate(its):
        for _ in range(index):
            _ = next(it, None)
    yield from zip(*its)

k = 2

kmers = ["".join(p) for p in product("ATGC", repeat=k)]

counters = []
for col in zip(*[chunkwise(seq, n=k) for seq in seqs]):
    counters.append(Counter("".join(chunk) for chunk in col))

for kmer in kmers:
    print("{}   {}".format(kmer, "    ".join("{:.2f}".format(c[kmer]/sum(c.values())) for c in counters)))

输出：

AA   0.00    0.14    0.00    0.00    0.00    0.14    0.00    0.00    0.00
AT   0.00    0.00    0.00    0.00    0.00    0.29    0.00    0.00    0.29
AG   0.14    0.14    0.14    0.00    0.00    0.14    0.29    0.00    0.14
AC   0.14    0.00    0.00    0.00    0.00    0.00    0.14    0.00    0.00
TA   0.00    0.00    0.00    0.00    0.57    0.00    0.00    0.00    0.00
TT   0.14    0.00    0.00    0.00    0.14    0.14    0.00    0.00    0.00
TG   0.00    0.29    0.14    0.00    0.14    0.00    0.43    0.00    0.00
TC   0.00    0.00    0.00    0.00    0.14    0.00    0.00    0.00    0.00
GA   0.00    0.00    0.00    0.00    0.00    0.14    0.00    0.43    0.00
GT   0.00    0.00    0.00    1.00    0.00    0.00    0.00    0.00    0.00
GG   0.00    0.14    0.71    0.00    0.00    0.00    0.00    0.14    0.00
GC   0.00    0.00    0.00    0.00    0.00    0.00    0.14    0.14    0.14
CA   0.29    0.00    0.00    0.00    0.00    0.14    0.00    0.00    0.43
CT   0.14    0.14    0.00    0.00    0.00    0.00    0.00    0.00    0.00
CG   0.00    0.14    0.00    0.00    0.00    0.00    0.00    0.00    0.00
CC   0.14    0.00    0.00    0.00    0.00    0.00    0.00    0.29    0.00
>>>

Answer 2

看看这样的行：

for (fst, sec) in zip(line[0::2], line[1::2]):

这将您锁定为 2 聚体：

>>> l = "CAGGTAGCCA"
>>> for (fst, sec) in zip(l[0::2], l[1::2]): print('{}{}'.format(fst, sec))
... 
CA
GG
TA
GC
CA

如果您想概括任何特定值 k 的 kmer 频率，请查找此类行以及您用于迭代的 start/stop 索引序列，以检索用于计数目的的键。

也许改用 slice 方法来获得更长的 kmers，使用正确的起始偏移和长度值，例如：

>>> k = 3
>>> for x in range(len(l)-k+1): print(l[x:x+k])
... 
CAG
AGG
GGT
GTA
TAG
AGC
GCC
CCA

或：

>>> k = 4
>>> for x in range(len(l)-k+1): print(l[x:x+k])
... 
CAGG
AGGT
GGTA
GTAG
TAGC
AGCC
GCCA

分析或计算生成用作密钥的 kmers 所需的时间可能很有用。一种方法是使用列表和 append ，另一种方法是构造一个 list comprehension.

一种使用 timeit 快速分析的方法，首先指定两个用于构建 kmer 列表的函数：

>>> def f1(l, k):
...     a=[]
...     for x in range(len(l)-k+1): a.append(l[x:x+k])
...     return a

>>> def f2(l, k):
...     return [l[x:x+k] for x in range(len(l)-k+1)]

然后运行来自 timeit class 的 timeit 函数。以下是我系统的一些临时结果：

>>> import timeit
>>> timeit.timeit('f1("CAGGTAGCCA", 3)', globals=globals())
1.3822405610000033
>>> timeit.timeit('f2("CAGGTAGCCA", 3)', globals=globals())
1.51058234300001

我很惊讶列表理解比附加到列表慢，但这就是分析很有用的原因。

计算每列的 kmer 核苷酸频率

calculating kmer nucleotide frequency per column

python

string

string-literals