Monday, November 14, 2011

Files Handling

the clearest and most efficient way to read FASTA file:

fh=open('seqA.fa')
name=fh.readline()[1:-1] #between '>' and '\n' is name of the sequence.
sequence=''
for line in fh:
    sequence+=fh.readline()[:-1]
fh.close()

It is better to use readline() than read(). read() will read the whole file at once and store it as a string. This i a potential problem if there is not enough memory available to store the file's content. (Genome sequence files are usually in MB or GB size.)

Some modules in Biopython package allow us to do more complicated operations on sequence file if we need more information than the plain sequences, like the name, id, description, and cross references and annotations.
Here are some important modules:
SeqRecord
SeqIO , which is common interface to input and output sequence file formats.
fh=open('seqA.fa')
for record in SeqIO.parse(fh,'fasta'):
     id=record.id
     seq=record.seq
     print id,'size:',len(seq)
fh.close()


Data types

New in Python 3: Dictionary
In python 3, keys(),values() and items() methods don't return a list as in previous versions. They return an iterable object. i.e. IUPAC={'A': 'Ala', 'C': 'Cys','E': 'Glu'}
>IUPAC.values() # under python 3.0, you will get:
<dict_values object at 0xb7x8s7s2>
>for x in IUPAC.values():
        print x
Ala
Cys
Glu

This is more memory efficient than Python 2.x, since Python doesn't store the whole list in memory.
If you need to get the list anyway, use list()
>list(IUPAC.values())
['Ala', 'Cys', 'Glu']