From 4095522fdf5142d3c0d5046f5a301479495d918c Mon Sep 17 00:00:00 2001 From: Yutaka Saito Date: Wed, 18 Mar 2015 15:43:31 +0900 Subject: [PATCH] update for bsf-call --- bsf-call/bsf-call | 4 ++-- bsf-call/bsf-call.txt | 2 +- bsf-call/bsfcall.py | 20 ++++++++++++++------ 3 files changed, 17 insertions(+), 9 deletions(-) diff --git a/bsf-call/bsf-call b/bsf-call/bsf-call index 3bdb0bc..3b7cec2 100755 --- 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 index 409959f..7114b7d 100644 --- 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 , --lover-bound= - Specifies the threshold for valid mC rate (default: 0.01). + Specifies the threshold for valid mC rate (default: 0.05). -s Specifies the threshold for valid alignment score (default: 150) diff --git a/bsf-call/bsfcall.py b/bsf-call/bsfcall.py index 0d57140..448cca5 100755 --- 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 '.'