Mckay/plot amino acids

CRISPResso2/BLOSUM62

@@ -0,0 +1,27 @@
+#
+# This image is the BLOSUM substitution matrix, which indicates the substitution score for replacing one amino acid with any other amino acid. 
+#
+#
+# Lowest score = -4, Highest score = 11
+#
+ A R N D C Q E G H I L K M F P S T W Y V
+A 4 -1 -2 -2 0 -1 -1 0 -2 -1 -1 -1 -1 -2 -1 1 0 -3 -2 0
+R -1 5 0 -2 -3 1 0 -2 0 -3 -2 2 -1 -3 -2 -1 -1 -3 -2 -3
+N -2 0 6 1 -3 0 0 0 1 -3 -3 0 -2 -3 -2 1 0 -4 -2 -3
+D -2 -2 1 6 -3 0 2 -1 -1 -3 -4 -1 -3 -3 -1 0 -1 -4 -3 -3
+C 0 -3 -3 -3 9 -3 -4 -3 -3 -1 -1 -3 -1 -2 -3 -1 -1 -2 -2 -1
+Q -1 1 0 0 -3 5 2 -2 0 -3 -2 1 0 -3 -1 0 -1 -2 -1 -2
+E -1 0 0 2 -4 2 5 -2 0 -3 -3 1 -2 -3 -1 0 -1 -3 -2 -2
+G 0 -2 0 -1 -3 -2 -2 6 -2 -4 -4 -2 -3 -3 -2 0 -2 -2 -3 -3
+H -2 0 1 -1 -3 0 0 -2 8 -3 -3 -1 -2 -1 -2 -1 -2 -2 2 -3
+I -1 -3 -3 -3 -1 -3 -3 -4 -3 4 2 -3 1 0 -3 -2 -1 -3 -1 3
+L -1 -2 -3 -4 -1 -2 -3 -4 -3 2 4 -2 2 0 -3 -2 -1 -2 -1 1
+K -1 2 0 -1 -3 1 1 -2 -1 -3 -2 5 -1 -3 -1 0 -1 -3 -2 -2
+M -1 -1 -2 -3 -1 0 -2 -3 -2 1 2 -1 5 0 -2 -1 -1 -1 -1 1
+F -2 -3 -3 -3 -2 -3 -3 -3 -1 0 0 -3 0 6 -4 -2 -2 1 3 -1
+P -1 -2 -2 -1 -3 -1 -1 -2 -2 -3 -3 -1 -2 -4 7 -1 -1 -4 -3 -2
+S 1 -1 1 0 -1 0 0 0 -1 -2 -2 0 -1 -2 -1 4 1 -3 -2 -2
+T 0 -1 0 -1 -1 -1 -1 -2 -2 -1 -1 -1 -1 -2 -1 1 5 -2 -2 0
+W -3 -3 -4 -4 -2 -2 -3 -2 -2 -3 -2 -3 -1 1 -4 -3 -2 11 2 -3
+Y -2 -2 -2 -3 -2 -1 -2 -3 2 -1 -1 -2 -1 3 -3 -2 -2 2 7 -1
+V 0 -3 -3 -3 -1 -2 -2 -3 -3 3 1 -2 1 -1 -2 -2 0 -3 -1 4

CRISPResso2/CRISPRessoAggregateCORE.py

@@ -71,6 +71,7 @@ def main():
 
         parser.add_argument('--debug', help='Show debug messages', action='store_true')
         parser.add_argument('-v', '--verbosity', type=int, help='Verbosity level of output to the console (1-4), 4 is the most verbose', default=3)
+        parser.add_argument('--halt_on_plot_fail', action="store_true", help="Halt execution if a plot fails to generate")
 
         # CRISPRessoPro params
         parser.add_argument('--use_matplotlib', action='store_true',
@@ -131,6 +132,7 @@ def main():
             num_processes=n_processes,
             process_pool=process_pool,
             process_futures=process_futures,
+            halt_on_plot_fail=args.halt_on_plot_fail,
         )
 
         #glob returns paths including the original prefix

