# HG changeset patch # User greg # Date 1679424402 0 # Node ID 03240ffe969a9536dab16c8566a4d61ef7cc845f # Parent 70073df30a0633f30f4308327eb3b32244c79282 Uploaded diff -r 70073df30a06 -r 03240ffe969a draw_amr_matrix.py --- a/draw_amr_matrix.py Tue Mar 21 14:06:11 2023 +0000 +++ b/draw_amr_matrix.py Tue Mar 21 18:46:42 2023 +0000 @@ -106,47 +106,70 @@ mutation_regions = pandas.read_csv(amr_mutation_regions_bed_file, header=0, sep='\t', index_col=False) # Validate mutation regions. if mutation_regions.shape[1] != 7: - efh.write("\nThe selected mutations regions BED file is invalid, it should be a six column file.\n") + efh.write("The selected mutations regions BED file is invalid, it should be a six column file.\n") elif mutation_regions.shape[0] == 0: - efh.write("\nThere are no rows in the selected mutation regions file.\n") + efh.write("There are no rows in the selected mutation regions file.\n") else: for region_i in range(mutation_regions.shape[0]): region = mutation_regions.iloc[region_i, :] if region[0] not in reference: - efh.write("\nMutation region '%s' not found in reference genome.\n" % str(region)) + efh.write("Mutation region '%s' not found in reference genome.\n" % str(region)) break if not isinstance(region[1], numpy.int64): - efh.write("\nNon-integer found in mutation region start (column 2): %s.\n" % str(region[1])) + efh.write("Non-integer found in mutation region start (column 2): %s.\n" % str(region[1])) break if not isinstance(region[2], numpy.int64): - efh.write("\nNon-integer found in mutation region start (column 3): %s.\n" % str(region[2])) + efh.write("Non-integer found in mutation region start (column 3): %s.\n" % str(region[2])) break if region[1] <= 0 or region[2] <= 0: - efh.write("\nMutation region '%s' starts before the reference sequence.\n" % str(region)) + efh.write("Mutation region '%s' starts before the reference sequence.\n" % str(region)) if region[1] > len(reference[region[0]].seq) or region[2] > len(reference[region[0]].seq): - efh.write("\nMutation region '%s' ends after the reference sequence.\n" % str(region)) + efh.write("Mutation region '%s' ends after the reference sequence.\n" % str(region)) if not region.get('type', default='No Type') in ['snp', 'small-indel', 'any']: 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')))) continue ofh.write("\nFinding AMR mutations for %s.\n" % str(region['name'])) region_bed = 'region_%s.bed' % region_i + ofh.write("region_bed: %s\n" % str(region_bed)) mutation_regions.loc[[region_i], ].to_csv(path_or_buf=region_bed, sep='\t', header=False, index=False) + ofh.write("mutation_regions.loc[[region_i], ]:\n%s\n" % str(mutation_regions.loc[[region_i], ])) region_mutations_tsv = os.path.join(mutation_regions_dir, 'region_%s_mutations.tsv' % region_i) - cmd = ' '.join(['bedtools intersect -nonamecheck -wb -a', - region_bed, - '-b', - varscan_vcf_file, + ofh.write("region_mutations_tsv: %s\n" % str(region_mutations_tsv)) + cmd = ' '.join(['bedtools intersect', + '-nonamecheck', + '-wb', + '-a', region_bed, + '-b', varscan_vcf_file, ' | awk \'BEGIN{getline < "' + amr_mutation_regions_bed_file + '";printf $0"\\t";', 'getline < "' + varscan_vcf_file + '"; getline < "' + varscan_vcf_file + '";print $0}{print}\'', '1>' + region_mutations_tsv]) ofh.write("\ncmd:\n%s\n" % cmd) run_command(cmd) + try: + ofh.write("After running command, os.path.getsize((region_mutations_tsv): %s\n" % str(os.path.getsize(region_mutations_tsv))) + region_mutations = pandas.read_csv(region_mutations_tsv, sep='\t', header=0, index_col=False) + except Exception: + continue + # Figure out what kind of mutations are in this region. + region_mutation_types = pandas.Series(['snp'] * region_mutations.shape[0], name='TYPE', index=region_mutations.index) + region_mutation_types[region_mutations['REF'].str.len() != region_mutations['ALT'].str.len()] = 'small-indel' + region_mutation_drugs = pandas.Series(region['drug'] * region_mutations.shape[0], name='DRUG', index=region_mutations.index) + region_notes = pandas.Series(region['note'] * region_mutations.shape[0], name='NOTE', index=region_mutations.index) + region_mutations = pandas.concat([region_mutations, region_mutation_types, region_mutation_drugs, region_notes], axis=1) + region_mutations = region_mutations[['#CHROM', 'POS', 'TYPE', 'REF', 'ALT', 'DRUG', 'NOTE']] + amr_mutations[region['name']] = region_mutations else: - ofh.write("\nMutation region BED not received.\n") + ofh.write("\nMutation region BED file not received.\n") # Roll up potentially resistance conferring mutations. + ofh.write("\n##### Rolling up potentially resistance conferring mutations..\n") for mutation_region, mutation_hits in amr_mutations.iteritems(): + ofh.write("mutation_region: %s\n" % str(mutation_region)) + ofh.write("mutation_hits: %s\n" % str(mutation_hits)) for mutation_idx, mutation_hit in mutation_hits.iterrows(): + ofh.write("mutation_idx: %s\n" % str(mutation_idx)) + ofh.write("mutation_hit: %s\n" % str(mutation_hit)) mutation_name = mutation_region + ' ' + mutation_hit['REF'] + '->' + mutation_hit['ALT'] + ofh.write("mutation_name: %s\n" % str(mutation_name)) 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)) if amr_deletions_file not in [None, 'None'] and os.path.getsize(amr_deletions_file) > 0: