annotate remove_tail.py @ 58:bbbae1ee87e0 draft default tip

fix for flexbar with small data issue
author rnateam
date Tue, 16 Feb 2016 10:08:58 -0500
parents 0b9aab6aaebf
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
2
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
1 #!/usr/bin/env python
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
2
50
0b9aab6aaebf Uploaded 16cfcafe8b42055c5dd64e62c42b82b455027a40
rnateam
parents: 2
diff changeset
3 import argparse
0b9aab6aaebf Uploaded 16cfcafe8b42055c5dd64e62c42b82b455027a40
rnateam
parents: 2
diff changeset
4 import logging
0b9aab6aaebf Uploaded 16cfcafe8b42055c5dd64e62c42b82b455027a40
rnateam
parents: 2
diff changeset
5 from sys import stdout
0b9aab6aaebf Uploaded 16cfcafe8b42055c5dd64e62c42b82b455027a40
rnateam
parents: 2
diff changeset
6 from Bio.SeqIO.QualityIO import FastqGeneralIterator
0b9aab6aaebf Uploaded 16cfcafe8b42055c5dd64e62c42b82b455027a40
rnateam
parents: 2
diff changeset
7 # avoid ugly python IOError when stdout output is piped into another program
0b9aab6aaebf Uploaded 16cfcafe8b42055c5dd64e62c42b82b455027a40
rnateam
parents: 2
diff changeset
8 # and then truncated (such as piping to head)
0b9aab6aaebf Uploaded 16cfcafe8b42055c5dd64e62c42b82b455027a40
rnateam
parents: 2
diff changeset
9 from signal import signal, SIGPIPE, SIG_DFL
0b9aab6aaebf Uploaded 16cfcafe8b42055c5dd64e62c42b82b455027a40
rnateam
parents: 2
diff changeset
10 signal(SIGPIPE, SIG_DFL)
0b9aab6aaebf Uploaded 16cfcafe8b42055c5dd64e62c42b82b455027a40
rnateam
parents: 2
diff changeset
11
2
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
12 tool_description = """
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
13 Remove a certain number of nucleotides from the 3'-tails of sequences in fastq
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
14 format.
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
15
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
16 Example usage:
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
17 - remove the last 7 nucleotides from file input.fastq, write result to file
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
18 output.fastq:
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
19 remove_tail.py input.fastq 7 --out output.fastq
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
20 """
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
21
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
22 epilog = """
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
23 Author: Daniel Maticzka
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
24 Copyright: 2015
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
25 License: Apache
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
26 Email: maticzkd@informatik.uni-freiburg.de
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
27 Status: Testing
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
28 """
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
29
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
30 # parse command line arguments
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
31 parser = argparse.ArgumentParser(description=tool_description,
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
32 epilog=epilog,
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
33 formatter_class=argparse.RawDescriptionHelpFormatter)
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
34 # positional arguments
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
35 parser.add_argument(
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
36 "infile",
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
37 help="Path to fastq file.")
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
38 parser.add_argument(
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
39 "length",
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
40 type=int,
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
41 help="Remove this many nts.")
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
42 # optional arguments
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
43 parser.add_argument(
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
44 "-o", "--outfile",
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
45 help="Write results to this file.")
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
46 parser.add_argument(
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
47 "-v", "--verbose",
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
48 help="Be verbose.",
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
49 action="store_true")
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
50 parser.add_argument(
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
51 "-d", "--debug",
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
52 help="Print lots of debugging information",
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
53 action="store_true")
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
54 parser.add_argument(
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
55 '--version',
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
56 action='version',
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
57 version='0.1.0')
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
58
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
59 args = parser.parse_args()
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
60 if args.debug:
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
61 logging.basicConfig(level=logging.DEBUG, format="%(asctime)s - %(filename)s - %(levelname)s - %(message)s")
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
62 elif args.verbose:
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
63 logging.basicConfig(level=logging.INFO, format="%(filename)s - %(levelname)s - %(message)s")
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
64 else:
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
65 logging.basicConfig(format="%(filename)s - %(levelname)s - %(message)s")
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
66 logging.info("Parsed arguments:")
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
67 logging.info(" infile: '{}'".format(args.infile))
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
68 logging.info(" length: '{}'".format(args.length))
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
69 if args.outfile:
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
70 logging.info(" outfile: enabled writing to file")
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
71 logging.info(" outfile: '{}'".format(args.outfile))
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
72 logging.info("")
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
73
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
74 # check length parameter
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
75 if args.length < 0:
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
76 raise ValueError("Length must be a positive integer, is '{}'.".format(args.length))
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
77
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
78 # remove tail
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
79 with (open(args.outfile, "w") if args.outfile is not None else stdout) as samout:
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
80 for header, seq, qual in FastqGeneralIterator(open(args.infile)):
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
81
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
82 # if removing tail would lead to an empty sequence,
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
83 # set sequence to a single N to keep fastq synchronized
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
84 if len(seq) <= args.length:
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
85 seq = "N"
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
86 qual = "B"
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
87 logging.debug("read '{}' was too short to remove full tail".format(header))
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
88 logging.debug("seq: {}".format(seq))
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
89 logging.debug("len(seq): {}".format(len(seq)))
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
90 else:
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
91 seq = seq[0:-args.length]
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
92 qual = qual[0:-args.length]
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
93
de4ea3aa1090 Uploaded
rnateam
parents:
diff changeset
94 samout.write("@%s\n%s\n+\n%s\n" % (header, seq, qual))