Replication: Zeighami et al, 2019¶

import livingpark_utils

utils = livingpark_utils.LivingParkUtils()
utils.notebook_init()

required_files = [
    FILENAME_PARTICIPANT_STATUS,
    FILENAME_PD_HISTORY,
    FILENAME_DEMOGRAPHICS,
    FILENAME_AGE,
    FILENAME_MOCA,
    FILENAME_UPDRS1A,
    FILENAME_UPDRS1B,
    FILENAME_UPDRS2,
    FILENAME_UPDRS3,
    FILENAME_ADL,
]

downloader = livingpark_utils.download.ppmi.Downloader(utils.study_files_dir)
utils.get_study_files(required_files, default=downloader)

df_status = load_ppmi_csv(utils, FILENAME_PARTICIPANT_STATUS)

cohort_t1_map = {}
for cohort_name in [MAIN_COHORT, VALIDATION_COHORT]:

    print(f"=============== {cohort_name.capitalize()} cohort ===============")

    df_t1_subset = get_t1_cohort(
        utils,
        cohort_name=cohort_name,
        filename=FILENAME_T1_INFO,
        sagittal_only=True,
    )
    cohort_t1_map[cohort_name] = df_t1_subset

    # cohort composition: number of PD patients/healthy controls
    print(
        df_status.loc[
            df_status[COL_PAT_ID].isin(df_t1_subset[COL_PAT_ID]), COL_STATUS
        ].value_counts()
    )

col_date_diff = "date_diff"
max_date_followup = MAX_DATES[VALIDATION_COHORT]

re_followup_compiled = re.compile(REGEX_VISIT_FOLLOWUP)

cohort_assessments_map: dict[str, pd.DataFrame] = {}
for cohort_name in cohort_assessments_map_orig:

    print(f"========== {cohort_name.upper()} COHORT ==========")

    date_diffs = []

    df_assessments_cohort: pd.DataFrame = cohort_assessments_map_orig[cohort_name]
    df_assessments_baseline = df_assessments_cohort.loc[
        df_assessments_cohort[COL_VISIT_TYPE] == VISIT_BASELINE
    ]
    df_assessments_screening = df_assessments_cohort.loc[
        df_assessments_cohort[COL_VISIT_TYPE] == VISIT_SCREENING
    ]

    # try to fill in missing baseline data
    for idx_row_baseline, row_baseline in df_assessments_baseline.iterrows():

        subject = row_baseline[COL_PAT_ID]
        date_baseline = row_baseline[COL_DATE_INFO]

        # for each score columns
        for col in [COL_MOCA]:

            # fill missing values with screening data
            if pd.isna(row_baseline[col]):

                df_screening_subject = df_assessments_screening.loc[
                    df_assessments_screening[COL_PAT_ID] == subject
                ]

                # some subjects in validation set had multiple screening visits
                # in this case we sort them by how close they are to the baseline visit
                n_screening = len(df_screening_subject)
                if n_screening > 1:
                    df_screening_subject[col_date_diff] = (
                        date_baseline - df_screening_subject[COL_DATE_INFO]
                    )
                    df_screening_subject = df_screening_subject.sort_values(
                        col_date_diff, ascending=True
                    )

                # find corresponding assessment score in screening visits
                for idx_row_screening, row_screening in df_screening_subject.iterrows():
                    new_value = row_screening[col]
                    date_diff = date_baseline - row_screening[COL_DATE_INFO]
                    if not pd.isna(new_value):
                        break

                # replace
                if not pd.isna(new_value):
                    df_assessments_baseline.loc[idx_row_baseline, col] = new_value
                    date_diffs.append(date_diff.days)  # for plotting

    # plot baseline-screening date difference
    date_diffs = pd.Series(date_diffs)
    ax = date_diffs.hist(
        bins=30, alpha=0.7, label=f"{cohort_name} (mean={date_diffs.mean():.2f})"
    )

    # pick latest follow-up visit
    df_assessments_followup = df_assessments_cohort.loc[
        (df_assessments_cohort[COL_PAT_ID].isin(df_assessments_baseline[COL_PAT_ID]))
        & (
            df_assessments_cohort[COL_VISIT_TYPE].apply(
                lambda v: re_followup_compiled.match(v)
            )
        )
    ]

    # keep most recent follow-up data up to retrieval date in original paper
    df_assessments_followup = filter_date(
        df_assessments_followup, max_date=max_date_followup
    )
    df_assessments_followup = (
        df_assessments_followup.groupby(COL_PAT_ID)
        .apply(lambda df: df.sort_values(COL_DATE_INFO, ascending=False).head(1))
        .reset_index(drop=True)
    )

    # we keep subjects who have both baseline and follow-up data
    subjects_common = set(df_assessments_cohort[COL_PAT_ID])
    for df_visit in [df_assessments_baseline, df_assessments_followup]:
        subjects_common = subjects_common & set(df_visit[COL_PAT_ID])

    # print cohort composition
    print(
        df_status.loc[
            df_status[COL_PAT_ID].isin(subjects_common), COL_STATUS
        ].value_counts()
    )

    # combine baseline and follow-up
    df_assessments_baseline[COL_FOLLOWUP] = False
    df_assessments_followup[COL_FOLLOWUP] = True
    df_assessments_merged = pd.concat(
        [df_assessments_baseline, df_assessments_followup], axis="index"
    )
    df_assessments_merged = df_assessments_merged.loc[
        df_assessments_merged[COL_PAT_ID].isin(subjects_common)
    ]

    # impute remaining missing MoCA values
    df_assessments_merged.loc[
        df_assessments_merged[COL_MOCA].isna(), COL_MOCA
    ] = mean_moca

    cohort_assessments_map[cohort_name] = df_assessments_merged

