# HG changeset patch # User peterjc # Date 1307481870 14400 # Node ID 5844f6a450ed30279ed1a3cba5a37beec81463e2 Migrated tool version 0.0.1 from old tool shed archive to new tool shed repository diff -r 000000000000 -r 5844f6a450ed tools/filters/seq_filter_by_id.py --- /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) diff -r 000000000000 -r 5844f6a450ed tools/filters/seq_filter_by_id.txt --- /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: + + + +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. diff -r 000000000000 -r 5844f6a450ed tools/filters/seq_filter_by_id.xml --- /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 @@ + + from a tabular file + +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 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + output_choice_cond["output_choice"] != "neg" + + + + + + + + + + + + output_choice_cond["output_choice"] != "pos" + + + + + + Bio + + + +**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. + + +