<jupyter_start><jupyter_text>Lesion ratios<jupyter_code>def lesions_in_df(df, lesion_masks, brain_masks, only_with_lesions=False, list_not_sum=False):
    """
    Return all lesion masks in a single DataFrame as list or sum.
    """
    
    # TODO its broken
    if not list_not_sum:
        lesion_sum = np.zeros(shape=shape)
        lesion_sum_hc = np.zeros(shape=shape)
        lesion_sum_ms = np.zeros(shape=shape)
    else:
        lesion_sum = []
        lesion_sum_hc = []
        lesion_sum_ms = []

    for idx, row in df.iterrows():
        lm, bm = lesion_masks[idx], brain_masks[idx]
        # catch scans without lesion masks i.e. contain random info
        if np.max(lm) != 1 and np.max(lm) != 0 and np.min(lm) != 0:
            print("No lesion mask in DataFrame. Label: {}".format(row["label"]))
        else:
            if not list_not_sum:
                lesion_sum += np.sum(lm)
            else:
                lesion_sum.append(lm)
            if row["label"] == "MS":
                if not list_not_sum:
                    lesion_sum_ms += np.sum(lm)
                else:
                    lesion_sum_ms.append(lm)
            else:
                if not list_not_sum:
                    lesion_sum_hc += np.sum(lm)
                else:
                    lesion_sum_hc.append(lm)
    return lesion_sum, lesion_sum_hc, lesion_sum_ms
# calculate ratio of all subjects WITH LESIONS
lesion_ratio = []
lesion_sum = []
lesion_sum_hc = []
lesion_sum_ms = []

for idx, row in train_df.iterrows():
    lm, bm = lesion_masks_train[idx], brain_masks_train[idx]
    if np.sum(lm) > 0:
        # catch scans without lesion masks i.e. contain random info
        if np.max(lm) != 1 and np.max(lm) != 0 and np.min(lm) != 0:
            print("No lesion mask")
        else:
            ratio = np.sum(lm)/np.sum(bm)
            lesion_ratio.append(ratio)
            lesion_sum.append(np.sum(lm))
            if row["label"] == "MS":
                lesion_sum_ms.append(np.sum(lm))
            else:
                lesion_sum_hc.append(np.sum(lm))

for idx, row in holdout_df.iterrows():
    lm, bm = lesion_masks_holdout[idx], brain_masks_holdout[idx]
    if np.sum(lm) > 0:
        # catch scans without lesion masks i.e. contain random info
        if np.max(lm) != 1 and np.max(lm) != 0 and np.min(lm) != 0:
            print("No lesion mask")
        else:
            ratio = np.sum(lm)/np.sum(bm)
            lesion_ratio.append(ratio)
            lesion_sum.append(np.sum(lm))
            if row["label"] == "MS":
                lesion_sum_ms.append(np.sum(lm))
            else:
                lesion_sum_hc.append(np.sum(lm))
    
print(np.mean(lesion_ratio))
np.max(lesion_sum_ms)
np.max(lesion_ratio)
np.min(lesion_ratio)
plt.hist(lesion_ratio, bins=60)
plt.show()
print("Mean lesion ratio {:.4f} %".format(np.mean(lesion_ratio) * 100))
print("Std lesion ratio {:.4f} %".format(np.std(lesion_ratio) * 100))
print("Max lesion ratio {:.4f} %".format(np.max(lesion_ratio) * 100))
print("Min lesion ratio {:.4f} %".format(np.min(lesion_ratio) * 100))
# calculate sums of ALL subjects
lesion_sum = []
lesion_sum_hc = []
lesion_sum_ms = []

for idx, row in train_df.iterrows():
    lm, bm = lesion_masks_train[idx], brain_masks_train[idx]
    # catch scans without lesion masks i.e. contain random info
    if np.max(lm) != 1 and np.max(lm) != 0 and np.min(lm) != 0:
        print("No lesion mask in training. abel: {}".format(row["label"]))
    else:
        lesion_sum.append(np.sum(lm))
        if row["label"] == "MS":
            lesion_sum_ms.append(np.sum(lm))
        else:
            lesion_sum_hc.append(np.sum(lm))

for idx, row in holdout_df.iterrows():
    lm, bm = lesion_masks_holdout[idx], brain_masks_holdout[idx]
    # catch scans without lesion masks i.e. contain random info
    if np.max(lm) != 1 and np.max(lm) != 0 and np.min(lm) != 0:
        print("No lesion mask in holdout")
    else:
        lesion_sum.append(np.sum(lm))
        if row["label"] == "MS":
            lesion_sum_ms.append(np.sum(lm))
        else:
            lesion_sum_hc.append(np.sum(lm))
