bugfix: handle the case where mv tags are not added to bam file

Closes #9

Author: krish8484

bin/plot_signal.py

@@ -31,6 +31,10 @@ def main():
     args = parser.parse_args()
 
     # Load the data generated by pod5_to_df.py
+    if not os.path.exists(args.csv):
+        print(f"Error: CSV file {args.csv} not found.")
+        return
+        
     raw_signal_df = pd.read_csv(args.csv)
 
     # Plotting configuration
@@ -39,37 +43,44 @@ def main():
     fig, ax = plt.subplots(1, 1, figsize=(args.figwidth, args.figheight))
     x_feature, y_feature = 'time', 'signal'
 
-    # 1. Unannotated part (ann == -2)
-    df_neg2 = raw_signal_df[raw_signal_df['ann'] == -2]
-    if not df_neg2.empty:
-        sns.scatterplot(data=df_neg2, x=x_feature, y=y_feature, color='blue', 
-                        label='unannotated part', s=50, zorder=4, ax=ax)
+    
+    # Check if 'ann' column exists (requires Dorado move tags)
+    if 'ann' in raw_signal_df.columns:
+        # 1. Unannotated part (ann == -2)
+        df_neg2 = raw_signal_df[raw_signal_df['ann'] == -2]
+        if not df_neg2.empty:
+            sns.scatterplot(data=df_neg2, x=x_feature, y=y_feature, color='blue', 
+                            label='unannotated part', s=50, zorder=4, ax=ax)
 
-    # 2. Trimmed primer/adapter (ann == -1)
-    df_neg1 = raw_signal_df[raw_signal_df['ann'] == -1]
-    if not df_neg1.empty:
-        sns.lineplot(data=df_neg1, x=x_feature, y=y_feature, color='green', 
-                     label='trimmed primer and adapter', zorder=2, ax=ax)
+        # 2. Trimmed primer/adapter (ann == -1)
+        df_neg1 = raw_signal_df[raw_signal_df['ann'] == -1]
+        if not df_neg1.empty:
+            sns.lineplot(data=df_neg1, x=x_feature, y=y_feature, color='green', 
+                         label='trimmed primer and adapter', zorder=2, ax=ax)
 
-    # 3. Basecalled region (ann is 0 or 1)
-    df_base = raw_signal_df[raw_signal_df['ann'].isin([0, 1])]
-    if not df_base.empty:
-        sns.lineplot(data=df_base, x=x_feature, y=y_feature, color='orange', 
-                     label='basecalled region', zorder=3, ax=ax)
+        # 3. Basecalled region (ann is 0 or 1)
+        df_base = raw_signal_df[raw_signal_df['ann'].isin([0, 1])]
+        if not df_base.empty:
+            sns.lineplot(data=df_base, x=x_feature, y=y_feature, color='orange', 
+                         label='basecalled region', zorder=3, ax=ax)
 
-    # 4. Poly-A Tail Region (polyA > -1)
-    # We highlight the region identified by the 'pa' tag anchors
-    df_polya = raw_signal_df[raw_signal_df['polyA'] > -1]
-    if not df_polya.empty:
-        # We use a distinct color (magenta) to show the estimated tail
-        sns.lineplot(data=df_polya, x=x_feature, y=y_feature, color='magenta', 
-                     label='polyA-tail region', zorder=5, linewidth=2, ax=ax)
+        # 5. Samples that emit bases (ann == 1) - Red Circles
+        df_emit = raw_signal_df[raw_signal_df['ann'] == 1]
+        if not df_emit.empty:
+            sns.scatterplot(data=df_emit, x=x_feature, y=y_feature, color='red', 
+                            label='samples that emit bases', s=50, fc="none", ec='red', zorder=6, ax=ax)
+    else:
+        # Fallback: Plot raw signal if 'ann' is missing (Move tags were likely missing in BAM)
+        print(f"Warning: 'ann' column missing in {args.csv}. Plotting raw signal only.")
+        sns.lineplot(data=raw_signal_df, x=x_feature, y=y_feature, color='grey', 
+                     alpha=0.5, label='raw signal (unannotated)', ax=ax)
 
-    # 5. Samples that emit bases (ann == 1) - Red Circles
-    df_emit = raw_signal_df[raw_signal_df['ann'] == 1]
-    if not df_emit.empty:
-        sns.scatterplot(data=df_emit, x=x_feature, y=y_feature, color='red', 
-                        label='samples that emit bases', s=50, fc="none", ec='red', zorder=6, ax=ax)
+    # 4. Poly-A Tail Region (independent check for 'polyA' column)
+    if 'polyA' in raw_signal_df.columns:
+        df_polya = raw_signal_df[raw_signal_df['polyA'] > -1]
+        if not df_polya.empty:
+            sns.lineplot(data=df_polya, x=x_feature, y=y_feature, color='magenta', 
+                         label='polyA-tail region', zorder=5, linewidth=2, ax=ax)
 
     ax.set(title=args.title)
 
@@ -78,6 +89,7 @@ def main():
     plt.tight_layout()
     plt.savefig(args.output, dpi=300)
     plt.close()
+    print(f"Successfully saved plot to {args.output}")
 
 if __name__ == "__main__":
     main()