Mercurial > repos > greg > draw_amr_matrix

comparison draw_amr_matrix.py @ 8:7fe8ea50a81d draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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():