Mercurial > repos > weilong-guo > bs_seeker2
comparison BSseeker2/bs_align/output.py @ 1:8b26adf64adc draft default tip
V2.0.5
| author | weilong-guo |
|---|---|
| date | Tue, 05 Nov 2013 01:55:39 -0500 |
| parents | e6df770c0e58 |
| children |
comparison
equal
deleted
inserted
replaced
| 0:e6df770c0e58 | 1:8b26adf64adc |
|---|---|
| 79 ('NM', N_mismatch), | 79 ('NM', N_mismatch), |
| 80 ('XM', methy), | 80 ('XM', methy), |
| 81 ('XG', output_genome)) | 81 ('XG', output_genome)) |
| 82 | 82 |
| 83 self.f.write(a) | 83 self.f.write(a) |
| 84 | |
| 85 | |
| 86 def store2(self, qname, flag, N_mismatch, FR, refname, strand, pos, cigar, original_BS, methy, STEVE, rnext = -1, pnext = -1, qual = None, output_genome = None, | |
| 87 rrbs = False, my_region_serial = -1, my_region_start = 0, my_region_end = 0): | |
| 88 | |
| 89 if self.format == BS_SEEKER1: | |
| 90 | |
| 91 # remove the soft clipped bases from the read | |
| 92 # this is done for backwards compatibility with the old format | |
| 93 r_start, r_end, _ = get_read_start_end_and_genome_length(cigar) | |
| 94 original_BS = original_BS[r_start : r_end] | |
| 95 | |
| 96 if rrbs: | |
| 97 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)) | |
| 98 else: | |
| 99 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)) | |
| 100 | |
| 101 | |
| 102 elif self.format == BAM or self.format == SAM: | |
| 103 | |
| 104 a = pysam.AlignedRead() | |
| 105 a.qname = qname | |
| 106 a.seq = original_BS if strand == '+' else reverse_compl_seq(original_BS) | |
| 107 a.flag = flag | |
| 108 a.tid = self.chrom_ids[refname] | |
| 109 a.pos = pos | |
| 110 a.mapq = 255 | |
| 111 a.cigar = cigar if strand == '+' else list(reversed(cigar)) | |
| 112 a.rnext = rnext if rnext == -1 else self.chrom_ids[rnext] | |
| 113 a.pnext = pnext | |
| 114 a.qual= qual | |
| 115 if rrbs: | |
| 116 a.tags = (('XO', FR), | |
| 117 ('XS', STEVE), | |
| 118 ('NM', N_mismatch), | |
| 119 ('XM', methy), | |
| 120 ('XG', output_genome), | |
| 121 ('YR', my_region_serial), | |
| 122 ('YS', my_region_start), | |
| 123 ('YE', my_region_end) | |
| 124 ) | |
| 125 | |
| 126 else: | |
| 127 a.tags = (('XO', FR), | |
| 128 ('XS', STEVE), | |
| 129 ('NM', N_mismatch), | |
| 130 ('XM', methy), | |
| 131 ('XG', output_genome)) | |
| 132 | |
| 133 self.f.write(a) |