annotate draw_amr_matrix.py @ 3:7d7884f2d921 draft

Uploaded
author greg
date Fri, 10 Feb 2023 20:04:31 +0000
parents 5c923c77cf5f
children 33a0ea992043
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
1 #!/usr/bin/env python
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
2
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
3 import argparse
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
4 import csv
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
5 import os
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
6
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
7 import numpy
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
8 import pandas
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
9 import matplotlib.pyplot as pyplot
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
10
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
11
1
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
12 def get_amr_in_feature_hits(amr_feature_hits):
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
13 for k in amr_feature_hits.keys():
0
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
14 if k.lower().find('amr') >= 0:
1
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
15 return amr_feature_hits[k]
0
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
16 return None
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
17
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
18
1
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
19 def draw_amr_matrix(amr_feature_hits_files, amr_deletions_file, amr_mutations_file, amr_gene_drug_file, output_dir):
0
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
20 ofh = open('process_log', 'w')
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
21
1
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
22 # Read amr_feature_hits_files.
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
23 amr_feature_hits = pandas.Series(dtype=object)
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
24 for amr_feature_hits_file in amr_feature_hits_files:
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
25 feature_name = os.path.basename(amr_feature_hits_file)
0
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
26 # Make sure the file is not empty.
1
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
27 if os.path.isfile(amr_feature_hits_file) and os.path.getsize(amr_feature_hits_file) > 0:
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
28 best_hits = pandas.read_csv(filepath_or_buffer=amr_feature_hits_file, sep='\t', header=None)
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
29 ofh.write("\nFeature file %s will be processed\n" % os.path.basename(amr_feature_hits_file))
0
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
30 else:
1
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
31 ofh.write("\nEmpty feature file %s will NOT be processed\n" % os.path.basename(amr_feature_hits_file))
0
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
32 best_hits = None
1
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
33 amr_feature_hits[feature_name] = best_hits
0
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
34
1
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
35 amr_hits = get_amr_in_feature_hits(amr_feature_hits)
0
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
36 ofh.write("\namr_hits:\n%s\n" % str(amr_hits))
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
37 if amr_hits is not None:
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
38 amr_to_draw = pandas.DataFrame(columns=['gene', 'drug'])
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
39 ofh.write("\namr_to_draw:\n%s\n" % str(amr_to_draw))
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
40 # Read amr_drug_gene_file.
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
41 amr_gene_drug = pandas.read_csv(amr_gene_drug_file, index_col=None, sep='\t', quoting=csv.QUOTE_NONE, header=None)
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
42 ofh.write("\namr_gene_drug:\n%s\n" % str(amr_gene_drug))
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
43
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
44 # Roll up AMR gene hits.
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
45 ofh.write("\namr_hits.shape[0]:%s\n" % str(amr_hits.shape[0]))
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
46 if amr_hits.shape[0] > 0:
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
47 for gene_idx, gene in amr_hits.iterrows():
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
48 ofh.write("gene_idx:%s\n" % str(gene_idx))
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
49 ofh.write("gene:%s\n" % str(gene))
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
50 gene_name = gene[3]
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
51 ofh.write("gene_name: %s\n" % str(gene_name))
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
52 ofh.write("amr_gene_drug[0]: %s\n" % str(amr_gene_drug[0]))
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
53 drugs = amr_gene_drug.loc[amr_gene_drug[0] == gene_name, :][1]
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
54 ofh.write("drugs:%s\n" % str(drugs))
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
55 for drug in drugs:
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
56 amr_to_draw = amr_to_draw.append(pandas.Series([gene_name, drug], name=amr_to_draw.shape[0], index=amr_to_draw.columns))
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
57 ofh.write("\amr_to_draw:%s\n" % str(amr_to_draw))
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
58
1
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
59 if amr_mutations_file is not None:
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
60 # TODO: So far, no samples have produced mutations, so we haven't been able
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
61 # to produce a populated VarScan VCF file of mutations - https://github.com/appliedbinf/pima_md/blob/main/pima.py#L2923.
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
62 # The call_amr_mutations Galaxy tool will currently produce this VarScan VCF file, but we need a sample that
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
63 # will produce a populated file. After we find one, we'll need to figure out how to implement this loop
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
64 # https://github.com/appliedbinf/pima_md/blob/main/pima.py#L2925 in a Galaxy tool so that the VarScan VCF
3
7d7884f2d921 Uploaded
greg
parents: 1
diff changeset
65 # file will be converted to the TSV amr_mutations_file that thsi tool expects.
7d7884f2d921 Uploaded
greg
parents: 1
diff changeset
66 amr_mutations = pandas.DataFrame()
0
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
67 # Roll up potentially resistance conferring mutations.
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
68 for mutation_region, mutation_hits in amr_mutations.iteritems():
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
69 for mutation_idx, mutation_hit in mutation_hits.iterrows():
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
70 mutation_name = mutation_region + ' ' + mutation_hit['REF'] + '->' + mutation_hit['ALT']
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
71 amr_to_draw = amr_to_draw.append(pandas.Series([mutation_name, mutation_hit['DRUG']], name=amr_to_draw.shape[0], index=amr_to_draw.columns))
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
72
1
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
73 if amr_deletions_file is not None:
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
74 # TODO: So far, no samples have produced deletions, but we do have all the pices in place
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
75 # within the workflow to receive the amr_deletions_file here, although it is currently
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
76 # always empty...
0
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
77 # Roll up deletions that might confer resistance.
3
7d7884f2d921 Uploaded
greg
parents: 1
diff changeset
78 try:
7d7884f2d921 Uploaded
greg
parents: 1
diff changeset
79 amr_deletions = pandas.read_csv(filepath_or_buffer=amr_deletions_file, sep='\t', header=None)
7d7884f2d921 Uploaded
greg
parents: 1
diff changeset
80 except Exception:
7d7884f2d921 Uploaded
greg
parents: 1
diff changeset
81 amr_deletions = pandas.DataFrame()
1
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
82 if amr_deletions.shape[0] > 0:
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
83 amr_deletions.columns = ['contig', 'start', 'stop', 'name', 'type', 'drug', 'note']
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
84 amr_deletions = amr_deletions.loc[amr_deletions['type'].isin(['large-deletion', 'any']), :]
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
85 for deletion_idx, deleted_gene in amr_deletions.iterrows():
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
86 amr_to_draw = amr_to_draw.append(pandas.Series(['\u0394' + deleted_gene[3], deleted_gene[5]], name=amr_to_draw.shape[0], index=amr_to_draw.columns))
0
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
87
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
88 if amr_to_draw.shape[0] > 1:
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
89 ofh.write("\nDrawing AMR matrix...\n")
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
90 present_genes = amr_to_draw['gene'].unique()
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
91 present_drugs = amr_to_draw['drug'].unique()
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
92 amr_matrix = pandas.DataFrame(0, index=present_genes, columns=present_drugs)
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
93 for hit_idx, hit in amr_to_draw.iterrows():
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
94 amr_matrix.loc[hit[0], hit[1]] = 1
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
95 amr_matrix_png = os.path.join(output_dir, 'amr_matrix.png')
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
96 int_matrix = amr_matrix[amr_matrix.columns].astype(int)
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
97 figure, axis = pyplot.subplots()
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
98 axis.invert_yaxis()
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
99 axis.set_yticks(numpy.arange(0.5, len(amr_matrix.index)), minor=False)
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
100 axis.set_yticklabels(int_matrix.index.values)
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
101 axis.set_xticks(numpy.arange(0.5, len(amr_matrix.columns)), minor=False)
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
102 axis.set_xticklabels(amr_matrix.columns.values, rotation=90)
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
103 axis.xaxis.tick_top()
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
104 axis.xaxis.set_label_position('top')
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
105 pyplot.tight_layout()
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
106 pyplot.savefig(amr_matrix_png, dpi=300)
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
107 else:
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
108 ofh.write("\nEmpty AMR matrix, nothing to draw...\n")
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
109 ofh.close()
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
110
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
111
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
112 if __name__ == '__main__':
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
113 parser = argparse.ArgumentParser()
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
114
1
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
115 parser.add_argument('--amr_feature_hits_dir', action='store', dest='amr_feature_hits_dir', help='Directory of tabular files containing feature hits')
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
116 parser.add_argument('--amr_deletions_file', action='store', dest='amr_deletions_file', default=None, help='AMR deletions BED file')
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
117 parser.add_argument('--amr_mutations_file', action='store', dest='amr_mutations_file', default=None, help='AMR mutations TSV file')
0
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
118 parser.add_argument('--amr_gene_drug_file', action='store', dest='amr_gene_drug_file', help='AMR_gene_drugs tsv file')
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
119 parser.add_argument('--output_dir', action='store', dest='output_dir', help='Output directory')
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
120
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
121 args = parser.parse_args()
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
122
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
123 # Get thge collection of feature hits files. The collection
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
124 # will be sorted alphabetically and will contain 2 files
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
125 # named something like AMR_CDS_311_2022_12_20.fasta and
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
126 # Incompatibility_Groups_2023_01_01.fasta.
1
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
127 amr_feature_hits_files = []
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
128 for file_name in sorted(os.listdir(args.amr_feature_hits_dir)):
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
129 file_path = os.path.abspath(os.path.join(args.amr_feature_hits_dir, file_name))
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
130 amr_feature_hits_files.append(file_path)
0
6044ca44e9f6 Uploaded
greg
parents:
diff changeset
131
1
5c923c77cf5f Uploaded
greg
parents: 0
diff changeset
132 draw_amr_matrix(amr_feature_hits_files, args.amr_deletions_file, args.amr_mutations_file, args.amr_gene_drug_file, args.output_dir)