poretools-basenames.py

#!/usr/bin/env python

version_info=\
"""
poretools-basenames.py: part of the Personal Identification Pipeline
https://github.com/TeamErlich/personal-identification-pipeline

Copyright (C) 2016 Yaniv Erlich (yaniv@cs.columbia.edu)
All Rights Reserved.
This program is licensed under GPL version 3.
See LICENSE file for full details.

Version 0.3
"""


# In the future, perhaps poretools will add an option to generate
# clearer IDs in the fastq/times file (note: they need to match,
# for later usage of 'minion-merge-fastq-times.py').
#
# TODO: improve I/O error handling (especially with sigpipe on stdout).

import sys,argparse, errno
from os.path import basename

version_info=\
"""
Poretools Basename  - version 0.1
"""

def parse_command_line():
    # Define parameters
    parser = argparse.ArgumentParser(
        formatter_class=argparse.RawDescriptionHelpFormatter,
        description="Poretools basename manipulator",
        version=version_info,
        epilog="""

This script reads fastq/times files generated by poretools,
and removes the path of each file (leaving only the 'basename').

If neither --fastq,--times is specified, an attempt is made
to guess based on the file extension (not recommended).

Example:

    # Generate FASTQ,TIMES by poretools.
    # the files will contain the full path of each FAST5 file.

    $ poretools fastq /server/data/myruns/test1/ > test1.fastq
    $ poretools times /server/data/myruns/test1/ > test1.times

    # Remove the extraneous dirnames:

    $ %(prog)s --fastq test1.fastq > test1.basename.fastq
    $ %(prog)s --times test1.times > test1.basename.times
        """)


    parser.add_argument('--fastq-id', dest='mode', action='store_const',
                        const="fastq-id",
                        help="process FASTQ file, replace each sequence ID" \
                             "with the filename (from 3rd part of seq-ID)")

    parser.add_argument('--fastq', dest='mode', action='store_const',
                        const="fastq", help="process FASTQ file, remove " \
                        "fullpath from 3rd part of sequence ID")

    parser.add_argument('--times', dest='mode', action='store_const',
                        const="times", help="process TIMES file")

    # Positional parameter
    parser.add_argument('filename', metavar='FILE', help='file to process');
    args = parser.parse_args()

    if not args.mode:
        # Guess type if needed
        if args.filename.endswith(".fastq") or args.filename.endswith(".fq"):
            args.mode = "fastq-id"
        elif args.filename.endswith(".times"):
            args.mode = "times"
        else:
            sys.exit("failed to guess '--fastq-id' or '--times' based on input"\
                     " filename (%s)" % (args.filename))

    return args


def process_fastq(filename,change_id_to_filename=False):
    """
    remove the path (dirname) from the last part of the sequence ID
    in the FASTQ file.
    Optionally, set the ID to the filename.
    """
    try:
        f = open(filename,'r')

        seen_ids = set()

        linenum = 0
        seqnum = 0
        while True:
            ##
            ## Read 4 lines from the FASTQ file, validate format
            ##
            id  = f.readline().strip()
            if id == '':
                ## Ugly hack: assume EOF
                break

            linenum = linenum + 1
            # Validate standard FASTQ sequence ID
            if not id.startswith("@"):
                err = "input error in '%s' line %d: invalid sequence ID (%s)" \
                      % (filename, linenum, id)
                sys.exit(err)

            ## Detect and remove path from the Sequence ID
            ## (NOTE: this assumes new poretools format, no ':')
            id_parts = id.split(' ')
            if len(id_parts) != 3:
                err = "input error in '%s' line %d: unrecognized sequence ID " \
                      "format (%s) - is this the output of 'poretools fastq' ?"\
                      % (filename, linenum, id)
                sys.exit(err)

            # TODO: should really validate the sequence as well...
            seq = f.readline().strip()
            linenum = linenum + 1

            id2 = f.readline().strip()
            linenum = linenum + 1
            if not id2.startswith("+"):
                err = "input error in '%s' line %d: invalid 2nd sequence ID (%s)" \
                      % (filename, linenum, id)
                sys.exit(err)

            qual= f.readline().strip()
            linenum = linenum + 1

            seqnum = seqnum + 1

            ##
            ## Valid sequence from poretools fastq, trim the filename.
            ##
            fullpath = id_parts[2]
            b = basename(fullpath)

            if change_id_to_filename:
                if b in seen_ids:
                    err = "input error in '%s' line %d: filename '%s' appears"\
                          "multiple sequences, can not use --fastq-id (will " \
                          "result in multiple sequences with same seq-id" \
                          % (filename, linenum-3, b)
                    sys.exit(err)

                seen_ids.add(b)
                newid = '@' + b
            else:
                newid = ' '.join([id_parts[0], id_parts[1], b])

            # Print a FASTQ with the modified ID line
            print newid
            print seq
            print id2
            print qual

    except IOError as e:
        if e.errno == errno.EPIPE:
            sys.exit(0) # exit silently

        # TODO: this is a cop-out, hard to tell what's the exact error
        #       and give informative,useful error message to the user.
        sys.exit("I/O error: %s" % (str(e)))


def process_times(filename):
    """
    Remove the path (dirname) from the second field of the 'times' file.
    """
    try:
        f=open(filename,'r')
        for linenum,line in enumerate(f):
            err = "input error in '%s' line %d: " % (filename, linenum+1)

            line = line.strip()
            flds = line.split('\t')
            if len(flds) != 11:
                sys.exit(err + "expecting 11 fields, found %d - is this a " \
                         "'poretools times' file?" % (len(flds)))

            # First line: validate, then print as-is
            if linenum==0:
                if flds[0] != "channel":
                    sys.exit(err + "expecting header line (first word: " \
                             "'channel') - is this a 'poretools times' file?")
                print '\t'.join(flds)
                continue


            # other lines - remove path from second field.
            flds[1] = basename(flds[1])
            print '\t'.join(flds)


    except IOError as e:
        if e.errno == errno.EPIPE:
            sys.exit(0) # exit silently

        # TODO: this is a cop-out, hard to tell what's the exact error
        #       and give informative,useful error message to the user.
        sys.exit("I/O error: %s" % (str(e)))



if __name__ == "__main__":
    args = parse_command_line()

    if args.mode == "fastq":
        process_fastq(args.filename, False)
    elif args.mode == "fastq-id":
        process_fastq(args.filename, True)
    elif args.mode == "times":
        process_times(args.filename)