Mercurial > repos > greg > draw_amr_matrix
changeset 4:33a0ea992043 draft
Uploaded
author | greg |
---|---|
date | Tue, 21 Feb 2023 19:55:42 +0000 (22 months ago) |
parents | 7d7884f2d921 |
children | 46d31ec5d24c |
files | draw_amr_matrix.py draw_amr_matrix.xml macros.xml |
diffstat | 3 files changed, 161 insertions(+), 24 deletions(-) [+] |
line wrap: on
line diff
--- a/draw_amr_matrix.py Fri Feb 10 20:04:31 2023 +0000 +++ b/draw_amr_matrix.py Tue Feb 21 19:55:42 2023 +0000 @@ -3,7 +3,11 @@ import argparse import csv import os +import subprocess +import sys +import tempfile +import Bio.SeqIO import numpy import pandas import matplotlib.pyplot as pyplot @@ -16,7 +20,45 @@ return None -def draw_amr_matrix(amr_feature_hits_files, amr_deletions_file, amr_mutations_file, amr_gene_drug_file, output_dir): +def load_fasta(fasta_file): + sequence = pandas.Series(dtype=object) + for contig in Bio.SeqIO.parse(fasta_file, 'fasta'): + sequence[contig.id] = contig + return sequence + + +def run_command(cmd): + try: + tmp_name = tempfile.NamedTemporaryFile(dir=".").name + tmp_stderr = open(tmp_name, 'wb') + proc = subprocess.Popen(args=cmd, shell=True, stderr=tmp_stderr.fileno()) + returncode = proc.wait() + tmp_stderr.close() + if returncode != 0: + # Get stderr, allowing for case where it's very large. + tmp_stderr = open(tmp_name, 'rb') + stderr = '' + buffsize = 1048576 + try: + while True: + stderr += tmp_stderr.read(buffsize) + if not stderr or len(stderr) % buffsize != 0: + break + except OverflowError: + pass + tmp_stderr.close() + os.remove(tmp_name) + stop_err(stderr) + except Exception as e: + stop_err('Command:\n%s\n\nended with error:\n%s\n\n' % (cmd, str(e))) + + +def stop_err(msg): + sys.stderr.write(msg) + sys.exit(1) + + +def draw_amr_matrix(amr_feature_hits_files, amr_deletions_file, amr_mutations_file, amr_mutation_regions_file, amr_gene_drug_file, reference, reference_size, region_mutations_output_file, mutations_dir, output_dir): ofh = open('process_log', 'w') # Read amr_feature_hits_files. @@ -42,35 +84,80 @@ ofh.write("\namr_gene_drug:\n%s\n" % str(amr_gene_drug)) # Roll up AMR gene hits. - ofh.write("\namr_hits.shape[0]:%s\n" % str(amr_hits.shape[0])) + ofh.write("\namr_hits.shape[0]: %s\n" % str(amr_hits.shape[0])) if amr_hits.shape[0] > 0: for gene_idx, gene in amr_hits.iterrows(): - ofh.write("gene_idx:%s\n" % str(gene_idx)) - ofh.write("gene:%s\n" % str(gene)) + ofh.write("gene_idx: %s\n" % str(gene_idx)) + ofh.write("gene: %s\n" % str(gene)) gene_name = gene[3] ofh.write("gene_name: %s\n" % str(gene_name)) ofh.write("amr_gene_drug[0]: %s\n" % str(amr_gene_drug[0])) drugs = amr_gene_drug.loc[amr_gene_drug[0] == gene_name, :][1] - ofh.write("drugs:%s\n" % str(drugs)) + ofh.write("drugs: %s\n" % str(drugs)) for drug in drugs: amr_to_draw = amr_to_draw.append(pandas.Series([gene_name, drug], name=amr_to_draw.shape[0], index=amr_to_draw.columns)) - ofh.write("\amr_to_draw:%s\n" % str(amr_to_draw)) + ofh.write("\amr_to_draw: %s\n" % str(amr_to_draw)) - if amr_mutations_file is not None: - # TODO: So far, no samples have produced mutations, so we haven't been able - # to produce a populated VarScan VCF file of mutations - https://github.com/appliedbinf/pima_md/blob/main/pima.py#L2923. - # The call_amr_mutations Galaxy tool will currently produce this VarScan VCF file, but we need a sample that - # will produce a populated file. After we find one, we'll need to figure out how to implement this loop - # https://github.com/appliedbinf/pima_md/blob/main/pima.py#L2925 in a Galaxy tool so that the VarScan VCF - # file will be converted to the TSV amr_mutations_file that thsi tool expects. - amr_mutations = pandas.DataFrame() + ofh.write("\namr_mutations_file si None: %s\n" % str(amr_mutations_file == 'None')) + if amr_mutations_file not in [None, 'None'] and os.path.getsize(amr_mutations_file) > 0: + amr_mutations = pandas.Series(dtype=object) + if amr_mutation_regions_file is not None: + mutation_regions = pandas.read_csv(amr_mutation_regions_file, header=0, sep='\t', index_col=False) + if mutation_regions.shape[1] != 7: + ofh.write("\nMutation regions should be a six column file.\n") + elif mutation_regions.shape[0] == 0: + ofh.write("\nNo rows in mutation regions file.\n") + else: + # Make sure the positions in the BED file fall within the chromosomes provided in the reference sequence. + for mutation_region in range(mutation_regions.shape[0]): + mutation_region = mutation_regions.iloc[mutation_region, :] + if not (mutation_region[0] in reference): + ofh.write("\nMutation region :%s not found in reference genome.\n" % ' '.join(mutation_region.astype(str))) + continue + if not isinstance(mutation_region[1], numpy.int64): + ofh.write("\nNon-integer found in mutation region start (column 2): %s.\n" % str(mutation_region[1])) + break + elif not isinstance(mutation_region[2], numpy.int64): + ofh.write("\nNon-integer found in mutation region start (column 3): %s.\n" % str(mutation_region[2])) + break + if mutation_region[1] <= 0 or mutation_region[2] <= 0: + ofh.write("\nMutation region %s starts before the reference sequence.\n" % ' '.join(mutation_region.astype(str))) + if mutation_region[1] > len(reference[mutation_region[0]].seq) or mutation_region[2] > len(reference[mutation_region[0]].seq): + ofh.write("\nMutation region %s ends after the reference sequence.\n" % ' '.join(mutation_region.astype(str))) + for region_i in range(mutation_regions.shape[0]): + region = mutation_regions.iloc[region_i, :] + if not region.get('type', default='No Type') in ['snp', 'small-indel', 'any']: + continue + ofh.write("\nFinding AMR mutations for %s.\n" % str(region['name'])) + region_dir = os.path.join(mutations_dir, 'region_' + str(region_i)) + os.mkdir(region_dir) + region_bed = os.path.join(region_dir, 'region.bed') + mutation_regions.loc[[region_i], ].to_csv(path_or_buf=region_bed, sep='\t', header=False, index=False) + cmd = "bedtools intersect -nonamecheck -wb -a '%s' -b '%s' | awk '{BEGIN{getline < \"%s\" ;printf $0\"\t\";getline < \"%s\"; getline < \"%s\";print $0}{print}' > %s" % (region_bed, amr_mutations_file, amr_mutation_regions_file, amr_mutations_file, amr_mutations_file, region_mutations_output_file) + run_command(cmd) + try: + region_mutations = pandas.read_csv(region_mutations_output_file, sep='\t', header=0, index_col=False) + except Exception: + region_mutations = pandas.DataFrame() + if region_mutations.shape[0] == 0: + 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") # Roll up potentially resistance conferring mutations. for mutation_region, mutation_hits in amr_mutations.iteritems(): for mutation_idx, mutation_hit in mutation_hits.iterrows(): mutation_name = mutation_region + ' ' + mutation_hit['REF'] + '->' + mutation_hit['ALT'] 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 is not None: + if amr_deletions_file not in [None, 'None'] and os.path.getsize(amr_deletions_file) > 0: # TODO: So far, no samples have produced deletions, but we do have all the pices in place # within the workflow to receive the amr_deletions_file here, although it is currently # always empty... @@ -115,12 +202,16 @@ parser.add_argument('--amr_feature_hits_dir', action='store', dest='amr_feature_hits_dir', help='Directory of tabular files containing feature hits') parser.add_argument('--amr_deletions_file', action='store', dest='amr_deletions_file', default=None, help='AMR deletions BED file') parser.add_argument('--amr_mutations_file', action='store', dest='amr_mutations_file', default=None, help='AMR mutations TSV file') + parser.add_argument('--amr_mutation_regions_file', action='store', dest='amr_mutation_regions_file', default=None, help='AMR mutation regions BED file') parser.add_argument('--amr_gene_drug_file', action='store', dest='amr_gene_drug_file', help='AMR_gene_drugs tsv file') + parser.add_argument('--reference_genome', action='store', dest='reference_genome', help='Reference genome fasta file') + parser.add_argument('--region_mutations_output_file', action='store', dest='region_mutations_output_file', default=None, help='Region mutations TSV output file') + parser.add_argument('--mutations_dir', action='store', dest='mutations_dir', help='Mutations directory') parser.add_argument('--output_dir', action='store', dest='output_dir', help='Output directory') args = parser.parse_args() - # Get thge collection of feature hits files. The collection + # Get the collection of feature hits files. The collection # will be sorted alphabetically and will contain 2 files # named something like AMR_CDS_311_2022_12_20.fasta and # Incompatibility_Groups_2023_01_01.fasta. @@ -129,4 +220,10 @@ file_path = os.path.abspath(os.path.join(args.amr_feature_hits_dir, file_name)) amr_feature_hits_files.append(file_path) - draw_amr_matrix(amr_feature_hits_files, args.amr_deletions_file, args.amr_mutations_file, args.amr_gene_drug_file, args.output_dir) + # Load the reference genome into memory. + reference = load_fasta(args.reference_genome) + reference_size = 0 + for i in reference: + reference_size += len(i.seq) + + draw_amr_matrix(amr_feature_hits_files, args.amr_deletions_file, args.amr_mutations_file, args.amr_mutation_regions_file, args.amr_gene_drug_file, reference, reference_size, args.region_mutations_output_file, args.mutations_dir, args.output_dir)
--- a/draw_amr_matrix.xml Fri Feb 10 20:04:31 2023 +0000 +++ b/draw_amr_matrix.xml Tue Feb 21 19:55:42 2023 +0000 @@ -8,8 +8,15 @@ #import re mkdir amr_feature_hits_dir && +mkdir mutations_dir && mkdir output_dir && +#if $reference_source.reference_source_selector == 'history': + ln -f -s '$reference_source.ref_file' reference.fa && +#else: + ln -f -s '$reference_source.ref_file.fields.path' reference.fa && +#end if + #for $i in $amr_feature_hits: #set file_name = $i.file_name #set identifier = re.sub('[^\s\w\-\\.]', '_', str($i.element_identifier)) @@ -17,18 +24,47 @@ #end for python '$__tool_directory__/draw_amr_matrix.py' +--amr_feature_hits_dir 'amr_feature_hits_dir' +#if str($amr_deletions_file) != 'None': + --amr_deletions_file '$amr_deletions_file' +#end if +#if str($amr_mutations_file) != 'None': + --amr_mutations_file '$amr_mutations_file' +#end if +#if str($amr_mutation_regions_file) != 'None': + --amr_mutation_regions_file '$amr_mutation_regions_file' + --region_mutations_output_file '$region_mutations_output_file' +#end if --amr_gene_drug_file '$amr_gene_drug_file' ---amr_deletions_file '$amr_deletions_file' ---amr_feature_hits_dir 'amr_feature_hits_dir' +--reference_genome reference.fa +--mutations_dir 'mutations_dir' --output_dir 'output_dir' #if str($output_process_log) == 'yes': && mv 'process_log' '$process_log' #end if ]]></command> <inputs> + <conditional name="reference_source"> + <param name="reference_source_selector" type="select" label="Select a reference genome from your history or use a cached genome index?"> + <option value="cached">Use a cached genome index</option> + <option value="history">Select a genome from the history and build the index</option> + </param> + <when value="cached"> + <param name="ref_file" type="select" label="Using reference genome" help="Select reference genome"> + <options from_data_table="all_fasta"> + <filter type="sort_by" column="2"/> + <validator type="no_options" message="No reference genomes are available"/> + </options> + </param> + </when> + <when value="history"> + <param name="ref_file" type="data" format="fasta,fastq" label="Select the reference sequence" help="You can upload a FASTA file and use it as reference"/> + </when> + </conditional> <param argument="--amr_feature_hits" format="bed" type="data_collection" collection_type="list" label="Collection of feature hits BED files"/> <param argument="--amr_deletions_file" type="data" format="bed" optional="true" label="AMR deletions file" help="Optional, leave blank to ignore"/> <param argument="--amr_mutations_file" type="data" format="tabular,tsv" optional="true" label="AMR mutations file" help="Optional, leave blank to ignore"/> + <param argument="--amr_mutation_regions_file" type="data" format="bed" optional="true" label="AMR mutation regions BED file" help="Optional, leave blank to ignore"/> <param argument="--amr_gene_drug_file" type="data" format="tabular,tsv" label="AMR gene drugs file"/> <param name="output_process_log" type="select" display="radio" label="Output process log file?"> <option value="no" selected="true">No</option> @@ -36,18 +72,23 @@ </param> </inputs> <outputs> - <data name="process_log" format="txt" label="${tool.name} on ${on_string} (Process log)"> + <data name="process_log" format="txt" label="${tool.name} on ${on_string} (process log)"> <filter>output_process_log == 'yes'</filter> </data> + <data name="region_mutations_output_file" format="tsv" label="${tool.name} on ${on_string} (region mutations)"> + <filter>amr_mutation_regions_file not in [None, 'None']</filter> + </data> <collection name="amr_matrix_png" type="list" format="png"> <discover_datasets pattern="(?P<designation>.+)\.(?P<ext>png)" directory="output_dir"/> </collection> </outputs> <tests> <test> + <param name="reference_source_selector" value="history"/> + <param name="ref_file" ftype="fasta" value="ref_genome.fasta"/> <param name="amr_feature_hits"> <collection type="list"> - <element name="amr_cds.bed" value="amr_cds.bed"/> + <element name="amr_pima_md.bed" value="amr_pima_md.bed"/> </collection> </param> <param name="amr_gene_drug_file" value="amr_gene_drug.tsv" ftype="tsv"/>
--- a/macros.xml Fri Feb 10 20:04:31 2023 +0000 +++ b/macros.xml Tue Feb 21 19:55:42 2023 +0000 @@ -4,13 +4,12 @@ <token name="@PROFILE@">21.01</token> <xml name="requirements"> <requirements> - <requirement type="package" version="9.1">coreutils</requirement> + <requirement type="package" version="2.30.0">bedtools</requirement> + <requirement type="package" version="1.79">biopython</requirement> <requirement type="package" version="5.1.0">gawk</requirement> - <requirement type="package" version="3.4">grep</requirement> <requirement type="package" version="3.6.2">matplotlib</requirement> <requirement type="package" version="1.23.5">numpy</requirement> <requirement type="package" version="1.5.2">pandas</requirement> - <requirement type="package" version="1.16.1">samtools</requirement> </requirements> </xml> <xml name="citations">