print("Mean lesion sum {:.2f}".format(np.mean(lesion_sum)))
print("Std lesion sum {:.2f}".format(np.std(lesion_sum)))
print("Max lesion sum {:.2f}".format(np.max(lesion_sum)))
print("Min lesion sum {:.2f}".format(np.min(lesion_sum)))
print("Mean lesion sum HC {:.2f}".format(np.mean(lesion_sum_hc)))
print("Std lesion sum HC {:.2f}".format(np.std(lesion_sum_hc)))
print("Max lesion sum HC {:.2f}".format(np.max(lesion_sum_hc)))
print("Min lesion sum HC {:.2f}".format(np.min(lesion_sum_hc)))
print("Mean lesion sum MS {:.2f}".format(np.mean(lesion_sum_ms)))
print("Std lesion sum MS {:.2f}".format(np.std(lesion_sum_ms)))
print("Max lesion sum MS {:.2f}".format(np.max(lesion_sum_ms)))
print("Min lesion sum MS {:.2f}".format(np.min(lesion_sum_ms)))
len(lesion_sum_ms)
len(lesion_sum_hc)<jupyter_output><empty_output><jupyter_text>Lesion area distribution<jupyter_code>from utils import load_nifti
# load original sized version!
wm_atlas_path = "/analysis/ritter/misc/masks/FSL_atlases/JHU/JHU-ICBM-labels-1mm.nii.gz"
wm_atlas = load_nifti(wm_atlas_path)
plt.imshow(wm_atlas[:,:,96])
plt.show()
shape = (182, 218, 182)
def get_lesion_sum(list_not_sum=False, include_train=True, include_holdout=True, include_brain_masks=False):
    if not list_not_sum:
        # calculate sums of all subjects over both datasets!
        lesion_sum = np.zeros(shape=shape)
        lesion_sum_hc = np.zeros(shape=shape)
        lesion_sum_ms = np.zeros(shape=shape)
    else:
        lesion_sum = []
        lesion_sum_hc = []
        lesion_sum_ms = []
        
    if include_brain_masks:
        brain_avg = np.zeros(shape=shape)
        brain_avg_hc = np.zeros(shape=shape)
        brain_avg_ms = np.zeros(shape=shape)
        

    if include_train:
        for idx, row in train_df.iterrows():
            lm, bm = lesion_masks_train[idx], brain_masks_train[idx]
            # catch scans without lesion masks i.e. contain random info
            if np.max(lm) != 1 and np.max(lm) != 0 and np.min(lm) != 0:
                print("No lesion mask in training. Label: {}".format(row["label"]))
            else:
                if not list_not_sum:
                    lesion_sum += lm
                else:
                    lesion_sum.append(lm)
                if include_brain_masks:
                    brain_avg += bm
                if row["label"] == "MS":
                    if not list_not_sum:
                        lesion_sum_ms += lm
                    else:
                        lesion_sum_ms.append(lm)
                    if include_brain_masks:
                        brain_avg_ms += bm
                else:
                    if not list_not_sum:
                        lesion_sum_hc += lm
                    else:
                        lesion_sum_hc.append(lm)
                    if include_brain_masks:
                        brain_avg_hc += bm

    if include_holdout:
        for idx, row in holdout_df.iterrows():
            lm, bm = lesion_masks_holdout[idx], brain_masks_holdout[idx]
            # catch scans without lesion masks i.e. contain random info
            if np.max(lm) != 1 and np.max(lm) != 0 and np.min(lm) != 0:
                print("No lesion mask in holdout")
            else:
                if not list_not_sum:
                    lesion_sum += lm
                else:
                    lesion_sum.append(lm)
                if include_brain_masks:
                    brain_avg += bm
                if row["label"] == "MS":
                    if not list_not_sum:
                        lesion_sum_ms += lm
                    else:
                        lesion_sum_ms.append(lm)
                    if include_brain_masks:
                        brain_avg_ms += bm
                else:
                    if not list_not_sum:
                        lesion_sum_hc += lm
                    else:
                        lesion_sum_hc.append(lm)
                    if include_brain_masks:
                        brain_avg_hc += bm
    if not include_brain_masks:
        return lesion_sum, lesion_sum_ms, lesion_sum_hc
    else:
        return lesion_sum, lesion_sum_ms, lesion_sum_hc, brain_avg, brain_avg_ms, brain_avg_hc
lesion_sum, lesion_sum_ms, lesion_sum_hc, brain_avg, brain_avg_ms, brain_avg_hc = get_lesion_sum(include_holdout=False, include_brain_masks=True)
# Divide by number of samples per class
num_train_hc = 61 - 2
num_train_ms = 63 - 2

