changeset 0:9cff9a1176ea

Uploaded v0.0.1
author peterjc
date Thu, 19 Jan 2012 10:17:10 -0500
parents
children 922d69bd5258
files tools/filters/get_orfs_or_cdss.py tools/filters/get_orfs_or_cdss.txt tools/filters/get_orfs_or_cdss.xml
diffstat 3 files changed, 441 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/filters/get_orfs_or_cdss.py	Thu Jan 19 10:17:10 2012 -0500
@@ -0,0 +1,219 @@
+#!/usr/bin/env python
+"""Find ORFs in a nucleotide sequence file.
+
+get_orfs_or_cdss.py $input_fasta $input_format $table $ftype $ends $mode $min_len $strand $out_nuc_file $out_prot_file
+
+Takes ten command line options, input sequence filename, format, genetic
+code, CDS vs ORF, end type (open, closed), selection mode (all, top, one),
+minimum length (in amino acids), strand (both, forward, reverse), output
+nucleotide filename, and output protein filename.
+
+This tool is a short Python script which requires Biopython. If you use
+this tool in scientific work leading to a publication, please cite the
+Biopython application note:
+
+Cock et al 2009. Biopython: freely available Python tools for computational
+molecular biology and bioinformatics. Bioinformatics 25(11) 1422-3.
+http://dx.doi.org/10.1093/bioinformatics/btp163 pmid:19304878.
+
+This script is copyright 2011 by Peter Cock, The James Hutton Institute
+(formerly SCRI), Dundee, UK. All rights reserved.
+
+See accompanying text file for licence details (MIT/BSD style).
+
+This is version 0.0.1 of the script.
+"""
+import sys
+import re
+
+def stop_err(msg, err=1):
+    sys.stderr.write(msg.rstrip() + "\n")
+    sys.exit(err)
+
+try:
+    from Bio.Seq import Seq, reverse_complement, translate
+    from Bio.SeqRecord import SeqRecord
+    from Bio import SeqIO
+    from Bio.Data import CodonTable
+except ImportError:
+    stop_err("Missing Biopython library")
+
+#Parse Command Line
+try:
+    input_file, seq_format, table, ftype, ends, mode, min_len, strand, out_nuc_file, out_prot_file = sys.argv[1:]
+except ValueError:
+    stop_err("Expected ten arguments, got %i:\n%s" % (len(sys.argv)-1, " ".join(sys.argv)))
+
+try:
+    table = int(table)
+except ValueError:
+    stop_err("Expected integer for genetic code table, got %s" % table)
+
+try:
+    table_obj = CodonTable.ambiguous_generic_by_id[table]
+except KeyError:
+    stop_err("Unknown codon table %i" % table)
+
+if ftype not in ["CDS", "ORF"]:
+    stop_err("Expected CDS or ORF, got %s" % ftype)
+
+if ends not in ["open", "closed"]:
+    stop_err("Expected open or closed for end treatment, got %s" % ends)
+
+try:
+    min_len = int(min_len)
+except ValueError:
+    stop_err("Expected integer for min_len, got %s" % min_len)
+
+if seq_format.lower()=="sff":
+    seq_format = "sff-trim"
+elif seq_format.lower()=="fasta":
+    seq_format = "fasta"
+elif seq_format.lower().startswith("fastq"):
+    seq_format = "fastq"
+else:
+    stop_err("Unsupported file type %r" % seq_format)
+
+print "Genetic code table %i" % table
+print "Minimum length %i aa" % min_len
+#print "Taking %s ORF(s) from %s strand(s)" % (mode, strand)
+
+starts = sorted(table_obj.start_codons)
+assert "NNN" not in starts
+re_starts = re.compile("|".join(starts))
+
+stops = sorted(table_obj.stop_codons)
+assert "NNN" not in stops
+re_stops = re.compile("|".join(stops))
+
+def start_chop_and_trans(s, strict=True):
+    """Returns offset, trimmed nuc, protein."""
+    if strict:
+        assert s[-3:] in stops, s
+    assert len(s) % 3 == 0
+    for match in re_starts.finditer(s):
+        #Must check the start is in frame
+        start = match.start()
+        if start % 3 == 0:
+            n = s[start:]
+            assert len(n) % 3 == 0, "%s is len %i" % (n, len(n))
+            if strict:
+                t = translate(n, table, cds=True)
+            else:
+                #Use when missing stop codon,
+                t = "M" + translate(n[3:], table, to_stop=True)
+            return start, n, t
+    return None, None, None
+
+def break_up_frame(s):
+    """Returns offset, nuc, protein."""
+    start = 0
+    for match in re_stops.finditer(s):
+        index = match.start() + 3
+        if index % 3 != 0:
+            continue
+        n = s[start:index]
+        if ftype=="CDS":
+            offset, n, t = start_chop_and_trans(n)
+        else:
+            offset = 0
+            t = translate(n, table, to_stop=True)
+        if n and len(t) >= min_len:
+            yield start + offset, n, t
+        start = index
+    if ends == "open":
+        #No stop codon, Biopython's strict CDS translate will fail
+        n = s[start:]
+        #Ensure we have whole codons
+        #TODO - Try appending N instead?
+        #TODO - Do the next four lines more elegantly
+        if len(n) % 3:
+            n = n[:-1]
+        if len(n) % 3:
+            n = n[:-1]
+        if ftype=="CDS":
+            offset, n, t = start_chop_and_trans(n, strict=False)
+        else:
+            offset = 0
+            t = translate(n, table, to_stop=True)
+        if n and len(t) >= min_len:
+            yield start + offset, n, t
+                        
+
+def get_all_peptides(nuc_seq):
+    """Returns start, end, strand, nucleotides, protein.
+
+    Co-ordinates are Python style zero-based.
+    """
+    #TODO - Refactor to use a generator function (in start order)
+    #rather than making a list and sorting?
+    answer = []
+    full_len = len(nuc_seq)
+    if strand != "reverse":
+        for frame in range(0,3):
+            for offset, n, t in break_up_frame(nuc_seq[frame:]):
+                start = frame + offset #zero based
+                answer.append((start, start + len(n), +1, n, t))
+    if strand != "forward":
+        rc = reverse_complement(nuc_seq)
+        for frame in range(0,3) :
+            for offset, n, t in break_up_frame(rc[frame:]):
+                start = full_len - frame - offset #zero based
+                answer.append((start, start + len(n), -1, n ,t))
+    answer.sort()
+    return answer
+
+def get_top_peptides(nuc_seq):
+    """Returns all peptides of max length."""
+    values = list(get_all_peptides(nuc_seq))
+    if not values:
+        raise StopIteration
+    max_len = max(len(x[-1]) for x in values)
+    for x in values:
+        if len(x[-1]) == max_len:
+            yield x
+
+def get_one_peptide(nuc_seq):
+    """Returns first (left most) peptide with max length."""
+    values = list(get_top_peptides(nuc_seq))
+    if not values:
+        raise StopIteration
+    yield values[0]
+
+if mode == "all":
+    get_peptides = get_all_peptides
+elif mode == "top":
+    get_peptides = get_top_peptides
+elif mode == "one":
+    get_peptides = get_one_peptide
+
+in_count = 0
+out_count = 0
+if out_nuc_file == "-":
+    out_nuc = sys.stdout
+else:
+    out_nuc = open(out_nuc_file, "w")
+if out_prot_file == "-":
+    out_prot = sys.stdout
+else:
+    out_prot = open(out_prot_file, "w")
+for record in SeqIO.parse(input_file, seq_format):
+    for i, (f_start, f_end, f_strand, n, t) in enumerate(get_peptides(str(record.seq).upper())):
+        out_count += 1
+        if f_strand == +1:
+            loc = "%i..%i" % (f_start+1, f_end)
+        else:
+            loc = "complement(%i..%i)" % (f_start+1, f_end)
+        descr = "length %i aa, %i bp, from %s of %s" \
+                % (len(t), len(n), loc, record.description)
+        r = SeqRecord(Seq(n), id = record.id + "|%s%i" % (ftype, i+1), name = "", description= descr)
+        t = SeqRecord(Seq(t), id = record.id + "|%s%i" % (ftype, i+1), name = "", description= descr)
+        SeqIO.write(r, out_nuc, "fasta")
+        SeqIO.write(t, out_prot, "fasta")
+    in_count += 1
+if out_nuc is not sys.stdout:
+    out_nuc.close()
+if out_prot is not sys.stdout:
+    out_prot.close()
+
+print "Found %i %ss in %i sequences" % (out_count, ftype, in_count)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/filters/get_orfs_or_cdss.txt	Thu Jan 19 10:17:10 2012 -0500
@@ -0,0 +1,75 @@
+Galaxy tool to find ORFs or simple CDSs
+=======================================
+
+This tool is copyright 2011 by Peter Cock, The James Hutton Institute
+(formerly SCRI, Scottish Crop Research Institute), UK. All rights reserved.
+See the licence text below.
+
+This tool is a short Python script (using Biopython library functions)
+to search nucleotide sequences for open reading frames (ORFs) or coding
+sequences (CDSs) where the first potential start codon is used. See the
+help text in the XML file for more information.
+
+There are just two files to install:
+
+* get_orfs_or_cdss.py (the Python script)
+* get_orfs_or_cdss.xml (the Galaxy tool definition)
+
+The suggested location is in the Galaxy folder tools/filters next to the tool
+for calling sff_extract.py for converting SFF to FASTQ or FASTA + QUAL.
+
+You will also need to modify the tools_conf.xml file to tell Galaxy to offer the
+tool. One suggested location is in the filters section. Simply add the line:
+
+<tool file="filters/get_orfs_or_cdss.xml" />
+
+You will also need to install Biopython 1.54 or later. If you want to run
+the unit tests, include this line in tools_conf.xml.sample and the sample
+FASTA files under the test-data directory. That's it.
+
+
+History
+=======
+
+v0.0.1 - Initial version.
+
+
+Developers
+==========
+
+This script and related tools are being developed on the following hg branch:
+http://bitbucket.org/peterjc/galaxy-central/src/tools
+
+For making the "Galaxy Tool Shed" http://community.g2.bx.psu.edu/ tarball use
+the following command from the Galaxy root folder:
+
+tar -czf get_orfs_or_cdss.tar.gz tools/filters/get_orfs_or_cdss.*
+
+Check this worked:
+
+$ tar -tzf get_orfs_or_cdss.tar.gz
+filter/get_orfs_or_cdss.py
+filter/get_orfs_or_cdss.txt
+filter/get_orfs_or_cdss.xml
+
+
+Licence (MIT/BSD style)
+=======================
+
+Permission to use, copy, modify, and distribute this software and its
+documentation with or without modifications and for any purpose and
+without fee is hereby granted, provided that any copyright notices
+appear in all copies and that both those copyright notices and this
+permission notice appear in supporting documentation, and that the
+names of the contributors or copyright holders not be used in
+advertising or publicity pertaining to distribution of the software
+without specific prior permission.
+
+THE CONTRIBUTORS AND COPYRIGHT HOLDERS OF THIS SOFTWARE DISCLAIM ALL
+WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING ALL IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL THE
+CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY SPECIAL, INDIRECT
+OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
+OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE
+OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE
+OR PERFORMANCE OF THIS SOFTWARE.
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/filters/get_orfs_or_cdss.xml	Thu Jan 19 10:17:10 2012 -0500
@@ -0,0 +1,147 @@
+<tool id="get_orfs_or_cdss" name="Get open reading frames (ORFs) or coding sequences (CDSs)" version="0.0.1">
+	<description>e.g. to get peptides from ESTs</description>
+	<command interpreter="python">
+get_orfs_or_cdss.py $input_file $input_file.ext $table $ftype $ends $mode $min_len $strand $out_nuc_file $out_prot_file
+	</command>
+	<inputs>
+		<param name="input_file" type="data" format="fasta,fastq,sff" label="Sequence file (nucleotides)" help="FASTA, FASTQ, or SFF format." />
+		<param name="table" type="select" label="Genetic code" help="Tables from the NCBI, these determine the start and stop codons">
+			<option value="1">1. Standard</option>
+			<option value="2">2. Vertebrate Mitochondrial</option>
+			<option value="3">3. Yeast Mitochondrial</option>
+			<option value="4">4. Mold, Protozoan, Coelenterate Mitochondrial and Mycoplasma/Spiroplasma</option>
+			<option value="5">5. Invertebrate Mitochondrial</option>
+			<option value="6">6. Ciliate Macronuclear and Dasycladacean</option>
+			<option value="9">9. Echinoderm Mitochondrial</option>
+			<option value="10">10. Euplotid Nuclear</option>
+			<option value="11">11. Bacterial</option>
+			<option value="12">12. Alternative Yeast Nuclear</option>
+			<option value="13">13. Ascidian Mitochondrial</option>
+			<option value="14">14. Flatworm Mitochondrial</option>
+			<option value="15">15. Blepharisma Macronuclear</option>
+			<option value="16">16. Chlorophycean Mitochondrial</option>
+			<option value="21">21. Trematode Mitochondrial</option>
+			<option value="22">22. Scenedesmus obliquus</option>
+			<option value="23">23. Thraustochytrium Mitochondrial</option>
+		</param>
+		<param name="ftype" type="select" value="True" label="Look for ORFs or CDSs">
+                        <option value="ORF">Look for ORFs (check for stop codons only, ignore start codons)</option>
+                        <option value="CDS">Look for CDSs (with start and stop codons)</option>
+		</param>
+                <param name="ends" type="select" value="open" label="Sequence end treatment">
+			<option value="open">Open ended (will allow missing start/stop codons at the ends)</option>
+                        <option value="closed">Complete (will check for start/stop codons at the ends)</option>
+                        <!-- TODO? Circular, for using this on finished bacteria etc -->
+                </param>
+
+		<param name="mode" type="select" label="Selection criteria" help="Suppose a sequence has ORFs/CDSs of lengths 100, 102 and 102 -- which should be taken? These options would return 3, 2 or 1 ORF.">
+                    <option value="all">All ORFs/CDSs from each sequence</option>
+                    <option value="top">All ORFs/CDSs from each sequence with the maximum length</option>
+                    <option value="one">First ORF/CDS from each sequence with the maximum length</option>
+		</param>
+                <param name="min_len" type="integer" size="5" value="30" label="Minimum length ORF/CDS (in amino acids, e.g. 30 aa = 90 bp plus any stop codon)">
+                </param>
+                <param name="strand" type="select" label="Strand to search" help="Use the forward only option if your sequence directionality is known (e.g. from poly-A tails, or strand specific RNA sequencing.">
+                    <option value="both">Search both the forward and reverse strand</option>
+                    <option value="forward">Only search the forward strand</option>
+                    <option value="reverse">Only search the reverse strand</option>
+                </param>
+	</inputs>
+	<outputs>
+		<data name="out_nuc_file" format="fasta" label="${ftype.value}s (nucleotides)" />
+		<data name="out_prot_file" format="fasta" label="${ftype.value}s (amino acids)" />
+	</outputs>
+	<tests>
+                <test>
+                        <param name="input_file" value="get_orf_input.fasta" />
+                        <param name="table" value="1" />
+                        <param name="ftype" value="CDS" />
+                        <param name="ends" value="open" />
+                        <param name="mode" value="all" />
+                        <param name="min_len" value="10" />
+                        <param name="strand" value="forward" />
+                        <output name="out_nuc_file" file="get_orf_input.t1_nuc_out.fasta" />
+                        <output name="out_prot_file" file="get_orf_input.t1_prot_out.fasta" />
+                </test>
+		<test>
+			<param name="input_file" value="get_orf_input.fasta" />
+			<param name="table" value="11" />
+			<param name="ftype" value="CDS" />
+			<param name="ends" value="closed" />
+			<param name="mode" value="all" />
+			<param name="min_len" value="10" />
+			<param name="strand" value="forward" />
+			<output name="out_nuc_file" file="get_orf_input.t11_nuc_out.fasta" />
+			<output	name="out_prot_file" file="get_orf_input.t11_prot_out.fasta" />
+		</test>
+		<test>
+                        <param name="input_file" value="get_orf_input.fasta" />
+                        <param name="table" value="11" />
+                        <param name="ftype" value="CDS" />
+                        <param name="ends" value="open" />
+                        <param name="mode" value="all" />
+                        <param name="min_len" value="10" />
+                        <param name="strand" value="forward" />
+                        <output name="out_nuc_file" file="get_orf_input.t11_open_nuc_out.fasta" />
+                        <output name="out_prot_file" file="get_orf_input.t11_open_prot_out.fasta" />
+		</test>
+	</tests>
+	<requirements>
+		<requirement type="python-module">Bio</requirement>
+	</requirements>
+	<help>
+
+**What it does**
+
+Takes an input file of nucleotide sequences (typically FASTA, but also FASTQ
+and Standard Flowgram Format (SFF) are supported), and searches each sequence
+for open reading frames (ORFs) or potential coding sequences (CDSs) of the
+given minimum length. These are returned as FASTA files of nucleotides and
+protein sequences.
+
+You can choose to have all the ORFs/CDSs above the minimum length for each
+sequence (similar to the EMBOSS getorf tool), those with the longest length
+equal, or the first ORF/CDS with the longest length (in the special case
+where a sequence encodes two or more long ORFs/CDSs of the same length). The
+last option is a reasonable choice when the input sequences represent EST or
+mRNA sequences, where only one ORF/CDS is expected.
+
+Note that if no ORFs/CDSs in a sequence match the criteria, there will be no
+output for that sequence.
+
+Also note that the ORFs/CDSs are assigned modified identifiers to distinguish
+them from the original full length sequences, by appending a suffix.
+
+The start and stop codons are taken from the `NCBI Genetic Codes
+&lt;http://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi&gt;`_.
+When searching for ORFs, the sequences will run from stop codon to stop
+codon, and any start codons are ignored. When searching for CDSs, the first
+potential start codon will be used, giving the longest possible CDS within
+each ORF, and thus the longest possible protein sequence. This is useful
+for things like BLAST or domain searching, but since this may not be the
+correct start codon may not be appropriate for signal peptide detection
+etc.
+
+**Example Usage**
+
+Given some EST sequences (Sanger capillary reads) assembled into unigenes,
+or a transcriptome assembly from some RNA-Seq, each of your nucleotide
+sequences should (barring sequencing, assembly errors, frame-shifts etc)
+encode one protein as a single ORF/CDS, which you wish to extract (and
+perhaps translate into amino acids).
+
+If your RNS-Seq data was strand specific, and assembled taking this into
+account, you should only search for ORFs/CDSs on the forward strand.
+
+**Citation**
+
+This tool uses Biopython. If you use this tool in scientific work leading
+to a publication, please cite the Biopython application note (and Galaxy
+too of course):
+
+Cock et al 2009. Biopython: freely available Python tools for computational
+molecular biology and bioinformatics. Bioinformatics 25(11) 1422-3.
+http://dx.doi.org/10.1093/bioinformatics/btp163 pmid:19304878.
+
+	</help>
+</tool>