· 8 years ago · Mar 17, 2018, 10:10 AM
1import re
2import csv
3import textwrap
4import matplotlib.pyplot as plt
5
6from collections import Counter
7
8class ReadError(LookupError):
9 """Raise this when a property of a read is called for and no data exists"""
10 pass
11
12
13class PileupError(LookupError):
14 """Raise this when a Pileup file is corrupt or invalid"""
15 pass
16
17
18class Read(object):
19 def __init__(self):
20 self._read_bases = []
21
22 @property
23 def strand(self):
24 if '.' in self._read_bases:
25 return '+'
26 elif ',' in self._read_bases:
27 return '-'
28 else:
29 raise ReadError('No sequence set')
30
31 @property
32 def sequence(self):
33 if self.strand == '+':
34 return ''.join(self._read_bases)
35 else:
36 return ''.join(self._read_bases[::-1])
37
38 @property
39 def read_length(self):
40 return len(self._read_bases)
41
42 def __repr__(self):
43 return('Strand: {}\nSequence:\n{}'.format(
44 self.strand,
45 textwrap.fill(self.sequence, width=80)))
46
47
48class Pileup(object):
49 def __init__(self):
50 self.Ns_by_position = Counter()
51 self.indels_by_position = Counter()
52 self.mutations_by_position = Counter()
53 self.sequenced = {base: Counter() for base in dna_bases}
54
55 def read(self, infile, start=0, end=None):
56 reads = []
57 self.infile = infile
58
59 with open(self.infile) as pileup, open(self.outfile, 'w') as mutpos:
60 out = csv.writer(mutpos, delimiter='\t')
61 for i, line in enumerate(pileup):
62 # if i >= 2000: break
63 # Strip newlines split on tabs, then unpack and change type.
64 chrom, pos, ref, depth, read_base, _ = line.strip().split('\t')
65 pos, depth = int(pos), int(depth)
66
67 # If an end location is defined and we have reached it, break.
68 if self.end is not None and pos > self.end:
69 break
70
71 # Create new Read and append to list of current reads for all
72 # new reads starting at this location.
73 for m in re.finditer(r'\^.', read_base):
74 read_base = read_base.replace(m.group(), '', 1)
75 reads.append(Read())
76
77 # Loop through indel lengths existing at this position.
78 # For each of the indels at this position:
79 # 1. Record them and remove them from the read_base
80 # 2. Capitalize and check to see if previously found
81 indel_positions, inserts, deletes = [], [], []
82 for length in [int(s) for s in re.findall(r'\d+', read_base)]:
83 regex = r'[+-]{}[acgtnACGTN]{{{}}}'.format(length)
84 m = re.search(regex, read_base)
85 indel_positions.append(
86 m.start() - 1 - read_base[:m.start()].count('$'))
87 read_base = read_base.replace(m.group(), '', 1)
88
89 indel = m.group().upper()
90 if '+' in indel:
91 inserts.append(indel)
92 elif '-' in indel:
93 deletes.append(indel)
94 else:
95 raise PileupError
96
97 # The only double spaced grapheme remaining is the $-symbol
98 # which is placed after a read_base symbol to signal the end
99 # of a read. We must first find all of the $-symbols and their
100 # respective locations and then subtract cumulatively from
101 # their locations as we find more of them. These positions will
102 # then reflect the location of the read_base in the string
103 # after removal of the dollar symbols from the string.
104 closing = []
105 for i, match in enumerate(re.finditer(r'.\$', read_base)):
106 closing.append(match.start() - i)
107
108 # Redefine reads as those that do not need to be closed
109 # Get sequences from all reads to be closed
110 sequences = []
111 for read in reads:
112 if reads.index(read) in closing:
113 sequences.append(read.sequence)
114 reads = [r for r in reads if reads.index(r) not in closing]
115
116 # Continue to reconstruct each read sequence.
117 for i, base in enumerate(read_base.replace('$', '')):
118 # The asterisk is simply a deletion placeholder and
119 # shouldn't be counted.
120 if base is '*':
121 continue
122 elif i in indel_positions:
123 reads[i]._read_bases.append('X')
124 else:
125 reads[i]._read_bases.append(base)
126
127 # Filter out records that do not match filter criteria
128 if pos < self.start:
129 continue
130
131 # Count all times these bases are observed in the read_base
132 A, C, G, T, N = [read_base.upper().count(_) for _ in 'ACGTN']
133
134 # Correct depth at this location by subtracting Ns
135 depth = 1 if depth - N == 0 else depth - N
136
137 # Tally all sums, nucleotides sequenced, and update Counters.
138 for base, amount in zip(dna_bases, (A, C, G, T)):
139 self.sequenced[ref].update(base * amount)
140 self.sequenced[ref].update(ref * depth)
141
142 for i, base in enumerate(list(zip(*sequences))):
143 self.Ns_by_position.update([i + 1] * base.count('N'))
144 self.indels_by_position.update([i + 1] * base.count('X'))
145 self.mutations_by_position.update(
146 [i + 1] * sum([base.count(_) for _ in dna_bases]))
147
148 # Write position summary output to .mutpos file
149 out.writerow([chrom, ref, pos, depth, A, C, G, T, N,
150 ':'.join(inserts), ':'.join(deletes)])
151
152 @property
153 def total_nucleotides(self):
154 return sum(self.sequenced[base][base] for base in dna_bases)
155
156 @property
157 def outfile(self):
158 if self.infile.endswith('.pileup'):
159 return self.infile.replace('.pileup', '.mutpos')
160 else:
161 return ''.join([self.infile, '.mutpos'])
162
163 def plot_statistics(self, title, **kwargs):
164
165 titles = ['Base Substitutions', 'Undeciphered (Ns)', 'Indels']
166
167 xlabels = ['', '', 'Position in Read']
168
169 counters = [self.mutations_by_position,
170 self.Ns_by_position,
171 self.indels_by_position,
172 self.depth_by_reference]
173
174 fig, axes = plt.subplots(3, figsize=(6, 8))
175 fig.tight_layout(h_pad=4)
176
177 plot_data = zip(titles, xlabels, axes, counters)
178
179 for title, xlabel, ax, counter in plot_data:
180 ax = plot_numeric_counter(counter,
181 ax=ax,
182 read_length=self.read_length,
183 num_reads=self.total,
184 title=title,
185 kwargs=kwargs)
186 ax.set_xlabel(xlabel)
187
188 title_on = kwargs.pop('title_on', True)
189
190 if self.id is not '' and title_on is True:
191 fig.suptitle(self.id,
192 y=1.04,
193 fontsize=kwargs.pop('title_fontize', 20))
194
195 return fig, axes
196
197 @property
198 def _summary(self):
199 from statsmodels.stats.proportion import proportion_confint
200
201 header_style = '<div style="text-align:left;"><b>{}<br></b></div>' # noqa
202
203 total_deletions = 0
204 total_insertions = 0
205 total_substitutions = 0
206
207 table = [
208 ['ID:', self.id, '', '', ''],
209 ['Name:', self.name, '', '', ''],
210 ['Description:', self.description, '', '', ''],
211 ['Minimum Depth:', self.min_depth, '', '', ''],
212 ['Clonality:', '[{},{}]'.format(self.c, self.C), '', '', ''],
213 ['', '', '', '', ''],
214
215 ['CONVERSION', 'COUNT', 'FREQUENCY', 'LOWER_CONF', 'UPPER_CONF']]
216
217 # Base substitutions
218 for base in dna_bases:
219 n = self.sequenced[base][base]
220 table.append(['{:,}\'s Sequenced:'.format(base), n, '', '', ''])
221
222 for substitution in dna_bases:
223 if substitution == base:
224 continue
225
226 count = self.sequenced[base][substitution]
227 table.append(['{} to {}:'.format(base, substitution),
228 count,
229 count / n,
230 *proportion_confint(count, n, method='jeffrey')])
231 total_substitutions += count
232
233 table.append(['', '', '', '', ''])
234
235 table.extend([
236 ['Nucleotides Sequenced', total_substitutions, '', '', ''],
237 ['Point Substitutions',
238 total_substitutions,
239 total_substitutions / self.total,
240 *proportion_confint(total_substitutions, self.total)]])
241
242 for insertion, count in sorted(self.insertions.items()):
243 table.append(['Total Insertions',
244 count,
245 count / self.total,
246 *proportion_confint(count,
247 self.total,
248 method='jeffrey')])
249 total_insertions += count
250
251 for deletion, count in sorted(self.deletions.items()):
252 table.append(['Total Deletions',
253 count,
254 count / self.total,
255 *proportion_confint(count,
256 self.total,
257 method='jeffrey')])
258 total_deletions += count
259
260 for name, total in zip(['Insertions', 'Deletions'],
261 [total_insertions, total_deletions]):
262
263 table.append([name,
264 total,
265 total / self.total,
266 *proportion_confint(total,
267 self.total,
268 method='jeffrey')])
269 return table
270
271 def _repr(self):
272 from tabulate import tabulate
273
274 return(tabulate(
275 self._summary,
276 stralign='right',
277 numalign='right',
278 tablefmt='plain'))
279
280 def _repr_html(self):
281 from tabulate import tabulate
282
283 return(tabulate(
284 self._summary,
285 stralign='right',
286 numalign='right',
287 tablefmt='html'))