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()
No comments:
Post a Comment