# NOTE: date resolution is months
ax.set_title("Number of days between baseline and screening visits")
ax.legend()

dfs_duration = []
for df_assessments in cohort_assessments_map.values():
    dfs_duration.append(
        df_assessments.groupby(COL_PAT_ID).apply(
            lambda df: (
                df.loc[df[COL_FOLLOWUP], COL_DATE_INFO].item()
                - df.loc[~df[COL_FOLLOWUP], COL_DATE_INFO].item()
            )
        )
    )

df_duration = pd.concat(dfs_duration)
print(
    f"{df_duration.mean().days / 365:.3f} "
    f"(\u00B1 {df_duration.std().days / 365:.3f}) years"
)

def to_1_decimal_str(f):
    return str(round(f, 1))


df_age = load_ppmi_csv(utils, FILENAME_AGE)
df_demographics = load_ppmi_csv(utils, FILENAME_DEMOGRAPHICS)
df_pd_history = load_ppmi_csv(utils, FILENAME_PD_HISTORY)

col_male = "is_male"
col_cohort = "cohort"
col_age_diag = "age_diag"

dfs_summary = []
for cohort_name, df_assessments in cohort_assessments_map.items():

    subjects = df_assessments[COL_PAT_ID].drop_duplicates()
    subjects_pd = df_status.loc[
        (df_status[COL_PAT_ID].isin(subjects)) & (df_status[COL_STATUS] == STATUS_PD),
        COL_PAT_ID,
    ]

    # the paper only report measures for the PD patients
    df_assessments = df_assessments.loc[df_assessments[COL_PAT_ID].isin(subjects_pd)]

    df_summary = df_demographics.merge(
        df_pd_history,
        on=[COL_PAT_ID],
    )
    df_summary[col_age_diag] = (
        df_summary[COL_DATE_PD] - df_summary[COL_DATE_BIRTH]
    ) / np.timedelta64(1, "Y")
    df_summary[col_male] = (df_summary[COL_SEX] == SEX_MALE).apply(
        lambda v: 100 if v else 0
    )
    df_summary = df_summary.loc[
        df_summary[COL_PAT_ID].isin(subjects_pd),
        [COL_PAT_ID, col_age_diag, col_male],
    ]
    df_summary[COL_FOLLOWUP] = False

    # combine with clinical/cognitive scores
    df_summary = df_summary.merge(
        df_assessments, on=[COL_PAT_ID, COL_FOLLOWUP], how="outer"
    )
    df_summary = df_summary.drop(columns=COL_PAT_ID)

    # append
    df_summary[col_cohort] = f"{cohort_name} (n = {len(subjects_pd)})"
    dfs_summary.append(df_summary)

