Use local script to plot SV calls and mark complex regions.

Snakefile

@@ -250,19 +250,19 @@ rule plot_SV_calls:
     input:
         counts = "counts/{sample}/{windows}.txt.gz",
         calls  = "sv_calls/{sample}/{windows}.{bpdens}/{method}.txt",
-        strand = "strand_states/{sample}/final.txt",
+        complex = "sv_calls/{sample}/{windows}.{bpdens}/{method}.complex.tsv",
+        strand = "strand_states/{sample}/{windows}.{bpdens}/final.txt",
         segments = "segmentation2/{sample}/{windows}.{bpdens}.txt"
     output:
         "sv_calls/{sample}/{windows}.{bpdens,selected_j[0-9\\.]+_s[0-9\\.]+}/plots/sv_calls/{method}.{chrom}.pdf"
     log:
         "log/plot_SV_calls/{sample}/{windows}.{bpdens}.{method}.{chrom}.log"
-    params:
-        sv_plot_script = config["sv_plot_script"]
     shell:
         """
-        Rscript {params.sv_plot_script} \
+        Rscript utils/plot-sv-calls.R \
             segments={input.segments} \
             strand={input.strand} \
+            complex={input.complex} \
             calls={input.calls} \
             {input.counts} \
             {wildcards.chrom} \
@@ -280,11 +280,9 @@ rule plot_SV_calls_simulated:
         "sv_calls/simulation{seed}-{window}/{window}_fixed.{bpdens,selected_j[0-9\\.]+_s[0-9\\.]+}/plots/sv_calls/{method}.{chrom}.pdf"
     log:
         "log/plot_SV_calls_simulated/simulation{seed}-{window}/{window}_fixed.{bpdens}.{method}.{chrom}.log"
-    params:
-        sv_plot_script = config["sv_plot_script"]
     shell:
         """
-        Rscript {params.sv_plot_script} \
+        Rscript utils/plot-sv-calls.R \
             segments={input.segments} \
             strand={input.strand} \
             truth={input.truth} \
@@ -575,6 +573,15 @@ rule mosaiClassifier_make_call_biallelic:
     script:
         "utils/mosaiClassifier_call_biallelic.snakemake.R"
 
+rule call_complex_regions:
+    input:
+        calls  = "sv_calls/{sample}/{windows}.{bpdens}/{method}.txt",
+    output:
+        complex = "sv_calls/{sample}/{windows}.{bpdens}/{method}.complex.tsv",
+    log:
+        "log/call_complex_regions/{sample}/{windows}.{bpdens}.{method}.log"
+    shell:
+        "utils/call-complex-regions.py {input.calls} > {output.complex} 2>{log}"
 
 
 ################################################################################

utils/call-complex-regions.py

+#!/usr/bin/env python
+
+import sys
+from argparse import ArgumentParser
+import pandas as pd
+import numpy as np
+
+def is_complex(x):
+	s = set(x)
+	if (len(s) > 1) or ('complex' in s):
+		return True
+	else:
+		return False
+
+
+def main():
+	parser = ArgumentParser(prog='call-complex-regions.py', description=__doc__)
+	parser.add_argument('--window_size', default=5000000, type=int, 
+		help='Window size in bp.')
+
+	parser.add_argument('callset', metavar='CALLSET', help='Callset file (tsv) as output by MosaiClassifier')
+
+	args = parser.parse_args()
+
+	print('Reading', args.callset, file=sys.stderr)
+	calls = pd.read_csv(args.callset, sep='\t')
+
+	# Identify "complex" intervals
+	segments = calls.groupby(by=['chrom','start','end']).sv_call_name.agg({'is_complex':is_complex}).reset_index()
+	
+	complex_segments = segments[segments.is_complex]
+	total_complex = sum(complex_segments.end - complex_segments.start)
+	
+	print('Total amount of complex sequence: {}Mbp'.format(total_complex/1000000), file=sys.stderr)
+	complex_segments[['chrom','start','end']].to_csv(sys.stdout, index=False, sep='\t')
+	#print(complex_segments, file=sys.stderr)
+
+
+if __name__ == '__main__':
+	main()