Merge pull request #88 from icgc-argo/prep-mutect2-qc.0.1.2.0

Prep mutect2 qc.0.1.2.0

Author: junjun-zhang

tests/data/filter-mutect-calls.filtering-stats

@@ -0,0 +1,37 @@
+#<METADATA>Ln prior of deletion of length 10=-20.72326583694641
+#<METADATA>Ln prior of deletion of length 9=-20.72326583694641
+#<METADATA>Ln prior of deletion of length 8=-20.72326583694641
+#<METADATA>Ln prior of deletion of length 7=-20.72326583694641
+#<METADATA>Ln prior of deletion of length 6=-20.72326583694641
+#<METADATA>Ln prior of deletion of length 5=-20.72326583694641
+#<METADATA>Ln prior of deletion of length 4=-20.72326583694641
+#<METADATA>Ln prior of deletion of length 3=-18.120246963748748
+#<METADATA>Ln prior of deletion of length 2=-20.72326583694641
+#<METADATA>Ln prior of deletion of length 1=-18.286701229849623
+#<METADATA>Ln prior of SNV=-14.821551779841153
+#<METADATA>Ln prior of insertion of length 1=-20.72326583694641
+#<METADATA>Ln prior of insertion of length 2=-18.254241200832894
+#<METADATA>Ln prior of insertion of length 3=-20.72326583694641
+#<METADATA>Ln prior of insertion of length 4=-20.72326583694641
+#<METADATA>Ln prior of insertion of length 5=-20.72326583694641
+#<METADATA>Ln prior of insertion of length 6=-20.72326583694641
+#<METADATA>Ln prior of insertion of length 7=-20.72326583694641
+#<METADATA>Ln prior of insertion of length 8=-20.72326583694641
+#<METADATA>Ln prior of insertion of length 9=-20.72326583694641
+#<METADATA>Ln prior of insertion of length 10=-20.72326583694641
+#<METADATA>Background beta-binomial cluster=weight = 0.2243, alpha = 1.02, beta = 4.78
+#<METADATA>High-AF beta-binomial cluster=weight = 0.0245, alpha = 9.99, beta = 1.04
+#<METADATA>Binomial cluster=weight = 0.5275, mean = 0.304
+#<METADATA>Binomial cluster=weight = 0.2236, mean = 0.079
+#<METADATA>threshold=0.577
+#<METADATA>fdr=0.178
+#<METADATA>sensitivity=0.768
+filter	FP	FDR	FN	FNR
+weak_evidence	1.77	0.05	3.7	0.1
+strand_bias	3.0	0.09	3.02	0.08
+contamination	0.02	0.0	0.0	0.0
+normal_artifact	0.56	0.02	0.67	0.02
+orientation	0.49	0.01	0.13	0.0
+slippage	0.0	0.0	0.0	0.0
+haplotype	0.45	0.01	0.12	0.0
+germline	0.48	0.01	0.43	0.01

tests/data/merged-mutect-stats.stats

@@ -0,0 +1,2 @@
+statistic	value
+callable	2.834460576E9

tools/prep-mutect2-qc/prep-mutect2-qc.nf

@@ -23,7 +23,7 @@
  */
 
 nextflow.enable.dsl = 2
-version = '0.1.1.0'
+version = '0.1.2.0'
 
 params.qc_files = ""
 params.container_version = ""

tools/prep-mutect2-qc/prep-mutect2-qc.py

@@ -33,9 +33,47 @@ def run_cmd(cmd):
     return
 
 
+def get_filtering_stats_extra_info(file_path):
+    extra_info = {
+        'filtering_stats': {}
+    }
+    # just to make it simpler, go through the file twice
+    with open(file_path, 'r') as fp:
+        for row in fp.readlines():
+            row.strip()
+            row = row.replace('#<METADATA>', '')
+            if row.startswith('threshold=') or row.startswith('fdr=') or row.startswith('sensitivity='):
+                key, value = row.split('=')
+                extra_info[key] = float(value)
+
+    with open(file_path, 'r') as fp:
+        rdr = csv.DictReader(filter(lambda row: row[0] != '#', fp), delimiter='\t')
+        for row in rdr:
+            filter_name = row.pop('filter')
+            values = row
+            for k, v in values.items():
+                values[k] = float(v)
+            extra_info['filtering_stats'][filter_name] = values
+
+    return extra_info
+
+
+def get_callable_stats_extra_info(file_path):
+    extra_info = {}
+    with open(file_path, 'r') as fp:
+        for row in fp.readlines():
+            if row.startswith('callable'):
+                cols = row.strip().split()
+                extra_info.update({
+                    'callable': int(float(cols[1]))
+                })
+                break
+
+    return extra_info
+
+
 def get_contamination_extra_info(file_path):
     extra_info = {}
