annotate tools/get_orfs_or_cdss/get_orfs_or_cdss.py @ 8:09a8be9247ca draft

v0.2.0 with GFF3 output
author peterjc
date Sat, 09 Jan 2016 23:42:32 -0500
parents 705a2e2df7fb
children a06ad07431ba
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
1 #!/usr/bin/env python
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
2 """Find ORFs in a nucleotide sequence file.
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
3
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
4 For more details, see the help text and argument descriptions in the
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
5 accompanying get_orfs_or_cdss.xml file which defines a Galaxy interface.
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
6
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
7 This tool is a short Python script which requires Biopython. If you use
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
8 this tool in scientific work leading to a publication, please cite the
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
9 Biopython application note:
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
10
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
11 Cock et al 2009. Biopython: freely available Python tools for computational
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
12 molecular biology and bioinformatics. Bioinformatics 25(11) 1422-3.
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
13 http://dx.doi.org/10.1093/bioinformatics/btp163 pmid:19304878.
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
14
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
15 This script is copyright 2011-2013 by Peter Cock, The James Hutton Institute
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
16 (formerly SCRI), Dundee, UK. All rights reserved.
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
17
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
18 See accompanying text file for licence details (MIT licence).
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
19 """
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
20 import sys
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
21 import re
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
22 from optparse import OptionParser
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
23
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
24 usage = """Use as follows:
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
25
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
26 $ python get_orfs_or_cdss.py -i genome.fa -f fasta --table 11 -t CDS -e open -m all -s both --on cds.nuc.fa --op cds.protein.fa --ob cds.bed --og cds.gff3
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
27 """
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
28
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
29 try:
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
30 from Bio.Seq import Seq, reverse_complement, translate
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
31 from Bio.SeqRecord import SeqRecord
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
32 from Bio import SeqIO
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
33 from Bio.Data import CodonTable
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
34 except ImportError:
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
35 sys.exit("Missing Biopython library")
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
36
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
37
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
38 parser = OptionParser(usage=usage)
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
39 parser.add_option('-i', '--input', dest='input_file',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
40 default=None, help='Input fasta file',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
41 metavar='FILE')
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
42 parser.add_option('-f', '--format', dest='seq_format',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
43 default='fasta', help='Sequence format (e.g. fasta, fastq, sff)')
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
44 parser.add_option('--table', dest='table',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
45 default=1, help='NCBI Translation table', type='int')
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
46 parser.add_option('-t', '--ftype', dest='ftype', type='choice',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
47 choices=['CDS', 'ORF'], default='ORF',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
48 help='Find ORF or CDSs')
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
49 parser.add_option('-e', '--ends', dest='ends', type='choice',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
50 choices=['open', 'closed'], default='closed',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
51 help='Open or closed. Closed ensures start/stop codons are present')
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
52 parser.add_option('-m', '--mode', dest='mode', type='choice',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
53 choices=['all', 'top', 'one'], default='all',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
54 help='Output all ORFs/CDSs from sequence, all ORFs/CDSs '
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
55 'with max length, or first with maximum length')
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
56 parser.add_option('--min_len', dest='min_len',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
57 default=10, help='Minimum ORF/CDS length', type='int')
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
58 parser.add_option('-s', '--strand', dest='strand', type='choice',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
59 choices=['forward', 'reverse', 'both'], default='both',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
60 help='Strand to search for features on')
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
61 parser.add_option('--on', dest='out_nuc_file',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
62 default=None, help='Output nucleotide sequences, or - for STDOUT',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
63 metavar='FILE')
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
64 parser.add_option('--op', dest='out_prot_file',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
65 default=None, help='Output protein sequences, or - for STDOUT',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
66 metavar='FILE')
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
67 parser.add_option('--ob', dest='out_bed_file',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
68 default=None, help='Output BED file, or - for STDOUT',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
69 metavar='FILE')
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
70 parser.add_option('--og', dest='out_gff3_file',
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
71 default=None, help='Output GFF3 file, or - for STDOUT',
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
72 metavar='FILE')
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
73 parser.add_option('-v', '--version', dest='version',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
74 default=False, action='store_true',
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
75 help='Show version and quit')
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
76
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
77 options, args = parser.parse_args()
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
78
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
79 if options.version:
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
80 print("v0.2.0")
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
81 sys.exit(0)
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
82
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
83 if not options.input_file:
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
84 sys.exit("Input file is required")
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
85
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
86 if not any((options.out_nuc_file, options.out_prot_file, options.out_bed_file, options.out_gff3_file)):
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
87 sys.exit("At least one output file is required")
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
88
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
89 try:
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
90 table_obj = CodonTable.ambiguous_generic_by_id[options.table]
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
91 except KeyError:
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
92 sys.exit("Unknown codon table %i" % options.table)
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
93
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
94 if options.seq_format.lower() == "sff":
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
95 seq_format = "sff-trim"
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
96 elif options.seq_format.lower() == "fasta":
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
97 seq_format = "fasta"
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
98 elif options.seq_format.lower().startswith("fastq"):
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
99 seq_format = "fastq"
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
100 else:
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
101 sys.exit("Unsupported file type %r" % options.seq_format)
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
102
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
103 print "Genetic code table %i" % options.table
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
104 print "Minimum length %i aa" % options.min_len
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
105 # print "Taking %s ORF(s) from %s strand(s)" % (mode, strand)
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
106
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
107 starts = sorted(table_obj.start_codons)
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
108 assert "NNN" not in starts
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
109 re_starts = re.compile("|".join(starts))
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
110
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
111 stops = sorted(table_obj.stop_codons)
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
112 assert "NNN" not in stops
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
113 re_stops = re.compile("|".join(stops))
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
114
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
115
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
116 def start_chop_and_trans(s, strict=True):
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
117 """Returns offset, trimmed nuc, protein."""
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
118 if strict:
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
119 assert s[-3:] in stops, s
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
120 assert len(s) % 3 == 0
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
121 for match in re_starts.finditer(s):
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
122 # Must check the start is in frame
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
123 start = match.start()
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
124 if start % 3 == 0:
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
125 n = s[start:]
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
126 assert len(n) % 3 == 0, "%s is len %i" % (n, len(n))
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
127 if strict:
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
128 t = translate(n, options.table, cds=True)
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
129 else:
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
130 # Use when missing stop codon,
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
131 t = "M" + translate(n[3:], options.table, to_stop=True)
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
132 return start, n, t
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
133 return None, None, None
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
134
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
135
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
136 def break_up_frame(s):
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
137 """Returns offset, nuc, protein."""
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
138 start = 0
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
139 for match in re_stops.finditer(s):
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
140 index = match.start() + 3
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
141 if index % 3 != 0:
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
142 continue
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
143 n = s[start:index]
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
144 if options.ftype == "CDS":
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
145 offset, n, t = start_chop_and_trans(n)
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
146 else:
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
147 offset = 0
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
148 t = translate(n, options.table, to_stop=True)
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
149 if n and len(t) >= options.min_len:
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
150 yield start + offset, n, t
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
151 start = index
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
152 if options.ends == "open":
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
153 # No stop codon, Biopython's strict CDS translate will fail
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
154 n = s[start:]
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
155 # Ensure we have whole codons
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
156 # TODO - Try appending N instead?
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
157 # TODO - Do the next four lines more elegantly
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
158 if len(n) % 3:
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
159 n = n[:-1]
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
160 if len(n) % 3:
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
161 n = n[:-1]
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
162 if options.ftype == "CDS":
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
163 offset, n, t = start_chop_and_trans(n, strict=False)
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
164 else:
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
165 offset = 0
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
166 t = translate(n, options.table, to_stop=True)
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
167 if n and len(t) >= options.min_len:
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
168 yield start + offset, n, t
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
169
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
170
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
171 def get_all_peptides(nuc_seq):
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
172 """Returns start, end, strand, nucleotides, protein.
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
173
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
174 Co-ordinates are Python style zero-based.
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
175 """
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
176 # TODO - Refactor to use a generator function (in start order)
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
177 # rather than making a list and sorting?
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
178 answer = []
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
179 full_len = len(nuc_seq)
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
180 if options.strand != "reverse":
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
181 for frame in range(0, 3):
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
182 for offset, n, t in break_up_frame(nuc_seq[frame:]):
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
183 start = frame + offset # zero based
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
184 answer.append((start, start + len(n), +1, n, t))
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
185 if options.strand != "forward":
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
186 rc = reverse_complement(nuc_seq)
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
187 for frame in range(0, 3):
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
188 for offset, n, t in break_up_frame(rc[frame:]):
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
189 start = full_len - frame - offset # zero based
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
190 answer.append((start - len(n), start, -1, n, t))
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
191 answer.sort()
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
192 return answer
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
193
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
194
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
195 def get_top_peptides(nuc_seq):
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
196 """Returns all peptides of max length."""
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
197 values = list(get_all_peptides(nuc_seq))
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
198 if not values:
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
199 raise StopIteration
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
200 max_len = max(len(x[-1]) for x in values)
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
201 for x in values:
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
202 if len(x[-1]) == max_len:
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
203 yield x
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
204
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
205
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
206 def get_one_peptide(nuc_seq):
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
207 """Returns first (left most) peptide with max length."""
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
208 values = list(get_top_peptides(nuc_seq))
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
209 if not values:
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
210 raise StopIteration
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
211 yield values[0]
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
212
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
213 if options.mode == "all":
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
214 get_peptides = get_all_peptides
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
215 elif options.mode == "top":
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
216 get_peptides = get_top_peptides
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
217 elif options.mode == "one":
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
218 get_peptides = get_one_peptide
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
219
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
220 in_count = 0
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
221 out_count = 0
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
222 if options.out_nuc_file == "-":
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
223 out_nuc = sys.stdout
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
224 elif options.out_nuc_file:
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
225 out_nuc = open(options.out_nuc_file, "w")
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
226 else:
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
227 out_nuc = None
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
228
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
229 if options.out_prot_file == "-":
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
230 out_prot = sys.stdout
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
231 elif options.out_prot_file:
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
232 out_prot = open(options.out_prot_file, "w")
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
233 else:
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
234 out_prot = None
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
235
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
236 if options.out_bed_file == "-":
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
237 out_bed = sys.stdout
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
238 elif options.out_bed_file:
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
239 out_bed = open(options.out_bed_file, "w")
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
240 else:
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
241 out_bed = None
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
242
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
243 if options.out_gff3_file == "-":
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
244 out_gff3 = sys.stdout
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
245 elif options.out_gff3_file:
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
246 out_gff3 = open(options.out_gff3_file, "w")
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
247 else:
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
248 out_gff3 = None
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
249
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
250 if out_gff3:
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
251 out_gff3.write('##gff-version 3\n')
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
252
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
253 for record in SeqIO.parse(options.input_file, seq_format):
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
254 for i, (f_start, f_end, f_strand, n, t) in enumerate(get_peptides(str(record.seq).upper())):
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
255 out_count += 1
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
256 if f_strand == +1:
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
257 loc = "%i..%i" % (f_start + 1, f_end)
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
258 else:
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
259 loc = "complement(%i..%i)" % (f_start + 1, f_end)
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
260 descr = "length %i aa, %i bp, from %s of %s" \
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
261 % (len(t), len(n), loc, record.description)
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
262 fid = record.id + "|%s%i" % (options.ftype, i + 1)
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
263 r = SeqRecord(Seq(n), id=fid, name="", description=descr)
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
264 t = SeqRecord(Seq(t), id=fid, name="", description=descr)
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
265 if out_nuc:
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
266 SeqIO.write(r, out_nuc, "fasta")
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
267 if out_prot:
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
268 SeqIO.write(t, out_prot, "fasta")
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
269 nice_strand = '+' if f_strand == +1 else '-'
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
270 if out_bed:
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
271 out_bed.write('\t'.join(map(str, [record.id, f_start, f_end, fid, 0, nice_strand])) + '\n')
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
272 if out_gff3:
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
273 out_gff3.write('\t'.join(map(str, [record.id, 'getOrfsOrCds', 'CDS', f_start + 1, f_end, '.',
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
274 nice_strand, 0, 'ID=%s%s' % (options.ftype, i + 1)])) + '\n')
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
275 in_count += 1
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
276 if out_nuc and out_nuc is not sys.stdout:
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
277 out_nuc.close()
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
278 if out_prot and out_prot is not sys.stdout:
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
279 out_prot.close()
8
09a8be9247ca v0.2.0 with GFF3 output
peterjc
parents: 7
diff changeset
280 if out_bed and out_bed is not sys.stdout:
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
281 out_bed.close()
5
5208c15805ec Uploaded v0.0.5 dependant on Biopython 1.62
peterjc
parents:
diff changeset
282
7
705a2e2df7fb v0.1.1 fix typo; v0.1.0 BED output (Eric Rasche), NCBI genetic code 24; v0.0.7 embeds citation
peterjc
parents: 5
diff changeset
283 print "Found %i %ss in %i sequences" % (out_count, options.ftype, in_count)