Mercurial > repos > weilong-guo > bs_seeker2
comparison BSseeker2/bs_align/output.py @ 0:e6df770c0e58 draft
Initial upload
author | weilong-guo |
---|---|
date | Fri, 12 Jul 2013 18:47:28 -0400 |
parents | |
children | 8b26adf64adc |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:e6df770c0e58 |
---|---|
1 try : | |
2 import pysam | |
3 except ImportError : | |
4 print "[Error] It seems that you haven't install \"pysam\" package.. Please do it before you run this script." | |
5 exit(-1) | |
6 | |
7 import sys | |
8 from bs_align_utils import * | |
9 | |
10 BAM = 'bam' | |
11 SAM = 'sam' | |
12 BS_SEEKER1 = 'bs_seeker1' | |
13 | |
14 formats = [BAM, SAM, BS_SEEKER1] | |
15 | |
16 class outfile: | |
17 def __init__(self, filename, format, chrom_len, cmd_line, suppress_SAM_header): | |
18 self.filename = filename | |
19 self.format = format | |
20 self.chrom_ids = dict((k, i) for i, k in enumerate(sorted(chrom_len))) | |
21 | |
22 if format == BS_SEEKER1: | |
23 self.f = open(filename, 'w') | |
24 elif format in [SAM, BAM]: | |
25 header = { 'HD' : { 'VN': '1.0'}, | |
26 'SQ' : [ {'LN' : chrom_len[c], 'SN' : c} for c in sorted(chrom_len) ], | |
27 'PG' : [ { 'ID' : 1, 'PN' : 'BS Seeker 2', 'CL' : cmd_line} ] | |
28 } | |
29 self.f = pysam.Samfile(filename, 'w' + ('b' if format == BAM else ('' if suppress_SAM_header else 'h')), header = header) | |
30 | |
31 | |
32 | |
33 def close(self): | |
34 self.f.close() | |
35 | |
36 def store(self, qname, N_mismatch, FR, refname, strand, pos, cigar, original_BS, methy, STEVE, rnext = -1, pnext = -1, qual = None, output_genome = None, | |
37 rrbs = False, my_region_serial = -1, my_region_start = 0, my_region_end = 0): | |
38 | |
39 if self.format == BS_SEEKER1: | |
40 | |
41 # remove the soft clipped bases from the read | |
42 # this is done for backwards compatibility with the old format | |
43 r_start, r_end, _ = get_read_start_end_and_genome_length(cigar) | |
44 original_BS = original_BS[r_start : r_end] | |
45 | |
46 if rrbs: | |
47 self.f.write('%s\t%2d\t%s\t%s%s%s\t%s\t%s\t%s\t%d\n' % (qname, N_mismatch, FR, refname, strand, str(pos+1).zfill(10), output_genome, original_BS, methy, STEVE)) | |
48 else: | |
49 self.f.write('%s\t%2d\t%s\t%s%s%s\t%s\t%s\t%s\t%d\t%d\t%d\t%d\n' % (qname, N_mismatch, FR, refname, strand, str(pos+1).zfill(10), output_genome, original_BS, methy, my_region_serial, my_region_start, my_region_end, STEVE)) | |
50 | |
51 | |
52 elif self.format == BAM or self.format == SAM: | |
53 | |
54 a = pysam.AlignedRead() | |
55 a.qname = qname | |
56 a.seq = original_BS if strand == '+' else reverse_compl_seq(original_BS) | |
57 a.flag = 0x10 if strand == '-' else 0 | |
58 a.tid = self.chrom_ids[refname] | |
59 a.pos = pos | |
60 a.mapq = 255 | |
61 a.cigar = cigar if strand == '+' else list(reversed(cigar)) | |
62 a.rnext = rnext if rnext == -1 else self.chrom_ids[rnext] | |
63 a.pnext = pnext | |
64 a.qual= qual | |
65 if rrbs: | |
66 a.tags = (('XO', FR), | |
67 ('XS', STEVE), | |
68 ('NM', N_mismatch), | |
69 ('XM', methy), | |
70 ('XG', output_genome), | |
71 ('YR', my_region_serial), | |
72 ('YS', my_region_start), | |
73 ('YE', my_region_end) | |
74 ) | |
75 | |
76 else: | |
77 a.tags = (('XO', FR), | |
78 ('XS', STEVE), | |
79 ('NM', N_mismatch), | |
80 ('XM', methy), | |
81 ('XG', output_genome)) | |
82 | |
83 self.f.write(a) |