df_summary = pd.concat(dfs_summary)
df_summary_means = (
    df_summary.groupby([col_cohort, COL_FOLLOWUP]).mean().applymap(to_1_decimal_str)
)
df_summary_stds = (
    df_summary.groupby([col_cohort, COL_FOLLOWUP]).std().applymap(to_1_decimal_str)
)
df_summary_stds = " (" + df_summary_stds + ")"
df_summary_stds.loc[:, col_male] = ""
df_summary_combined = (df_summary_means + df_summary_stds).T
df_summary_combined = df_summary_combined.applymap(lambda x: "-" if "nan" in x else x)
df_summary_combined = df_summary_combined.rename(
    index={
        col_age_diag: "Age at diagnosis",
        col_male: "Male (%)",
        COL_UPDRS2: "UPDRS Part II",
        COL_UPDRS3: "UPDRS Part III",
        COL_ADL: "S&E ADL",
        COL_PIGD: "PIGD",
        COL_UPDRS1: "UPDRS Part I",
        COL_MOCA: "MoCA",
        COL_GCO: "Global composite outcome",
    }
)
df_summary_combined

cohort_ids = {}
fpaths_cohort = {}
for cohort_name, df_assessments in cohort_assessments_map.items():

    # get cohort ID
    cohort_id = livingpark_utils.dataset.ppmi.cohort_id(
        df_assessments.drop_duplicates(COL_PAT_ID)
    )
    print(f"{cohort_name} cohort id: {cohort_id}")

    df_t1 = cohort_t1_map[cohort_name]

    # save patient ids
    filename = f"zeighami-etal-2019-cohort-{cohort_id}.csv"
    df_to_save = df_t1.loc[df_t1[COL_PAT_ID].isin(df_assessments[COL_PAT_ID])]
    df_to_save[COL_PAT_ID] = df_to_save[COL_PAT_ID].astype(np.int64)
    df_to_save.drop_duplicates().sort_values(by=COL_PAT_ID).to_csv(
        filename, index=False
    )
    print(f"\tCohort patient ids saved in {filename}")

    cohort_ids[cohort_name] = cohort_id
    fpaths_cohort[cohort_name] = filename

print("Do not share these files publicly!")

_, dpath_outputs, _, _ = check_result_availability(
    cohort_id=cohort_ids[MAIN_COHORT],
    cohort_file=f"zeighami-etal-2019-cohort-{cohort_ids[MAIN_COHORT]}.csv",
    dpath_repo="/data/origami/livingpark/zeighami-etal-2019/",
)

df_scores = pd.read_csv(dpath_outputs / "scores.csv", index_col=0)
df_scores.index = df_status.set_index(COL_PAT_ID).loc[
    df_scores.index.astype(np.int64), COL_STATUS
]
df_scores.index.value_counts()

from scipy.stats import ttest_ind

correct = True
p_values_data = {}

for component in df_scores.columns:
    t_stat, p = ttest_ind(
        df_scores.loc[STATUS_PD, component],
        df_scores.loc[STATUS_HC, component],
    )
    p_values_data[component] = p

p_values = pd.Series(p_values_data)

if correct:
    p_values = p_values * len(p_values)  # Bonferroni correction

i_components_sorted = np.argsort(list(p_values))
i_component_top = i_components_sorted[0]

p_values.sort_values()

import seaborn as sns

x_label = "IC"
y_label = "Mean DBM value"

df_scores_to_plot = pd.melt(
    df_scores.reset_index(),
    id_vars=[COL_STATUS],
    var_name=x_label,
    value_name=y_label,
)

ax = sns.stripplot(
    df_scores_to_plot,
    x=x_label,
    y=y_label,
    hue=COL_STATUS,
    order=p_values.sort_values().index,
    alpha=0.4,
)
ax.tick_params(axis="x", rotation=90)

fig_IC = plot_IC(
    dpath_outputs / "melodic_results" / "melodic_IC.nii.gz",
    i_component=i_components_sorted[0],
)

%%bash
if [ ! -f ${FNAME_PD_ICA} ]
then
    wget "https://neurovault.org/media/images/860/melodic_IC_0007.nii.gz"
else
    echo 'ICA file already exists, skipping download'
fi

fig_IC_orig = plot_IC("melodic_IC_0007.nii.gz")