diff --git a/Snakefile b/Snakefile index 85b7ea398bee9a277364d6d248ebec866ebd90c3..4936aa0e9e0d78ed18f7b0876650f300d2235a97 100644 --- a/Snakefile +++ b/Snakefile @@ -508,7 +508,7 @@ rule segmentation_selection: "log/segmentation_selection/{sample}/{window}_{file_name}.log" params: cellnames = lambda wc: ",".join(cell for cell in CELL_PER_SAMPLE[wc.sample]), - sce_min_distance = 400000, + sce_min_distance = 500000, shell: "./utils/detect_strand_states.py --sce_min_distance {params.sce_min_distance} --output_jointseg {output.jointseg} --output_strand_states {output.strand_states} --samplename {wildcards.sample} --cellnames {params.cellnames} {input.info} {input.counts} {input.jointseg} {input.singleseg} > {log} 2>&1" diff --git a/utils/detect_strand_states.py b/utils/detect_strand_states.py index 0690a905f2cf99faba65d5fdb7cf62dd6ae15900..9e08274a5f767be980524f4d0d975305e4627faa 100755 --- a/utils/detect_strand_states.py +++ b/utils/detect_strand_states.py @@ -168,9 +168,9 @@ def get_strand_state(w, c): if (w is None) or (c is None) or (w+c == 0): return (0,0) r = w/(w+c) - if r < 0.1: + if r < 0.2: return (0,2) - elif r > 0.9: + elif r > 0.8: return (2,0) else: return (1,1) @@ -256,7 +256,7 @@ def main(): print(' ... done.', file=sys.stderr) jointseg = Segmentation(args.jointseg) - jointseg.select_k() + jointseg.select_k(max_diff = 0.5) print('Selected breakpoint numbers for joint segmentation:', file=sys.stderr) for chromosome in sorted(jointseg.selected_k.keys()): print(chromosome, jointseg.selected_k[chromosome], file=sys.stderr) @@ -272,7 +272,7 @@ def main(): print('='*100, filename, file=sys.stderr) print('Processing', filename, file=sys.stderr) singleseg = Segmentation(filename) - singleseg.select_k() + singleseg.select_k(max_diff = 1) for chromosome in singleseg.chromosomes: print(' -- chromosome', chromosome, file=sys.stderr) breaks = singleseg.get_selected_segmentation(chromosome)