Summary improvements

ariba/__init__.py

@@ -39,6 +39,7 @@
     'sequence_variant',
     'summary',
     'summary_cluster',
+    'summary_cluster_variant',
     'summary_sample',
     'tasks',
     'versions',

ariba/summary_cluster.py

@@ -1,4 +1,4 @@
-from ariba import flag, report
+from ariba import flag, report, summary_cluster_variant
 
 class Error (Exception): pass
 
@@ -56,7 +56,7 @@ def line2dict(cls, line):
             d['var_group'] = '.'
         else:
             try:
-                d['var_group'] = d['var_description'].split(':')[3]
+                d['var_group'] = d['var_description'].split(':')[4]
             except:
                 raise Error('Error getting variant group from the following line:\n' + line)
 
@@ -118,6 +118,8 @@ def _to_cluster_summary_assembled(self):
                 return 'yes'
             else:
                 return 'yes_nonunique'
+        elif self.flag.has('assembled_into_one_contig'):
+            return 'interrupted'
         else:
             return 'fragmented'
 
@@ -218,6 +220,7 @@ def _get_known_noncoding_het_snp(data_dict):
             return None
 
 
+
     @staticmethod
     def _get_nonsynonymous_var(data_dict):
         '''if data_dict has a non synonymous variant, return string:
@@ -295,5 +298,26 @@ def known_noncoding_het_snps(self):
         for d in self.data:
             snp_tuple = self._get_known_noncoding_het_snp(d)
             if snp_tuple is not None:
-                snps[snp_tuple[0]] = snp_tuple[1]
+                snp_id = d['var_description'].split(':')[4]
+                if snp_id not in snps:
+                    snps[snp_id] = {}
+                snps[snp_id][snp_tuple[0]] = snp_tuple[1]
         return snps
+
+
+    @classmethod
+    def _get_all_nonsynon_variants_set(cls, data_dicts):
+        variants = set()
+
+        for data_dict in data_dicts:
+            cluster_var = summary_cluster_variant.SummaryClusterVariant(data_dict)
+            if cluster_var.has_nonsynon:
+                variants.add(cluster_var)
+
+        return variants
+
+
+    def gather_data(self):
+        self.summary = self.column_summary_data()
+        self.variants = self._get_all_nonsynon_variants_set(self.data)
+

ariba/summary_cluster_variant.py

@@ -0,0 +1,83 @@
+from ariba import flag, report
+
+class Error (Exception): pass
+
+class SummaryClusterVariant:
+    def __init__(self, data_dict):
+        self._get_nonsynon_variant_data(data_dict)
+
+
+    def __eq__(self, other):
+       return type(other) is type(self) and self.__dict__ == other.__dict__
+
+
+    def __hash__(self):
+        return hash(tuple([self.__dict__[x] for x in sorted(self.__dict__.keys())]))
+
+
+    def __str__(self):
+        if self.has_nonsynon:
+            return ', '.join((str(self.known), self.var_group, str(self.coding), self.var_string, str(self.het_percent)))
+        else:
+            return 'None'
+
+
+    @classmethod
+    def _has_nonsynonymous(cls, data_dict):
+        return data_dict['ref_ctg_effect'] != 'SYN' and \
+          (
+              data_dict['has_known_var'] == '1' or \
+              (data_dict['known_var'] != '1' and (data_dict['ref_ctg_change'] != '.' or data_dict['ref_ctg_effect'] != '.'))
+          )
+
+
+    @classmethod
+    def _get_het_percent(cls, data_dict):
+        if data_dict['gene'] == '1' or data_dict['ref_ctg_effect'] != 'SNP' or data_dict['smtls_alt_nt'] == '.' or ';' in data_dict['smtls_alt_nt']:
+            return None
+        else:
+            nucleotides = [data_dict['ctg_nt']] + data_dict['smtls_alt_nt'].split(',')
+            depths = data_dict['smtls_alt_depth'].split(',')
+
+            if len(nucleotides) != len(depths):
+                raise Error('Mismatch in number of inferred nucleotides from ctg_nt, smtls_alt_nt, smtls_alt_depth columns. Cannot continue\n' + str(data_dict))
+
+            try:
+                var_nucleotide = data_dict['known_var_change'][-1] if data_dict['known_var_change'] != '.' else data_dict['ref_ctg_change'][-1]
+                if var_nucleotide == '.':
+                    return None
+                depths = [int(x) for x in depths]
+                nuc_to_depth = dict(zip(nucleotides, depths))
+                total_depth = sum(depths)
+                var_depth = nuc_to_depth.get(var_nucleotide, 0)
+                return round(100 * var_depth / total_depth, 1)
+            except:
+                return None
+
+
+    def _get_nonsynon_variant_data(self, data_dict):
+        if not SummaryClusterVariant._has_nonsynonymous(data_dict):
+            self.has_nonsynon = False
+            return
+
+        self.has_nonsynon = True
+
+        if data_dict['known_var_change'] == data_dict['ref_ctg_change'] == '.' == data_dict['ref_ctg_effect']:
+            raise Error('Unexpected data in ariba summary... \n' + str(data_dict) + '\n... known_var_change, ref_ctg_change, ref_ctg_effect all equal to ".", but has a non synonymous change. Something is inconsistent. Cannot continue')
+        elif '.' not in [data_dict['known_var_change'], data_dict['ref_ctg_change']] and \
+          data_dict['known_var_change'] != data_dict['ref_ctg_change']:
+            raise Error('Unexpected data in ariba summary... \n' + str(data_dict) + '\n... known_var_change != ref_ctg_change. Cannot continue')
+
+        self.known = data_dict['known_var'] == '1'
+        self.var_group = data_dict['var_group']
+        self.coding = data_dict['gene'] == '1'
+
+        if data_dict['known_var'] == '1' and data_dict['known_var_change'] != '.':
+            self.var_string = data_dict['known_var_change']
+        elif data_dict['ref_ctg_change'] != '.':
+            self.var_string = data_dict['ref_ctg_change']
+        else:
+            self.var_string = data_dict['ref_ctg_effect']
+
+        self.het_percent = SummaryClusterVariant._get_het_percent(data_dict)
+

ariba/summary_sample.py

@@ -4,9 +4,10 @@
 class Error (Exception): pass
 
 class SummarySample:
-    def __init__(self, report_tsv, min_pc_id=90):
+    def __init__(self, report_tsv, min_pc_id=90, only_clusters=None):
         self.report_tsv = report_tsv
         self.min_pc_id = min_pc_id
+        self.only_clusters = only_clusters
         self.clusters = {}
 
 
@@ -15,7 +16,7 @@ def __eq__(self, other):
 
 
     @staticmethod
-    def _load_file(filename, min_pc_id):
+    def _load_file(filename, min_pc_id, only_clusters=None):
         f = pyfastaq.utils.open_file_read(filename)
         clusters = {}
 
@@ -28,11 +29,25 @@ def _load_file(filename, min_pc_id):
 
             data_dict = summary_cluster.SummaryCluster.line2dict(line)
             cluster = data_dict['cluster']
+            if only_clusters is not None and cluster not in only_clusters:
+                continue
+
             if cluster not in clusters:
                 clusters[cluster] = summary_cluster.SummaryCluster(min_pc_id=min_pc_id)
             clusters[cluster].add_data_dict(data_dict)
 
         pyfastaq.utils.close(f)
+
+        to_delete = set()
+
+        for cluster_name, cluster in clusters.items():
+            cluster.gather_data()
+            if cluster.name is None:
+                to_delete.add(cluster_name)
+
+        for name in to_delete:
+            del clusters[name]
+
         return clusters
 
 
@@ -58,7 +73,7 @@ def _variant_column_names_tuples_and_het_snps(self):
 
 
     def run(self):
-        self.clusters = self._load_file(self.report_tsv, self.min_pc_id)
+        self.clusters = self._load_file(self.report_tsv, self.min_pc_id, only_clusters=self.only_clusters)
         self.column_summary_data = self._column_summary_data()
         self.variant_column_names_tuples, self.het_snps = self._variant_column_names_tuples_and_het_snps()
         self.var_groups = self._var_groups()

ariba/tasks/summary.py

@@ -9,66 +9,33 @@ def use_preset(options):
     preset_to_vals = {
         'minimal': {
             'cluster_cols': 'match',
-            'variant_cols': '',
             'col_filter': 'y',
             'row_filter': 'y',
-            'var_groups': 'n',
-            'known_vars': 'n',
-            'novel_vars': 'n'
         },
         'cluster_small': {
             'cluster_cols': 'assembled,match,ref_seq,known_var',
-            'variant_cols': '',
             'col_filter': 'y',
             'row_filter': 'y',
-            'var_groups': 'n',
-            'known_vars': 'n',
-            'novel_vars': 'n'
         },
         'cluster_all': {
             'cluster_cols': 'assembled,match,ref_seq,pct_id,known_var,novel_var',
-            'variant_cols': '',
             'col_filter': 'y',
             'row_filter': 'y',
-            'var_groups': 'n',
-            'known_vars': 'n',
-            'novel_vars': 'n'
         },
         'cluster_var_groups': {
             'cluster_cols': 'assembled,match,ref_seq,pct_id,known_var,novel_var',
-            'variant_cols': 'groups',
             'col_filter': 'y',
             'row_filter': 'y',
-            'var_groups': 'y',
-            'known_vars': 'n',
-            'novel_vars': 'n'
-        },
-        'cluster_known_vars': {
-            'cluster_cols': 'assembled,match,ref_seq,pct_id,known_var,novel_var',
-            'variant_cols': 'groups,grouped,ungrouped',
-            'col_filter': 'y',
-            'row_filter': 'y',
-            'var_groups': 'y',
-            'known_vars': 'y',
-            'novel_vars': 'n'
         },
         'all': {
             'cluster_cols': 'assembled,match,ref_seq,pct_id,known_var,novel_var',
-            'variant_cols': 'groups,grouped,ungrouped,novel',
             'col_filter': 'y',
             'row_filter': 'y',
-            'var_groups': 'y',
-            'known_vars': 'y',
-            'novel_vars': 'y'
         },
         'all_no_filter': {
             'cluster_cols': 'assembled,match,ref_seq,pct_id,known_var,novel_var',
-            'variant_cols': 'groups,grouped,ungrouped,novel',
             'col_filter': 'n',
             'row_filter': 'n',
-            'var_groups': 'y',
-            'known_vars': 'y',
-            'novel_vars': 'y'
         },
     }
 
@@ -77,6 +44,12 @@ def use_preset(options):
     for key, val in preset_to_vals[options.preset].items():
         exec('options.' + key + ' = "' + val + '"')
 
+    if options.preset in {'cluster_var_groups', 'all', 'all_no_filter'}:
+        options.v_groups = True
+
+    if options.preset in {'all', 'all_no_filter'}:
+        options.variants = True
+
     return options
 
 
@@ -93,9 +66,11 @@ def run(options):
         filter_rows=options.col_filter == 'y',
         filter_columns=options.row_filter == 'y',
         min_id=options.min_id,
-        show_known_het=options.het,
         cluster_cols=options.cluster_cols,
-        variant_cols=options.var_cols,
+        make_phandango_tree=(not options.no_tree),
+        only_clusters=None if options.only_cluster is None else {options.only_cluster},
+        show_var_groups=options.v_groups,
+        show_vars=options.variants,
         verbose=options.verbose
     )
     s.run()

ariba/tests/data/summary_gather_unfiltered_output_data.in.1.tsv

@@ -0,0 +1,5 @@
+#ref_name	gene	var_only	flag	reads	cluster	ref_len	ref_base_assembled	pc_ident	ctg	ctg_len	ctg_cov	known_var	var_type	var_seq_type	known_var_change	has_known_var	ref_ctg_change	ref_ctg_effect	ref_start	ref_end	ref_nt	ctg_start	ctg_end	ctg_nt	smtls_total_depth	smtls_alt_nt	smtls_alt_depth	var_description	free_text
+noncoding_ref1	0	0	19	78	noncoding1	120	120	98.33	noncoding1.scaffold.1	279	10.0	1	SNP	n	A14T	1	A14T	SNP	13	13	A	84	84	T	17	.	17	noncoding_ref1:0:0:A14T:id1:ref has wild type, reads have variant so should report	generic description of noncoding1
+noncoding_ref2	0	0	19	78	noncoding2	120	120	98.33	noncoding2.scaffold.1	279	10.0	1	SNP	n	A42T	1	A42T	SNP	42	42	A	84	84	T	17	.	17	noncoding_ref2:0:0:A42T:id2:ref has wild type, reads have variant so should report	generic description of noncoding1
+noncoding_ref2	0	0	19	78	noncoding2	120	120	98.33	noncoding2.scaffold.1	279	10.0	1	SNP	n	A52T	1	A52T	SNP	42	42	A	84	84	T	17	G	20,30	noncoding_ref2:0:0:A42T:id2:ref has wild type, reads have variant so should report	generic description of noncoding1
+presence_absence_ref1	1	0	27	88	presence_absence1	96	96	98.96	presence_absence1.scaffold.1	267	20.1	0	SNP	p	A10V	.	A10V	NONSYN	28	28	C	113	113	T	29	.	29	presence_absence_ref1:1:0:A10V:.:Ref has wild, reads have variant so report	Generic description of presence_absence1

ariba/tests/data/summary_gather_unfiltered_output_data.in.2.tsv

@@ -0,0 +1,6 @@
+#ref_name	gene	var_only	flag	reads	cluster	ref_len	ref_base_assembled	pc_ident	ctg	ctg_len	ctg_cov	known_var	var_type	var_seq_type	known_var_change	has_known_var	ref_ctg_change	ref_ctg_effect	ref_start	ref_end	ref_nt	ctg_start	ctg_end	ctg_nt	smtls_total_depth	smtls_alt_nt	smtls_alt_depth	var_description	free_text
+noncoding_ref1	0	0	19	78	noncoding1	120	120	98.33	noncoding1.scaffold.1	279	50.1	1	SNP	n	A14T	1	A14T	SNP	13	13	A	84	84	T	50	G	40,10	noncoding1:0:0:A14T:id1:ref has wild type, reads have variant so should report	generic description of noncoding1
+noncoding_ref1	0	0	19	78	noncoding1	120	120	98.33	noncoding1.scaffold.1	279	50.1	1	SNP	n	A6G	1	.	.	6	6	G	77	77	G	18	.	18	noncoding1:0:0:A6G:id3:variant in ref and reads so should report	generic description of noncoding1
+noncoding_ref2	0	0	19	78	noncoding2	120	120	98.33	noncoding2.scaffold.1	279	10.0	1	SNP	n	A52T	1	A52T	SNP	42	42	A	84	84	T	17	G	20,30	noncoding_ref2:0:0:A42T:id2:ref has wild type, reads have variant so should report	generic description of noncoding1
+presence_absence1	1	0	27	88	presence_absence1	96	96	98.96	presence_absence1.scaffold.1	267	51.1	0	SNP	p	A10V	.	A10V	NONSYN	28	28	C	113	113	T	29	.	29	presence_absence1:1:0:A10V:.:Ref has wild, reads have variant so report	Generic description of presence_absence1
+variants_only1	1	1	64	12	variants_only1	90	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.

ariba/tests/data/summary_sample_test_column_names_tuples_and_het_snps.tsv

@@ -1,8 +1,8 @@
 #ref_name	gene	var_only	flag	reads	cluster	ref_len	ref_base_assembled	pc_ident	ctg	ctg_len	ctg_cov	known_var	var_type	var_seq_type	known_var_change	has_known_var	ref_ctg_change	ref_ctg_effect	ref_start	ref_end	ref_nt	ctg_start	ctg_end	ctg_nt	smtls_total_depth	smtls_alt_nt	smtls_alt_depth	var_description	free_text
-noncoding1	0	0	19	78	cluster.n	120	120	98.33	noncoding1.scaffold.1	279	35.4	1	SNP	n	A14T	1	A14T	SNP	13	13	A	84	84	T	50	G	40,10	noncoding1:n:A14T:.:ref has wild type, reads have variant so should report	generic description of noncoding1
+noncoding1	0	0	19	78	cluster.n	120	120	98.33	noncoding1.scaffold.1	279	35.4	1	SNP	n	A14T	1	A14T	SNP	13	13	A	84	84	T	50	G	40,10	noncoding1:0:0:A14T:.:ref has wild type, reads have variant so should report	generic description of noncoding1
 noncoding1	0	0	19	78	cluster.n	120	120	98.33	noncoding1.scaffold.1	279	35.4	0	SNP	.	.	.	G15T	SNP	15	15	G	85	85	T	17	.	17	.	generic description of noncoding1
-noncoding1	0	0	19	78	cluster.n	120	120	98.33	noncoding1.scaffold.1	279	35.4	1	SNP	n	A6G	1	.	.	6	6	G	77	77	G	18	.	18	noncoding1:n:A6G:id2:variant in ref and reads so should report	generic description of noncoding1
+noncoding1	0	0	19	78	cluster.n	120	120	98.33	noncoding1.scaffold.1	279	35.4	1	SNP	n	A6G	1	.	.	6	6	G	77	77	G	18	.	18	noncoding1:0:0:A6G:id2:variant in ref and reads so should report	generic description of noncoding1
 noncoding1	0	0	19	78	cluster.n	120	120	98.33	noncoding2.scaffold.1	279	35.4	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	generic description of noncoding2
-presence_absence1	1	0	27	88	cluster.p	96	96	98.96	presence_absence1.scaffold.1	267	35.1	1	SNP	p	A10V	1	A10V	NONSYN	28	28	C	113	113	T	29	.	29	presence_absence1:p:A10V:id3:Ref has wild, reads have variant so report	Generic description of presence_absence1
+presence_absence1	1	0	27	88	cluster.p	96	96	98.96	presence_absence1.scaffold.1	267	35.1	1	SNP	p	A10V	1	A10V	NONSYN	28	28	C	113	113	T	29	.	29	presence_absence1:1:0:A10V:id3:Ref has wild, reads have variant so report	Generic description of presence_absence1
 presence_absence1	1	0	27	88	cluster.p	96	96	98.96	presence_absence1.scaffold.2	267	35.1	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	Generic description of presence_absence2
-variants_only1	1	1	27	64	cluster.v	90	90	100.0	variants_only1.scaffold.1	260	42.4	1	SNP	p	S5T	1	.	.	13	15	A;C;C	96	98	A;C;C	12;13;13	.;.;.	12;13;13	variants_only1:p:S5T:.:Ref and reads have variant so report	Generic description of variants_only1
+variants_only1	1	1	27	64	cluster.v	90	90	100.0	variants_only1.scaffold.1	260	42.4	1	SNP	p	S5T	1	.	.	13	15	A;C;C	96	98	A;C;C	12;13;13	.;.;.	12;13;13	variants_only1:1:0:S5T:.:Ref and reads have variant so report	Generic description of variants_only1

ariba/tests/data/summary_sample_test_column_summary_data.tsv

@@ -1,8 +1,8 @@
 #ref_name	gene	var_only	flag	reads	cluster	ref_len	ref_base_assembled	pc_ident	ctg	ctg_len	ctg_cov	known_var	var_type	var_seq_type	known_var_change	has_known_var	ref_ctg_change	ref_ctg_effect	ref_start	ref_end	ref_nt	ctg_start	ctg_end	ctg_nt	smtls_total_depth	smtls_alt_nt	smtls_alt_depth	var_description	free_text
-noncoding1	0	0	19	78	cluster.n	120	120	98.33	noncoding1.scaffold.1	279	35.4	1	SNP	n	A14T	1	A14T	SNP	13	13	A	84	84	T	17	.	17	noncoding1_n:A14T:id1:ref has wild type, reads have variant so should report	generic description of noncoding1
+noncoding1	0	0	19	78	cluster.n	120	120	98.33	noncoding1.scaffold.1	279	35.4	1	SNP	n	A14T	1	A14T	SNP	13	13	A	84	84	T	17	.	17	noncoding1:0:0:A14T:id1:ref has wild type, reads have variant so should report	generic description of noncoding1
 noncoding1	0	0	19	78	cluster.n	120	120	98.33	noncoding1.scaffold.1	279	35.4	0	SNP	.	.	.	G15T	SNP	15	15	G	85	85	T	17	.	17	.	generic description of noncoding1
-noncoding1	0	0	19	78	cluster.n	120	120	98.33	noncoding1.scaffold.1	279	35.4	1	SNP	n	A6G	1	.	.	6	6	G	77	77	G	18	.	18	noncoding1_n:A6G:id2:variant in ref and reads so should report	generic description of noncoding1
+noncoding1	0	0	19	78	cluster.n	120	120	98.33	noncoding1.scaffold.1	279	35.4	1	SNP	n	A6G	1	.	.	6	6	G	77	77	G	18	.	18	noncoding1:0:0:A6G:id2:variant in ref and reads so should report	generic description of noncoding1
 noncoding1	0	0	19	78	cluster.n	120	120	98.33	noncoding2.scaffold.1	279	35.4	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	generic description of noncoding2
-presence_absence1	1	0	27	88	cluster.p	96	96	98.96	presence_absence1.scaffold.1	267	35.1	1	SNP	p	A10V	1	A10V	NONSYN	28	28	C	113	113	T	29	.	29	presence_absence1:p:A10V:id3:Ref has wild, reads have variant so report	Generic description of presence_absence1
+presence_absence1	1	0	27	88	cluster.p	96	96	98.96	presence_absence1.scaffold.1	267	35.1	1	SNP	p	A10V	1	A10V	NONSYN	28	28	C	113	113	T	29	.	29	presence_absence1:1:0:A10V:id3:Ref has wild, reads have variant so report	Generic description of presence_absence1
 presence_absence1	1	0	27	88	cluster.p	96	96	98.96	presence_absence1.scaffold.2	267	35.1	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	Generic description of presence_absence2
-variants_only1	1	1	27	64	cluster.v	90	90	100.0	variants_only1.scaffold.1	260	42.4	1	SNP	p	S5T	1	.	.	13	15	A;C;C	96	98	A;C;C	12;13;13	.;.;.	12;13;13	variants_only1:p:S5T:id4:Ref and reads have variant so report	Generic description of variants_only1
+variants_only1	1	1	27	64	cluster.v	90	90	100.0	variants_only1.scaffold.1	260	42.4	1	SNP	p	S5T	1	.	.	13	15	A;C;C	96	98	A;C;C	12;13;13	.;.;.	12;13;13	variants_only1:1:0:S5T:id4:Ref and reads have variant so report	Generic description of variants_only1