update for bsf-call

yutaka-saito · Mar 18, 2015 · 4095522 · 4095522
1 parent 0f0eda2
commit 4095522
Show file tree

Hide file tree

Showing 3 changed files with 17 additions and 9 deletions.
diff --git a/bsf-call/bsf-call b/bsf-call/bsf-call
@@ -8,7 +8,7 @@ is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Un
 http://creativecommons.org/licenses/by-nc-sa/3.0/
 """
 
-__version__= "1.2"
+__version__= "1.3"
 
 from optparse import OptionParser
 import os
@@ -28,7 +28,7 @@ op.add_option("-c", "--coverage", type="int", default=5, metavar="C",
 # op.add_option("-d", "--pe-direction", type="string", default="ff",
 #               help="direction of paired-end probes: ff, fr, rf, rr (default: %default)")
 # 
-op.add_option("-l", "--lower-bound", type="float", default=0.005, metavar="L",
+op.add_option("-l", "--lower-bound", type="float", default=0.05, metavar="L",
               help="threshold of mC ratio (default: %default)")
 
 op.add_option("-p", "--multi-thread", type="int", default=1, metavar="P",

diff --git a/bsf-call/bsf-call.txt b/bsf-call/bsf-call.txt
@@ -11,7 +11,7 @@ OPTIONS
 		Specifies the threshold for valid read coverage for mC detection (default: 5).
 
 	-l <float>, --lover-bound=<float>
-		Specifies the threshold for valid mC rate (default: 0.01).
+		Specifies the threshold for valid mC rate (default: 0.05).
 
 	-s <integer>
 		Specifies the threshold for valid alignment score (default: 150)

diff --git a/bsf-call/bsfcall.py b/bsf-call/bsfcall.py
@@ -8,7 +8,7 @@
 http://creativecommons.org/licenses/by-nc-sa/3.0/
 """
 
-__version__= "1.2"
+__version__= "1.3"
 
 import sys
 import os
@@ -1450,13 +1450,19 @@ def processMafFile(self, resultFile):
                             read_seq = block[2].split()[-1]
                             read_len = len(read_seq)
                             read_qual = block[3].split()[-1]
-                            # strand = self.findStrandFromAlignment(chr, ref_start, read_seq, read_len)
-                            strand = block[2].split()[4]
+                            strand = self.findStrandFromAlignment(chr, ref_start, read_seq, read_len)
+                            # strand = block[2].split()[4]
                             if strand == '+' or strand == '-':
                                 lines = self.extractMcContextsByOneRead(chr, strand, mismap_prob, ref_seq, ref_start, read_seq, read_qual, read_len)
                                 for line in lines:
                                     fout.write(line)
-                                logging.debug("processMafFile: a maf block is successfully captured.")
+                                logging.debug("processMafFile: a maf block(%s) is successfully captured." % strand)
+                            elif strand == '+-':
+                                for st in ('+', '-'):
+                                    lines = self.extractMcContextsByOneRead(chr, st, mismap_prob, ref_seq, ref_start, read_seq, read_qual, read_len)
+                                    for line in lines:
+                                        fout.write(line)
+                                logging.debug("processMafFile: a maf block(%s) is successfully captured." % strand)
                             else:
                                 logging.debug("processMafFile: a maf block does not show strand-specific info.")
                         else:
@@ -1592,10 +1598,10 @@ def findStrandFromAlignment(self, chr, ref_start, read_seq, read_len):
         ref_seq = self.refGenomeBuf[chr]
         for i in range(read_len):
             j = ref_start + i
-            if ref_seq[j]=='C' and (read_seq[i]=='C' or read_seq[i]=='T'):
+            if ref_seq[j]=='C' and read_seq[i]=='T':
                 plus_sup += 1
                 base_size += 1
-            elif ref_seq[j]=='G' and (read_seq[i]=='G' or read_seq[i]=='A'):
+            elif ref_seq[j]=='G' and read_seq[i]=='A':
                 minus_sup += 1
                 base_size += 1
             elif (ref_seq[j]!='-' and read_seq[i]!='-') and ref_seq[j]!=read_seq[i]:
@@ -1612,6 +1618,8 @@ def findStrandFromAlignment(self, chr, ref_start, read_seq, read_len):
         # if minus_sup > plus_sup and mismatch_rate < 0.05:
         if minus_sup > plus_sup:
             return '-'
+        if plus_sup == minus_sup:
+            return '+-'
         return '.'