Mercurial > repos > greg > draw_amr_matrix

comparison draw_amr_matrix.py @ 10:03240ffe969a draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Tue, 21 Mar 2023 18:46:42 +0000
parents 70073df30a06
children da1c9c1be421
comparison
equal deleted inserted replaced
9:70073df30a06 10:03240ffe969a
104 amr_mutations = pandas.Series(dtype=object) 104 amr_mutations = pandas.Series(dtype=object)
105 if amr_mutation_regions_bed_file is not None: 105 if amr_mutation_regions_bed_file is not None:
106 mutation_regions = pandas.read_csv(amr_mutation_regions_bed_file, header=0, sep='\t', index_col=False) 106 mutation_regions = pandas.read_csv(amr_mutation_regions_bed_file, header=0, sep='\t', index_col=False)
107 # Validate mutation regions. 107 # Validate mutation regions.
108 if mutation_regions.shape[1] != 7: 108 if mutation_regions.shape[1] != 7:
109 efh.write("\nThe selected mutations regions BED file is invalid, it should be a six column file.\n") 109 efh.write("The selected mutations regions BED file is invalid, it should be a six column file.\n")
110 elif mutation_regions.shape[0] == 0: 110 elif mutation_regions.shape[0] == 0:
111 efh.write("\nThere are no rows in the selected mutation regions file.\n") 111 efh.write("There are no rows in the selected mutation regions file.\n")
112 else: 112 else:
113 for region_i in range(mutation_regions.shape[0]): 113 for region_i in range(mutation_regions.shape[0]):
114 region = mutation_regions.iloc[region_i, :] 114 region = mutation_regions.iloc[region_i, :]
115 if region[0] not in reference: 115 if region[0] not in reference:
116 efh.write("\nMutation region '%s' not found in reference genome.\n" % str(region)) 116 efh.write("Mutation region '%s' not found in reference genome.\n" % str(region))
117 break 117 break
118 if not isinstance(region[1], numpy.int64): 118 if not isinstance(region[1], numpy.int64):
119 efh.write("\nNon-integer found in mutation region start (column 2): %s.\n" % str(region[1])) 119 efh.write("Non-integer found in mutation region start (column 2): %s.\n" % str(region[1]))
120 break 120 break
121 if not isinstance(region[2], numpy.int64): 121 if not isinstance(region[2], numpy.int64):
122 efh.write("\nNon-integer found in mutation region start (column 3): %s.\n" % str(region[2])) 122 efh.write("Non-integer found in mutation region start (column 3): %s.\n" % str(region[2]))
123 break 123 break
124 if region[1] <= 0 or region[2] <= 0: 124 if region[1] <= 0 or region[2] <= 0:
125 efh.write("\nMutation region '%s' starts before the reference sequence.\n" % str(region)) 125 efh.write("Mutation region '%s' starts before the reference sequence.\n" % str(region))
126 if region[1] > len(reference[region[0]].seq) or region[2] > len(reference[region[0]].seq): 126 if region[1] > len(reference[region[0]].seq) or region[2] > len(reference[region[0]].seq):
127 efh.write("\nMutation region '%s' ends after the reference sequence.\n" % str(region)) 127 efh.write("Mutation region '%s' ends after the reference sequence.\n" % str(region))
128 if not region.get('type', default='No Type') in ['snp', 'small-indel', 'any']: 128 if not region.get('type', default='No Type') in ['snp', 'small-indel', 'any']:
129 ofh.write("\n\nSkipping mutation region '%s' with invalid type '%s', valid types are 'snp', 'small-indel', 'any'.\n\n" % (str(region), str(region.get('type', default='No Type')))) 129 ofh.write("\n\nSkipping mutation region '%s' with invalid type '%s', valid types are 'snp', 'small-indel', 'any'.\n\n" % (str(region), str(region.get('type', default='No Type'))))
130 continue 130 continue
131 ofh.write("\nFinding AMR mutations for %s.\n" % str(region['name'])) 131 ofh.write("\nFinding AMR mutations for %s.\n" % str(region['name']))
132 region_bed = 'region_%s.bed' % region_i 132 region_bed = 'region_%s.bed' % region_i
133 ofh.write("region_bed: %s\n" % str(region_bed))
133 mutation_regions.loc[[region_i], ].to_csv(path_or_buf=region_bed, sep='\t', header=False, index=False) 134 mutation_regions.loc[[region_i], ].to_csv(path_or_buf=region_bed, sep='\t', header=False, index=False)
135 ofh.write("mutation_regions.loc[[region_i], ]:\n%s\n" % str(mutation_regions.loc[[region_i], ]))
134 region_mutations_tsv = os.path.join(mutation_regions_dir, 'region_%s_mutations.tsv' % region_i) 136 region_mutations_tsv = os.path.join(mutation_regions_dir, 'region_%s_mutations.tsv' % region_i)
135 cmd = ' '.join(['bedtools intersect -nonamecheck -wb -a', 137 ofh.write("region_mutations_tsv: %s\n" % str(region_mutations_tsv))
136 region_bed, 138 cmd = ' '.join(['bedtools intersect',
137 '-b', 139 '-nonamecheck',
138 varscan_vcf_file, 140 '-wb',
141 '-a', region_bed,
142 '-b', varscan_vcf_file,
139 ' | awk \'BEGIN{getline < "' + amr_mutation_regions_bed_file + '";printf $0"\\t";', 143 ' | awk \'BEGIN{getline < "' + amr_mutation_regions_bed_file + '";printf $0"\\t";',
140 'getline < "' + varscan_vcf_file + '"; getline < "' + varscan_vcf_file + '";print $0}{print}\'', 144 'getline < "' + varscan_vcf_file + '"; getline < "' + varscan_vcf_file + '";print $0}{print}\'',
141 '1>' + region_mutations_tsv]) 145 '1>' + region_mutations_tsv])
142 ofh.write("\ncmd:\n%s\n" % cmd) 146 ofh.write("\ncmd:\n%s\n" % cmd)
143 run_command(cmd) 147 run_command(cmd)
148 try:
149 ofh.write("After running command, os.path.getsize((region_mutations_tsv): %s\n" % str(os.path.getsize(region_mutations_tsv)))
150 region_mutations = pandas.read_csv(region_mutations_tsv, sep='\t', header=0, index_col=False)
151 except Exception:
152 continue
153 # Figure out what kind of mutations are in this region.
154 region_mutation_types = pandas.Series(['snp'] * region_mutations.shape[0], name='TYPE', index=region_mutations.index)
155 region_mutation_types[region_mutations['REF'].str.len() != region_mutations['ALT'].str.len()] = 'small-indel'
156 region_mutation_drugs = pandas.Series(region['drug'] * region_mutations.shape[0], name='DRUG', index=region_mutations.index)
157 region_notes = pandas.Series(region['note'] * region_mutations.shape[0], name='NOTE', index=region_mutations.index)
158 region_mutations = pandas.concat([region_mutations, region_mutation_types, region_mutation_drugs, region_notes], axis=1)
159 region_mutations = region_mutations[['#CHROM', 'POS', 'TYPE', 'REF', 'ALT', 'DRUG', 'NOTE']]
160 amr_mutations[region['name']] = region_mutations
144 else: 161 else:
145 ofh.write("\nMutation region BED not received.\n") 162 ofh.write("\nMutation region BED file not received.\n")
146 # Roll up potentially resistance conferring mutations. 163 # Roll up potentially resistance conferring mutations.
164 ofh.write("\n##### Rolling up potentially resistance conferring mutations..\n")
147 for mutation_region, mutation_hits in amr_mutations.iteritems(): 165 for mutation_region, mutation_hits in amr_mutations.iteritems():
166 ofh.write("mutation_region: %s\n" % str(mutation_region))
167 ofh.write("mutation_hits: %s\n" % str(mutation_hits))
148 for mutation_idx, mutation_hit in mutation_hits.iterrows(): 168 for mutation_idx, mutation_hit in mutation_hits.iterrows():
169 ofh.write("mutation_idx: %s\n" % str(mutation_idx))
170 ofh.write("mutation_hit: %s\n" % str(mutation_hit))
149 mutation_name = mutation_region + ' ' + mutation_hit['REF'] + '->' + mutation_hit['ALT'] 171 mutation_name = mutation_region + ' ' + mutation_hit['REF'] + '->' + mutation_hit['ALT']
172 ofh.write("mutation_name: %s\n" % str(mutation_name))
150 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)) 173 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))
151 174
152 if amr_deletions_file not in [None, 'None'] and os.path.getsize(amr_deletions_file) > 0: 175 if amr_deletions_file not in [None, 'None'] and os.path.getsize(amr_deletions_file) > 0:
153 # Roll up deletions that might confer resistance. 176 # Roll up deletions that might confer resistance.
154 try: 177 try: