Update analysis scripts

Author: alexander-ngo

src/analyses/03_atrophyAssociation.py

@@ -25,14 +25,19 @@ def load_multi_tle():
     """
     return util.load_data(
         "../../data/processed/multi_tle_data.npz",
-        ["age", "sex", "group", "focus", "ct"],
+        ["age", "sex", "focus", "ct"],
     )
 
 
 # Main analysis
 def main():
     load_dataset = {"local": load_local_tle, "multi": load_multi_tle}
     for key in load_dataset:
+        
+        print("------------------------")
+        print(f"Dataset: {key}")
+        print("------------------------")
+        
         # Load LTE data
         age, sex, focus, ct = load_dataset[key]()
 

src/analyses/04_epicentreAssociation.py

@@ -9,7 +9,7 @@ def main():
         "network"
     )
 
-    for site in ["local", "multi"]:
+    for site in ["multi"]: #["local", "multi"]:
         print("------------------------")
         print(f"Dataset: {site}")
         print("------------------------")

src/analyses/s01_subcorticalCorrelation.py

@@ -83,7 +83,7 @@ def main():
         print(f"{sctx}: r = {r[i]}, r-squared = {r2[i]}, p = {p[i]}")
 
     global_association = pd.DataFrame(
-        {"subcortical": subcortical, "r": r, "r2": r2, "p": p}
+        {"subcortical": subcortical.keys(), "r": r, "r2": r2, "p": p}
     )
 
     print()

src/analyses/s02_igeSpecificity.py

@@ -37,7 +37,7 @@ def main():
     x = util.zscore_flip(ct, group, "HC", "_")
     covar = pd.DataFrame({"age": age, "sex": sex})
 
-    slm = util.casecontrol_difference(x, covar, group, "C", "X")
+    slm = util.casecontrol_difference(x, covar, group, "HC", "IGE")
     atrophy["ige"] = slm
 
     print()

src/analyses/s03_psychiatryAtrophySpecificity.py

@@ -19,20 +19,12 @@ def load_psychiatry():
           the cohen's D effect size from the summary statistics for that condition.
     """
     atrophy = {
-        "adhd": load_summary_stats("adhd")["CortThick_case_vs_controls_adult"]["d_icv"],
-        "asd": load_summary_stats("asd")["CortThick_case_vs_controls_meta_analysis"][
-            "d_icv"
-        ],
-        "bd": load_summary_stats("bipolar")["CortThick_case_vs_controls_adult"][
-            "d_icv"
-        ],
-        "mdd": load_summary_stats("depression")["CortThick_case_vs_controls_adult"][
-            "d_icv"
-        ],
-        "ocd": load_summary_stats("ocd")["CortThick_case_vs_controls_adult"]["d_icv"],
-        "scz": load_summary_stats("schizophrenia")["CortThick_case_vs_controls"][
-            "d_icv"
-        ],
+        "adhd": load_summary_stats("adhd")["CortThick_case_vs_controls_adult"]["d_icv"].to_numpy(),
+        "asd": load_summary_stats("asd")["CortThick_case_vs_controls_meta_analysis"]["d_icv"].to_numpy(),
+        "bd": load_summary_stats("bipolar")["CortThick_case_vs_controls_adult"]["d_icv"].to_numpy(),
+        "mdd": load_summary_stats("depression")["CortThick_case_vs_controls_adult"]["d_icv"].to_numpy(),
+        "ocd": load_summary_stats("ocd")["CortThick_case_vs_controls_adult"]["d_icv"].to_numpy(),
+        "scz": load_summary_stats("schizophrenia")["CortThick_case_vs_controls"]["d_icv"].to_numpy(),
     }
 
     return atrophy
@@ -44,7 +36,7 @@ def main():
     atrophy = load_psychiatry()
 
     # Load imaging genetic result
-    imaging_genetics = util.load_imaging_genetic("regional")
+    imaging_genetics = util.load_imaging_genetic("regional").t
 
     association = {}
     print("----------------------")
@@ -54,7 +46,7 @@ def main():
         print()
         print(f"{psy}")
         print("----------------------")
-        r, p, null = util.spatial_correlation(atrophy[psy], imaging_genetics.t)
+        r, p, null = util.spatial_correlation(atrophy[psy], imaging_genetics)
 
         association[psy] = {"r": r, "p": p, "null": null}
 

src/analyses/s04_psychiatryEpicentreSpecificity.py

@@ -82,12 +82,12 @@ def main():
     print("Save results")
     print("------------")
     util.save_to_pickle(
-        "../../data/results/s04_psychiatricEpicentreSpecificity/psychiatric_epicentre.pkl",
+        "../../data/results/s04_psychiatryEpicentreSpecificity/psychiatric_epicentre.pkl",
         {"epicentre": epicentre},
     )
 
     util.save_to_pickle(
-        "../../data/results/s04_psychiatricsEpicentreSpecificity/psychiatric_association.pkl",
+        "../../data/results/s04_psychiatryEpicentreSpecificity/psychiatric_association.pkl",
         {"association": association},
     )
 

src/analyses/s05_thresholdConsistency.py

@@ -21,19 +21,33 @@ def load_abcd():
     )
 
 
-def load_local_tle():
+def load_multi_tle():
     """
-    Load local TLE (EpiC x MICs x NKG) data
+    Load multisite TLE (ENIGMA) data
 
     Returns:
-    numpy.ndarray: The loaded data containing age, sex, dataset, focus, and ct information.
+    numpy.ndarray: The loaded data containing age, sex, dataset, group, focus, and ct information.
     """
     return util.load_data(
-        "../../data/processed/local_tle_data.npz",
-        ["age", "sex", "dataset", "focus", "ct"],
+        "../../data/processed/multi_tle_data.npz",
+        ["age", "sex", "focus", "ct"],
+    )
+    
+    
+def load_multi_ige():
+    """
+    Load multisite IGE (ENIGMA) data
+
+    Returns:
+    numpy.ndarray: The loaded data containing age, sex, dataset, group, focus, and ct information.
+    """
+    return util.load_data(
+        "../../data/processed/multi_ige_data.npz",
+        ["age", "sex", "group", "ct"],
     )
 
 
+
 # Main analysis
 def main():
     # Load data
@@ -72,30 +86,7 @@ def main():
 
     # Lobe indices
     lobe_names = ["whole", "frontal", "limbic", "occipital", "parietal", "temporal"]
-    frontal = [
-        17,
-        18,
-        25,
-        16,
-        22,
-        2,
-        26,
-        10,
-        12,
-        30,
-        15,
-        46,
-        64,
-        60,
-        49,
-        44,
-        59,
-        52,
-        50,
-        56,
-        36,
-        51,
-    ]
+    frontal = [17, 18, 25, 16, 22, 2, 26, 10, 12, 30, 15, 46, 64, 60, 49, 44, 59, 52, 50, 56, 36, 51,]
     limbic = [8, 21, 1, 24, 33, 58, 35, 55, 42, 67]
     occipital = [11, 19, 9, 3, 45, 53, 37, 43]
     parietal = [20, 27, 29, 6, 23, 54, 61, 40, 63, 57]
@@ -182,49 +173,27 @@ def main():
     print("Comparison to case-control atrophy")
     print("----------------------------------")
     # Generate atrophy for each site:
