Mercurial > repos > peterjc > seq_filter_by_id

changeset 0:5844f6a450ed

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Migrated tool version 0.0.1 from old tool shed archive to new tool shed repository
author peterjc
date Tue, 07 Jun 2011 17:24:30 -0400
parents
children 262f08104540
files tools/filters/seq_filter_by_id.py tools/filters/seq_filter_by_id.txt tools/filters/seq_filter_by_id.xml
diffstat 3 files changed, 358 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/filters/seq_filter_by_id.py	Tue Jun 07 17:24:30 2011 -0400
@@ -0,0 +1,182 @@
+#!/usr/bin/env python
+"""Filter a FASTA, FASTQ or SSF file with IDs from a tabular file.
+
+Takes six command line options, tabular filename, ID column numbers (comma
+separated list using one based counting), input filename, input type (e.g.
+FASTA or SFF) and two output filenames (for records with and without the
+given IDs, same format as input sequence file).
+
+If either output filename is just a minus sign, that file is not created.
+This is intended to allow output for just the matched (or just the non-matched)
+records.
+
+When filtering an SFF file, any Roche XML manifest in the input file is
+preserved in both output files.
+
+Note in the default NCBI BLAST+ tabular output, the query sequence ID is
+in column one, and the ID of the match from the database is in column two.
+Here sensible values for the column numbers would therefore be "1" or "2".
+
+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 2010 by Peter Cock, SCRI, 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, cols_arg, in_file, seq_format, out_positive_file, out_negative_file = sys.argv[1:]
+except ValueError:
+    stop_err("Expected six arguments, got %i:\n%s" % (len(sys.argv)-1, " ".join(sys.argv)))
+try:
+    columns = [int(arg)-1 for arg in cols_arg.split(",")]
+except ValueError:
+    stop_err("Expected list of columns (comma separated integers), got %s" % cols_arg)
+
+if out_positive_file == "-" and out_negative_file == "-":
+    stop_err("Neither output file requested")
+
+
+#Read tabular file and record all specified identifiers
+ids = set()
+handle = open(tabular_file, "rU")
+if len(columns)>1:
+    #General case of many columns
+    for line in handle:
+        if line.startswith("#"):
+            #Ignore comments
+            continue
+        parts = line.rstrip("\n").split("\t")
+        for col in columns:
+            ids.add(parts[col])
+    print "Using %i IDs from %i columns of tabular file" % (len(ids), len(columns))
+else:
+    #Single column, special case speed up
+    col = columns[0]
+    for line in handle:
+        if not line.startswith("#"):
+            ids.add(line.rstrip("\n").split("\t")[col])
+    print "Using %i IDs from tabular file" % (len(ids))
+handle.close()
+
+
+if seq_format.lower()=="sff":
+    #Now write filtered SFF file based on IDs from BLAST file
+    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
+    #This makes two passes though the SFF file with isn't so efficient,
+    #but this makes the code simple.
+    if out_positive_file != "-":
+        out_handle = open(out_positive_file, "wb")
+        writer = SffWriter(out_handle, xml=manifest)
+        in_handle.seek(0) #start again after getting manifest
+        pos_count = writer.write_file(rec for rec in SffIterator(in_handle) if rec.id in ids)
+        out_handle.close()
+    if out_negative_file != "-":
+        out_handle = open(out_negative_file, "wb")
+        writer = SffWriter(out_handle, xml=manifest)
+        in_handle.seek(0) #start again
+        neg_count = writer.write_file(rec for rec in SffIterator(in_handle) if rec.id not in ids)
+        out_handle.close()
+    #And we're done
+    in_handle.close()
+    #At the time of writing, Galaxy doesn't show SFF file read counts,
+    #so it is useful to put them in stdout and thus shown in job info.
+    if out_positive_file != "-" and out_negative_file != "-":
+        print "%i with and %i without specified IDs" % (pos_count, neg_count)
+    elif out_positive_file != "-":
+        print "%i with specified IDs" % pos_count
+    elif out_negative_file != "-":
+        print "%i without specified IDs" % neg_count
+elif seq_format.lower()=="fasta":
+    #Write filtered FASTA file based on IDs from tabular file
+    from galaxy_utils.sequence.fasta import fastaReader, fastaWriter
+    reader = fastaReader(open(in_file, "rU"))
+    if out_positive_file != "-" and out_negative_file != "-":
+        print "Generating two FASTA files"
+        positive_writer = fastaWriter(open(out_positive_file, "w"))
+        negative_writer = fastaWriter(open(out_negative_file, "w"))
+        for record in reader:
+            #The [1:] is because the fastaReader leaves the > on the identifer.
+            if record.identifier and record.identifier.split()[0][1:] in ids:
+                positive_writer.write(record)
+            else:
+                negative_writer.write(record)
+        positive_writer.close()
+        negative_writer.close()
+    elif out_positive_file != "-":
+        print "Generating matching FASTA file"
+        positive_writer = fastaWriter(open(out_positive_file, "w"))
+        for record in reader:
+            #The [1:] is because the fastaReader leaves the > on the identifer.
+            if record.identifier and record.identifier.split()[0][1:] in ids:
+                positive_writer.write(record)
+        positive_writer.close()
+    elif out_negative_file != "-":
+        print "Generating non-matching FASTA file"
+        negative_writer = fastaWriter(open(out_negative_file, "w"))
+        for record in reader:
+            #The [1:] is because the fastaReader leaves the > on the identifer.
+            if not record.identifier or record.identifier.split()[0][1:] not in ids:
+                negative_writer.write(record)
+        negative_writer.close()
+elif seq_format.lower().startswith("fastq"):
+    #Write filtered FASTQ file based on IDs from tabular file
+    from galaxy_utils.sequence.fastq import fastqReader, fastqWriter
+    reader = fastqReader(open(in_file, "rU"))
+    if out_positive_file != "-" and out_negative_file != "-":
+        print "Generating two FASTQ files"
+        positive_writer = fastqWriter(open(out_positive_file, "w"))
+        negative_writer = fastqWriter(open(out_negative_file, "w"))
+        for record in reader:
+            #The [1:] is because the fastaReader leaves the @ on the identifer.
+            if record.identifier and record.identifier.split()[0][1:] in ids:
+                positive_writer.write(record)
+            else:
+                negative_writer.write(record)
+        positive_writer.close()
+        negative_writer.close()
+    elif out_positive_file != "-":
+        print "Generating matching FASTQ file"
+        positive_writer = fastqWriter(open(out_positive_file, "w"))
+        for record in reader:
+            #The [1:] is because the fastaReader leaves the @ on the identifer.
+            if record.identifier and record.identifier.split()[0][1:] in ids:
+                positive_writer.write(record)
+        positive_writer.close()
+    elif out_negative_file != "-":
+        print "Generating non-matching FASTQ file"
+        negative_writer = fastqWriter(open(out_negative_file, "w"))
+        for record in reader:
+            #The [1:] is because the fastaReader leaves the @ on the identifer.
+            if not record.identifier or record.identifier.split()[0][1:] not in ids:
+                negative_writer.write(record)
+        negative_writer.close()
+else:
+    stop_err("Unsupported file type %r" % seq_format)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/filters/seq_filter_by_id.txt	Tue Jun 07 17:24:30 2011 -0400
@@ -0,0 +1,76 @@
+Galaxy tool to filter FASTA, FASTQ or SFF sequences by ID
+=========================================================
+
+This tool is copyright 2010 by Peter Cock, SCRI, UK. All rights reserved.
+See the licence text below.
+
+This tool is a short Python script (using both the Galaxy and Biopython library
+functions) which divides a FASTA, FASTQ, or SFF file in two, those sequences with
+or without an ID present in the specified column(s) of a tabular file. Example uses
+include filtering based on search results from a tool like NCBI BLAST before
+assembly.
+
+There are just two files to install:
+
+* seq_filter_by_id.py (the Python script)
+* seq_filter_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_filter_by_id.xml" />
+
+You will also need to install Biopython 1.54 or later. That's it.
+
+
+History
+=======
+
+v0.0.1 - Initial version, combining three separate scripts for each file format.
+
+
+Developers
+==========
+
+This script and related tools are being developed on the following hg branch:
+http://bitbucket.org/peterjc/galaxy-central/src/tools
+
+This incorporates the previously used hg branch:
+http://bitbucket.org/peterjc/galaxy-central/src/fasta_filter
+
+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_filter_by_id.tar.gz tools/filters/seq_filter_by_id.*
+
+Check this worked:
+
+$ tar -tzf seq_filter_by_id.tar.gz
+filter/seq_filter_by_id.py
+filter/seq_filter_by_id.txt
+filter/seq_filter_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_filter_by_id.xml	Tue Jun 07 17:24:30 2011 -0400
@@ -0,0 +1,100 @@
+<tool id="seq_filter_by_id" name="Filter sequences by ID" version="0.0.1">
+	<description>from a tabular file</description>
+	<command interpreter="python">
+seq_filter_by_id.py $input_tabular $columns $input_file $input_file.ext
+#if $output_choice_cond.output_choice=="both"
+ $output_pos $output_neg
+#elif $output_choice_cond.output_choice=="pos"
+ $output_pos -
+#elif $output_choice_cond.output_choice=="neg"
+ - $output_neg
+#end if
+	</command>
+	<inputs>
+		<param name="input_file" type="data" format="fasta,fastq,sff" label="Sequence file to filter on the identifiers" description="FASTA, FASTQ, or SFF format." />
+		<param name="input_tabular" type="data" format="tabular" label="Tabular file containing sequence identifiers"/>
+		<param name="columns" type="data_column" data_ref="input_tabular" multiple="True" numerical="False" label="Column(s) containing sequence identifiers" help="Multi-select list - hold the appropriate key while clicking to select multiple columns">
+			<validator type="no_options" message="Pick at least one column"/>
+		</param>
+		<conditional name="output_choice_cond">
+			<param name="output_choice" type="select" label="Output positive matches, negative matches, or both?">
+				<option value="both">Both positive matches (ID on list) and negative matches (ID not on list), as two files</option>
+				<option value="pos">Just positive matches (ID on list), as a single file</option>
+				<option value="neg">Just negative matches (ID not on list), as a single file</option>
+			</param>
+			<!-- Seems need these dummy entries here, compare this to indels/indel_sam2interval.xml -->
+			<when value="both" />
+			<when value="pos" />
+			<when value="neg" />
+		</conditional>
+	</inputs>
+	<outputs>
+		<data name="output_pos" format="fasta" label="With matched ID">
+            <!-- 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>
+			<filter>output_choice_cond["output_choice"] != "neg"</filter>
+		</data>
+		<data name="output_neg" format="fasta" label="Without matched ID">
+            <!-- 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>
+			<filter>output_choice_cond["output_choice"] != "pos"</filter>
+		</data>
+	</outputs>
+	<tests>
+	</tests>
+	<requirements>
+		<requirement type="python-module">Bio</requirement>
+	</requirements>
+	<help>
+
+**What it does**
+
+By default it divides a FASTA, FASTQ or Standard Flowgram Format (SFF) file in
+two, those sequences with or without an ID present in the tabular file column(s)
+specified. You can opt to have a single output file of just the matching records,
+or just the non-matching ones.
+
+Note that the order of sequences in the original sequence file is preserved, as
+is any Roche XML Manifest in an SFF file. Also, if any sequences share an
+identifier (which would be very unusual in SFF files), duplicates are not removed.
+
+**Example Usage**
+
+You may have performed some kind of contamination search, for example running
+BLASTN against a database of cloning vectors or bacteria, giving you a tabular
+file containing read identifiers. You could use this tool to extract only the
+reads without BLAST matches (i.e. those which do not match your contaminant
+database).
+
+You may have a file of FASTA sequences which has been run some some analysis
+tool giving tabular output, which has then been filtered on some criteria.
+You can then use this tool to divide the original FASTA file into those entries
+matching or not matching your criteria (those with or without their identifier
+in the filtered tabular file).
+
+**Citation**
+
+This tool uses Biopython to read and write SFF 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>