Getting the frequency of letters in each position

Question

I have a text file like this example:

>chr12:86512-86521
CGGCCAAAG
>chr16:96990-96999
CTTTCATTT
>chr16:97016-97025
TTTTGATTA
>chr16:97068-97077
ATTTAGGGA

This file is divided into different parts, and every part has 2 lines. The line which starts with > is ID and the 2nd line is a sequence of letters and the letters are A, T, C or G and also the length of each sequence is 9 so, for every sequence of letters there are 9 positions. I want to get the frequency of the 4 mentioned letters in every position (we have 9 positions).

Here is the expected output for the small example:

one = {'T': 1, 'A': 1, 'C': 2, 'G': 0}
two = {'T': 3, 'A': 0, 'C': 0, 'G': 1}
three = {'T': 3, 'A': 0, 'C': 0, 'G': 1}
four = {'T': 3, 'A': 0, 'C': 1, 'G': 0}
five = {'T': 0, 'A': 1, 'C': 2, 'G': 1}
six  = {'T': 0, 'A': 3, 'C': 0, 'G': 1}
seven = {'T': 2, 'A': 1, 'C': 0, 'G': 1}
eight = {'T': 2, 'A': 1, 'C': 0, 'G': 1}
nine = {'T': 1, 'A': 2, 'C': 0, 'G': 1}

I am doing that in Python using the following command. This command has 3 steps. Steps 1 and 2 work fine, but would you help me to improve step 3, which made this pipeline so slow for big files?

Step 1: to parse the file into a comma-separated file

def fasta_to_textfile(filename, outfile): 
    with open(filename) as f, open(outfile, 'w') as outfile:
        header = sequence = None
        out = csv.writer(outfile, delimiter=',')
        for line in f:
            if line.startswith('>'):
                if header:
                    entry = header + [''.join(sequence)]
                    out.writerow(entry)
                header = line.strip('>\n').split('|')
                sequence = []
            else:
                sequence.append(line.strip())
        if header:
            entry = header + [''.join(sequence)]
            out.writerow(entry)

Step 2: comma-separated file to a Python dictionary

def file_to_dict(filename):
    f = open(filename, 'r')
    answer = {}
    for line in f:
        k, v = line.strip().split(',')
        answer[k.strip()] = v.strip()
    return answer

To print functions from steps 1 and 2:

a = fasta_to_textfile('infile.txt', 'out.txt')
d = file_to_dict('out.txt')

Step 3: to get the frequency

one=[]
two=[]
three=[]
four=[]
five=[]
six=[]
seven=[]
eight=[]
nine=[]
mylist = d.values()
for seq in mylist:
    one.append(seq[0])
    two.append(seq[1])
    se.append(seq[2])
    four.append(seq[3])
    five.append(seq[4])
    six.append(seq[5])
    seven.append(seq[6])
    eight.append(seq[7])
    nine.append(seq[8])

from collections import Counter
one=Counter(one)    
two=Counter(two)
three=Counter(three)
four=Counter(four)
five=Counter(five)

Answer 1

That's a crazy amount of code for such a simple task. There's no need to generate any intermediate files. There's no need to create nine separate variables for the nine columns.

To process an input file (whether from sys.stdin or as a named argument on the command line), use fileinput. That way, you don't have to hard-code 'infile.txt'. Then, simply ignore the lines that start with > and strip off the newlines.

To work on columns rather than rows, use zip().

Suggested solution

These seven lines could replace your entire code.

import fileinput
from collections import Counter

def nucleotide_lines(f):
    for line in f:
        if not line.startswith('>'):
            yield line.rstrip()

print([Counter(col) for col in zip(*nucleotide_lines(fileinput.input()))])

Answer 2

You forgot to add import csv and from collections import Counter. Probably missed it while copy pasting. Also, your = signs are inconsistent in step 3. Try to follow PEP8. Also, a is useless in this line:

a = fasta_to_textfile('infile.txt', 'out.txt')

Since you've programmed a void function, a = None because it returns nothing.

Is the conversion to the CSV file really necessary? This would be an example of the pipeline:

1. Read the file.
2. Extract the sequence and load it into a N*9 table, where N is the number of sequences
3. Swap the rows and columns (numpy can help you out here)
4. A simple for loop that uses the Counter function on each row (but really column), refactored into less lines. Sadly I don't have time right to rewrite bits of your code right now.

One last thing - are you sure your example is correct? I tried loading it but got: ValueError: not enough values to unpack (expected 2, got 1)...