diff BSseeker2/bs_align/output.py @ 1:8b26adf64adc draft default tip

V2.0.5

--- a/BSseeker2/bs_align/output.py	Fri Jul 12 18:47:28 2013 -0400
+++ b/BSseeker2/bs_align/output.py	Tue Nov 05 01:55:39 2013 -0500
@@ -81,3 +81,53 @@
                           ('XG', output_genome))
 
             self.f.write(a)
+
+
+    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,
+              rrbs = False, my_region_serial = -1, my_region_start = 0, my_region_end = 0):
+
+        if self.format == BS_SEEKER1:
+
+            # remove the soft clipped bases from the read
+            # this is done for backwards compatibility with the old format
+            r_start, r_end, _ = get_read_start_end_and_genome_length(cigar)
+            original_BS = original_BS[r_start : r_end]
+
+            if rrbs:
+                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))
+            else:
+                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))
+
+
+        elif self.format == BAM or self.format == SAM:
+
+            a = pysam.AlignedRead()
+            a.qname = qname
+            a.seq = original_BS if strand == '+' else reverse_compl_seq(original_BS)
+            a.flag = flag
+            a.tid = self.chrom_ids[refname]
+            a.pos = pos
+            a.mapq = 255
+            a.cigar = cigar if strand == '+' else list(reversed(cigar))
+            a.rnext = rnext if rnext == -1 else self.chrom_ids[rnext]
+            a.pnext = pnext
+            a.qual= qual
+            if rrbs:
+                a.tags = (('XO', FR),
+                          ('XS', STEVE),
+                          ('NM', N_mismatch),
+                          ('XM', methy),
+                          ('XG', output_genome),
+                          ('YR', my_region_serial),
+                          ('YS', my_region_start),
+                          ('YE', my_region_end)
+                          )
+
+            else:
+                a.tags = (('XO', FR),
+                          ('XS', STEVE),
+                          ('NM', N_mismatch),
+                          ('XM', methy),
+                          ('XG', output_genome))
+
+            self.f.write(a)