<p>I better try and explain the issues a bit better. A dataset of 0.5 million at 1 GB will struggle to load into <code>Biopython</code>. <code>Biopython</code> 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 &quot;RAM bottleneck&quot;. A generator is an easy workaround.</p>
<p>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 <em>ad hoc</em> 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.</p>
<ol>
<li>Reformatting the <em>ad hoc</em> alignment to permit import into <code>Biopython</code> on a <em>pilot data set</em> (not 0.5 million sequences). This is close to @terdon's answer.</li>
<li>Code that will process very large files in <code>Biopython</code>. This is the answer below.</li>
<li>An alignment format via <code>Biopython</code> (probably output - could be input) that integrates linear sequence with secondary structure.</li>
<li>Plotting a histogram for 0.5 million sequences.</li>
</ol>
<p>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.</p>
<p>If this is a RAM bottleneck then solution is below.</p>
<hr />
<p>The easiest approach is to use a generator (please see 'format').</p>
<pre><code>example = (record for record in SeqIO.parse(filename,'fasta'))
print (next(example)) # example.id ?
</code></pre>
<p>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).</p>
<p><strong>Note</strong> I always used <code>yield</code> commands but @Steve introduced the <code>()</code> 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).</p>
<hr />
<p><strong>Format</strong>
The format supplied isn't fasta so the <code>SeqIO</code> will not work in <code>fasta</code> mode and should throw an error/exception. The format you need to leverage <code>Biopython</code> is very different. I would recommend a standard 'open' command and passing this to a custom dictionary. Also the output of <code>print(next(record))</code> should produce an object symbol, not the actual data e.g. <code>record.id</code> would be needed to see data: it does however demonstrate the generator is working.</p>

CC BY-SA 4.0
2023-02-07T14:34:12.377Z
2023-02-08T08:58:51.587Z
2023-02-08T08:58:51.587Z
M__
<p>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:</p>
<pre><code>$ awk 'NR%4==2' pdb.fasta_qual.dataset 
247
</code></pre>
<p><code>NR</code> is the current line number, and <code>%</code> is the <a href="https://en.wikipedia.org/wiki/Modulo" rel="nofollow noreferrer">modulo</a> operator. Therefore, <code>NR % 4</code> will equal <code>2</code> on the second line of each group of 4 lines in the file. In <code>awk</code>, 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).</p>
<p>If you want to add some sort of message to be printed, you could do:</p>
<pre><code>$ awk 'NR%4==2{ print &quot;The length of protein&quot;,++c,&quot;is:&quot;,$0}' pdb.fasta_qual.dataset 
The length of protein 1 is: 247
</code></pre>
<p>If you really must do this in Python, you could do something like this:</p>
<pre class="lang-py prettyprint-override"><code>#!/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(&quot;The length of protein %d is: %d&quot; % (line_number,int(line.strip())))

</code></pre>
<p>And then you would run it like this (assuming you named the script <code>foo.py</code>):</p>
<pre><code>$ python3 foo.py pdb.fasta_qual.dataset 
The length of protein 2 is: 247
</code></pre>
<hr />
<p>Out of curiosity, I created a test file by repeating your example 10000000 times, which left me with a 2.4G file:</p>
<pre><code>$ perl -e 'open(A,&quot;pdb.fasta_qual.dataset&quot;); @lines=&lt;A&gt;; print &quot;@lines&quot; x 10000000' &gt; testFile
$ ls -lh testFile 
-rw-r--r-- 1 terdon terdon 2.4G Feb  7 17:06 testFile
</code></pre>
<p>I then tested both approaches for speed. Unsurprisingly, <code>awk</code> was  much faster:</p>
<pre class="lang-none prettyprint-override"><code>$ time awk 'NR%4==2' testFile &gt; /dev/null 

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

$ time foo.py testFile &gt; /dev/null 

real    0m10.856s
user    0m10.111s
sys     0m0.523s
</code></pre>

CC BY-SA 4.0
2023-02-07T17:00:58.623Z
2023-02-07T17:08:16.537Z
2023-02-07T17:08:16.537Z
terdon
<p>I have a big <code>fasta.dataset</code> file containing half a million proteins (1.0 GB). I have four lines for each protein code:</p>
<ul>
<li>line 1:the protein code</li>
<li>line 2: protein length in amino acids</li>
<li>line 3: amino acid sequence</li>
<li>line 4: secondary structure</li>
</ul>
<p>Now, I am trying to open and read it in python (Biopython), and it does not work:</p>
<pre><code>filename = 'pdb.fasta_qual.dataset'
sequences = SeqIO.parse ( filename,'fasta')
for record in sequences:
    example = record
    break
