|
2
|
1 #!/usr/bin/env python
|
|
|
2
|
|
50
|
3 import argparse
|
|
|
4 import logging
|
|
|
5 from sys import stdout
|
|
|
6 import pandas as pd
|
|
|
7 from subprocess import check_call
|
|
|
8 from shutil import rmtree
|
|
|
9 from tempfile import mkdtemp
|
|
|
10 from os.path import isfile
|
|
|
11 # avoid ugly python IOError when stdout output is piped into another program
|
|
|
12 # and then truncated (such as piping to head)
|
|
|
13 from signal import signal, SIGPIPE, SIG_DFL
|
|
|
14 signal(SIGPIPE, SIG_DFL)
|
|
|
15
|
|
2
|
16 tool_description = """
|
|
|
17 Merge PCR duplicates according to random barcode library.
|
|
|
18
|
|
|
19 Barcodes containing uncalled base 'N' are removed. By default output is written
|
|
|
20 to stdout.
|
|
|
21
|
|
|
22 Input:
|
|
|
23 * bed6 file containing alignments with fastq read-id in name field
|
|
8
|
24 * fastq library of random barcodes
|
|
2
|
25
|
|
|
26 Output:
|
|
|
27 * bed6 file with random barcode in name field and number of PCR duplicates as
|
|
|
28 score, sorted by fields chrom, start, stop, strand, name
|
|
|
29
|
|
|
30 Example usage:
|
|
|
31 - read PCR duplicates from file duplicates.bed and write merged results to file
|
|
|
32 merged.bed:
|
|
|
33 merge_pcr_duplicates.py duplicates.bed bclibrary.fa --out merged.bed
|
|
|
34 """
|
|
|
35
|
|
|
36 epilog = """
|
|
|
37 Author: Daniel Maticzka
|
|
|
38 Copyright: 2015
|
|
|
39 License: Apache
|
|
|
40 Email: maticzkd@informatik.uni-freiburg.de
|
|
|
41 Status: Testing
|
|
|
42 """
|
|
|
43
|
|
|
44 # parse command line arguments
|
|
|
45 parser = argparse.ArgumentParser(description=tool_description,
|
|
|
46 epilog=epilog,
|
|
|
47 formatter_class=argparse.RawDescriptionHelpFormatter)
|
|
|
48 # positional arguments
|
|
|
49 parser.add_argument(
|
|
|
50 "alignments",
|
|
|
51 help="Path to bed6 file containing alignments.")
|
|
|
52 parser.add_argument(
|
|
|
53 "bclib",
|
|
50
|
54 help="Path to fastq barcode library.")
|
|
2
|
55 # optional arguments
|
|
|
56 parser.add_argument(
|
|
|
57 "-o", "--outfile",
|
|
|
58 help="Write results to this file.")
|
|
|
59 # misc arguments
|
|
|
60 parser.add_argument(
|
|
|
61 "-v", "--verbose",
|
|
|
62 help="Be verbose.",
|
|
|
63 action="store_true")
|
|
|
64 parser.add_argument(
|
|
|
65 "-d", "--debug",
|
|
|
66 help="Print lots of debugging information",
|
|
|
67 action="store_true")
|
|
|
68 parser.add_argument(
|
|
|
69 '--version',
|
|
|
70 action='version',
|
|
50
|
71 version='0.2.0')
|
|
2
|
72
|
|
|
73 args = parser.parse_args()
|
|
|
74
|
|
|
75 if args.debug:
|
|
|
76 logging.basicConfig(level=logging.DEBUG, format="%(asctime)s - %(filename)s - %(levelname)s - %(message)s")
|
|
|
77 elif args.verbose:
|
|
|
78 logging.basicConfig(level=logging.INFO, format="%(filename)s - %(levelname)s - %(message)s")
|
|
|
79 else:
|
|
|
80 logging.basicConfig(format="%(filename)s - %(levelname)s - %(message)s")
|
|
|
81 logging.info("Parsed arguments:")
|
|
|
82 logging.info(" alignments: '{}'".format(args.alignments))
|
|
|
83 logging.info(" bclib: '{}'".format(args.bclib))
|
|
|
84 if args.outfile:
|
|
|
85 logging.info(" outfile: enabled writing to file")
|
|
|
86 logging.info(" outfile: '{}'".format(args.outfile))
|
|
|
87 logging.info("")
|
|
|
88
|
|
50
|
89 # see if alignments are empty and the tool can quit
|
|
|
90 n_alns = sum(1 for line in open(args.alignments))
|
|
|
91 if n_alns == 0:
|
|
|
92 logging.warning("WARNING: Working on empty set of alignments, writing empty output.")
|
|
|
93 eventalnout = (open(args.outfile, "w") if args.outfile is not None else stdout)
|
|
|
94 eventalnout.close()
|
|
|
95 exit(0)
|
|
|
96
|
|
|
97 # check input filenames
|
|
|
98 if not isfile(args.bclib):
|
|
|
99 raise Exception("ERROR: barcode library '{}' not found.")
|
|
|
100 if not isfile(args.alignments):
|
|
|
101 raise Exception("ERROR: alignments '{}' not found.")
|
|
|
102
|
|
|
103 try:
|
|
|
104 tmpdir = mkdtemp()
|
|
|
105 logging.debug("tmpdir: " + tmpdir)
|
|
2
|
106
|
|
50
|
107 # prepare barcode library
|
|
|
108 syscall1 = "cat " + args.bclib + " | awk 'BEGIN{OFS=\"\\t\"}NR%4==1{gsub(/^@/,\"\"); id=$1}NR%4==2{bc=$1}NR%4==3{print id,bc}' | sort -k1,1 > " + tmpdir + "/bclib.csv"
|
|
|
109 check_call(syscall1, shell=True)
|
|
|
110
|
|
|
111 # prepare alinments
|
|
|
112 syscall2 = "cat " + args.alignments + " | awk -F \"\\t\" 'BEGIN{OFS=\"\\t\"}{split($4, a, \" \"); $4 = a[1]; print}'| sort -k4,4 > " + tmpdir + "/alns.csv"
|
|
|
113 check_call(syscall2, shell=True)
|
|
2
|
114
|
|
50
|
115 # join barcode library and alignments
|
|
|
116 syscall3 = "join -1 1 -2 4 " + tmpdir + "/bclib.csv " + tmpdir + "/alns.csv " + " | awk 'BEGIN{OFS=\"\\t\"}{print $3,$4,$5,$2,$6,$7}' > " + tmpdir + "/bcalib.csv"
|
|
|
117 check_call(syscall3, shell=True)
|
|
14
|
118
|
|
50
|
119 # get alignments combined with barcodes
|
|
|
120 bcalib = pd.read_csv(
|
|
|
121 tmpdir + "/bcalib.csv",
|
|
|
122 sep="\t",
|
|
|
123 names=["chrom", "start", "stop", "bc", "score", "strand"])
|
|
|
124 finally:
|
|
|
125 logging.debug("removed tmpdir: " + tmpdir)
|
|
|
126 rmtree(tmpdir)
|
|
|
127
|
|
|
128 # fail if alignments given but combined library is empty
|
|
2
|
129 if bcalib.empty:
|
|
|
130 raise Exception("ERROR: no common entries for alignments and barcode library found. Please check your input files.")
|
|
50
|
131
|
|
|
132 # warn if not all alignments could be assigned a barcode
|
|
2
|
133 n_bcalib = len(bcalib.index)
|
|
|
134 if n_bcalib < n_alns:
|
|
|
135 logging.warning(
|
|
50
|
136 "{} of {} alignments could not be associated with a random barcode.".format(n_alns - n_bcalib, n_alns))
|
|
2
|
137
|
|
|
138 # remove entries with barcodes that has uncalled base N
|
|
|
139 bcalib_cleaned = bcalib.drop(bcalib[bcalib.bc.str.contains("N")].index)
|
|
|
140 n_bcalib_cleaned = len(bcalib_cleaned)
|
|
50
|
141 # if n_bcalib_cleaned < n_bcalib:
|
|
|
142 # msg = "{} of {} alignments had random barcodes containing uncalled bases and were dropped.".format(
|
|
|
143 # n_bcalib - n_bcalib_cleaned, n_bcalib)
|
|
|
144 # if n_bcalib_cleaned < (0.8 * n_bcalib):
|
|
|
145 # logging.warning(msg)
|
|
|
146 # else:
|
|
|
147 # logging.info(msg)
|
|
2
|
148
|
|
|
149 # count and merge pcr duplicates
|
|
|
150 # grouping sorts by keys, so the ouput will be properly sorted
|
|
|
151 merged = bcalib_cleaned.groupby(['chrom', 'start', 'stop', 'strand', 'bc']).size().reset_index()
|
|
|
152 merged.rename(columns={0: 'ndupes'}, copy=False, inplace=True)
|
|
|
153
|
|
|
154 # write coordinates of crosslinking event alignments
|
|
|
155 eventalnout = (open(args.outfile, "w") if args.outfile is not None else stdout)
|
|
|
156 merged.to_csv(
|
|
|
157 eventalnout,
|
|
|
158 columns=['chrom', 'start', 'stop', 'bc', 'ndupes', 'strand'],
|
|
|
159 sep="\t", index=False, header=False)
|
|
|
160 eventalnout.close()
|
