changeset 7:389c98d344ce draft

Uploaded
author greg
date Fri, 03 Mar 2023 22:03:09 +0000 (22 months ago)
parents caf554e039b2
children 7fe8ea50a81d
files draw_amr_matrix.py draw_amr_matrix.xml
diffstat 2 files changed, 40 insertions(+), 36 deletions(-) [+]
line wrap: on
line diff
--- a/draw_amr_matrix.py	Wed Feb 22 21:23:25 2023 +0000
+++ b/draw_amr_matrix.py	Fri Mar 03 22:03:09 2023 +0000
@@ -58,7 +58,7 @@
     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):
+def draw_amr_matrix(amr_feature_hits_files, amr_deletions_file, varscan_vcf_file, amr_mutation_regions_bed_file, amr_gene_drug_file, reference, reference_size, mutation_regions_dir, amr_matrix_png_dir):
     ofh = open('process_log', 'w')
 
     # Read amr_feature_hits_files.
@@ -98,21 +98,23 @@
                     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("\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:
+        ofh.write("\nvarscan_vcf_file si None: %s\n" % str(varscan_vcf_file == 'None'))
+        if varscan_vcf_file not in [None, 'None'] and os.path.getsize(varscan_vcf_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 amr_mutation_regions_bed_file is not None:
+                mutation_regions = pandas.read_csv(amr_mutation_regions_bed_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:
+                    """
+                    TODO: move this coder to the pima_report tool...
                     # 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)))
+                            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]))
@@ -124,26 +126,28 @@
                             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')
+                        region_bed = 'region_%s.bed' % region_i
                         mutation_regions.loc[[region_i], ].to_csv(path_or_buf=region_bed, sep='\t', header=False, index=False)
+                        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',
-                                        amr_mutations_file,
-                                        ' | awk \'BEGIN{getline < "' + amr_mutation_regions_file + '";printf $0"\\t";',
-                                        'getline < "' + amr_mutations_file + '"; getline < "' + amr_mutations_file + '";print $0}{print}\'',
-                                        '1>' + region_mutations_output_file])
+                                        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)
+                        """
+                        TODO: move this coder to the pima_report tool...
                         try:
-                            region_mutations = pandas.read_csv(region_mutations_output_file, sep='\t', header=0, index_col=False)
+                            region_mutations = pandas.read_csv(region_mutations_tsv, sep='\t', header=0, index_col=False)
                         except Exception:
                             region_mutations = pandas.DataFrame()
                         if region_mutations.shape[0] == 0:
@@ -156,6 +160,7 @@
                         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.
@@ -183,7 +188,7 @@
             amr_matrix = pandas.DataFrame(0, index=present_genes, columns=present_drugs)
             for hit_idx, hit in amr_to_draw.iterrows():
                 amr_matrix.loc[hit[0], hit[1]] = 1
-            amr_matrix_png = os.path.join(output_dir, 'amr_matrix.png')
+            amr_matrix_png = os.path.join(amr_matrix_png_dir, 'amr_matrix.png')
             int_matrix = amr_matrix[amr_matrix.columns].astype(int)
             figure, axis = pyplot.subplots()
             heatmap = axis.pcolor(int_matrix, cmap=pyplot.cm.Blues, linewidth=0)
@@ -206,13 +211,12 @@
 
     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('--varscan_vcf_file', action='store', dest='varscan_vcf_file', default=None, help='Varscan VCF file produced by the call_amr_mutations tool')
+    parser.add_argument('--amr_mutation_regions_bed_file', action='store', dest='amr_mutation_regions_bed_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')
+    parser.add_argument('--mutation_regions_dir', action='store', dest='mutation_regions_dir', help='Directory for mutation regions TSV files produced by this tool')
+    parser.add_argument('--amr_matrix_png_dir', action='store', dest='amr_matrix_png_dir', help='Directory for PNG files produced by this tool')
 
     args = parser.parse_args()
 
@@ -231,4 +235,4 @@
     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)
+    draw_amr_matrix(amr_feature_hits_files, args.amr_deletions_file, args.varscan_vcf_file, args.amr_mutation_regions_bed_file, args.amr_gene_drug_file, reference, reference_size, args.mutation_regions_dir, args.amr_matrix_png_dir)
--- a/draw_amr_matrix.xml	Wed Feb 22 21:23:25 2023 +0000
+++ b/draw_amr_matrix.xml	Fri Mar 03 22:03:09 2023 +0000
@@ -8,8 +8,8 @@
 #import re
 
 mkdir amr_feature_hits_dir &&
-mkdir mutations_dir &&
-mkdir output_dir &&
+mkdir mutation_regions_dir &&
+mkdir amr_matrix_png_dir &&
 
 #if $reference_source.reference_source_selector == 'history':
     ln -f -s '$reference_source.ref_file' reference.fa &&
@@ -28,17 +28,16 @@
 #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'
+#if str($varscan_vcf_file) != 'None':
+    --varscan_vcf_file '$varscan_vcf_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'
+#if str($amr_mutation_regions_bed_file) != 'None':
+    --amr_mutation_regions_bed_file '$amr_mutation_regions_bed_file'
+    --mutation_regions_dir 'mutation_regions_dir'
 #end if
 --amr_gene_drug_file '$amr_gene_drug_file'
 --reference_genome reference.fa
---mutations_dir 'mutations_dir'
---output_dir 'output_dir'
+--amr_matrix_png_dir 'amr_matrix_png_dir'
 #if str($output_process_log) == 'yes':
     && mv 'process_log' '$process_log'
 #end if
@@ -63,8 +62,8 @@
         </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="--varscan_vcf_file" type="data" format="vcf" optional="true" label="Varscan VCF file" help="Optional, leave blank to ignore"/>
+        <param argument="--amr_mutation_regions_bed_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>
@@ -75,11 +74,12 @@
         <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="mutation_regions_tsv" type="list" format="tsv" label="${tool.name} on ${on_string} (mutation regions)">
+            <filter>amr_mutation_regions_bed_file not in [None, 'None']</filter>
+            <discover_datasets pattern="(?P&lt;designation&gt;.+)\.(?P&lt;ext&gt;tsv)" directory="mutation_regions_dir"/>
+        </collection>
         <collection name="amr_matrix_png" type="list" format="png">
-            <discover_datasets pattern="(?P&lt;designation&gt;.+)\.(?P&lt;ext&gt;png)" directory="output_dir"/>
+            <discover_datasets pattern="(?P&lt;designation&gt;.+)\.(?P&lt;ext&gt;png)" directory="amr_matrix_png_dir"/>
         </collection>
     </outputs>
     <tests>