|
| 1 | +diff --git a/svtyper/classic.py b/svtyper/classic.py |
| 2 | +index 76dc97c..97561f9 100755 |
| 3 | +--- a/svtyper/classic.py |
| 4 | ++++ b/svtyper/classic.py |
| 5 | +@@ -13,6 +13,10 @@ from svtyper.parsers import Vcf, Variant, Sample, confidence_interval |
| 6 | + from svtyper.utils import * |
| 7 | + from svtyper.statistics import bayes_gt |
| 8 | + |
| 9 | ++try: |
| 10 | ++ xrange |
| 11 | ++except NameError: |
| 12 | ++ xrange = range |
| 13 | + # -------------------------------------- |
| 14 | + # define functions |
| 15 | + |
| 16 | +@@ -126,7 +130,7 @@ def sv_genotype(bam_string, |
| 17 | + bam_list.append(pysam.AlignmentFile(b, mode='rb')) |
| 18 | + elif b.endswith('.cram'): |
| 19 | + bam_list.append(pysam.AlignmentFile(b, |
| 20 | +- mode='rc',reference_filename=ref_fasta,format_options=["required_fields=7167"])) |
| 21 | ++ mode='rc',reference_filename=ref_fasta,format_options=[b"required_fields=7167"])) |
| 22 | + else: |
| 23 | + sys.stderr.write('Error: %s is not a valid alignment file (*.bam or *.cram)\n' % b) |
| 24 | + exit(1) |
| 25 | +@@ -413,12 +417,12 @@ def sv_genotype(bam_string, |
| 26 | + out_bam_written_reads = write_alignment(read, out_bam, out_bam_written_reads) |
| 27 | + |
| 28 | + if debug: |
| 29 | +- print '--------------------------' |
| 30 | +- print 'ref_span:', ref_span |
| 31 | +- print 'alt_span:', alt_span |
| 32 | +- print 'ref_seq:', ref_seq |
| 33 | +- print 'alt_seq:', alt_seq |
| 34 | +- print 'alt_clip:', alt_clip |
| 35 | ++ print('--------------------------') |
| 36 | ++ print('ref_span:', ref_span) |
| 37 | ++ print('alt_span:', alt_span) |
| 38 | ++ print('ref_seq:', ref_seq) |
| 39 | ++ print('alt_seq:', alt_seq) |
| 40 | ++ print('alt_clip:', alt_clip) |
| 41 | + |
| 42 | + # in the absence of evidence for a particular type, ignore the reference |
| 43 | + # support for that type as well |
| 44 | +@@ -443,12 +447,12 @@ def sv_genotype(bam_string, |
| 45 | + QR = int(split_weight * ref_seq) + int(disc_weight * ref_span) |
| 46 | + QA = int(split_weight * alt_splitters) + int(disc_weight * alt_span) |
| 47 | + gt_lplist = bayes_gt(QR, QA, is_dup) |
| 48 | +- best, second_best = sorted([ (i, e) for i, e in enumerate(gt_lplist) ], key=lambda(x): x[1], reverse=True)[0:2] |
| 49 | ++ best, second_best = sorted([ (i, e) for i, e in enumerate(gt_lplist) ], key=lambda x: x[1], reverse=True)[0:2] |
| 50 | + gt_idx = best[0] |
| 51 | + |
| 52 | + # print log probabilities of homref, het, homalt |
| 53 | + if debug: |
| 54 | +- print gt_lplist |
| 55 | ++ print(gt_lplist) |
| 56 | + |
| 57 | + # set the overall variant QUAL score and sample specific fields |
| 58 | + var.genotype(sample.name).set_format('GL', ','.join(['%.0f' % x for x in gt_lplist])) |
| 59 | +@@ -573,7 +577,7 @@ def main(): |
| 60 | + def cli(): |
| 61 | + try: |
| 62 | + sys.exit(main()) |
| 63 | +- except IOError, e: |
| 64 | ++ except IOError as e: |
| 65 | + if e.errno != 32: # ignore SIGPIPE |
| 66 | + raise |
| 67 | + |
| 68 | +diff --git a/svtyper/parsers.py b/svtyper/parsers.py |
| 69 | +index 1357b8c..6f116ec 100644 |
| 70 | +--- a/svtyper/parsers.py |
| 71 | ++++ b/svtyper/parsers.py |
| 72 | +@@ -3,15 +3,16 @@ from __future__ import print_function |
| 73 | + import time, re, json, sys |
| 74 | + from collections import Counter |
| 75 | + |
| 76 | +-from svtyper.statistics import mean, stdev, median, upper_mad |
| 77 | ++from svtyper.statistics import mean, stdev, median, upper_mad, xrange |
| 78 | ++ |
| 79 | + |
| 80 | + # ================================================== |
| 81 | + # VCF parsing tools |
| 82 | + # ================================================== |
| 83 | + def confidence_interval(var, tag, alt_tag, max_ci_dist): |
| 84 | +- ci = map(int, var.info[tag].split(',')) |
| 85 | ++ ci = list(map(int, var.info[tag].split(','))) |
| 86 | + if ci[1] - ci[0] > max_ci_dist: |
| 87 | +- return map(int, var.info[alt_tag].split(',')) |
| 88 | ++ return list(map(int, var.info[alt_tag].split(','))) |
| 89 | + return ci |
| 90 | + |
| 91 | + |
| 92 | +@@ -480,7 +481,7 @@ class Library(object): |
| 93 | + for r in bam.header['RG']: |
| 94 | + try: |
| 95 | + in_lib = r['LB'] == lib_name |
| 96 | +- except KeyError, e: |
| 97 | ++ except KeyError as e: |
| 98 | + in_lib = lib_name == '' |
| 99 | + |
| 100 | + if in_lib: |
| 101 | +@@ -520,8 +521,10 @@ class Library(object): |
| 102 | + for read in self.bam.fetch(): |
| 103 | + if read.get_tag('RG') not in self.readgroups: |
| 104 | + continue |
| 105 | +- if read.infer_query_length() > max_rl: |
| 106 | +- max_rl = read.infer_query_length() |
| 107 | ++ rl = read.infer_query_length() |
| 108 | ++ if rl is None: continue |
| 109 | ++ if rl > max_rl: |
| 110 | ++ max_rl = rl |
| 111 | + if counter == num_samp: |
| 112 | + break |
| 113 | + counter += 1 |
| 114 | +@@ -669,7 +672,7 @@ class Sample(object): |
| 115 | + for r in bam.header['RG']: |
| 116 | + try: |
| 117 | + lib_name=r['LB'] |
| 118 | +- except KeyError, e: |
| 119 | ++ except KeyError as e: |
| 120 | + lib_name='' |
| 121 | + |
| 122 | + # add the new library |
| 123 | +@@ -877,7 +880,7 @@ class SamFragment(object): |
| 124 | + try: |
| 125 | + p = float(self.lib.dens[ospan_length]) * conc_prior / (conc_prior * self.lib.dens[ospan_length] + disc_prior * (self.lib.dens[ospan_length - var_length])) |
| 126 | + except ZeroDivisionError: |
| 127 | +- p = None |
| 128 | ++ return False |
| 129 | + |
| 130 | + return p > 0.5 |
| 131 | + |
| 132 | +diff --git a/svtyper/singlesample.py b/svtyper/singlesample.py |
| 133 | +index 4acb845..d3b16b0 100644 |
| 134 | +--- a/svtyper/singlesample.py |
| 135 | ++++ b/svtyper/singlesample.py |
| 136 | +@@ -47,17 +47,33 @@ description: Compute genotype of structural variants based on breakpoint depth o |
| 137 | + return args |
| 138 | + |
| 139 | + def ensure_valid_alignment_file(afile): |
| 140 | +- if not (afile.endswith('.bam') or afile.endswith('.cram')): |
| 141 | +- die('Error: %s is not a valid alignment file (*.bam or *.cram)\n' % afile) |
| 142 | ++ print(afile.__class__.__name__, file=sys.stderr) |
| 143 | ++ if not isinstance(afile, str): |
| 144 | ++ afile = afile.decode() |
| 145 | ++ |
| 146 | ++ try: |
| 147 | ++ if not (afile.endswith('.bam') or afile.endswith('.cram')): |
| 148 | ++ die('Error: %s is not a valid alignment file (*.bam or *.cram)\n' % afile) |
| 149 | ++ except TypeError: |
| 150 | ++ if not (afile.endswith(b'.bam') or afile.endswith(b'.cram')): |
| 151 | ++ die('Error: %s is not a valid alignment file (*.bam or *.cram)\n' % afile) |
| 152 | + |
| 153 | + def open_alignment_file(afile, reference_fasta): |
| 154 | + fd = None |
| 155 | +- if afile.endswith('.bam'): |
| 156 | +- fd = pysam.AlignmentFile(afile, mode='rb') |
| 157 | +- elif afile.endswith('.cram'): |
| 158 | +- fd = pysam.AlignmentFile(afile, mode='rc', reference_filename=reference_fasta) |
| 159 | +- else: |
| 160 | +- die('Error: %s is not a valid alignment file (*.bam or *.cram)\n' % afile) |
| 161 | ++ try: |
| 162 | ++ if afile.endswith('.bam'): |
| 163 | ++ fd = pysam.AlignmentFile(afile, mode='rb') |
| 164 | ++ elif afile.endswith('.cram'): |
| 165 | ++ fd = pysam.AlignmentFile(afile, mode='rc', reference_filename=reference_fasta) |
| 166 | ++ else: |
| 167 | ++ die('Error: %s is not a valid alignment file (*.bam or *.cram)\n' % afile) |
| 168 | ++ except TypeError: |
| 169 | ++ if afile.endswith(b'.bam'): |
| 170 | ++ fd = pysam.AlignmentFile(afile, mode='rb') |
| 171 | ++ elif afile.endswith(b'.cram'): |
| 172 | ++ fd = pysam.AlignmentFile(afile, mode='rc', reference_filename=reference_fasta) |
| 173 | ++ else: |
| 174 | ++ die('Error: %s is not a valid alignment file (*.bam or *.cram)\n' % afile) |
| 175 | + |
| 176 | + return fd |
| 177 | + |
| 178 | +@@ -417,7 +433,7 @@ def bayesian_genotype(breakpoint, counts, split_weight, disc_weight, debug): |
| 179 | + |
| 180 | + # the actual bayesian calculation and decision |
| 181 | + gt_lplist = bayes_gt(QR, QA, is_dup) |
| 182 | +- best, second_best = sorted([ (i, e) for i, e in enumerate(gt_lplist) ], key=lambda(x): x[1], reverse=True)[0:2] |
| 183 | ++ best, second_best = sorted([ (i, e) for i, e in enumerate(gt_lplist) ], key=lambda x: x[1], reverse=True)[0:2] |
| 184 | + gt_idx = best[0] |
| 185 | + |
| 186 | + # print log probabilities of homref, het, homalt |
| 187 | +@@ -846,7 +862,7 @@ def main(): |
| 188 | + def cli(): |
| 189 | + try: |
| 190 | + sys.exit(main()) |
| 191 | +- except IOError, e: |
| 192 | ++ except IOError as e: |
| 193 | + if e.errno != 32: # ignore SIGPIPE |
| 194 | + raise |
| 195 | + |
| 196 | +diff --git a/svtyper/statistics.py b/svtyper/statistics.py |
| 197 | +index 51cc3cd..1e005a2 100644 |
| 198 | +--- a/svtyper/statistics.py |
| 199 | ++++ b/svtyper/statistics.py |
| 200 | +@@ -1,6 +1,11 @@ |
| 201 | + import math |
| 202 | + from collections import Counter |
| 203 | + |
| 204 | ++try: |
| 205 | ++ xrange |
| 206 | ++except NameError: |
| 207 | ++ xrange = range |
| 208 | ++ |
| 209 | + # ================================================== |
| 210 | + # Statistical tools |
| 211 | + # ================================================== |
| 212 | +diff --git a/svtyper/utils.py b/svtyper/utils.py |
| 213 | +index e80793e..f3fac57 100644 |
| 214 | +--- a/svtyper/utils.py |
| 215 | ++++ b/svtyper/utils.py |
| 216 | +@@ -47,6 +47,10 @@ def write_sample_json(sample_list, lib_info_file): |
| 217 | + lib_info[sample.name] = s |
| 218 | + |
| 219 | + # write the json file |
| 220 | ++ for k, sample in lib_info.items(): |
| 221 | ++ if isinstance(sample["bam"], bytes): |
| 222 | ++ lib_info[k]["bam"] = sample["bam"].decode() |
| 223 | ++ |
| 224 | + json.dump(lib_info, lib_info_file, indent=4) |
| 225 | + lib_info_file.close() |
| 226 | + |
0 commit comments