Mercurial > repos > greg > draw_amr_matrix
comparison draw_amr_matrix.py @ 8:7fe8ea50a81d draft
Uploaded
author | greg |
---|---|
date | Wed, 15 Mar 2023 13:31:18 +0000 |
parents | 389c98d344ce |
children | 70073df30a06 |
comparison
equal
deleted
inserted
replaced
7:389c98d344ce | 8:7fe8ea50a81d |
---|---|
106 if mutation_regions.shape[1] != 7: | 106 if mutation_regions.shape[1] != 7: |
107 ofh.write("\nMutation regions should be a six column file.\n") | 107 ofh.write("\nMutation regions should be a six column file.\n") |
108 elif mutation_regions.shape[0] == 0: | 108 elif mutation_regions.shape[0] == 0: |
109 ofh.write("\nNo rows in mutation regions file.\n") | 109 ofh.write("\nNo rows in mutation regions file.\n") |
110 else: | 110 else: |
111 """ | |
112 TODO: move this coder to the pima_report tool... | |
113 # Make sure the positions in the BED file fall within the chromosomes provided in the reference sequence. | |
114 for mutation_region in range(mutation_regions.shape[0]): | |
115 mutation_region = mutation_regions.iloc[mutation_region, :] | |
116 if not (mutation_region[0] in reference): | |
117 ofh.write("\nMutation region: %s not found in reference genome.\n" % ' '.join(mutation_region.astype(str))) | |
118 continue | |
119 if not isinstance(mutation_region[1], numpy.int64): | |
120 ofh.write("\nNon-integer found in mutation region start (column 2): %s.\n" % str(mutation_region[1])) | |
121 break | |
122 elif not isinstance(mutation_region[2], numpy.int64): | |
123 ofh.write("\nNon-integer found in mutation region start (column 3): %s.\n" % str(mutation_region[2])) | |
124 break | |
125 if mutation_region[1] <= 0 or mutation_region[2] <= 0: | |
126 ofh.write("\nMutation region %s starts before the reference sequence.\n" % ' '.join(mutation_region.astype(str))) | |
127 if mutation_region[1] > len(reference[mutation_region[0]].seq) or mutation_region[2] > len(reference[mutation_region[0]].seq): | |
128 ofh.write("\nMutation region %s ends after the reference sequence.\n" % ' '.join(mutation_region.astype(str))) | |
129 """ | |
130 for region_i in range(mutation_regions.shape[0]): | 111 for region_i in range(mutation_regions.shape[0]): |
131 region = mutation_regions.iloc[region_i, :] | 112 region = mutation_regions.iloc[region_i, :] |
132 if not region.get('type', default='No Type') in ['snp', 'small-indel', 'any']: | 113 if not region.get('type', default='No Type') in ['snp', 'small-indel', 'any']: |
133 continue | 114 continue |
134 ofh.write("\nFinding AMR mutations for %s.\n" % str(region['name'])) | 115 ofh.write("\nFinding AMR mutations for %s.\n" % str(region['name'])) |
142 ' | awk \'BEGIN{getline < "' + amr_mutation_regions_bed_file + '";printf $0"\\t";', | 123 ' | awk \'BEGIN{getline < "' + amr_mutation_regions_bed_file + '";printf $0"\\t";', |
143 'getline < "' + varscan_vcf_file + '"; getline < "' + varscan_vcf_file + '";print $0}{print}\'', | 124 'getline < "' + varscan_vcf_file + '"; getline < "' + varscan_vcf_file + '";print $0}{print}\'', |
144 '1>' + region_mutations_tsv]) | 125 '1>' + region_mutations_tsv]) |
145 ofh.write("\ncmd:\n%s\n" % cmd) | 126 ofh.write("\ncmd:\n%s\n" % cmd) |
146 run_command(cmd) | 127 run_command(cmd) |
147 """ | |
148 TODO: move this coder to the pima_report tool... | |
149 try: | |
150 region_mutations = pandas.read_csv(region_mutations_tsv, sep='\t', header=0, index_col=False) | |
151 except Exception: | |
152 region_mutations = pandas.DataFrame() | |
153 if region_mutations.shape[0] == 0: | |
154 continue | |
155 # Figure out what kind of mutations are in this region. | |
156 region_mutation_types = pandas.Series(['snp'] * region_mutations.shape[0], name='TYPE', index=region_mutations.index) | |
157 region_mutation_types[region_mutations['REF'].str.len() != region_mutations['ALT'].str.len()] = 'small-indel' | |
158 region_mutation_drugs = pandas.Series(region['drug'] * region_mutations.shape[0], name='DRUG', index=region_mutations.index) | |
159 region_notes = pandas.Series(region['note'] * region_mutations.shape[0], name='NOTE', index=region_mutations.index) | |
160 region_mutations = pandas.concat([region_mutations, region_mutation_types, region_mutation_drugs, region_notes], axis=1) | |
161 region_mutations = region_mutations[['#CHROM', 'POS', 'TYPE', 'REF', 'ALT', 'DRUG', 'NOTE']] | |
162 amr_mutations[region['name']] = region_mutations | |
163 """ | |
164 else: | 128 else: |
165 ofh.write("\nMutation region BED not received.\n") | 129 ofh.write("\nMutation region BED not received.\n") |
166 # Roll up potentially resistance conferring mutations. | 130 # Roll up potentially resistance conferring mutations. |
167 for mutation_region, mutation_hits in amr_mutations.iteritems(): | 131 for mutation_region, mutation_hits in amr_mutations.iteritems(): |
168 for mutation_idx, mutation_hit in mutation_hits.iterrows(): | 132 for mutation_idx, mutation_hit in mutation_hits.iterrows(): |