-    age, sex, dataset, focus, ct = load_local_tle()
-    local_tle_atrophy = util.load_result(
-        "../../data/results/03_atrophyAssociation/local_atrophy.pkl", ["atrophy"]
+    epilepsy_atrophy = {}
+    
+    multi_tle_atrophy = util.load_result(
+        "../../data/results/03_atrophyAssociation/multi_atrophy.pkl", ["atrophy"]
     )[0]
-    # Combined
-    tle_atrophy = {}
-    tle_atrophy["lh_all"] = local_tle_atrophy["ltle"]
-    tle_atrophy["rh_all"] = local_tle_atrophy["rtle"]
-
-    # Across different datasets
-    for d in ["EpiC", "MICs"]:
-        dataset_idx = dataset == d
-        dataset_ct = ct[dataset_idx, :]
-        dataset_age = age[dataset_idx]
-        dataset_sex = sex[dataset_idx]
-        dataset_focus = focus[dataset_idx]
-
-        x = util.zscore_flip(dataset_ct, dataset_focus, "C", "R")
-        for hemi in ["L", "R"]:
-            covar = pd.DataFrame({"age": dataset_age, "sex": dataset_sex})
-            tle_atrophy[f"{hemi.lower()}h_{d.lower()}"] = util.casecontrol_difference(
-                x, covar, dataset_focus, "C", hemi
-            )
-
-        print(f"{d} atrophy done")
-
-    site_labels = [
-        "lh_all",
-        "lh_epic",
-        "lh_mics",
-        "rh_all",
-        "rh_epic",
-        "rh_mics",
-    ]
-    atrophy_similarity_r = np.zeros((len(thresholds), len(site_labels)))
-    atrophy_similarity_p = np.zeros((len(thresholds), len(site_labels)))
+    epilepsy_atrophy["ltle"] = multi_tle_atrophy["ltle"]
+    epilepsy_atrophy["rtle"] = multi_tle_atrophy["rtle"]
+    
+    multi_ige_atrophy = util.load_result(
+        "../../data/results/s02_igeSpecificity/ige_atrophy.pkl", ["atrophy"]
+    )[0]
+    epilepsy_atrophy["ige"] = multi_ige_atrophy["ige"]
+    
+    atrophy_similarity_r = np.zeros((len(thresholds), len(epilepsy_atrophy)))
+    atrophy_similarity_p = np.zeros((len(thresholds), len(epilepsy_atrophy)))
     for i in range(len(thresholds)):
         map1 = regional_association[thresholds[i]]
-        for site in site_labels:
-            map2 = tle_atrophy[site]
-            atrophy_similarity_r[i, j], atrophy_similarity_p[i, j], _ = (
-                util.spatial_correlation(map1.t, map2.t, n_rot=5000)
-            )
+        for j, subtype in enumerate(epilepsy_atrophy):
+            map2 = epilepsy_atrophy[subtype]
+            atrophy_similarity_r[i, j], atrophy_similarity_p[i, j], _ = util.spatial_correlation(map1.t, map2.t, n_rot=5000)
+
 
     print()
     print("Save results")

src/analyses/s06_epicentreConsistency.py

@@ -2,20 +2,6 @@
 import utilities as util
 
 
-# Helper functions
-def load_abcd():
-    """
-    Load abcd data
-
-    Returns:
-    numpy.ndarray: The loaded data containing age, sex, dataset, focus, and ct information.
-    """
-    return util.load_data(
-        "../../data/processed/abcd_data.npz",
-        ["age", "sex", "site", "pc10", "thresh", "prs_all", "ct_vertex", "ct_aparc"],
-    )
-
-
 # Main analysis
 def main():
     print("------------------------------------------")
@@ -33,9 +19,9 @@ def main():
     ]
 
     regional_association = util.load_result(
-        "../../data/results/s01_thresholdConsistency/threshold_regional_association.pkl",
+        "../../data/results/s05_thresholdConsistency/threshold_regional_association.pkl",
         ["regional_association"],
-    )
+    )[0]
 
     fc_ctx, fc_sctx, sc_ctx, sc_sctx = util.load_connectomes()
 

src/datasets.py

@@ -254,7 +254,7 @@ def main():
     print("Thickness data")
     print("-------------------------")
     multi_ige_ct = pd.read_csv("../data/raw/multi_ige_ct.csv", index_col="participant")
-    ct = np.transpose(multi_ige_ct.to_numpy())
+    ct = multi_ige_ct.to_numpy()
 
     print()
     print("Save data")