0
$\begingroup$

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.

$\endgroup$
10
  • 1
    $\begingroup$ Hi and welcome to the site! We need more detail to be able to help you, so please edit your question and i) add a few lines of your fasta.dataset file; maybe 3-4 protein entries so we can see what we are working with since what you describe is not fasta; ii) show us the exact output you would want to see from that example input; iii) show us the code you wrote so we don't need to start from 0; iv) explain exactly how it "did not work"; v) tell us what operating system you are working with. Note that since this is not fasta, I doubt there is any reason to use biopython. $\endgroup$
    – terdon
    Feb 7, 2023 at 13:57
  • $\begingroup$ Oh, and please use the formatting tools to format your input and output examples and the code you show as code (the {} button). $\endgroup$
    – terdon
    Feb 7, 2023 at 14:01
  • $\begingroup$ Thanks for the edit, but please don't post images of text: we can't use those to test our solutions. We need a few lines of your actual file and the output you want to see from those lines. Both as text so we can copy/paste them. $\endgroup$
    – terdon
    Feb 7, 2023 at 15:09
  • $\begingroup$ And what output do you want from this? Just the number 247? I ask you again, for the third time, please show us the output you are expecting. And does it have to be in python? I mean, all you want here is awk 'NR%4==2' file, is that enough or do you absolutely need a whole python script? $\endgroup$
    – terdon
    Feb 7, 2023 at 16:11
  • $\begingroup$ awk 'NR%4==3' file is calculated just the first sequence while I have half a million of sequences in the same dataset $\endgroup$
    – Amal
    Feb 7, 2023 at 16:39

2 Answers 2

0
$\begingroup$

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.

$\endgroup$
4
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.