use get_row_around_cut_asymmetrical for allele plot (#117) (pinellolab#527)

* use get_row_around_cut_asymmetrical for allele plot

* tests branch

* removed commented code

* update to use asym function

* white space

* spaces

* integration tests

* Move test branch back to master

---------

Co-authored-by: mbowcut2 <55161542+mbowcut2@users.noreply.github.com>

Author: Colelyman

.github/workflows/integration_tests.yml

@@ -133,3 +133,18 @@ jobs:
         if: success() || failure()
         run: |
           make aggregate test print
+
+      - name: Run Asym Allele Plot Left
+        if: success() || failure()
+        run: |
+          make asym-left test print
+
+      - name: Run Asym Allele Plot Right
+        if: success() || failure()
+        run: |
+          make asym-right test print
+
+      - name: Run Asym Allele Plot Both
+        if: success() || failure()
+        run: |
+          make asym-both test print

CRISPResso2/CRISPRessoCORE.py

@@ -4581,8 +4581,30 @@ def count_alternate_alleles(sub_base_vectors, ref_name, ref_sequence, ref_total_
                     sgRNA_legend = sgRNA_name + " (" + sgRNA +")"
                 sgRNA_label = CRISPRessoShared.slugify(sgRNA_label)
 
-                plot_half_window = max(1, args.plot_window_size)
-                df_alleles_around_cut=CRISPRessoShared.get_dataframe_around_cut(df_alleles.loc[df_alleles['Reference_Name'] == ref_name], cut_point, plot_half_window)
+                pass_cut_point = False
+
+                #Set left window size
+                if cut_point - args.plot_window_size + 1 >= 0:
+                    plot_half_window_left = args.plot_window_size
+                else:
+                    plot_half_window_left = cut_point + 1
+                    pass_cut_point = True
+                    warn(f'sgRNA {0} is too close to the start of the amplicon to plot the full window. Truncating the window.')
+                
+                #Set right window size
+                if cut_point + args.plot_window_size < ref_len:
+                    plot_half_window_right = args.plot_window_size
+                else:
+                    plot_half_window_right = ref_len - cut_point - 1
+                    pass_cut_point = True
+                    warn(f'sgRNA {0} is too close to the end of the amplicon to plot the full window. Truncating the window.')
+                    
+                df_alleles_around_cut = CRISPRessoShared.get_dataframe_around_cut_asymmetrical(
+                    df_alleles.loc[df_alleles['Reference_Name'] == ref_name],
+                    cut_point,
+                    plot_half_window_left,
+                    plot_half_window_right,
+                )
                 count_total = counts_total[ref_name]
                 if args.allele_plot_pcts_only_for_assigned_reference:
                     df_alleles_around_cut['%AllReads']=df_alleles_around_cut['%Reads']
@@ -4593,7 +4615,7 @@ def count_alternate_alleles(sub_base_vectors, ref_name, ref_sequence, ref_total_
                 df_alleles_around_cut.to_csv(allele_filename, sep='\t', header=True)
                 crispresso2_info['results']['refs'][ref_name]['allele_frequency_files'].append(os.path.basename(allele_filename))
 
-                ref_seq_around_cut=refs[ref_name]['sequence'][cut_point-plot_half_window+1:cut_point+plot_half_window+1]
+                ref_seq_around_cut=refs[ref_name]['sequence'][cut_point-plot_half_window_left+1:cut_point+plot_half_window_right+1]
                 fig_filename_root = _jp('9.'+ref_plot_name+'Alleles_frequency_table_around_'+sgRNA_label)
                 n_good = df_alleles_around_cut[df_alleles_around_cut['%Reads']>=args.min_frequency_alleles_around_cut_to_plot].shape[0]
                 if not args.suppress_plots and n_good > 0:
@@ -4605,11 +4627,12 @@ def count_alternate_alleles(sub_base_vectors, ref_name, ref_sequence, ref_total_
 
                     new_sgRNA_intervals = []
                     #adjust coordinates of sgRNAs
-                    new_sel_cols_start = cut_point - plot_half_window
+                    new_sel_cols_start = cut_point - plot_half_window_left
                     for (int_start, int_end) in refs[ref_name]['sgRNA_intervals']:
                         new_sgRNA_intervals += [(int_start - new_sel_cols_start - 1, int_end - new_sel_cols_start - 1)]
-
-
+                        if int_start <= cut_point and cut_point <= int_end:
+                            new_cut_point = cut_point - new_sel_cols_start - 1
+         
                     prepped_df_alleles, annotations, y_labels, insertion_dict, per_element_annot_kws, is_reference = CRISPRessoPlot.prep_alleles_table(
                         df_to_plot,
                         ref_seq_around_cut,
@@ -4628,6 +4651,7 @@ def count_alternate_alleles(sub_base_vectors, ref_name, ref_sequence, ref_total_
                         'custom_colors': custom_config["colors"],
                         'SAVE_ALSO_PNG': save_png,
                         'plot_cut_point': plot_cut_point,
+                        'cut_point_ind': new_cut_point if pass_cut_point else None,
                         'sgRNA_intervals': new_sgRNA_intervals,
                         'sgRNA_names': sgRNA_names,
                         'sgRNA_mismatches': sgRNA_mismatches,

CRISPResso2/CRISPRessoShared.py

@@ -1309,32 +1309,21 @@ def split_interleaved_fastq(fastq_filename, output_filename_r1, output_filename_
 # allele modification functions
 ######
 
-def get_row_around_cut(row, cut_point, offset):
-    cut_idx = row['ref_positions'].index(cut_point)
-    return row['Aligned_Sequence'][cut_idx - offset + 1:cut_idx + offset + 1], row['Reference_Sequence'][
-                                                                               cut_idx - offset + 1:cut_idx + offset + 1], \
-           row['Read_Status'] == 'UNMODIFIED', row['n_deleted'], row['n_inserted'], row['n_mutated'], row['#Reads'], \
-           row['%Reads']
+def get_row_around_cut_asymmetrical(row,cut_point,plot_left,plot_right):
+    cut_idx=row['ref_positions'].index(cut_point)
+    return row['Aligned_Sequence'][cut_idx-plot_left+1:cut_idx+plot_right+1],row['Reference_Sequence'][cut_idx-plot_left+1:cut_idx+plot_right+1],row['Read_Status']=='UNMODIFIED',row['n_deleted'],row['n_inserted'],row['n_mutated'],row['#Reads'], row['%Reads']
 
 
-def get_dataframe_around_cut(df_alleles, cut_point, offset, collapse_by_sequence=True):
+def get_dataframe_around_cut_asymmetrical(df_alleles, cut_point,plot_left,plot_right,collapse_by_sequence=True):
     if df_alleles.shape[0] == 0:
         return df_alleles
     ref1 = df_alleles['Reference_Sequence'].iloc[0]
-    ref1 = ref1.replace('-', '')
-    if (cut_point + offset + 1 > len(ref1)):
-        raise (BadParameterException(
-            'The plotting window cannot extend past the end of the amplicon. Amplicon length is ' + str(
-                len(ref1)) + ' but plot extends to ' + str(cut_point + offset + 1)))
-
-    df_alleles_around_cut = pd.DataFrame(
-        list(df_alleles.apply(lambda row: get_row_around_cut(row, cut_point, offset), axis=1).values),
-        columns=['Aligned_Sequence', 'Reference_Sequence', 'Unedited', 'n_deleted', 'n_inserted', 'n_mutated', '#Reads',
-                 '%Reads'])
+    ref1 = ref1.replace('-','')
+    
+    df_alleles_around_cut=pd.DataFrame(list(df_alleles.apply(lambda row: get_row_around_cut_asymmetrical(row,cut_point,plot_left,plot_right),axis=1).values),
+                    columns=['Aligned_Sequence','Reference_Sequence','Unedited','n_deleted','n_inserted','n_mutated','#Reads','%Reads'])
 
-    df_alleles_around_cut = df_alleles_around_cut.groupby(
-        ['Aligned_Sequence', 'Reference_Sequence', 'Unedited', 'n_deleted', 'n_inserted',
-         'n_mutated']).sum().reset_index().set_index('Aligned_Sequence')
+    df_alleles_around_cut=df_alleles_around_cut.groupby(['Aligned_Sequence','Reference_Sequence','Unedited','n_deleted','n_inserted','n_mutated']).sum().reset_index().set_index('Aligned_Sequence')
 
     df_alleles_around_cut.sort_values(by=['#Reads', 'Aligned_Sequence', 'Reference_Sequence'], inplace=True, ascending=[False, True, True])
     df_alleles_around_cut['Unedited'] = df_alleles_around_cut['Unedited'] > 0
@@ -1631,13 +1620,11 @@ def get_amplicon_info_for_guides(ref_seq, guides, guide_mismatches, guide_names,
     if this_sgRNA_cut_points and plot_window_size > 0:
         for cut_p in this_sgRNA_cut_points:
             if cut_p - window_around_cut + 1 < 0:
-                raise BadParameterException(
-                    'Offset around cut would extend to the left of the amplicon. Please decrease plot_window_size parameter. Cut point: ' + str(
-                        cut_p) + ' window: ' + str(window_around_cut) + ' reference: ' + str(ref_seq_length))
+                logging.warning('Offset around cut would extend to the left of the amplicon. Window will be truncated.')
+                window_around_cut = cut_p + 1
             if cut_p + window_around_cut > ref_seq_length - 1:
-                raise BadParameterException(
-                    'Offset around cut would be greater than reference sequence length. Please decrease plot_window_size parameter. Cut point: ' + str(
-                        cut_p) + ' window: ' + str(window_around_cut) + ' reference: ' + str(ref_seq_length))
+                logging.warning('Offset around cut would be greater than reference sequence length. Window will be truncated.')
+                window_around_cut = ref_seq_length - cut_p - 1
             st = max(0, cut_p - window_around_cut + 1)
             en = min(ref_seq_length - 1, cut_p + window_around_cut + 1)
             this_sgRNA_plot_idxs.append(sorted(list(range(st, en))))

scripts/plotAmbiguous.py

@@ -104,7 +104,7 @@ def plot_ambiguous_alleles_tables_from_folder(crispresso_output_folder,fig_filen
             ref_seq_around_cut=refs[ref_name]['sequence'][cut_point-plot_half_window+1:cut_point+plot_half_window+1]
 
             ambiguous_ref_name = "AMBIGUOUS_"+ref_name
-            df_alleles_around_cut=CRISPRessoShared.get_dataframe_around_cut(df_alleles.loc[df_alleles['Reference_Name'] == ambiguous_ref_name],cut_point,plot_half_window)
+            df_alleles_around_cut=CRISPRessoShared.get_dataframe_around_cut_asymmetrical(df_alleles.loc[df_alleles['Reference_Name'] == ambiguous_ref_name],cut_point,plot_half_window, plot_half_window)
             this_ambig_allele_count = len(df_alleles_around_cut.index)
             if this_ambig_allele_count < 1:
                 print('No ambiguous reads found for ' + ref_name)