lesion_sum_hc /= num_train_hc
lesion_sum_ms /= num_train_ms
sl = 92
fig, axes = plt.subplots(ncols=2, figsize=(7, 6),
                         gridspec_kw = {'wspace':0.1, 'hspace':0})

for ax in axes.ravel():
    ax.axis('off')
    
axes[0].imshow(np.rot90(brain_avg_ms[:,:,sl], axes=(0, 1)), cmap='gray')
axes[0].imshow(np.rot90(lesion_sum_ms[:,:,sl], axes=(0, 1)), cmap='hot', alpha=0.6)
axes[0].set_title("MS")

axes[1].imshow(np.rot90(brain_avg_hc[:,:,sl], axes=(0, 1)), cmap='gray')
axes[1].imshow(np.rot90(lesion_sum_hc[:,:,sl], axes=(0, 1)), cmap='hot', alpha=0.6)
axes[1].set_title("Healthy Controls")

plt.subplots_adjust(left=0.0, right=1, top=1, bottom=0.0)

file_path = os.path.join(image_file_dir, "lesion_distribution_on_average_brain.png")
#plt.savefig(file_path, format='png', transparent=True, bbox_inches='tight')
plt.show()

lesion_sum, lesion_sum_ms, lesion_sum_hc = get_lesion_sum()
# Summary of JHU white matter areas
# Name: (min, max)
all_areas= {
    "Middle cerebellar peduncle": (1, 2),
    "Corpus callosum": (3, 5),
    "Fornix": (6, 6),
    "Corticospinal tract": (7, 8),
    "Medial lemniscus": (9, 10),
    "Inferior cerebellar peduncle": (11, 12),
    "Superior cerebellar peduncle": (13, 14),
    "Cerebral peduncle": (15, 16),
    "Anterior limb of internal capsule": (17, 18),
    "Posterior limb of internal capsule": (19, 20),
    "Retrolenticular part of internal capsule": (21, 22),
    "Anterior corona radiata": (23, 24),
    "Superior corona radiata": (25, 26),
    "Posterior corona radiata": (27, 28),
    "Posterior thalamic radiation": (29, 30),
    "Sagittal stratum": (31, 32),
    "External capsule": (33, 34),
    "Cingulum": (35, 38),
    "Superior longitudinal fasciculus": (41, 42),
    "Superior fronto-occipital fasciculus": (43, 44),
    "Uncinate fasciculus": (45, 46),
    "Tapetum": (47, 48),   
}<jupyter_output><empty_output><jupyter_text>Absolute values<jupyter_code>regional_lm = heatmap_per_region(lesion_sum_ms, wm_atlas, positive=True, size_normalize=False)
area_lm = aggregate_regions(regional_lm, all_areas)
# sort values
area_lm_sorted_ms = sorted(area_lm.items(), key=lambda kv: kv[1])
keys_sorted_ms = [row[0] for row in area_lm_sorted_ms]
values_sorted_ms = [row[1] for row in area_lm_sorted_ms]

regional_lm = heatmap_per_region(lesion_sum_hc, wm_atlas, positive=True, size_normalize=False)
area_lm = aggregate_regions(regional_lm, all_areas)
# sort values
area_lm_sorted_hc = sorted(area_lm.items(), key=lambda kv: kv[1])
keys_sorted_hc = [row[0] for row in area_lm_sorted_hc]
values_sorted_hc = [row[1] for row in area_lm_sorted_hc]

plt.figure(figsize=(10, 6))
plt.plot(keys_sorted_ms, values_sorted_ms, 'o', color=ms_color, label="MS")
plt.plot(keys_sorted_hc, values_sorted_hc, 'o', color=hc_color,  label="HC")
plt.xticks(rotation='vertical')
plt.legend(loc="center left")
plt.title("Lesions within white matter areas; nor normalized")
plt.show()
dict(area_lm_sorted_ms)
pd.DataFrame(dict(area_lm_sorted_ms), index=["MS"])
df_ms = pd.DataFrame(dict(area_lm_sorted_ms), index=["MS"])
df_hc = pd.DataFrame(dict(area_lm_sorted_hc), index=["HC"])
df = df_ms.append(df_hc)
df
#plt.figure()
fig, axes = plt.subplots(figsize=(12, 10))
axes.violinplot(dataset=[df_ms[key].values.astype(np.float) for key in keys_sorted_ms], showmeans=True)
axes.violinplot(dataset=[df_hc[key].values.astype(np.float) for key in keys_sorted_ms], showmeans=True)