-    print(file_path)
     with open(file_path, 'r') as fp:
         for row in fp.readlines():
             if row.startswith('sample'):
@@ -53,12 +91,18 @@ def get_contamination_extra_info(file_path):
 def main(args):
     qc_file_patterns = {
         'tumour_contamination': '*.tumour.*_metrics',
-        'normal_contamination': '*.normal.*_metrics'
+        'normal_contamination': '*.normal.*_metrics',
+        'callable_stats': '*.stats',
+        'filtering_stats': '*.filtering-stats'
     }
 
     # only contamination metrics for now, may have more later
     description = {
-        'contamination': 'Cross sample contamination estimated by GATK CalculateContamination tool'
+        'contamination': 'Cross sample contamination estimated by GATK CalculateContamination tool',
+        'callable_stats': 'Number of sites that are considered callable for Mutect stats with read depth equals or '
+                          'is higher than callable-depth which we set to default 10',
+        'filtering_stats': 'Information on the probability threshold chosen to optimize the F score '
+                           'and the number of false positives and false negatives from each filter to be expected from this choice.'
     }
 
     for qc_file_pttn in qc_file_patterns:
@@ -68,21 +112,47 @@ def main(args):
         }
         metrics = {}
         tar_name = None
-        cont_metrics_file = None
-        qc_files = glob.glob(qc_file_patterns[qc_file_pttn])
-        for f in qc_files:
-            if f.endswith('contamination_metrics'):
-                cont_metrics_file = f
-                tar_name = f'{f}.tgz'
-            extra_info['files_in_tgz'].append(f)
+        if qc_file_pttn.endswith('_contamination'):
+            cont_metrics_file = None
+            qc_files = glob.glob(qc_file_patterns[qc_file_pttn])
+            for f in qc_files:
+                if f.endswith('contamination_metrics'):
+                    cont_metrics_file = f
+                    tar_name = f'{f}.tgz'
+                extra_info['files_in_tgz'].append(f)
+
+            extra_info['description'] = description['contamination']
+
+            metrics = get_contamination_extra_info(cont_metrics_file)
+
+        elif qc_file_pttn.endswith('filtering_stats'):
+            filtering_stats_file = None
+            qc_files = glob.glob(qc_file_patterns[qc_file_pttn])
+            for f in qc_files:
+                if f.endswith('.filtering-stats'):
+                    filtering_stats_file = f
+                    tar_name = f'{f}.filtering_metrics.tgz'
+                extra_info['files_in_tgz'].append(f)
+
+            extra_info['description'] = description['filtering_stats']
+
+            metrics = get_filtering_stats_extra_info(filtering_stats_file)
+
+        elif qc_file_pttn.endswith('callable_stats'):
+            callable_stats_file = None
+            qc_files = glob.glob(qc_file_patterns[qc_file_pttn])
+            for f in qc_files:
+                if f.endswith('.stats'):
+                    callable_stats_file = f
+                    tar_name = f'{f}.callable_metrics.tgz'
+                extra_info['files_in_tgz'].append(f)
 
-        extra_info['description'] = description['contamination']
+            extra_info['description'] = description['callable_stats']
 
-        # TODO: populate metrics info
-        metrics = get_contamination_extra_info(cont_metrics_file)
+            metrics = get_callable_stats_extra_info(callable_stats_file)
 
         extra_info.update({"metrics": metrics})
-        extra_json = f'{cont_metrics_file}.extra_info.json'
+        extra_json = f'{qc_file_pttn}.extra_info.json'
         extra_info['files_in_tgz'].append(extra_json)
         with open(extra_json, 'w') as j:
             j.write(json.dumps(extra_info, indent=2))

tools/prep-mutect2-qc/tests/test-1.nf.json

@@ -3,6 +3,8 @@
     "data/141f7e0b869a059953894fbf920e8134.tumour.contamination_metrics",
     "data/141f7e0b869a059953894fbf920e8134.tumour.segmentation_metrics",
     "data/82b4e5f5d6a28262cb3d6974f40d7f9a.normal.contamination_metrics",
-    "data/82b4e5f5d6a28262cb3d6974f40d7f9a.normal.segmentation_metrics"
+    "data/82b4e5f5d6a28262cb3d6974f40d7f9a.normal.segmentation_metrics",
+    "data/merged-mutect-stats.stats",
+    "data/filter-mutect-calls.filtering-stats"
   ]
 }