Посчитать число 20-ти в фасте по питону

Question

Обычный файл fasta с длиной чтения 120 nt: 'single_mapped.fa'

CSV-файл содержит 10000 20-метровых данных и число для каждого 20-метрового файла: "20frequent_20mers.txt", например:

AAAAAGTATAGGAGATAGAA    35
AAAAATAGGAGGACTATTCA    26
AAAAATAGGAGGACTATTTA    24
AAAAATAGGAGGCCTATTCA    62

Я хочу пройти через single_mapped.fa, рассчитать накопленные значения всех 20-метров в 20frequent_20mers.txt для каждого чтения, то есть для чтения:

AAAAAGTATAGGAGATAGAA AAAAATAGGAGGACTATTCA, я хочу 61 (35+26)

мой код:

file2 = open('20frequent_20mers.txt','r')
kmer_list = csv.reader(file2, delimiter='\t')

for seq_record in SeqIO.parse("single_mapped.fa", "fasta"):
    print(seq_record.id)
    score_fre = 0
    sequence_string = str(seq_record.seq)
    for i in range(0,101):
            seq = sequence_string[i:i+20]
            for row in kmer_list:
                if row[0] == seq:
                    score_fre = score_fre + int(row[1])            
    print(score_fre)

Каждый цикл работает хорошо, когда я запускаю их отдельно, но не работает, как указано выше, кто-нибудь может сказать мне, откуда ошибки? или если есть более умный и эффективный способ сделать это? Заранее спасибо!

Answer 1

Имея такой код, вам нужно будет перечитать файл kmer с самого начала для каждой последовательности и i значение. Это будет очень медленно и его следует избегать. Поскольку вы не перемещаете указатель файла обратно на начало, он будет работать только один раз.

Указатель файла можно переместить, добавив перед for row in kmer_list: линия:

file2.seek(0)

Гораздо лучшим подходом было бы сначала загрузить все ваши записи на kmer вместе с соответствующим количеством. Таким образом, их можно быстро найти:

import csv

kmers = {}

with open('20frequent_20mers.txt') as f_kmers:
    for kmer, count in csv.reader(f_kmers, delimiter='\t'):
        kmers[kmer] = int(count)

for seq_record in SeqIO.parse("single_mapped.fa", "fasta"):
    print(seq_record.id)
    score_fre = 0
    sequence_string = str(seq_record.seq)

    for i in range(0, 101):
        seq = sequence_string[i:i+20]
        score_fre += kmers.get(seq, 0)

    print(score_fre) 

Если seq не найден в словаре, значение по умолчанию 0 возвращается

Answer 2

Альтернативная реализация (не обязательно лучше и быстрее) со словарем @MartinEvans, но с использованием re.findall() генерировать Kmers для тестирования и использования map а также sum вместо (явного) внутреннего цикла:

from Bio import SeqIO
from re import findall
from itertools import repeat

kmers = {}

with open('20frequent_20mers.txt') as f_kmers:
    for line in f_kmers:
        kmer, count = line.strip().split('\t')
        kmers[kmer] = int(count)

for seq_record in SeqIO.parse("single_mapped.fa", "fasta"):
    print(seq_record.id)

    # use forward lookahead to make findall() find overlapping results;

    score_fre = sum(map(kmers.get, findall(r'(?=([ACTG]{20}))', str(seq_record.seq)), repeat(0)))

    print(score_fre)