print(example)
</code></pre>
<p>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?</p>
<p>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)</p>
<pre><code>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       .
</code></pre>
<p>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.
,,,</p>
<p>The output i am expecting is :</p>
<p>The len for the first seq is  =</p>
<p>The len for the second seq is  =</p>
<p>Until</p>
<p>the len of the last sequence  which is number (500.000)=</p>
<p>And then i can plot a histogram for the len of the sequences
,,,</p>
<p>NOTE: I have opened and read the file's info by Linux, but I could not by python.</p>

CC BY-SA 4.0
2023-02-07T13:57:42.503Z
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.
terdon
CC BY-SA 4.0
2023-02-07T14:01:04.507Z
Oh, and please use the [formatting tools](https://bioinformatics.stackexchange.com/help/formatting) to format your input and output examples and the code you show as code (the `{}` button).
terdon
CC BY-SA 4.0
2023-02-07T15:09:24.823Z
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.
terdon
CC BY-SA 4.0
2023-02-07T16:11:48.370Z
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?
terdon
CC BY-SA 4.0
2023-02-07T16:39:27.180Z
awk 'NR%4==3' file is calculated just the first sequence while I have half a million of sequences in the same dataset
Amal
CC BY-SA 4.0
2023-02-07T16:46:01.863Z
no, that will print every line in the file whose line number modulo 4 is equal to 3, so it will print the 2nd, 6th, 10th, 14th and so on. But are you now saying you want some string along with the length? Should the solutions print `The length for the Nth seq is` and then the length?
terdon
CC BY-SA 4.0
2023-02-07T17:05:07.050Z
Sorry before we continue we need to establish why it's not working. Is it the file size is too big and it just freezes or is it throwing a bug? There are bugs in the code, but passing this via Biopython will likely hit a RAM bottleneck
M__
CC BY-SA 4.0
2023-02-08T09:15:14.363Z
This ```awk 'NR%5==3{ print length ($0);}' file name ``` works and gives me the len of characters in the sequences line 3. Now how can I plot a histogram for it? Also calculate the frequency of the same line?
Amal
CC BY-SA 4.0
2023-02-19T21:28:42.107Z
Based on your description, your file is not a FASTA file. I suggest you modify your question. Here is the description of FASTA file> `https://www.ncbi.nlm.nih.gov/genbank/fastaformat/`. Perhaps that is why it's not working.
Supertech
CC BY-SA 4.0
2023-02-19T21:32:59.860Z
Just read every third line in the file and do sth with it.
Supertech
<fasta><proteins><protein-structure><biopython>
open and read Fasta file (raw data)
# open and read Fasta file (raw data)

(Asked by: username_0 on 2023-02-07)

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.

**Tags:** fasta, proteins, protein-structure, biopython
> **Comments:**
> **username_1**: 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.  
> **username_1**: Oh, and please use the [formatting tools](https://bioinformatics.stackexchange.com/help/formatting) to format your input and output examples and the code you show as code (the `{}` button).  
> **username_1**: 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.  
> **username_1**: 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?  
> **username_0**: awk 'NR%4==3' file is calculated just the first sequence while I have half a million of sequences in the same dataset  
> **username_1**: no, that will print every line in the file whose line number modulo 4 is equal to 3, so it will print the 2nd, 6th, 10th, 14th and so on. But are you now saying you want some string along with the length? Should the solutions print `The length for the Nth seq is` and then the length?  
> **username_2**: Sorry before we continue we need to establish why it's not working. Is it the file size is too big and it just freezes or is it throwing a bug? There are bugs in the code, but passing this via Biopython will likely hit a RAM bottleneck  
> **username_0**: This ```awk 'NR%5==3{ print length ($0);}' file name ``` works and gives me the len of characters in the sequences line 3. Now how can I plot a histogram for it? Also calculate the frequency of the same line?  
> **username_3**: Based on your description, your file is not a FASTA file. I suggest you modify your question. Here is the description of FASTA file> `https://www.ncbi.nlm.nih.gov/genbank/fastaformat/`. Perhaps that is why it's not working.  
> **username_3**: Just read every third line in the file and do sth with it.  

---

## 2 answers

### Answer by username_2 on 2023-02-07. Score: 0 | Accepted answer

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.

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.
Code that will process very large files in Biopython. This is the answer below.
An alignment format via Biopython (probably output - could be input) that integrates linear sequence with secondary structure.
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 by username_1 on 2023-02-07. Score: 4

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=; 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


---

