在for环中又称核苷酸对的“ Digram”计数

Question

我需要编写一个函数，该函数接受一个fasta文件并计算该文件中的图（AT，CG，TT，CC等）。

My for循环当前逐行读取文件，并生成该行的计数。然后，它在下一行重新开始计数。 （所有这些都整理成字典）

我想保持每一行的计数，所以我得到了整个文件的计数，而不仅仅是单个行。

这是我要修复的代码：

dinucleotides = ['AA','AT','AG','AC',
                 'TA','TT','TG','TC',
                 'GA','GT','GG','GC',
                 'CA','CT','CG','CT']

all_counts = {}

with open('short.fasta', 'r') as dna_file:
    dna_file.readline()

    for line in dna_file:
        my_line = line.strip()

        for pairs in dinucleotides:
            count = my_line.count(pairs)
            all_counts[pairs] = count

谢谢！

Answer 1

将它添加到您拥有的最后一个计数中，

all_counts[pairs] = all_counts.get(pairs, 0) + count

Answer 2

您可以将collections.defaultdict和int一起用作default_factory。并将all_counts[pairs] = count更改为all_counts[pairs] += count。

from collections import defaultdict


dinucleotides = ['AA','AT','AG','AC',
                 'TA','TT','TG','TC',
                 'GA','GT','GG','GC',
                 'CA','CT','CG','CT']

all_counts = defaultdict(int)

with open('short.fasta', 'r') as dna_file:
    dna_file.readline()

    for line in dna_file:
        my_line = line.strip()

        for pairs in dinucleotides:
            count = my_line.count(pairs)
            all_counts[pairs] += count

或者，使用dict.setdefault方法。

...
all_counts = {}
...
            all_counts.setdefault(pairs, 0) += count

Answer 3

一个想法是初始化一个Python字典，将每个2gram映射为零，然后根据每一行递增该字典。在这里，我将假设FASQ文件仅包含“ ATGC”中的基数。另外，对每条线的每个可能的对进行迭代需要每条线进行16次传递。可以通过以下方式避免这种情况：每行预读一次并保存每对。也许如下：

import random

def random_dnukes(lines=1000, cols=40):
  return [''.join(random.choices('ATGC', k=cols)) for _ in range(lines)]

# e.g.
# ['TGACTCGTCAAAGGTACGTTAATCCTTGGGCAGTTACGGG',
#  'ATTGTTCAATCGAACGTTCGCTACTCGACTCGCGCCCCCT',
#  'TCCCGTGGGACAGGTTCCCAATTGACCGGACGCCGGACTC',
#  'TCGTCGTGCCCCGACATTGCTTCACGGCGGTGCGCGCTGG',
#  'GGTCCGGTCTAGGCGATCCCTAATAGTCAAGCACCGATTA',
#  'CCGAGCCTTGTGTATACTCTGTAAACACTTCTTCCCATAC',
#  'CGGATAGCAGCTAGTGGTTCCCGCAGTACAGGATGACCAA',
#  'CTCGGACGAGAAATCAGGCCAACCTCCACTGGCGACAGAA',
#  'TCTGACCTGCAGTGCAGTCCAGTTATAGTGGAACACCAGC',
#  'GTCAGCCCTTATCCGTTAGCCCAGGTGCCTCAATAGGAGG']

fake_file_iterator = iter(random_dnukes(1000, 40))

from collections import defaultdict
total_counts = defaultdict(int)

for line in fake_file_iterator:
  line = line.strip()
  for i in range(len(line) - 1):
    total_counts[line[i:i+2]] += 1

for k, v in total_counts.items():
  print(k, v)

产生结果

GC 2497
CC 2382
CG 2444
GT 2422
TT 2508
TA 2373
AC 2466
GG 2408
TG 2473
CA 2462
AA 2412
CT 2448
AG 2454
GA 2470
TC 2400
AT 2381

我要补充一点，这可以优化到荒谬的程度。

---

顺便说一句，GC, TT, GG和TG似乎不适合拼写检查，因为它们不是美国的州，股票代码或原子元素（截至目前）。