посчитать число 20-ти в фасте по питону - PullRequest
Новая
       6

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

0 голосов
/ 15 мая 2018

Обычный файл 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, Iхочу иметь 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)

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

Ответы [ 2 ]

0 голосов
/ 28 июля 2018

Альтернативная реализация (не обязательно лучше и быстрее) со словарем @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)
0 голосов
/ 15 мая 2018

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

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

file2.seek(0)

Aгораздо лучше было бы сначала загрузить все ваши записи на 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.

...