open and read Fasta file (raw data)

Question

I have a big fasta.dataset file containing half a million proteins (1.0 GB). I have four lines for each protein code:

• line 1:the protein code
• line 2: protein length in amino acids
• line 3: amino acid sequence
• line 4: secondary structure

Now, I am trying to open and read it in python (Biopython), and it does not work:

filename = 'pdb.fasta_qual.dataset'
sequences = SeqIO.parse ( filename,'fasta')
for record in sequences:
    example = record
    break
print(example)

How can I read it in python and loop through the file to look at line 3 for each protein to count the sequence length and distribution?

here is the first 5 lines of my file: which my file contains 500,000 proteins for each one has a 4 lines (name ,len of protein in amino acid,the seq represents by letters which what I would to calculate,the secondary structure)

4LGTD
247
M       S       E       K       L       Q       K       V       L       A       R       A       G       H       G      T
.       .       E       E       H       H       H       H       H       H       H       T       T       S       S       .

I want to open and read the file and loop through line 3 for each protein to calculate the length of the sequences and plot a histogram ! to check the distribution. ,,,

The output i am expecting is :

The len for the first seq is =

The len for the second seq is =

Until

the len of the last sequence which is number (500.000)=

And then i can plot a histogram for the len of the sequences ,,,

NOTE: I have opened and read the file's info by Linux, but I could not by python.

Answer 1

I better try and explain the issues a bit better. A dataset of 0.5 million at 1 GB will struggle to load into Biopython. Biopython isn't great at handling large numbers of sequences - for a regular desktop/laptop - and what will happen is the parser will freeze or else take an enormous amount of time (many hours). I call this a "RAM bottleneck". A generator is an easy workaround.

The other issue with the code is there are lots of bugs combining format and code. There is a clarity issue as well it looks like you want to calculate the secondary structure - rather than import the secondary structure alongside the alignment. I could be wrong on this though, certainly secondary structure is present in the ad hoc alignment. There is an alignment format that handles both secondary structure and sequence data. You might think about restructuring the question and your approach into multiple questions.

1. Reformatting the ad hoc alignment to permit import into Biopython on a pilot data set (not 0.5 million sequences). This is close to @terdon's answer.
2. Code that will process very large files in Biopython. This is the answer below.
3. An alignment format via Biopython (probably output - could be input) that integrates linear sequence with secondary structure.
4. Plotting a histogram for 0.5 million sequences.

Points (Questions) 1. to 3. can be leveraged by Biopython. Point (Question) 4 can be done via Python too - graphing/histograms are pretty good. There seems an underlying issue of project management for coding based work.

If this is a RAM bottleneck then solution is below.

---

The easiest approach is to use a generator (please see 'format').

example = (record for record in SeqIO.parse(filename,'fasta'))
print (next(example)) # example.id ?

A generator will deal with the sequences one by one instead of trying to load the whole thing to memory. It may not be suitable for certain complex calculations, but it is by far the easiest approach to the problem (if the problem is the OP is struggling to load a million sequences into memory).

Note I always used yield commands but @Steve introduced the () format and I definitely prefer it. Also note that manipulating the output of a generator can be different to a list, tuple or dictionary (calculation depending).

---

Format The format supplied isn't fasta so the SeqIO will not work in fasta mode and should throw an error/exception. The format you need to leverage Biopython is very different. I would recommend a standard 'open' command and passing this to a custom dictionary. Also the output of print(next(record)) should produce an object symbol, not the actual data e.g. record.id would be needed to see data: it does however demonstrate the generator is working.

Answer 2

It sounds like all you really want is to print every 2nd line of a text file. If so, you don't need Python, let alone BioPython, you can do it with basic *nix tools:

$ awk 'NR%4==2' pdb.fasta_qual.dataset 
247

NR is the current line number, and % is the modulo operator. Therefore, NR % 4 will equal 2 on the second line of each group of 4 lines in the file. In awk, the default action when something evaluates to true is to print the current line, so this will print out all length lines from your input file (assuming you don't have blank lines anywhere, which we don't know since you only showed the first entry).

If you want to add some sort of message to be printed, you could do:

$ awk 'NR%4==2{ print "The length of protein",++c,"is:",$0}' pdb.fasta_qual.dataset 
The length of protein 1 is: 247

If you really must do this in Python, you could do something like this:

#!/usr/bin/env python3

import sys

#my_file= sys.argv[1]
line_number = 0
with open(sys.argv[1]) as f:
    for line in f:
        line_number += 1
        if (line_number % 4 == 2):
            print("The length of protein %d is: %d" % (line_number,int(line.strip())))

And then you would run it like this (assuming you named the script foo.py):

$ python3 foo.py pdb.fasta_qual.dataset 
The length of protein 2 is: 247

---

Out of curiosity, I created a test file by repeating your example 10000000 times, which left me with a 2.4G file:

$ perl -e 'open(A,"pdb.fasta_qual.dataset"); @lines=<A>; print "@lines" x 10000000' > testFile
$ ls -lh testFile 
-rw-r--r-- 1 terdon terdon 2.4G Feb  7 17:06 testFile

I then tested both approaches for speed. Unsurprisingly, awk was much faster:

$ time awk 'NR%4==2' testFile > /dev/null 

real    0m6.447s
user    0m5.917s
sys     0m0.492s

$ time foo.py testFile > /dev/null 

real    0m10.856s
user    0m10.111s
sys     0m0.523s