axes.set_xticks(range(1, len(keys_sorted_ms)+1))
axes.set_xticklabels([key for key in keys_sorted_ms])
plt.xticks(rotation='vertical')
ms_patch = mpatches.Patch(color=ms_color, alpha=0.5)
hc_patch = mpatches.Patch(color=hc_color, alpha=0.5)
axes.legend([ms_patch, hc_patch], ["MS", "HC"], loc="center left")
plt.title("Positive relevance; not normalized")
plt.show()<jupyter_output><empty_output><jupyter_text>Normalized by size<jupyter_code>regional_lm = heatmap_per_region(lesion_sum_ms, wm_atlas, positive=True, size_normalize=True)
area_lm = aggregate_regions(regional_lm, all_areas)
# sort values
area_lm_sorted_ms = sorted(area_lm.items(), key=lambda kv: kv[1])
keys_sorted_ms = [row[0] for row in area_lm_sorted_ms]
values_sorted_ms = [row[1] for row in area_lm_sorted_ms]

regional_lm = heatmap_per_region(lesion_sum_hc, wm_atlas, positive=True, size_normalize=True)
area_lm = aggregate_regions(regional_lm, all_areas)
# sort values
area_lm_sorted_hc = sorted(area_lm.items(), key=lambda kv: kv[1])
keys_sorted_hc = [row[0] for row in area_lm_sorted_hc]
values_sorted_hc = [row[1] for row in area_lm_sorted_hc]

plt.figure(figsize=(10, 6))
plt.plot(keys_sorted_ms, values_sorted_ms, 'o', color=ms_color, label="MS")
plt.plot(keys_sorted_hc, values_sorted_hc, 'o', color=hc_color,  label="HC")
plt.xticks(rotation='vertical')
plt.legend(loc="center left")
plt.title("Lesions within white matter areas; normalized by size")
plt.show()<jupyter_output><empty_output><jupyter_text>Violin plots<jupyter_code>lesion_sum, lesion_sum_ms, lesion_sum_hc = get_lesion_sum(list_not_sum=True)
# individual lesion mask
plt.imshow(lesion_sum_ms[0][:,:,96], cmap='gray')<jupyter_output><empty_output><jupyter_text>JHU ICBM Atlas<jupyter_code>keys_sorted_ms, values_sorted_ms = get_area_relevance(lesion_sum_ms, wm_atlas, wm_areas, positive=True, size_normalize=False)
keys_sorted_hc, values_sorted_hc = get_area_relevance(lesion_sum_hc, wm_atlas, wm_areas, positive=True, size_normalize=False)
hm_distrib_ms = {}
for key, ind_value in zip(keys_sorted_ms, values_sorted_ms):
    for area, value in zip(key, ind_value):
        if area in hm_distrib_ms.keys():
            hm_distrib_ms[area].append(value)
        else:
            hm_distrib_ms[area] = [value,]
            
hm_distrib_hc = {}
for key, ind_value in zip(keys_sorted_hc, values_sorted_hc):
    for area, value in zip(key, ind_value):
        if area in hm_distrib_hc.keys():
            hm_distrib_hc[area].append(value)
        else:
            hm_distrib_hc[area] = [value,]
df_ms = pd.DataFrame(hm_distrib_ms)
df_hc = pd.DataFrame(hm_distrib_hc)<jupyter_output><empty_output><jupyter_text>Not normalized<jupyter_code>#plt.figure()
fig, axes = plt.subplots(figsize=(10, 5))
axes.violinplot(dataset=[df_ms[key] for key in keys_sorted_ms[0]], showmeans=True)
axes.violinplot(dataset=[df_hc[key] for key in keys_sorted_ms[0]], showmeans=True)

axes.set_xticks(range(1, len(keys_sorted_ms[0])+1))
axes.set_xticklabels([key for key in keys_sorted_ms[0]])
plt.xticks(rotation='vertical')
ms_patch = mpatches.Patch(color=ms_color, alpha=0.5)
hc_patch = mpatches.Patch(color=hc_color, alpha=0.5)
axes.legend([ms_patch, hc_patch], ["MS", "HC"], loc="center left")
plt.title("Distribution of lesions according to JHU-ICBM-DTI atlas")

file_path = os.path.join(image_file_dir, "lesion_distribution_WM_absolute.png")
#plt.savefig(file_path, format='png', transparent=True, bbox_inches='tight')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13172721973965418, 'fraction_numerical': 0.015154458908101806, 'mean_word_length': 3.2053376906318083, 'pattern_counts': {'":': 29, '<': 29, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '426064', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5995, 'n_tokens_neox': 5582, 'n_words': 1114, 'path': 'Other investigations/Count lesion ratio of total tissue and area wise.ipynb', 'repo_name': 'derEitel/explainableMS', 'size': '15441'}