在for循环中计算“Digram”即核苷酸对

Question

我需要写一个函数，它接受一个fasta文件，并计算文件中的digram（AT、CG、TT、CC等）

我的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

您可以将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 2

一个想法是初始化一个Python dict，将每个2gram映射为零，并根据每一行递增。在这里，我假设FASQ文件只包含“ATGC”中的base。此外，对每一行的每一个可能对进行迭代需要对每一行进行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

Answer 3

把它加到你的最后一次计数中

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