计算每列的 kmer 核苷酸频率

Question

我有一个序列列表：

CAGGTAGCCA
CCGGTCAGAT
AGGGTTTGAT
TTGGTGAGGC
CAAGTATGAG
ACTGTATGCA
CTGGTAACCA

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

例如，第一两列将有这些碱基 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”，因为您不需要先构造一个全为零的字典。如果您尝试访问计数器中不存在的键值对，默认情况下它将返回零：

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-mers 中：

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

如果您想概括 k 的任何特定值的 kmer 频率，请查找该类型的行以及您用于迭代序列的开始/停止索引，以检索使用的键用于计数目的。

也许改为使用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 类运行timeit 函数。以下是我系统的一些临时结果：

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

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