Mercurial > repos > jasper > dna_protein_transle

comparison translate.py @ 0:a55eb751d2d8 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author jasper
date Tue, 09 May 2017 13:09:07 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:a55eb751d2d8
1 import sys
2 import warnings
3 import argparse
4 from Bio import SeqIO, SeqRecord, BiopythonWarning
5 from Bio.Seq import Seq
6 from Bio.Alphabet import IUPAC
7
8
9 warnings.simplefilter("ignore", BiopythonWarning)
10 parser = argparse.ArgumentParser(description="Translation")
11 parser.add_argument('-r', action='store', dest='read_frame', type=int)
12 parser.add_argument('-i', action='store', dest='input_file')
13 parser.add_argument('-o', action='store', dest='output_file')
14
15 inputs = parser.parse_args()
16 with open(inputs.output_file, 'w') as f:
17 for seq_record in SeqIO.parse(inputs.input_file, 'fasta'):
18 seqid = seq_record.id
19 name=seq_record.name
20 description = seq_record.description
21 if inputs.read_frame > 0:
22 seq_record = seq_record[(inputs.read_frame - 1):]
23 elif inputs.read_frame == -1:
24 seq_record = seq_record.reverse_complement()
25 elif inputs.read_frame == -2:
26 seq_record = seq_record.reverse_complement()
27 seq_record = seq_record[2:]
28 elif inputs.read_frame == -3:
29 seq_record = seq_record.reverse_complement()
30 seq_record = seq_record[1:]
31 dna = Seq(str(seq_record.seq), IUPAC.ambiguous_dna)
32 protein = SeqRecord.SeqRecord(seq=dna.translate(), id=seqid, name=name, description=description)
33 SeqIO.write(protein, f, 'fasta')
34 f.write('\n')
35 f.close()