diff --git a/Snakefile b/Snakefile
index 77c671b09f646b928a2a752dbc832a48749da40e..88eb6eb0c7010032e015dcaf097a7bd840b79f65 100644
--- a/Snakefile
+++ b/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}"
 
 
 ################################################################################
diff --git a/utils/call-complex-regions.py b/utils/call-complex-regions.py
new file mode 100755
index 0000000000000000000000000000000000000000..70012ba3f42fa00e204c820cbc16c2925df3fa7e
--- /dev/null
+++ b/utils/call-complex-regions.py
@@ -0,0 +1,40 @@
+#!/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()