<jupyter_start><jupyter_text>Binary Simplification in MetabolomicsNotebook to support the study on the application of **Bin**ary **Sim**plification as a competing form of pre-processing procedure for high-resolution metabolomics data.This is notebook `paper_binsim_permuts.ipynb`. Organization of the Notebook- Set up database of data sets- Application of different pre-treatments (including BinSim) to each data set- **Permutation tests generation and figure representation**Permutation tests are slow to generate.**Warning**: Permutation tests will write a json file into a folder called **permuts** (if it doesn't exist, it will cause an error). Permutation test will happen if Generate = True in the cell before application of each one. Thus, a folder named **permuts** has to exist to run the notebook. Since the models built depend on the folds created for cross-validation analysis, if the notebook is ran, the results may vary slightly from the exact results presented in the paper.  Needed Imports<jupyter_code>import itertools
from pathlib import Path

# json for persistence
import json
from time import perf_counter

import numpy as np
import pandas as pd
from pandas.testing import assert_frame_equal

import scipy.spatial.distance as dist
import scipy.cluster.hierarchy as hier
import scipy.stats as stats

import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.patches as mpatches
from matplotlib import ticker

import seaborn as sns

# Metabolinks package
import metabolinks as mtl
import metabolinks.transformations as transf

# Python files in the repository
import multianalysis as ma
from elips import plot_confidence_ellipse
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Description of dataset records`datasets` is the global dict that holds all data sets. It is a **dict of dict's**.Each data set is **represented as a dict**.Each record has the following fields (keys):- `name`: the table/figure name of the data set- `source`: the biological source for each dataset- `mode`: the aquisition mode- `alignment`: the alignment used to generate the data matrix- `data`: the data matrix- `target`: the sample labels, possibly already integer encoded- ``: transformed data matrix. These treatment names can be    - `original`: an alias to `data`    - `Ionly`: missing value imputed data, only    - `P`: Pareto scaled data    - `NP`: Pareto scaled and normalized    - `NGP`: normalized, glog transformed and Pareto scaled    - `BinSim`: binary simplified dataThe keys of `datasets` may be shared with dicts holding records resulting from comparison analysis.Here are the keys (and respective names) of datasets used in this study:- GD_neg_global2 (GDg2-)- GD_pos_global2 (GDg2+)- GD_neg_class2 (GDc2-)- GD_pos_class2 (GDc2+)- YD (YD 2/15)- YD2 (YD 6/15)- vitis_types (GD types)  Description of grapevine data setsGrapevine Datasets (Negative and Positive) - 33 samples belonging to 11 different grapevine varieties (3 samples per variety/biological group) of FT-ICR-MS metabolomics data obtained in negative and positive ionization mode.5 different _Vitis_ species (other than _V. vinifera_) varieties:- CAN - 3 Samples (14, 15, 16) of _V. candicans Engelmann_ (VIVC variety number: 13508)- RIP - 3 Samples (17, 18, 19) of _V. ri<NAME>_ (Riparia Gloire de Montpellier, VIVC variety number: 4824) - ROT - 3 Samples (20, 21, 22) of _V. rotundifolia_ (Muscadinia Rotundifolia Michaux cv. Rotundifolia, VIVC variety number: 13586)- RU - 3 Samples (35, 36, 37) of _V. rupestris Scheele_ (Rupestris du lot, VIVC variety number: 10389)- LAB - 3 Samples (8, 9, 10) of _V. labrusca_ (Isabella, VIVC variety number: 5560)6 different _V. vinifera_ cultivars varieties are:- SYL - 3 samples (11, 12, 13) of the subspecies _sylvestris_ (VIVC variety number: -)- CS - 3 Samples (29, 30, 31) of the subspecies _sativa_ cultivar Cabernet Sauvignon (VIVC variety number: 1929)- PN - 3 Samples (23, 24, 25) of the subspecies _sativa_ cultivar Pinot Noir (VIVC variety number: 9279)- REG - 3 Samples (38, 39, 40) of the subspecies _sativa_ cultivar Regent (VIVC variety number: 4572)- RL - 3 Samples (26, 27, 28) of the subspecies _sativa_ cultivar Riesling Weiss (VIVC variety number: 10077)- TRI - 3 Samples (32, 33, 34) of the subspecies _sativa_ cultivar Cabernet Sauvignon (VIVC variety number: 15685)Data acquired by Maia et al. (2020):- <NAME>, <NAME>, <NAME>, et al. Integrating metabolomics and targeted gene expression to uncover potential biomarkers of fungal / oomycetes ‑ associated disease susceptibility in grapevine. Sci Rep. Published online 2020:1-15. doi:10.1038/s41598-020-72781-2- <NAME>, <NAME>, <NAME>, <NAME>. Grapevine untargeted metabolomics to uncover potential biomarkers of fungal/oomycetes-associated diseases. 2020. doi:10.6084/m9.figshare.12357314.v1**Peak Alignment** and **Peak Filtering** were performed with function `metabolinks.peak_alignment.align()`. Human leucine enkephalin (Sigma Aldrich) was used as the reference feature (internal standard, [M+H]+ = 556.276575 Da or [M-H]- = 554.262022 Da).**4** data matrices were constructed from this data:- Data sets named `GD_pos_global2` (GDg2+) and `GD_neg_global2` (GDg2-) were generated after retaining only features that occur (globally) at least twice in all 33 samples of the data sets (filtering/alignment) for the **positive mode** data acquisition and the **negative mode** data acquisition, respectively.- Data sets named `GD_pos_class2` (GDc2+) and `GD_neg_class2` (GDc2-) were generated after retaining only features that occur at least twice in the three replicates of at least one _Vitis_ variety in the data sets (filtering/alignment) for the **positive mode** data acquisition and the **negative mode** data acquisition, respectively.For the purpose of assessing the performance of supervised methods each of these four datasets was used with target labels defining classes corresponding to replicates of each of the 11 Vitis species/cultivars.For the purpose of assessing the performance of supervised methods under a binary (two-class) problem, data set `GD_neg_class2` was also used with target labels defining two classes: Vitis vinifera cultivars and "wild", non-vinifera Vitis species. This is dataset `vitis_types` (GD types).  Description of the yeast data setYeast dataset - 15 samples belonging to 5 different yeast strains of _Saccharomyces cerevisiae_ (3 biological replicates per strain/biological group) of FT-ICR-MS metabolomics data obtained in positive ionization mode. The 5 strains were: the reference strain BY4741 (represented as BY) and 4 single-gene deletion mutants of this strain – ΔGLO1 (GLO1), ΔGLO2 (GLO2), ΔGRE3 (GRE3) and ΔENO1 (ENO1). These deleted genes are directly or indirectly related to methylglyoxal metabolism.Data acquired by Luz et al. (2021):- <NAME>, <NAME>, <NAME>, <NAME>. FT-ICR-MS based untargeted metabolomics for the discrimination of yeast mutants. 2021. doi:10.6084/m9.figshare.15173559.v1**Peak Alignment** and **Peak Filtering** was performed with MetaboScape 4.0 software (see reference above for details in sample preparation, pre-processing, formula assignment). In short, Yeast Dataset was obtained with Electrospray Ionization in Positive Mode and pre-processed by MetaboScape 4.0 (Bruker Daltonics). Human leucine enkephalin (Sigma Aldrich) was used as the reference feature (internal standard, [M+H]+ = 556.276575 Da or [M-H]- = 554.262022 Da).**2** data matrices were constructed from this data:- Data set named `YD` (YD 2/15) was generated after retaining only features that occur (globally) at least twice in all 15 samples (filtering/alignment).- Data set named `YD2` (YD 6/15) was generated after retaining only features that occur (globally) at least six times in all 15 samples (filtering/alignment).For the purpose of assessing the performance of supervised methods, this data set was used with target labels defining classes corresponding to replicates of each of the 4 yeast strains.  Load datasets<jupyter_code>path = Path.cwd() / "paperimages" / 'processed_data.json'
storepath = Path.cwd() / "paperimages" / 'processed_data.h5'
with pd.HDFStore(storepath) as store:

    with open(path, encoding='utf8') as read_file:
        datasets = json.load(read_file)
    
    for dskey, dataset in datasets.items():
        for key in dataset:
            value = dataset[key]
            if isinstance(value, str) and value.startswith("INSTORE"):
                storekey = value.split("_", 1)[1]
                dataset[key] = store[storekey]
            # convert colors to tuples, since they are read as lists from json file
            elif key == 'label_colors':
                dataset[key] = {lbl: tuple(c) for lbl, c in value.items()}
            elif key == 'sample_colors':
                dataset[key] = [tuple(c) for c in value]<jupyter_output><empty_output><jupyter_text>Supervised Statistical Analysis - Permutation TestsThe Supervised Statistical Analysis methods used will be Random Forest and PLS-DA.The performance of the classifiers will be evaluated by their predictive **accuracy** (which will always be estimated by internal stratified 3-fold cross-validation or 5-fold cross-validation in `vitis_types`).  Permutation Tests (Very Slow)Permutation tests is based on shuffling the labels of the different samples, shuflling the groups where they belong with the intent to see if the classifier tested, whether it is Random Forest or PLS-DA found a significant class structure in the data - assess the significance of the predictive accuracy results. For that a random k-fold cross-validation is performed on the original dataset (to serve as a comparation point) and on 500 permutations of datasets with labels randomly shuffled around. The models are evaluated by their predictive accuracies. The empirical p-value is given by (the number of times the permutation accuracy was bigger than the random k-fold cross-validation made with the original dataset + 1) / (number of permutations + 1) (source: Ojala and Garriga, 2010).<NAME>, Garriga GC. Permutation Tests for Studying Classifier Performance. In: 2009 Ninth IEEE International Conference on Data Mining. ; 2009:908-913. doi:10.1109/ICDM.2009.108Histograms with the prediction accuracy of the different permutations were plotted and compared to the accuracy got with the original dataset. Permutation Tests - Random ForestsUse of `permutation_RF` function from multianalysis.py. See details about the application of this function in the multianalysis.py file.<jupyter_code># json for persistence
import json
from time import perf_counter<jupyter_output><empty_output><jupyter_text>Use `GENERATE = True` to perform permutation tests and persist results in json<jupyter_code>GENERATE = True #False
if GENERATE:
    iter_num=500 # number of permutations


    permuts_RF = []

    to_permute = [name for name in datasets]# if 'global2' in name]
    for name in to_permute:
        for treatment in ('P', 'P_RF', 'NP', 'NP_RF', 'NGP', 'NGP_RF', 'BinSim'):
            dataset = datasets[name]
            print(f'{iter_num} permutations (Random Forest) for {name} with treatment {treatment}', end=' ...')
            n_fold = 5 if name == 'vitis_types' else 3
            start = perf_counter()
            permutations = ma.permutation_RF(dataset[treatment], dataset['target'], iter_num=iter_num, n_fold=n_fold, n_trees=100)
            res = {'dataset': name, 'treatment': treatment,
                   'non_permuted_CV': permutations[0],
                   'permutations': permutations[1],
                   'p-value': permutations[2]}       
            permuts_RF.append(res)
            end = perf_counter()
            pvalue = permutations[2]
            print(f'Done! took {(end - start):.3f} s, p-value = {pvalue}')
    
    # Store in json file
    fname = 'paperimages/permuts_rf.json'
    with open(fname, "w", encoding='utf8') as write_file:
        json.dump(permuts_RF, write_file)<jupyter_output>500 permutations (Random Forest) for GD_neg_global2 with treatment P ...Done! took 208.334 s, p-value = 0.001996007984031936
500 permutations (Random Forest) for GD_neg_global2 with treatment P_RF ...Done! took 217.220 s, p-value = 0.001996007984031936
500 permutations (Random Forest) for GD_neg_global2 with treatment NP ...Done! took 222.863 s, p-value = 0.001996007984031936
500 permutations (Random Forest) for GD_neg_global2 with treatment NP_RF ...Done! took 221.893 s, p-value = 0.001996007984031936
500 permutations (Random Forest) for GD_neg_global2 with treatment NGP ...Done! took 215.747 s, p-value = 0.001996007984031936
500 permutations (Random Forest) for GD_neg_global2 with treatment NGP_RF ...Done! took 215.858 s, p-value = 0.001996007984031936
500 permutations (Random Forest) for GD_neg_global2 with treatment BinSim ...Done! took 198.885 s, p-value = 0.001996007984031936
500 permutations (Random Forest) for GD_pos_global2 with treatment P ...Done! took 252.927 s, p-value = 0[...]<jupyter_text>Permutation Tests - PLS-DAUse of `permutation_PLSDA` function from multianalysis.py. See details about the application of this function in the multianalysis.py file.<jupyter_code>%%capture --no-stdout
if GENERATE:
    iter_num=500

    permuts_PLSDA = []

    to_permute = [name for name in datasets] # if 'global2' in name]
    for name in to_permute:
        for treatment in ('P', 'P_RF', 'NP', 'NP_RF', 'NGP', 'NGP_RF', 'BinSim'):
            dataset = datasets[name]
            print(f'Permutation test (PLS-DA) for {name} with treatment {treatment}', end=' ...')
            n_comp = 11 if name.startswith('GD') else 6
            n_fold = 5 if name == 'vitis_types' else 3
            start = perf_counter()
            permutations = ma.permutation_PLSDA(dataset[treatment], dataset['target'], n_comp=11,
                                                iter_num=iter_num, n_fold=n_fold)
            res = {'dataset': name, 'treatment': treatment,
                   'non_permuted_CV': permutations[0],
                   'permutations': permutations[1],
                   'p-value': permutations[2]}
            permuts_PLSDA.append(res)
            end = perf_counter()
            pvalue = permutations[2]
            print(f'Done! took {(end - start):.3f} s, p-value = {pvalue:10.5f}')
            
    # Store in json file
    fname = 'paperimages/permuts_plsda.json'
    with open(fname, "w", encoding='utf8') as write_file:
        json.dump(permuts_PLSDA, write_file)
# Get data from json file - random forests
fname = 'paperimages/permuts_rf.json'
with open(fname, "r", encoding='utf8') as read_file:
    permuts_RF = json.load(read_file)

for p in permuts_RF:
    print(f"{p['dataset']:<20}{p['treatment']:<8}{p['p-value']:10.5f}")
# Get data from json file - PLS-DA
fname = 'paperimages/permuts_plsda.json'
with open(fname, "r", encoding='utf8') as read_file:
    permuts_PLSDA = json.load(read_file)

for p in permuts_PLSDA:
    print(f"{p['dataset']:<20}{p['treatment']:<8}{p['p-value']:10.5f}")<jupyter_output>GD_neg_global2      P          0.00200
GD_neg_global2      P_RF       0.00200
GD_neg_global2      NP         0.00200
GD_neg_global2      NP_RF      0.00200
GD_neg_global2      NGP        0.00200
GD_neg_global2      NGP_RF     0.00200
GD_neg_global2      BinSim     0.00200
GD_pos_global2      P          0.00200
GD_pos_global2      P_RF       0.00200
GD_pos_global2      NP         0.00200
GD_pos_global2      NP_RF      0.00200
GD_pos_global2      NGP        0.00200
GD_pos_global2      NGP_RF     0.00200
GD_pos_global2      BinSim     0.00200
GD_neg_class2       P          0.00200
GD_neg_class2       P_RF       0.00200
GD_neg_class2       NP         0.00200
GD_neg_class2       NP_RF      0.00200
GD_neg_class2       NGP        0.00200
GD_neg_class2       NGP_RF     0.00200
GD_neg_class2       BinSim     0.00200
GD_pos_class2       P          0.00200
GD_pos_class2       P_RF       0.00200
GD_pos_class2       NP         0.00200
GD_pos_class2       NP_RF      0.00200
GD_pos_class2       NGP  [...]<jupyter_text>Plot the Permutations test results - HistogramsNº Occurences (of the permutations) vs CV prediction acuracy - The distribution of average prediction accuracy of 500 permutations1st Figure - Random Forest2nd Figure - PLS-DA<jupyter_code>with sns.axes_style("whitegrid"):
    with sns.plotting_context("notebook", font_scale=1.2):
        f, axes = plt.subplots(7, 4, figsize = (10,16), sharey='row', sharex='col')
        colors = sns.color_palette('tab10', 4)
        ylim = [0,200]
        treatments = ['P','NP','NGP','BinSim']
        
        for row, dskey in enumerate(datasets):

            to_plot = [p for p in permuts_RF if p['dataset'] == dskey]

            n_labels = len(datasets[dskey]['target'])
            n_bins = 34 if dskey == 'vitis_types' else 16

            for ax, p, tname, color in zip(axes[row].ravel(), to_plot, treatments, colors):
                ax.hist(np.array(p['permutations'])*100, range=(0, 100.01), label=name + ' Permutations',
                        bins=n_bins, edgecolor='black', color=color)
                #ax.axvline(p['non_permuted_CV']*100)

                ax.plot(2 * [p['non_permuted_CV'] * 100], ylim, '-', linewidth=3, color=color, #alpha = 0.5,
                         label=name + ' (pvalue %.5f)' % p['p-value'], solid_capstyle='round')
                #q.set(xlabel='Cross Validation Prediction Accuracy (%)', ylabel='Nº of occurrences', fontsize=12)
                #ax.set_xlabel('CV Prediction Accuracy (%)')
                #ax.set_ylabel('Nº of occurrences')
                props = dict(boxstyle='round', facecolor='white', alpha=1)
                ax.text(90, 190, 'p-value = %.3f' % p['p-value'], bbox=props, ha='right', fontsize='small')
                ax.set_title(f"{datasets[dskey]['name']}, {tname}", color=color)
                ax.set_ylim(0,250)

        f.text(0.5, 0.01, 'CV Prediction Accuracy (%)', ha='center', va='top')
        #f.text(0, 0.5, 'Counts', ha='center', va='top',rotation=90)
        
        
        f.suptitle(f'Permutation tests for Random Forests')
        plt.tight_layout()
        f.savefig('paperimages/permutations_RF.pdf', dpi=200)
        f.savefig('paperimages/permutations_RF.png', dpi=600)
        plt.show()
with sns.axes_style("whitegrid"):
    with sns.plotting_context("notebook", font_scale=1.2):
        f, axes = plt.subplots(7, 4, figsize = (10,16), sharey='row', sharex='col')
        colors = sns.color_palette('tab10', 4)
        ylim = [0,200]
        treatments = ['P','NP','NGP','BinSim']
        
        for row, dskey in enumerate(datasets):

            to_plot = [p for p in permuts_PLSDA if p['dataset'] == dskey]

            n_labels = len(datasets[dskey]['target'])
            
            n_bins = 34 if dskey == 'vitis_types' else 16

            for ax, p, tname, color in zip(axes[row].ravel(), to_plot, treatments, colors):
                ax.hist(np.array(p['permutations'])*100, range=(0, 100.01), label=name + ' Permutations',
                        bins=n_bins, edgecolor='black', color=color)
                #ax.axvline(p['non_permuted_CV']*100)

                ax.plot(2 * [p['non_permuted_CV'] * 100], ylim, '-', linewidth=3, color=color, #alpha = 0.5,
                         label=name + ' (pvalue %.5f)' % p['p-value'], solid_capstyle='round')
                #q.set(xlabel='Cross Validation Prediction Accuracy (%)', ylabel='Nº of occurrences', fontsize=12)
                #ax.set_xlabel('CV Prediction Accuracy (%)')
                #ax.set_ylabel('Nº of occurrences')
                props = dict(boxstyle='round', facecolor='white', alpha=1)
                ax.text(90, 190, 'p-value = %.3f' % p['p-value'], bbox=props, ha='right', fontsize='small')
                ax.set_title(f"{datasets[dskey]['name']}, {tname}", color=color)
                ax.set_ylim(0,250)

        f.text(0.5, 0.01, 'CV Prediction Accuracy (%)', ha='center', va='top')
        #f.text(0, 0.5, 'Counts', ha='center', va='top',rotation=90)
        
        
        f.suptitle(f'Permutation tests for PLS-DA')
        plt.tight_layout()
        f.savefig('paperimages/permutations_PLSDA.pdf', dpi=200)
        f.savefig('paperimages/permutations_PLSDA.png', dpi=600)
        plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10113742675661037, 'fraction_numerical': 0.04549707026441479, 'mean_word_length': 3.4549243255099804, 'pattern_counts': {'":': 0, '<': 40, '<?xml version=': 0, '>': 36, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '377750', 'license': 'MIT', 'n_tokens_mistral': 6995, 'n_tokens_neox': 6056, 'n_words': 2139, 'path': 'paper_binsim_permuts.ipynb', 'repo_name': 'aeferreira/binsim_paper', 'size': '20349'}