Skip to main content
Bioinformatics

Return to Answer

added 610 characters in body
Source Link
terdon ♦
terdon ♦
  • 10.6k
  • 5
  • 23
  • 48
#!/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())))


#!/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())))


$ 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

#!/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())))


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

#!/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())))


$ 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
Source Link
terdon ♦
terdon ♦
  • 10.6k
  • 5
  • 23
  • 48

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