comparison tools/seq_select_by_id/seq_select_by_id.py @ 6:91f55ee8fea5 draft

v0.0.11; more tests and assorting minor changes

5:1a83f5ab9e95

http://dx.doi.org/10.1093/bioinformatics/btp163 pmid:19304878.

This script is copyright 2011-2013 by Peter Cock, The James Hutton Institute UK.
All rights reserved. See accompanying text file for licence details (MIT
license).

This is version 0.0.6 of the script.
"""
import sys

def stop_err(msg, err=1):
    sys.stderr.write(msg.rstrip() + "\n")
    sys.exit(err)

if "-v" in sys.argv or "--version" in sys.argv:
    print "v0.0.6"
    sys.exit(0)

#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" % col_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 line.strip() and 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

    try:
        count = writer.write_file(iterator)
    except KeyError, err:
        out_handle.close()
        if name not in records:
            stop_err("Identifier %r not found in sequence file" % name)
        else:
            raise err
    out_handle.close()
else:
    #Avoid overhead of parsing into SeqRecord objects,

    for name in parse_ids(tabular_file, column):
        try:
            out_handle.write(records.get_raw(name))
        except KeyError:
            out_handle.close()
            stop_err("Identifier %r not found in sequence file" % name)
        count += 1
    out_handle.close()

print "Selected %i sequences by ID" % count

6:91f55ee8fea5

http://dx.doi.org/10.1093/bioinformatics/btp163 pmid:19304878.

This script is copyright 2011-2013 by Peter Cock, The James Hutton Institute UK.
All rights reserved. See accompanying text file for licence details (MIT
license).
"""
import sys

def sys_exit(msg, err=1):
    sys.stderr.write(msg.rstrip() + "\n")
    sys.exit(err)

if "-v" in sys.argv or "--version" in sys.argv:
    print "v0.0.9"
    sys.exit(0)

#Parse Command Line
try:
    tabular_file, col_arg, in_file, seq_format, out_file = sys.argv[1:]
except ValueError:
    sys_exit("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:
    sys_exit("Expected column number, got %s" % col_arg)

if seq_format == "fastqcssanger":
    sys_exit("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:
    sys_exit("Unrecognised file format %r" % seq_format)


try:
    from Bio import SeqIO
except ImportError:
    sys_exit("Biopython 1.54 or later is required")


def parse_ids(tabular_file, col):
    """Read tabular file and record all specified identifiers.

    Will print a single warning to stderr if any of the fields have
    non-trailing white space (only the first word will be used as
    the identifier).
    """
    handle = open(tabular_file, "rU")
    warn = False
    for line in handle:
        if line.strip() and not line.startswith("#"):
            field = line.rstrip("\n").split("\t")[col].strip()
            parts = field.split(None, 1)
            if len(parts) > 1 and not warn:
                warn = "WARNING: Some of your identifiers had white space in them, " + \
                       "using first word only. e.g.:\n%s\n" % field
            yield parts[0]
    handle.close()
    if warn:
        sys.stderr.write(warn)

#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:
        sys_exit("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

    try:
        count = writer.write_file(iterator)
    except KeyError, err:
        out_handle.close()
        if name not in records:
            sys_exit("Identifier %r not found in sequence file" % name)
        else:
            raise err
    out_handle.close()
else:
    #Avoid overhead of parsing into SeqRecord objects,

    for name in parse_ids(tabular_file, column):
        try:
            out_handle.write(records.get_raw(name))
        except KeyError:
            out_handle.close()
            sys_exit("Identifier %r not found in sequence file" % name)
        count += 1
    out_handle.close()

print "Selected %i sequences by ID" % count