CRISPResso2/CRISPRessoBatchCORE.py

@@ -400,6 +400,7 @@ def main():
             num_processes=n_processes_for_batch,
             process_futures=process_futures,
             process_pool=process_pool,
+            halt_on_plot_fail=args.halt_on_plot_fail,
         )
 
         window_nuc_pct_quilt_plot_names = []

CRISPResso2/CRISPRessoCORE.py

@@ -1606,6 +1606,9 @@ def rreplace(s, old, new):
                     raise CRISPRessoShared.NTException('The coding sequence contains bad characters:%s' % ' '.join(wrong_nt))
 
                 coding_seqs.append(exon_seq)
+
+        if len(coding_seqs) > 0:
+            crispresso2_info['running_info']['coding_seqs'] = coding_seqs
 
         ####SET REFERENCES TO COMPARE###
         ref_names = [] #ordered list of names
@@ -3746,6 +3749,7 @@ def count_alternate_alleles(sub_base_vectors, ref_name, ref_sequence, ref_total_
             num_processes=n_processes,
             process_pool=process_pool,
             process_futures=process_futures,
+            halt_on_plot_fail=args.halt_on_plot_fail,
         )
         ###############################################################################################################################################
         ### FIGURE 1: Alignment
@@ -4551,6 +4555,11 @@ def count_alternate_alleles(sub_base_vectors, ref_name, ref_sequence, ref_total_
             crispresso2_info['results']['refs'][ref_name]['plot_9_roots'] = []
             crispresso2_info['results']['refs'][ref_name]['plot_9_captions'] = []
             crispresso2_info['results']['refs'][ref_name]['plot_9_datas'] = []
+
+            crispresso2_info['results']['refs'][ref_name]['plot_9a_roots'] = []
+            crispresso2_info['results']['refs'][ref_name]['plot_9a_captions'] = []
+            crispresso2_info['results']['refs'][ref_name]['plot_9a_datas'] = []
+
             crispresso2_info['results']['refs'][ref_name]['allele_frequency_files'] = []
 
             crispresso2_info['results']['refs'][ref_name]['plot_10d_roots'] = []
@@ -4774,6 +4783,44 @@ def count_alternate_alleles(sub_base_vectors, ref_name, ref_sequence, ref_total_
                         crispresso2_info['results']['refs'][ref_name]['plot_10g_captions'].append("Figure 10g: Non-reference base counts. For target nucleotides in the plotting window, this plot shows the number of non-reference (non-" + args.conversion_nuc_from + ") bases. The number of each target base is annotated on the reference sequence at the bottom of the plot.")
                         crispresso2_info['results']['refs'][ref_name]['plot_10g_datas'].append([('Nucleotide frequencies at ' + args.conversion_nuc_from +'s', os.path.basename(quant_window_sel_nuc_freq_filename))])
 
+
+            if refs[ref_name]['contains_coding_seq']:
+                for i, coding_seq in enumerate(coding_seqs):
+                    fig_filename_root = _jp('9a.'+ref_plot_name+'amino_acid_table_around_'+coding_seq)
+                    coding_seq_amino_acids = CRISPRessoShared.get_amino_acids_from_nucs(coding_seq)
+                    amino_acid_cut_point = (cut_point - refs[ref_name]['exon_positions'][0] + 1)// 3  
+                    df_to_plot = CRISPRessoShared.get_amino_acid_dataframe(
+                        df_alleles.loc[df_alleles['Reference_Name'] == ref_name],
+                        refs[ref_name]['exon_intervals'][i][0], 
+                        len(coding_seq_amino_acids),
+                        os.path.join(_ROOT, "BLOSUM62"),
+                        amino_acid_cut_point)
+
+                    plot_9a_input = {
+                        'reference_seq': coding_seq_amino_acids,
+                        'df_alleles': df_to_plot,
+                        'fig_filename_root': fig_filename_root,
+                        'custom_colors': custom_config["colors"],
+                        'MIN_FREQUENCY': args.min_frequency_alleles_around_cut_to_plot,
+                        'MAX_N_ROWS': args.max_rows_alleles_around_cut_to_plot,
+                        'SAVE_ALSO_PNG': save_png,
+                        'plot_cut_point': plot_cut_point,
+                        'sgRNA_intervals': new_sgRNA_intervals,
+                        'sgRNA_names': sgRNA_names,
+                        'sgRNA_mismatches': sgRNA_mismatches,
+                        'annotate_wildtype_allele': args.annotate_wildtype_allele,
+                        'cut_point': amino_acid_cut_point,
+                    }
+
+                    amino_acid_filename = _jp(ref_plot_name+'amino_acid_table_for_'+coding_seq+'.txt')
+                    df_to_plot.to_csv(amino_acid_filename, sep='\t', header=True, index=True)
+
+                    debug('Plotting amino acids for {0}'.format(ref_name))
+                    plot(CRISPRessoPlot.plot_amino_acid_table, plot_9a_input)
+                    crispresso2_info['results']['refs'][ref_name]['plot_9a_roots'].append(os.path.basename(fig_filename_root))
+                    crispresso2_info['results']['refs'][ref_name]['plot_9a_captions'].append(
+                        "Figure 9a: Visualization of the distribution of identified amino acids based on the coding sequence (" + coding_seq+"). The vertical dashed line indicates the predicted cleavage site.")
+                    crispresso2_info['results']['refs'][ref_name]['plot_9a_datas'].append([('Amino Acid table', os.path.basename(amino_acid_filename))])
             info('Done!')
 
             #END GUIDE SPECIFIC PLOTS
@@ -5174,6 +5221,10 @@ def get_scaffold_len(row, scaffold_start_loc, scaffold_seq):
         print_stacktrace_if_debug()
         error('Filtering error, please check your input.\n\nERROR: %s' % e)
         sys.exit(13)
+    except CRISPRessoShared.PlotException as e:
+        print_stacktrace_if_debug()
+        error(e)
+        sys.exit(14)
     except Exception as e:
         print_stacktrace_if_debug()
         error('Unexpected error, please check your input.\n\nERROR: %s' % e)

CRISPResso2/CRISPRessoMultiProcessing.py

@@ -16,6 +16,9 @@
 import pandas as pd
 import traceback
 
+from CRISPResso2.CRISPRessoShared import PlotException
+
+
 def get_max_processes():
     return mp.cpu_count()
 
@@ -284,7 +287,7 @@ def run_parallel_commands(commands_arr, n_processes=1, descriptor='CRISPResso2',
     pool.join()
 
 
-def run_plot(plot_func, plot_args, num_processes, process_futures, process_pool):
+def run_plot(plot_func, plot_args, num_processes, process_futures, process_pool, halt_on_plot_fail):
     """Run a plot in parallel if num_processes > 1, otherwise in serial.
 
     Parameters
@@ -299,6 +302,8 @@ def run_plot(plot_func, plot_args, num_processes, process_futures, process_pool)
         The list of futures that submitting the parallel job will return.
     process_pool: ProcessPoolExecutor or ThreadPoolExecutor
         The pool to submit the job to.
+    halt_on_plot_fail: bool
+        If True, an exception will be raised if the plot fails
 
     Returns
     -------
@@ -311,5 +316,8 @@ def run_plot(plot_func, plot_args, num_processes, process_futures, process_pool)
         else:
             plot_func(**plot_args)
     except Exception as e:
+        if halt_on_plot_fail:
+            logger.critical(f"Plot error, halting execution \n")
+            raise PlotException(f'There was an error generating plot {plot_func.__name__}.')
         logger.warn(f"Plot error {e}, skipping plot \n")
         logger.debug(traceback.format_exc())