Apply a formula to a line of data where it crosses multiple lines

Question

I'm just making my first steps to try to learn a bit of Python. Currently working my way through the Rosalind online course which aims to teach bioinformatics python skills. (very good by the way, see: rosalind.info)

I am struggling with one particular problem. I have a file in FASTA format which has the form thus:

>Sequence_Header_1
ACGTACGTACGTACGTACGT
ACGTACGTACGTACGTACGT
>Sequence_Header_2
ACGTACGTACGTACGTACGT
ACGTACGTACGTACGTACGT

I need to calculate the percentage of G and C in each entry of the file (excluding the headers) and return this number, example:

>Sequence_Header_1
48.75%
>Sequence_header_2
52.43%

My code so far is:

file = open("input.txt" , "r")
for line in file:
    if line.startswith(">"):
        print(line.rstrip())        
    else:
        print ('%3.2f' % (line.count('G')+line.count('C')/len(line)*100))
file.close()

Which is doing almost what I need it to do. I am just having trouble where the sequence data crosses multiple lines. At the moment I get the % GC content for every line in the file rather than returning a single figure for each entry, example:

>Sequence_Header_1
48.75%
52.65%
>Sequence_header_2
52.43%
50.25%

How can I apply my formula to the data which crosses multiple lines?

Thanks in advance,

Answer 1

Not really a direct answer to your question, but I think it is a better way to go! If you are planning on doing more bioinformatics in python, have a look at biopython. It will handle fasta files and other common sequence manipulations (and a whole lot more!) for you.

so for example:

from Bio import SeqIO
from Bio.SeqUtils import GC

for rec in SeqIO.parse("input.txt", "fasta"):
    print rec.id,GC(rec.seq)

Answer 2

I think this is what you are looking for:

# Read the input file
with open("input.txt", "r") as f:
    s = f.read()

# Loop through the Sequences
for i, b in enumerate(s.split("Sequence_Header_")):
    if not b: continue # At least the first one is empty 
                       # because there is no data before the first sequence
    # Join the lines with data
    data = "".join(b.split("\n")[1:])

    # Print the result
    print("Sequence_Header_{i}\n{p:.2f}%".format(
        i=i, p=(data.count('G')+data.count('C'))/len(data)*100))

Note: I can't find the '>' sign in your example. If your header is starting with > then you can re-write the code to s.split(">") and the code should still be fine.

Answer 3

Try keeping a running tally, and then having this tally reset when a new header is found.

count = 0.0
line_length=0.0
seen_header = False
for line in open("input.txt" , "r"): #saves memory.
    if line.startswith('>'):
        if not seen_header:
            header = line.strip()
            seen_header=True
        if line_length > 0.0:
            print header,'\n',('%3.2f' % (100.0*count/line_length))
            count = 0.0
            line_length=0.0
            seen_header = False
    else:
        count += line.count('C') + line.count('C')
        line_length += len(line.strip())
print header,'\n',('%3.2f' % (100.0*count/line_length))

Also be careful about division in python, remember that the default is integer division. i.e. 5/2 = 2. You can avoid this by using decimals in your variables, or float().

edit: made nicer, also should have been line_length+=len(line.strip()), to avoid counting the newline "\n" as two characters.

Answer 4

It might not be possible to hold the entire file in memory. Assuming you cannot do s = f.read() all at once, you need to keep a running count of the number of letters and total letters until a new sequence starts. Something like this:

file = open("input.txt" , "r")
# for keeping count:
char_count = 0
char_total = 0
for line in file:
    if line.startswith(">"):
        if char_total != 0:
            # starting a new sequence, so calculate pct for last sequence
            char_pct = (char_count / char_total) * 100
            print ('%3.2f' % char_pct)
            # reset the count for the next sequence
            char_total = 0
            char_count = 0
        print(line.rstrip())        
    else:
        # continuing a sequence, so add to running counts
        char_count += line.count('G') + line.count('C')
        char_total += len(line)
file.close()

Answer 5

You can parse the fasta format and create a dictionary with the >ID as keys and the sequences as values, like so:

    from collections import defaultdict

    def parse_fasta(dataset):
        "Parses data in FASTA format, returning a dictionary of ID's and values"
        records = defaultdict(str)
        record_id = None
        for line in [l.strip() for l in dataset.splitlines()]:
            if line.startswith('>'):
                record_id = line[1:]
            else:
                records[record_id] += line
        return records

or you can rewrite this code a little and create a tuple / list. I prefer the dictionary, since it is already indexed. You can find me on the Rosalind site, if you still need help.