Mercurial > repos > peterjc > seq_select_by_id

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/filters/seq_select_by_id.py	Tue Jun 07 17:43:38 2011 -0400
@@ -0,0 +1,126 @@
+#!/usr/bin/env python
+"""Select FASTA, QUAL, FASTQ or SSF sequences by IDs from a tabular file.
+
+Takes five command line options, tabular filename, ID column number (using
+one based counting), input filename, input type (e.g. FASTA or SFF) and the
+output filename (same format as input sequence file).
+
+When selecting from an SFF file, any Roche XML manifest in the input file is
+preserved in both output files.
+
+This tool is a short Python script which requires Biopython 1.54 or later
+for SFF file support. 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 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
+
+def stop_err(msg, err=1):
+    sys.stderr.write(msg.rstrip() + "\n")
+    sys.exit(err)
+
+#Parse Command Line
+try:
+    tabular_file, col_arg, in_file, seq_format, out_file = sys.argv[1:]
+except ValueError:
+    stop_err("Expected five arguments, got %i:\n%s" % (len(sys.argv)-1, " ".join(sys.argv)))
+try:
+    if col_arg.startswith("c"):
+        column = int(col_arg[1:])-1
+    else:
+        column = int(col_arg)-1
+except ValueError:
+    stop_err("Expected column number, got %s" % cols_arg)
+
+if seq_format == "fastqcssanger":
+    stop_err("Colorspace FASTQ not supported.")
+elif seq_format.lower() in ["sff", "fastq", "qual", "fasta"]:
+    seq_format = seq_format.lower()
+elif seq_format.lower().startswith("fastq"):
+    #We don't care how the qualities are encoded
+    seq_format = "fastq"
+elif seq_format.lower().startswith("qual"):
+    #We don't care what the scores are
+    seq_format = "qual"
+else:
+    stop_err("Unrecognised file format %r" % seq_format)
+
+
+try:
+    from Bio import SeqIO
+except ImportError:
+    stop_err("Biopython 1.54 or later is required")
+
+
+def parse_ids(tabular_file, col):
+    """Read tabular file and record all specified identifiers."""
+    handle = open(tabular_file, "rU")
+    for line in handle:
+        if not line.startswith("#"):
+            yield line.rstrip("\n").split("\t")[col].strip()
+    handle.close()
+
+#Index the sequence file.
+#If very big, could use SeqIO.index_db() to avoid memory bottleneck...
+records = SeqIO.index(in_file, seq_format)
+print "Indexed %i sequences" % len(records)
+
+if seq_format.lower()=="sff":
+    #Special case to try to preserve the XML manifest
+    try:
+        from Bio.SeqIO.SffIO import SffIterator, SffWriter
+    except ImportError:
+        stop_err("Requires Biopython 1.54 or later")
+
+    try:
+        from Bio.SeqIO.SffIO import ReadRocheXmlManifest
+    except ImportError:
+        #Prior to Biopython 1.56 this was a private function
+        from Bio.SeqIO.SffIO import _sff_read_roche_index_xml as ReadRocheXmlManifest
+
+    in_handle = open(in_file, "rb") #must be binary mode!
+    try:
+        manifest = ReadRocheXmlManifest(in_handle)
+    except ValueError:
+        manifest = None
+    in_handle.close()
+
+    out_handle = open(out_file, "wb")
+    writer = SffWriter(out_handle, xml=manifest)
+    count = 0
+    #This does have the overhead of parsing into SeqRecord objects,
+    #but doing the header and index at the low level is too fidly.
+    iterator = (records[name] for name in parse_ids(tabular_file, column))
+    try:
+        count = writer.write_file(iterator)
+    except KeyError, err:
+        out_handle.close()
+        if name not in records:
+            stop_err("Identifier %s not found in sequence file" % name)
+        else:
+            raise err
+    out_handle.close()
+else:
+    #Avoid overhead of parsing into SeqRecord objects,
+    #just re-use the original formatting from the input file.
+    out_handle = open(out_file, "w")
+    count = 0
+    for name in parse_ids(tabular_file, column):
+        try:
+            out_handle.write(records.get_raw(name))
+        except KeyError:
+            out_handle.close()
+            stop_err("Identifier %s not found in sequence file" % name)
+        count += 1
+    out_handle.close()
+
+print "Selected %i sequences by ID" % count
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/filters/seq_select_by_id.txt	Tue Jun 07 17:43:38 2011 -0400
@@ -0,0 +1,78 @@
+Galaxy tool to select FASTA, FASTQ or SFF sequences by ID
+=========================================================
+
+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 extract
+sequences from a FASTA, QUAL, FASTQ, or SFF file based on the list of IDs given
+by a column of a tabular file. The output order follows that of the tabular file,
+and if there are duplicates in the tabular file, there will be duplicates in the
+output sequence file.
+
+See also the sister tool to filter sequence files according to IDs from column(s)
+of a tabular file, where the output order follows the sequence file, and any
+duplicate IDs are ignored.
+
+There are just two files to install:
+
+* seq_select_by_id.py (the Python script)
+* seq_select_by_id.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/sff_select_by_id.xml" />
+
+You will also need to install Biopython 1.54 or later. 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 seq_select_by_id.tar.gz tools/filters/seq_select_by_id.*
+
+Check this worked:
+
+$ tar -tzf seq_select_by_id.tar.gz
+filter/seq_select_by_id.py
+filter/seq_select_by_id.txt
+filter/seq_select_by_id.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/seq_select_by_id.xml	Tue Jun 07 17:43:38 2011 -0400
@@ -0,0 +1,51 @@
+<tool id="seq_select_by_id" name="Select sequences by ID" version="0.0.1">
+	<description>from a tabular file</description>
+	<command interpreter="python">
+seq_select_by_id.py $input_tabular $column $input_file $input_file.ext $output_file
+	</command>
+	<inputs>
+		<param name="input_file" type="data" format="fasta,qual,fastq,sff" label="Sequence file to select from" help="FASTA, QUAL, FASTQ, or SFF format." />
+		<param name="input_tabular" type="data" format="tabular" label="Tabular file containing sequence identifiers"/>
+		<param name="column" type="data_column" data_ref="input_tabular" multiple="False" numerical="False" label="Column containing sequence identifiers"/>
+	</inputs>
+	<outputs>
+		<data name="output_file" format="fasta" label="Selected sequences">
+			<!-- TODO - Replace this with format="input:input_fastq" if/when that works -->
+			<change_format>
+				<when input_dataset="input_file" attribute="extension" value="sff" format="sff" />
+				<when input_dataset="input_file" attribute="extension" value="fastq" format="fastq" />
+				<when input_dataset="input_file" attribute="extension" value="fastqsanger" format="fastqsanger" />
+				<when input_dataset="input_file" attribute="extension" value="fastqsolexa" format="fastqsolexa" />
+				<when input_dataset="input_file" attribute="extension" value="fastqillumina" format="fastqillumina" />
+				<when input_dataset="input_file" attribute="extension" value="fastqcssanger" format="fastqcssanger" />
+			</change_format>
+		</data>
+	</outputs>
+	<tests>
+	</tests>
+	<requirements>
+		<requirement type="python-module">Bio</requirement>
+	</requirements>
+	<help>
+
+**What it does**
+
+Takes a FASTA, QUAL, FASTQ or Standard Flowgram Format (SFF) file and produces a
+new sequence file (of the same format) containing only the records with identifiers
+in the tabular file (in the order from the tabular file).
+
+WARNING: If you have any duplicates in the tabular file identifiers, you will get
+duplicate sequences in the output.
+
+**Citation**
+
+This tool uses Biopython to read, write and index sequence files. 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>