Как вычислить одну частотную матрицу из всего файла генома?

Question

Итак, я просто пытаюсь рассчитать одиночные нуклеотидные частоты (A, T, C, G) в HUGE файле, который содержит шаблон, подобный этому: TTTGTATAAGAAAAAATAGG.

Это дало бы мне одну линию вывода всего файла, такие как:

The single nucleotide frequency matrix of T.volcanium Genome is: {'A': [234235], 'C': [234290], 'G': [32456], 'T': [346875]} 

вот мой код (без пути к файлу, открыть, закрыть и главный)

---

def freq_dict_of_lists_v1(dna_list): 
    n = max([len(dna) for dna in dna_list]) 
    frequency_matrix = { 
     'A': [0] * n, 
     'C': [0] * n, 
     'G': [0] * n, 
     'T': [0] * n, 
    } 
    for dna in dna_list: 
     for index, base in enumerate(dna): 
      frequency_matrix[base][index] += 1 

    return frequency_matrix 

for line in file: 
    dna_list = file.readline().rstrip("\n") 
    frequency_matrix = freq_dict_of_lists_v1(dna_list) 
    print("The single nucleotide frequency matrix of T.volcanium Genome is: ") 
    pprint.pprint(frequency_matrix) 

---

И это мой выход.

The single nucleotide frequency matrix of T.volcanium Genome is: 
{'A': [21], 'C': [10], 'G': [11], 'T': [18]} 
The single nucleotide frequency matrix of T.volcanium Genome is: 
{'A': [31], 'C': [6], 'G': [4], 'T': [19]} 
The single nucleotide frequency matrix of T.volcanium Genome is: 
{'A': [23], 'C': [9], 'G': [10], 'T': [18]} 
The single nucleotide frequency matrix of T.volcanium Genome is: 
{'A': [17], 'C': [8], 'G': [9], 'T': [26]} 
The single nucleotide frequency matrix of T.volcanium Genome is: 
{'A': [15], 'C': [13], 'G': [9], 'T': [23]} 
The single nucleotide frequency matrix of T.volcanium Genome is: 
{'A': [21], 'C': [12], 'G': [10], 'T': [17]} 
The single nucleotide frequency matrix of T.volcanium Genome is: 
{'A': [20], 'C': [9], 'G': [12], 'T': [19]} 
The single nucleotide frequency matrix of T.volcanium Genome is: 
{'A': [15], 'C': [15], 'G': [10], 'T': [20]} 
The single nucleotide frequency matrix of T.volcanium Genome is: 
{'A': [20], 'C': [11], 'G': [10], 'T': [19]} 
The single nucleotide frequency matrix of T.volcanium Genome is: 
{'A': [26], 'C': [13], 'G': [7], 'T': [14]} 
The single nucleotide frequency matrix of T.volcanium Genome is: 
{'A': [12], 'C': [13], 'G': [13], 'T': [22]} 
The single nucleotide frequency matrix of T.volcanium Genome is: 
{'A': [20], 'C': [16], 'G': [9], 'T': [15]} 
The single nucleotide frequency matrix of T.volcanium Genome is: 
{'A': [22], 'C': [12], 'G': [6], 'T': [20]} 

Поэтому он вычисляет его по строке. Я попытался вынуть цикл for или снять строки чтения, но тогда он даст мне только одну строку вывода только для одной строки в файле. не весь файл.

Я чувствую, что переусердствовал. Я уверен, что есть простой способ прочитать весь файл и распечатать одну строчку вывода с использованием всех частот ... Любое понимание очень ценится.

Answer 1

Я вижу две проблемы с вашим решением.

1. Вы отслеживаете базисы в положение, когда в вашем вопросе он говорит, вы хотите, чтобы следить за подсчеты по всем направлениям
2. Вы вызываете функцию только один раз в каждую строку.

Мое исправление ниже должно быть адресовано. Смотрите комментарии для объяснения

def freq_dict_of_lists_v1(dna_list): 
    frequency_matrix = { # We are only keeping one variable per base 
     'A': [0],   # so that we calculate counts across all lines 
     'C': [0], 
     'G': [0], 
     'T': [0], 
    } 
    for dna in dna_list: 
     for base in dna: # No longer need index, so I removed enumerate 
      frequency_matrix[base] += 1 # Change here since dict structure changed 

    return frequency_matrix 

# Unlike before, we are now appending all the lines into dna_list 
for line in file: 
    dna_list.append(file.readline().rstrip("\n")) 

# Calling freq_dict_of_lists_v1 on ALL the lines at once (it is now out of loop) 
frequency_matrix = freq_dict_of_lists_v1(dna_list) 
print("The single nucleotide frequency matrix of T.volcanium Genome is: ") 
pprint.pprint(frequency_matrix) 

Одно предостережение для этого решения является то, чтобы убедиться, что все базы в файле, верхний регистр. Кроме того, убедитесь, что нет символов, отличных от ACGT (некоторые последовательности имеют специальные символы пробела и т. Д.). Если это так, есть другие символы, вы можете обратиться к this thread, где ваша запись по умолчанию может быть чем-то вроде Gap.

Answer 2

Не знаете, что ОГРОМНОЕ означает МБ? GB ?, но это самое простое решение. Однако учтите, что он загружает весь файл в память.