diff --git a/statistics/generate_figures.py b/statistics/generate_figures.py
index 056b77e..2283a72 100644
--- a/statistics/generate_figures.py
+++ b/statistics/generate_figures.py
@@ -330,7 +330,7 @@ def create_lineplot_21_40_persex(df, path_out, show_cv=False):
     print('Figure saved: ' + path_filename)
 
 
-def create_lineplot(df, hue, path_out, show_cv=False):
+def create_lineplot(df, hue, path_out, show_cv=False, show_n=True):
     """
     Create lineplot for individual metrics per vertebral levels.
     Note: we are ploting slices not levels to avoid averaging across levels.
@@ -339,6 +339,8 @@ def create_lineplot(df, hue, path_out, show_cv=False):
         hue (str): column name of the dataframe to use for grouping; if None, no grouping is applied
         path_out (str): path to output directory
         show_cv (bool): if True, include coefficient of variation for each vertebral level to the plot
+        show_n (bool): if True, print the per-level subject/session count under each vertebral label.
+                       Set to False for wide coverage (e.g. C1-T12) where the counts would overlap.
     """
 
     mpl.rcParams['font.family'] = 'Arial'
@@ -443,11 +445,11 @@ def create_lineplot(df, hue, path_out, show_cv=False):
             if show_cv:
                 cv = compute_cv(df[(df['VertLevel'] == vert[x])], metric)
             n = n_per_level.get(vert[x], 0)
-            n_str = f"n={n}"
+            n_str = f"\nn={n}" if show_n else ""
             # Deal with T1 label (C8 -> T1)
             if vert[x] > 7:
                 level = 'T' + str(vert[x] - 7)
-                axs[index].text(df.loc[ind_vert_mid[idx], 'Slice (I->S)'], ymin, f'{level}\n{n_str}',
+                axs[index].text(df.loc[ind_vert_mid[idx], 'Slice (I->S)'], ymin, f'{level}{n_str}',
                                 horizontalalignment='center', verticalalignment='bottom', color='black',
                                 fontsize=TICKS_FONT_SIZE-2)
                 # Show CV
@@ -457,7 +459,7 @@ def create_lineplot(df, hue, path_out, show_cv=False):
                                     color='black')
             else:
                 level = 'C' + str(vert[x])
-                axs[index].text(df.loc[ind_vert_mid[idx], 'Slice (I->S)'], ymin, f'{level}\n{n_str}',
+                axs[index].text(df.loc[ind_vert_mid[idx], 'Slice (I->S)'], ymin, f'{level}{n_str}',
                                 horizontalalignment='center', verticalalignment='bottom', color='black',
                                 fontsize=TICKS_FONT_SIZE-2)
                 # Show CV
@@ -1293,19 +1295,43 @@ def compute_age_stats(df_participants):
           f'Max: {age_max}\n')
 
 
+def _discover_pam50_csvs(path_HC):
+    """
+    Find per-subject PAM50 CSV files under ``path_HC``.
+
+    Supports two filename conventions:
+      - flat layout used by this repo, e.g. ``whole-spine/sub-amuAL_T2w_PAM50.csv``
+      - nested BIDS-like, e.g. ``<sub>/<contrast>/sub-XXX_..._space-PAM50_desc-sct-morphometrics_stat.csv``
+
+    Returns a list of (basename, absolute_path) tuples.
+    """
+    # Flat layout first (existing behavior, fast path)
+    flat = [f for f in os.listdir(path_HC) if 'PAM50.csv' in f]
+    if flat:
+        return [(f, os.path.join(path_HC, f)) for f in flat]
+
+    # Fall back to recursive search for BIDS-like layout
+    found = []
+    for root, _, files in os.walk(path_HC):
+        for f in files:
+            if 'space-PAM50' in f and f.endswith('morphometrics_stat.csv'):
+                found.append((f, os.path.join(root, f)))
+    return found
+
+
 def read_csv_files(path_HC, participant_file=None, dataset_name=None):
     # Initialize pandas dataframe where data across all subjects will be stored
-    csv_files = [f for f in os.listdir(path_HC) if 'PAM50.csv' in f]
+    csv_files = _discover_pam50_csvs(path_HC)
     n_total = len(csv_files)
     print(f'Reading {path_HC} ({n_total} files)')
     dfs = []
-    for i, file in enumerate(csv_files, 1):
+    for i, (fname, fpath) in enumerate(csv_files, 1):
         # Read csv file as pandas dataframe for given subject
-        df_subject = pd.read_csv(os.path.join(path_HC, file), dtype=METRICS_DTYPE)
+        df_subject = pd.read_csv(fpath, dtype=METRICS_DTYPE)
         # Compute AP/RL ratio as MEAN(diameter_AP) / MEAN(diameter_RL)
         df_subject['MEAN(compression_ratio)'] = df_subject['MEAN(diameter_AP)'] / df_subject['MEAN(diameter_RL)']
         # Track source CSV filename to reliably extract participant_id
-        df_subject['source_file'] = file
+        df_subject['source_file'] = fname
         dfs.append(df_subject)
         if i % 100 == 0 or i == n_total:
             print(f'  {i}/{n_total} files read ({100 * i // n_total}%)', end='\r', flush=True)