<commit_before>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Instructions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MetAltas analysis workflow\n",
    "\n",
    "This notebook is for performing a targeted metabolomics analysis on your raw data files.  In order for this to work, you data files must be uploaded to the nersc raw data directory with appropriate permissions.  A fileconverter runs at nersc to convert as follows:\n",
    ".raw > .mzml > .h5 > pactolus.gz  This conversion must happen before you can proceed with data analysis.  A Northen Lab staff member can help you transfer and convert your files if they are not already at nersc.\n",
    "***\n",
    "/global/project/projectdirs/metatlas/raw_data/SUBFOLDER\n",
    "***\n",
    "\n",
    "A targeted analysis requires a list of compounds to search for in your sample files.  Thus, at a minimum, you must have some information on their m/z ratio in order to get started.  The inputs for this notebook, include the location of your raw files, and an atlas listing out these compounds.  \n",
    "\n",
    "Our lab has analyzed thousands of compounds using standardized LCMS methods.  You may use one of these libraries to create an atlas for your sample analysis, assuming you used one of the same standard LCMS methods.\n",
    "\n",
    "Analysis steps\n",
    "\n",
    "<ol>\n",
    "<li>Create your atlas csv external to jupyter: typically you will either have a custom atlas or generate an atlas from one of our EMA libraries.  These can be found in the shared google drive, under atlases.</li>\n",
    "<li>Run through with your initial atlas</li>\n",
    "<li>Use the RT adjuster to mark the quality of the peaks and MSMS matches, and also adjust your RT bounds</li>\n",
    "<li>Reload your atlas from the \"get atlas\" block</li>\n",
    "<li>Export your atlas, filter out the compounds marked under id notes as remove.  Save and upload your new atlas</li>\n",
    "<li>Repeat the analysis using your filtered atlas.</li>\n",
    "    <li>On your final export, make sure you are excluding system blanks and QC files in your output.  If you have extraction blanks you should include those.</li>\n",
    "    </ol>\n",
    "    \n",
    "There is a database that will store your atlas and registered run files.  During the analysis process you will pull these as local variables. These are used to pull out EIC and MSMS data which is stored in a variable called metatlas_dataset.  All output files are generated from that.  It is a local variable so is not stored in the db - thus if your kernel dies partway through, you need to rerun through the notebook to regenerate that variable.\n",
    "\n",
    "IMPORTANT: Anytime you adjust your RTs in the interactive plot, you need to retreive your atlas again before exporting any files.  The changes are stored in the db and the exports are made from locally stored variables/dataframes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# index\n",
    "1. Import python packages\n",
    "2. Create groups\n",
    "3. Select groups\n",
    "4. Create Atlases\n",
    "5. Select Atlas\n",
    "6. Annotate data\n",
    "7. Correct RT bounds\n",
    "8. Export results\n",
    "9. Tools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Helpful tips\n",
    "\n",
    "1. When entering search strings, use % for wildcard: %peas will return \"positive_peas\" and \"positive_greenpeas\", but not \"positive_greenpeas_HILIC\", for that use %peas%, you can also put multiple strings via %pos%peas%\n",
    "2. Set your kernel to mass_spec_cori.  If not sure how, check w Ben or Daniel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# optional: run this block to adjust the width of the notebook.  Change the width percent.\n",
    "\n",
    "from  IPython.core.display  import  display, HTML \n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 1. Import Python Packages\n",
    "## Instructions: \n",
    "\n",
    "1. On the 3rd line of the block below, add your directory where your most recent metatlas code is stored (to the sys.path.insert line).\n",
    "2. Only run this block once to load the metatlas python modules.  If the kernel dies, you will need to rerun this block.\n",
    "3. Run the block with the \"print fa.__file__\" code line to double check that the printed output matches the directory in the sys.path.insert line of the prior block.\n",
    "    If it does not match, then on the \"kernel\" dropdown, click restart and correct the path below before trying again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import sys, os\n",
    "\n",
    "# v edit this line\n",
    "sys.path.insert(0,'/global/homes/FIRST-INITIAL-OF-USERNAME/USERNAME/REPOFOLDER/metatlas-master-20200416/metatlas-master/')\n",
    "# ^ edit this line\n",
    "\n",
    "from metatlas.tools import fastanalysis as fa\n",
    "from metatlas.plots import dill2plots as dp\n",
    "from metatlas.io import metatlas_get_data_helper_fun as ma_data\n",
    "from metatlas.plots import chromplotplus as cpp\n",
    "from metatlas.datastructures import metatlas_objects as metob\n",
    "import qgrid\n",
    "from ipywidgets import interact, interactive, fixed\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display\n",
    "import time\n",
    "import pickle\n",
    "import dill\n",
    "import multiprocessing as mp\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plot\n",
    "import operator\n",
    "from importlib import reload\n",
    "\n",
    "pd.set_option('display.max_rows', 5000)\n",
    "pd.set_option('display.max_columns', 500)\n",
    "pd.set_option('display.max_colwidth', 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check that first block registered metatlas directory correctly.  If it does not look the same, seek help!\n",
    "print(fa.__file__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set atlas, project and output directories from your nersc home directory\n",
    "\n",
    "STEP 1\n",
    "1. Create a project folder name for this analysis by replacing the PROJECTDIRECTORY string text in red below.  Make sure to update the rest of the direcory to point to your home directory.  The pwd block will print out the directory where this jupyter notebook is stored.\n",
    "2. Create a subdirectory name for the output, each run through you may want to create a new output folder.\n",
    "3. When you run the block the folders will be created in your home directory.  If the directory already exists, the block will just set the path for use with future code blocks.\n",
    "\n",
    "STEP 2\n",
    "1.  Enter the nersc path where you have stored your atlas files using one of the two optional lines below:\n",
    "    1. Set a custom path\n",
    "    2. Use the project and output subdirectory from above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#STEP 1\n",
    "project_directory='/global/homes/FIRST-INITIAL-OF-USERNAME/USERNAME/PROJECTDIRECTORY/'  # <- edit this line, do not copy the path directly from NERSC (ex. the u1, or u2 directories)\n",
    "output_subfolder='HILIC_POS_20190830/'  # <- edit this as 'chromatography_polarity_yyyymmdd/'\n",
    "output_dir = os.path.join(project_directory,output_subfolder)\n",
    "\n",
    "if not os.path.exists(project_directory):\n",
    "    os.makedirs(project_directory)\n",
    "if not os.path.exists(output_dir):\n",
    "    os.makedirs(output_dir)\n",
    "\n",
    "#STEP 2\n",
    "pathtoatlas='/global/homes/FIRST-INITIAL-OF-USERNAME/USERNAME/DIRECTORYOFATLAS/'    # <- enter the directory where you have stored your atlas\n",
    "#pathtoatlas= '%s%s' % (project_directory,output_subfolder)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Create Groups (named variables that hold your replicates of each sample)\n",
    "\n",
    "### You must assign your raw files into experimental groups for analysis.  These are used for downstream statistics and for selection of specific groups for filtering to subsets of files for analysis (Ex. just pos or just neg).\n",
    "\n",
    "The groups are created from common file headers and the unique group names. The convention our lab group uses for filenames is as follows: \n",
    "***\n",
    "DATE_NORTHENLABINITIALS_COLLABINITIALS_PROJ_EXP_SAMPSET_SYSTEM_COLUMN-method_SERIAL_POL_ACQ_SAMPLENUMBER_ SAMPLEGROUP_REP_OPTIONAL_SEQ \n",
    "\n",
    "Ex.:20180105_SK_AD_ENIGMA_PseudoInt_R2ADec2017_QE119_50454_123456_POS_MSMS_001_Psyringae-R2A-30C-20hr_Rep01_NA_Seq001.raw\n",
    "***\n",
    "The common header consists of the fields 0-10: DATE_NORTHENLABINITIALS_COLLABINITIALS_PROJ_EXP_SAMPSET_SYSTEM_COLUMN-method_SERIAL_POL_ACQ \n",
    "\n",
    "The sample group name is commonly field # 12 (between underscore 11 and 12) -0 indexed-\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Find your files\n",
    "1. On the first line of the block below, set the 'experiment' and 'name' variables to find your files.  These fields require wildcards for partial string searches\n",
    "2. 'Experiment' is the folder name within global/project/projectdirs/metatlas/raw_data, that will be emailed to you when the files are uploaded to NERSC.  You can also look in the raw_data directory for the NERSC user who uploaded your files; your experiment folder should be in there.\n",
    "3. 'name' is string that will match a subset of your files within that folder.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "files = dp.get_metatlas_files(experiment = '%ENTERSTRING%',name = '%ENTERSTRING%',most_recent = True)\n",
    "# ^ edit the text string in experiment and name fields\n",
    "\n",
    "df = metob.to_dataframe(files)\n",
    "df[['experiment','name','username','acquisition_time']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(files)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# OPTION A: Automated Group Maker"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This will attempt to create groups in an automated fashion (rather than filling out a spreadsheet with a list of files and group names).  If your files are all in one folder at nersc, you can use this options.  If not, use option B below.\n",
    "\n",
    "A long group name consisting of the common header + either controlled vocab value or field #12 along with a short group name (just controlled vocab or field #12) will be stored in a local variable.  The short group names can be used on plots.\n",
    "\n",
    "\n",
    "1. STEP 1: View the groups\n",
    "    1. Pick an experiment folder to look for files in on the metob.retrieve function\n",
    "    2. Enter controlled vocabulary for control files to put select files into groups when control string may be in a different field (not #12) or as a randomly placed substring within a field (ex. if 'InjBl' is included in your controlled vocab list, files like _InjBl-MeOH_ and _StartInjBl_ will group together)\n",
    "    3. If your group name is not between _ 11 and 12 you can adjust those values in the split commands below.  All other (non-controlledvocab) groups will be created from that field.\n",
    "2. STEP 2: Create the groups variable after checking the output from STEP 1\n",
    "3. STEP 3: <br />\n",
    "    Option A: If everything looks fine the group names and short names, Store groups once you know you have files in correct groups by running and checking the output of STEPS 1 and 2.<br />\n",
    "    Option B (optional): If you would like to edit the groups, uncomment the options B-I and B-II. Run Option B-I to export a prefilled tab infosheet. Edit the file and then run Option B-II to import the new groups and save it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#STEP 1: View the groups\n",
    "\n",
    "files = metob.retrieve('lcmsruns',experiment='%ENTERSTRING%',username='*')\n",
    "controlled_vocab = ['QC','InjBl','ISTD'] #add _ to beginning. It will be stripped if at begining\n",
    "version_identifier = 'vs1'\n",
    "exclude_files = []  # Exclude files containing a substring (list) Eg., ['peas']\n",
    "file_dict = {}\n",
    "groups_dict = {}\n",
    "for f in files:\n",
    "    if not any(map(f.name.__contains__, exclude_files)):\n",
    "        k = f.name.split('.')[0]\n",
    "        #     get index if any controlled vocab in filename\n",
    "        indices = [i for i, s in enumerate(controlled_vocab) if s.lower() in k.lower()]\n",
    "        prefix = '_'.join(k.split('_')[:11])\n",
    "        if len(indices)>0:\n",
    "            short_name = controlled_vocab[indices[0]].lstrip('_')\n",
    "            group_name = '%s_%s_%s'%(prefix,version_identifier,short_name)\n",
    "            short_name = k.split('_')[9]+'_'+short_name # Prepending POL to short_name\n",
    "        else:\n",
    "            short_name = k.split('_')[12]\n",
    "            group_name = '%s_%s_%s'%(prefix,version_identifier,short_name)\n",
    "            short_name = k.split('_')[9]+'_'+k.split('_')[12]  # Prepending POL to short_name\n",
    "        file_dict[k] = {'file':f,'group':group_name,'short_name':short_name}\n",
    "        groups_dict[group_name] = {'items':[],'name':group_name,'short_name':short_name}\n",
    "df = pd.DataFrame(file_dict).T\n",
    "df.index.name = 'filename'\n",
    "df.reset_index(inplace=True)#['group'].unique()\n",
    "df.drop(columns=['file'],inplace=True)\n",
    "for ug in groups_dict.keys():\n",
    "    for file_key,file_value in file_dict.items():\n",
    "        if file_value['group'] == ug:\n",
    "            groups_dict[ug]['items'].append(file_value['file'])\n",
    "df.head(100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#STEP 2: create the groups variable, if the above looks OK\n",
    "\n",
    "groups = []\n",
    "for group_key,group_values in groups_dict.items():\n",
    "    g = metob.Group(name=group_key,items=group_values['items'],short_name=group_values['short_name'])\n",
    "    groups.append(g)        \n",
    "    for item in g.items:\n",
    "        print(g.name,g.short_name,item.name)\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# STEP 3 Option A: store the groups variable content in the DB (currently only the long group name is stored)\n",
    "metob.store(groups)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## STEP 3 Option B-I: OPTIONAL: Export groups to csv file for editing (filename, short_name, group, description)\n",
    "#dp.make_prefilled_fileinfo_sheet(groups,os.path.join(output_dir,'prefilled_fileinfo.tab'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## STEP 3 Option B-II: Import groups from csv file after editing the prefilled_fileinfo.tab\n",
    "#groups = dp.make_groups_from_fileinfo_sheet(os.path.join(output_dir,'prefilled_fileinfo.tab'), filetype='tab', store=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# OPTION B: Register LCMS Runs into categorical groups from a file.\n",
    "\n",
    "Typically, you will make one fileinfo sheet with all of your files (pos and neg) for this experiment.  At a minimum, group names MUST contain the first 11 underscore delimited fields (DATE_NORTHENLABINITIALS_COLLABINITIALS_PROJ_EXP_SAMPSET_SYSTEM_COLUMN-method_SERIAL_POL_ACQ) and the 'SAMPLEGROUP' field.\n",
    "\n",
    "Files can be from multiple folders at nersc.\n",
    "\n",
    "1. STEP 1: select files\n",
    "    1. Edit the experiment and name fields to find the files you want.\n",
    "2. STEP 2: create and save a .tab file to your project directory.\n",
    "    1. After running the block, find the .tab file in your project directory.\n",
    "    2. Open in excel or other spreadsheet editor.\n",
    "    3. Fill out the group names as per above in the editor.\n",
    "    4. Save the file as filled_fileinfo.txt\n",
    "3. STEP 3: Create groups from spreadsheet\n",
    "    1. Transfer the .txt. file back to your project directory \n",
    "    2. Run the make groups block using store=False\n",
    "4. STEP 4: CHECK groups\n",
    "    1. Run the next block 'metob.to_dataframe(g) and check that the information looks correct\n",
    "5. If it is correct, rerun the STEP 2 make groups block, using store=True.  If not, fix your file in excel and redo Steps 2&3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#STEP 1: Select files\n",
    "files = dp.get_metatlas_files(experiment = '%ENTERSTRING%',name = '%ENTERSTRING%',most_recent = True)\n",
    "# ^ edit the text string in experiment and name fields"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#STEP 2: Save spreadsheet file\n",
    "dp.make_empty_fileinfo_sheet('%s%s' % (output_dir,'empty_fileinfo.tab'),files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#STEP 3: create groups from file\n",
    "g = dp.make_groups_from_fileinfo_sheet('%s%s' % (output_dir,'filled_fileinfo.txt'),\n",
    "                                       filetype='tab',\n",
    "                                       store=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# STEP 4: check groups\n",
    "metob.to_dataframe(g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make data frame of short filenames and samplenames\n",
    "Uncomment the below 2 blocks to make short file names and smaple names.<br>\n",
    "This creates a dataframe and a csv file which can be edited, exported and imported. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make short_filename and short_samplename \n",
    "files = metob.retrieve('lcmsruns',experiment='%ENTERSTRING%',username='*')\n",
    "short_filename_delim_ids = [0,2,4,5,7,9,14]\n",
    "short_samplename_delim_ids = [9,12,13,14]\n",
    "short_names_df = pd.DataFrame(columns=['sample_treatment','short_filename','short_samplename'])\n",
    "ctr = 0\n",
    "for f in files:\n",
    "    short_filename = []\n",
    "    short_samplename = []\n",
    "    tokens = f.name.split('.')[0].split('_')\n",
    "    for id in short_filename_delim_ids:\n",
    "        short_filename.append(str(tokens[id]))\n",
    "    for id in short_samplename_delim_ids:\n",
    "        short_samplename.append(str(tokens[id]))\n",
    "    short_filename = \"_\".join(short_filename)\n",
    "    short_samplename = \"_\".join(short_samplename)\n",
    "    short_names_df.loc[ctr, 'full_filename'] = f.name.split('.')[0]\n",
    "    short_names_df.loc[ctr, 'sample_treatment'] = str(tokens[12]) # delim 12\n",
    "    short_names_df.loc[ctr, 'short_filename'] = short_filename\n",
    "    short_names_df.loc[ctr, 'short_samplename'] = short_samplename\n",
    "    short_names_df.loc[ctr, 'last_modified'] = pd.to_datetime(f.last_modified,unit='s')\n",
    "    ctr +=1\n",
    "short_names_df.sort_values(by='last_modified', inplace=True)\n",
    "short_names_df.drop(columns=['last_modified'], inplace=True)\n",
    "short_names_df.drop_duplicates(subset=['full_filename'], keep='last', inplace=True)\n",
    "short_names_df.set_index('full_filename', inplace=True)\n",
    "short_names_df.to_csv(os.path.join(output_dir, 'short_names.csv'), sep=',', index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optional import edited short_names.csv \n",
    "short_names_df = pd.read_csv(os.path.join(output_dir, 'short_names.csv'), sep=',', index_col='full_filename')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Select groups of files to operate on\n",
    "\n",
    "Here, you will assign your database groups to a local variable which will be used downstream in the notebook for analyzing your data with an atlas.\n",
    "\n",
    "1. in block below, fill out the fields for name, include_list and exclude_list using text strings from the group names you created in the previous step.  The include/exlcude lists do not need wildcards.  Name is a string unique to all of your groups (ex. fields 0-11 of your filenames)\n",
    "\n",
    "### Typically, you will run one polarity at a time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "polarity = 'POS'  #IMPORTANT: Please make sure you set the correct polarity for the analysis\n",
    "\n",
    "groups = dp.select_groups_for_analysis(name = '%ENTERSEARCHSTRING%',    # <- edit text search string here\n",
    "                                       most_recent = True,\n",
    "                                       remove_empty = True,\n",
    "                                       include_list = [], exclude_list = ['NEG','QC','InjBl'])# ex. ['QC','Blank'])\n",
    "print(\"sorted groups\")\n",
    "groups = sorted(groups, key=operator.attrgetter('name'))\n",
    "for i,a in enumerate(groups):\n",
    "    print(i, a.name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to view metadata about your groups, run the block below\n",
    "metob.to_dataframe(groups)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Create new Atlas entries in the Metatlas database from a csv file\n",
    "\n",
    "## QC, IS, and EMA template atlases are available on the google drive.\n",
    "\n",
    "1. Create your atlas as a csv file, check that it looks correct (has all the correct headers and no blank values in rows; all columns are the correct data type\n",
    "2. Save it with the type of atlas (EMA, QC or IS), your initials, the experiment name, the polarity, and the version or timestamp\n",
    "3. Upload it to your nersc project directory (the one you named above).  (If it doesn't work, double check your file permissions are set to at least rw-rw----).\n",
    "4. Run blocks below to create the DB entries for negative and positive mode atlases\n",
    "5. WARNING: Don't run this block over and over again - it will create multiple new DB entries with the same atlas name\n",
    "\n",
    "Required Atlas file headers:\n",
    "\n",
    "inchi_key,label,rt_min,rt_max,rt_peak,mz,mz_tolerance,adduct,polarity,identification_notes\n",
    "\n",
    "values in rows must be completed for all fields except inchi_key (leaving this blank will not allow you to perform MSMS matching below), and identification notes\n",
    "\n",
    "INFO: store=True will register your atlas in the database.  If you are not sure if your atlas structure is correct, set store=False for the first time your run the block to check if you get an error.  If there is no error, then rerun it with store=True."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NEGATIVE MODE ATLAS UPLOAD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "atlasfilename='%ENTERSTRING%'  # <- enter the exact name of your csv file without the file extension\n",
    "\n",
    "names = dp.make_atlas_from_spreadsheet('%s%s%s' % (pathtoatlas,atlasfilename,'.csv'),  # <- DO NOT EDIT THIS LINE\n",
    "                                       atlasfilename,\n",
    "                                       filetype='csv',\n",
    "                                       sheetname='',\n",
    "                                       polarity = 'negative',\n",
    "                                       store=True,\n",
    "                                       mz_tolerance = 12\n",
    "                                      )   \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# POSITIVE MODE ATLAS UPLOAD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "atlasfilename='%ENTERSTRING%'   # <- enter the exact name of your csv file without the file extension\n",
    "\n",
    "names = dp.make_atlas_from_spreadsheet('%s%s%s' % (pathtoatlas,atlasfilename,'.csv'),  # <- DO NOT EDIT THIS LINE\n",
    "                                       atlasfilename,\n",
    "                                       filetype='csv',\n",
    "                                       sheetname='',\n",
    "                                       polarity = 'positive',\n",
    "                                       store=True,\n",
    "                                       mz_tolerance = 12\n",
    "                                      )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Select Atlas to use\n",
    "\n",
    "1. The first block will retrieve a list of atlases matching the 'name' string that you enter.  Also, you must enter your username.\n",
    "2. The next block will select one from the list, using the index number.  Make sure to enter the index number for the atlas you want to use for your analysis by setting in this line: my_atlas = atlases[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "atlases = metob.retrieve('Atlas',name='%ENTERSTRING%',username='YOUR-NERSC-USERNAME')\n",
    "names = []\n",
    "for i,a in enumerate(atlases):\n",
    "    print(i,a.name,pd.to_datetime(a.last_modified,unit='s'))#len(a.compound_identifications)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_atlas = atlases[-1]\n",
    "atlas_df = ma_data.make_atlas_df(my_atlas)\n",
    "atlas_df['label'] = [cid.name for cid in my_atlas.compound_identifications]\n",
    "print(my_atlas.name)\n",
    "metob.to_dataframe([my_atlas])\n",
    "# the first line of the output will show the dimensions of the atlas dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# OPTIONAL: to view your atlas, run this block\n",
    "print(my_atlas.name)\n",
    "atlas_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. Get EICs and MSMS for all files in your groups, using all compounds in your atlas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### This block builds the metatlas_dataset variable.  This holds your EIC data (mz, rt, intensity values within your mz and rt ranges).\n",
    "\n",
    "The EIC data contains mz, intensity and RT values across your RT range.  There are two parameters that you will need to edit: extra_time and extra_mz.  Extra time will collect mz, intensity and RT values from outside of your atlas defined min and max rt values.  For example if your rt_min = 1.0, and rt_max = 2.0 and you set extra_time to 0.3, then your new rt range will be 0.7 to 2.3.  This is helpful for checking if you have nearby peaks at the same m/z.  Extra_mz should only be used for troubleshooting.  You should keep this at 0 unless you believe you have poor mass accuracy during your run.  Other ways to address this issue is by changing the mz_tolerance values in your atlas.  Before changing this value, you should check in with a metatlas experienced lab member to discuss when/how to use this value.\n",
    "\n",
    "1. Change the value in \"extra_time = 0.0\" to something like 0.5 to 1.0 for the first EMA runthrough on your files.  This will take longer but collect msms outside your retention windows which allows you to check the msms of nearby peaks before adjusting your rt bounds around the correct peak.\n",
    "2. extra_mz should almost always be set to 0.0   If you need to troubleshoot a low mz compound you could potentially use this value to run it back through with a larger mz error window than what was specified in your atlas (ppm tolerance).\n",
    "\n",
    ">On Your final runthrough, set extra_time to 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_files = []\n",
    "for my_group in groups:\n",
    "    for my_file in my_group.items:\n",
    "        extra_time = 0.75          # NOTE: 0.75 for the first run, 0.5 for final \n",
    "        extra_mz = 0.00\n",
    "        all_files.append((my_file,my_group,atlas_df,my_atlas,extra_time,extra_mz))\n",
    "pool = mp.Pool(processes=min(4, len(all_files)))\n",
    "t0 = time.time()\n",
    "metatlas_dataset = pool.map(ma_data.get_data_for_atlas_df_and_file, all_files)\n",
    "pool.close()\n",
    "pool.terminate()\n",
    "print(time.time() - t0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make data sources tables (atlas_metadata.tab, groups_metadata.tab, groups.tab and [atlasname]_originalatlas.tab within data_sources subfolder)\n",
    "ma_data.make_data_sources_tables(groups, my_atlas, output_dir) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6b Optional: Filter atlas for compounds with no or low signals\n",
    "\n",
    "Uncomment the below 3 blocks to filter the atlas.\n",
    "Please ensure that correct polarity is used for the atlases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dp = reload(dp)\n",
    "# num_data_points_passing = 5\n",
    "# peak_height_passing = 4e5\n",
    "# atlas_df_passing = dp.filter_atlas(atlas_df=atlas_df, input_dataset=metatlas_dataset, num_data_points_passing = num_data_points_passing, peak_height_passing = peak_height_passing)\n",
    "# print(\"# Compounds in Atlas: \"+str(len(atlas_df)))\n",
    "# print(\"# Compounds passing filter: \"+str(len(atlas_df_passing)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create new atlas and store in database\n",
    "This block creates a filtered atlas with a new name !!\n",
    "Automatically selects this atlas for processing. \n",
    "Make sure to use this atlas for downstream analyses. (NOTE: If you restart kernel or come back to the analysis, you need to reselect this newly created filtered atlas for processing)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# atlas_passing = my_atlas.name+'_filteredby-datapnts'+str(num_data_points_passing)+'-pkht'+str(peak_height_passing)\n",
    "# myAtlas_passing = dp.make_atlas_from_spreadsheet(atlas_df_passing,\n",
    "#                           atlas_passing,\n",
    "#                           filetype='dataframe',\n",
    "#                           sheetname='',\n",
    "#                           polarity = 'positive',\n",
    "#                           store=True,\n",
    "#                           mz_tolerance = 12)\n",
    "\n",
    "# atlases = dp.get_metatlas_atlas(name=atlas_passing,do_print = True, most_recent=True)\n",
    "\n",
    "# myAtlas = atlases[-1]\n",
    "# atlas_df = ma_data.make_atlas_df(myAtlas)\n",
    "# atlas_df['label'] = [cid.name for cid in myAtlas.compound_identifications]\n",
    "# print(myAtlas.name)\n",
    "# print(myAtlas.username)\n",
    "# metob.to_dataframe([myAtlas])# "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# all_files = []\n",
    "# for my_group in groups:\n",
    "#     for my_file in my_group.items:\n",
    "#         all_files.append((my_file,my_group,atlas_df,myAtlas))\n",
    "        \n",
    "# pool = mp.Pool(processes=min(4, len(all_files)))\n",
    "# t0 = time.time()\n",
    "# metatlas_dataset = pool.map(ma_data.get_data_for_atlas_df_and_file, all_files)\n",
    "# pool.close()\n",
    "# pool.terminate()\n",
    "# #If you're code crashes here, make sure to terminate any processes left open.\n",
    "#(print time.time() - t0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### One of the two blocks below builds the hits variable.  This holds your MSMS spectra (from within your mz, and rt ranges, and within the extra time indicated above).\n",
    "\n",
    "There are two options for generating the hits variable:\n",
    "1. block A: use when your files have msms. It create the hits variable and also saves a binary (pickled) serialized hits file to the output directory.\n",
    "2. block B: only run if your files were collected in MS1 mode\n",
    "3. If you have already run block A and then the kernel dies, you can skip block A and directly unplickle the binary hits file from the output directory. Skip block A, uncomment the Optional block and run it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##BLOCK A\n",
    "import warnings; warnings.simplefilter('ignore')\n",
    "t0 = time.time()\n",
    "\n",
    "hits=dp.get_msms_hits(metatlas_dataset,extra_time=True,keep_nonmatches=True, frag_mz_tolerance=0.01, ref_loc='/global/project/projectdirs/metatlas/projects/spectral_libraries/msms_refs_v3.tab')\n",
    "pickle.dump(hits, open(os.path.join(output_dir,polarity+'_hits.pkl'), \"wb\"))\n",
    "\n",
    "print(time.time() - t0)\n",
    "print('%s%s' % (len(hits),' <- total number of MSMS spectra found in your files'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## BLOCK B (uncomment lines below to run this.  Only use when all data files are MS1)\n",
    "#hits=pd.DataFrame([], columns=['database','id','file_name','msms_scan', u'score', u'num_matches', u'msv_query_aligned', u'msv_ref_aligned', u'name', u'adduct', u'inchi_key', u'precursor_mz', u'measured_precursor_mz'])\n",
    "#hits.set_index(['database','id','file_name','msms_scan'], inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optional: If you already have a pickled hits file and do not need to run get_msms_hits again, uncomment this block\n",
    "# hits = pickle.load(open(os.path.join(output_dir,polarity+'_hits.pkl'), \"rb\")) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. Adjust Retention Times. \n",
    "\n",
    "This block creates an interactive plot.  The top panel displays MSMS from within the two green RT bounds selected below (rt min and max, initially set in atlas).  When the database holds reference spectra, mirror plots are generated with the reference spectra inverted below the sample spectra.  The lower panel displays the EICs overlayed for all of the files in your selected groups.  You can highlight your groups different colors.  It is recommended that you do this, at least, for your extraction blank (or if not available, use a solvent injection blank).  This plot also displays radio buttons that can be interactively selected; the values will be exported in your final identifications table and in your atlas export.  Use these to mark peak/MSMS quality.\n",
    "\n",
    "How to use:\n",
    "1. STEP 1: Set peak flag radio buttons\n",
    "    1. OPTION A (custom flags): fill out the peak flags list (list of strings) \n",
    "        peak_flag_list = ('A','B') some recommendations are below.  \n",
    "    2. OPTION B (default flags): comment out the custom peak_flag_list line.  Uncomment the default peak_flags = \"\". \n",
    "        Flags default to: keep, remove, unresolvable isomers, check.\n",
    "2. STEP 2: Set EIC colors\n",
    "    1. Option A (custom EIC colors): fill out the colorlist in the format of below\n",
    "   > ***\n",
    "   > colorlist = [['color1nameorhexadec','partialgroupstring1'],\n",
    "   >         ['color2nameorhexadec','partialgroupstring2']]\n",
    "   > ***\n",
    "    <ul><li>You can add more comma delimited colors/groups as needed.</li>\n",
    "    <li>These are partial strings that match to you file names (not your group names).</li>\n",
    "    <li>The order they are listed in your list is the order they are displayed in the overlays (first is front, last is back)</li>\n",
    "    <li>Named colors available in matplotlib are here: https://matplotlib.org/3.1.0/gallery/color/named_colors.html\n",
    "          or use hexadecimal values '#000000'</li></ul>\n",
    "    B. Option B (default EIC colors): comment out the custom colorlist lines and uncomment the default colorlist = \"\". \n",
    "        Colors all default to black.\n",
    "3. User the right/left buttons on your keyboard to cycle through compounds in your atlas.\n",
    "4. Use the up/down buttons on your keyboard to cycle through MSMS spectra within the RT bounds of the lower plot.\n",
    "5. Use the horizontal rt min and rt max bars below the plots to adjust the rt bounds around your peak.  If there are multiple peaks, select one at a time and then click up/down to update the msms available in that new RT range.  If necessary evaluate your data in an external program such as mzmine to make sure you are selecting the correct peak.\n",
    "\n",
    "TIPS: use compound_idx = 0 in step 3 to change to a different compound in your atlas using the index number.  If your plot does not fit in your browser window, adjust height and width values.  Use alpha to change the transparency of the lines this is a value 0 (transparent) to 1 (opaque).\n",
    "\n",
    "DO NOT change your RT theoretical peak (the purple line).  It is locked from editing (unless you change a hidden parameter) and only to be changed in special cases.  The measured retention times of your peaks will be calculated and exported in your output files.  These will be compared with the RT theoreticals and used in your evidence of identification table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###STEP 1: Set the peak flag radio buttons using one of the two lines below, for custom flags or default flags\n",
    "import warnings; warnings.simplefilter('ignore')\n",
    "peak_flag_list=('','L1+ - 1 pk, good RT&MSMS','L1+ - known isomer overlap','L1+ - 1 pk, good RT, MSMS ok (coisolated mz/partial match/low int)','L1+ - 1 pk, good RT&MSMS from external library','L1 - 1 pk, correct RT, no MSMS or int too low for matching','L1 - 1 pk, good RT, very low intensity/poor pk shape','L2 put comp','L3 putative class','Remove - background/noise','Remove - bad EMA MSMS','Remove - bad MSMS NIST/MONA/Metlin')\n",
    "msms_flags_list = \"\" #\n",
    "#peak_flag_list =\"\"       # this will default to ('keep','remove','unresolvable isomers','poor peak shape')\n",
    "\n",
    "###STEP 2: Set the EIC line colors using on of the two lines below, for custom colors or default \n",
    "colorlist= [['red','peas'],                                                     \n",
    "           ['green','beets']]\n",
    "#colorlist=\"\"   # this will default to black\n",
    "\n",
    "###STEP 3\n",
    "a = dp.adjust_rt_for_selected_compound(metatlas_dataset, msms_hits=hits, peak_flags=peak_flag_list, msms_flags=msms_flags_list, color_me = colorlist, compound_idx=0,alpha=0.5,width=15,height=4.5)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Create filtered atlas excluding compounds marked removed\n",
    "\n",
    "Re-run the following before filtering atlas\n",
    "1. Get Groups (include InjBl)\n",
    "2. Get Atlas\n",
    "3. Get Data\n",
    "4. Get MSMS Hits"<commit_msg>Vs dev (#281)<commit_after>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Instructions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MetAltas analysis workflow\n",
    "\n",
    "This notebook is for performing a targeted metabolomics analysis on your raw data files.  In order for this to work, you data files must be uploaded to the nersc raw data directory with appropriate permissions.  A fileconverter runs at nersc to convert as follows:\n",
    ".raw > .mzml > .h5 > pactolus.gz  This conversion must happen before you can proceed with data analysis.  A Northen Lab staff member can help you transfer and convert your files if they are not already at nersc.\n",
    "***\n",
    "/global/project/projectdirs/metatlas/raw_data/SUBFOLDER\n",
    "***\n",
    "\n",
    "A targeted analysis requires a list of compounds to search for in your sample files.  Thus, at a minimum, you must have some information on their m/z ratio in order to get started.  The inputs for this notebook, include the location of your raw files, and an atlas listing out these compounds.  \n",
    "\n",
    "Our lab has analyzed thousands of compounds using standardized LCMS methods.  You may use one of these libraries to create an atlas for your sample analysis, assuming you used one of the same standard LCMS methods.\n",
    "\n",
    "Analysis steps\n",
    "\n",
    "<ol>\n",
    "<li>Create your atlas csv external to jupyter: typically you will either have a custom atlas or generate an atlas from one of our EMA libraries.  These can be found in the shared google drive, under atlases.</li>\n",
    "<li>Run through with your initial atlas</li>\n",
    "<li>Use the RT adjuster to mark the quality of the peaks and MSMS matches, and also adjust your RT bounds</li>\n",
    "<li>Reload your atlas from the \"get atlas\" block</li>\n",
    "<li>Export your atlas, filter out the compounds marked under id notes as remove.  Save and upload your new atlas</li>\n",
    "<li>Repeat the analysis using your filtered atlas.</li>\n",
    "    <li>On your final export, make sure you are excluding system blanks and QC files in your output.  If you have extraction blanks you should include those.</li>\n",
    "    </ol>\n",
    "    \n",
    "There is a database that will store your atlas and registered run files.  During the analysis process you will pull these as local variables. These are used to pull out EIC and MSMS data which is stored in a variable called metatlas_dataset.  All output files are generated from that.  It is a local variable so is not stored in the db - thus if your kernel dies partway through, you need to rerun through the notebook to regenerate that variable.\n",
    "\n",
    "IMPORTANT: Anytime you adjust your RTs in the interactive plot, you need to retreive your atlas again before exporting any files.  The changes are stored in the db and the exports are made from locally stored variables/dataframes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# index\n",
    "1. Import python packages\n",
    "2. Create groups\n",
    "3. Select groups\n",
    "4. Create Atlases\n",
    "5. Select Atlas\n",
    "6. Annotate data\n",
    "7. Correct RT bounds\n",
    "8. Export results\n",
    "9. Tools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Helpful tips\n",
    "\n",
    "1. When entering search strings, use % for wildcard: %peas will return \"positive_peas\" and \"positive_greenpeas\", but not \"positive_greenpeas_HILIC\", for that use %peas%, you can also put multiple strings via %pos%peas%\n",
    "2. Set your kernel to mass_spec_cori.  If not sure how, check w Ben or Daniel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# optional: run this block to adjust the width of the notebook.  Change the width percent.\n",
    "\n",
    "from  IPython.core.display  import  display, HTML \n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 1. Import Python Packages\n",
    "## Instructions: \n",
    "\n",
    "1. On the 3rd line of the block below, add your directory where your most recent metatlas code is stored (to the sys.path.insert line).\n",
    "2. Only run this block once to load the metatlas python modules.  If the kernel dies, you will need to rerun this block.\n",
    "3. Run the block with the \"print fa.__file__\" code line to double check that the printed output matches the directory in the sys.path.insert line of the prior block.\n",
    "    If it does not match, then on the \"kernel\" dropdown, click restart and correct the path below before trying again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import sys, os\n",
    "\n",
    "# v edit this line\n",
    "sys.path.insert(0,'/global/homes/FIRST-INITIAL-OF-USERNAME/USERNAME/REPOFOLDER/metatlas-master-20200416/metatlas-master/')\n",
    "# ^ edit this line\n",
    "\n",
    "from metatlas.tools import fastanalysis as fa\n",
    "from metatlas.plots import dill2plots as dp\n",
    "from metatlas.io import metatlas_get_data_helper_fun as ma_data\n",
    "from metatlas.plots import chromplotplus as cpp\n",
    "from metatlas.datastructures import metatlas_objects as metob\n",
    "import qgrid\n",
    "from ipywidgets import interact, interactive, fixed\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display\n",
    "import time\n",
    "import pickle\n",
    "import dill\n",
    "import multiprocessing as mp\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plot\n",
    "import operator\n",
    "from importlib import reload\n",
    "\n",
    "pd.set_option('display.max_rows', 5000)\n",
    "pd.set_option('display.max_columns', 500)\n",
    "pd.set_option('display.max_colwidth', 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check that first block registered metatlas directory correctly.  If it does not look the same, seek help!\n",
    "print(fa.__file__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set atlas, project and output directories from your nersc home directory\n",
    "\n",
    "STEP 1\n",
    "1. Create a project folder name for this analysis by replacing the PROJECTDIRECTORY string text in red below.  Make sure to update the rest of the direcory to point to your home directory.  The pwd block will print out the directory where this jupyter notebook is stored.\n",
    "2. Create a subdirectory name for the output, each run through you may want to create a new output folder.\n",
    "3. When you run the block the folders will be created in your home directory.  If the directory already exists, the block will just set the path for use with future code blocks.\n",
    "\n",
    "STEP 2\n",
    "1.  Enter the nersc path where you have stored your atlas files using one of the two optional lines below:\n",
    "    1. Set a custom path\n",
    "    2. Use the project and output subdirectory from above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#STEP 1\n",
    "project_directory='/global/homes/FIRST-INITIAL-OF-USERNAME/USERNAME/PROJECTDIRECTORY/'  # <- edit this line, do not copy the path directly from NERSC (ex. the u1, or u2 directories)\n",
    "output_subfolder='HILIC_POS_20190830/'  # <- edit this as 'chromatography_polarity_yyyymmdd/'\n",
    "output_dir = os.path.join(project_directory,output_subfolder)\n",
    "\n",
    "if not os.path.exists(project_directory):\n",
    "    os.makedirs(project_directory)\n",
    "if not os.path.exists(output_dir):\n",
    "    os.makedirs(output_dir)\n",
    "\n",
    "#STEP 2\n",
    "pathtoatlas='/global/homes/FIRST-INITIAL-OF-USERNAME/USERNAME/DIRECTORYOFATLAS/'    # <- enter the directory where you have stored your atlas\n",
    "#pathtoatlas= '%s%s' % (project_directory,output_subfolder)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Create Groups (named variables that hold your replicates of each sample)\n",
    "\n",
    "### You must assign your raw files into experimental groups for analysis.  These are used for downstream statistics and for selection of specific groups for filtering to subsets of files for analysis (Ex. just pos or just neg).\n",
    "\n",
    "The groups are created from common file headers and the unique group names. The convention our lab group uses for filenames is as follows: \n",
    "***\n",
    "DATE_NORTHENLABINITIALS_COLLABINITIALS_PROJ_EXP_SAMPSET_SYSTEM_COLUMN-method_SERIAL_POL_ACQ_SAMPLENUMBER_ SAMPLEGROUP_REP_OPTIONAL_SEQ \n",
    "\n",
    "Ex.:20180105_SK_AD_ENIGMA_PseudoInt_R2ADec2017_QE119_50454_123456_POS_MSMS_001_Psyringae-R2A-30C-20hr_Rep01_NA_Seq001.raw\n",
    "***\n",
    "The common header consists of the fields 0-10: DATE_NORTHENLABINITIALS_COLLABINITIALS_PROJ_EXP_SAMPSET_SYSTEM_COLUMN-method_SERIAL_POL_ACQ \n",
    "\n",
    "The sample group name is commonly field # 12 (between underscore 11 and 12) -0 indexed-\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Find your files\n",
    "1. On the first line of the block below, set the 'experiment' and 'name' variables to find your files.  These fields require wildcards for partial string searches\n",
    "2. 'Experiment' is the folder name within global/project/projectdirs/metatlas/raw_data, that will be emailed to you when the files are uploaded to NERSC.  You can also look in the raw_data directory for the NERSC user who uploaded your files; your experiment folder should be in there.\n",
    "3. 'name' is string that will match a subset of your files within that folder.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "files = dp.get_metatlas_files(experiment = '%ENTERSTRING%',name = '%ENTERSTRING%',most_recent = True)\n",
    "# ^ edit the text string in experiment and name fields\n",
    "\n",
    "df = metob.to_dataframe(files)\n",
    "df[['experiment','name','username','acquisition_time']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(files)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# OPTION A: Automated Group Maker"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This will attempt to create groups in an automated fashion (rather than filling out a spreadsheet with a list of files and group names).  If your files are all in one folder at nersc, you can use this options.  If not, use option B below.\n",
    "\n",
    "A long group name consisting of the common header + either controlled vocab value or field #12 along with a short group name (just controlled vocab or field #12) will be stored in a local variable.  The short group names can be used on plots.\n",
    "\n",
    "\n",
    "1. STEP 1: View the groups\n",
    "    1. Pick an experiment folder to look for files in on the metob.retrieve function\n",
    "    2. Enter controlled vocabulary for control files to put select files into groups when control string may be in a different field (not #12) or as a randomly placed substring within a field (ex. if 'InjBl' is included in your controlled vocab list, files like _InjBl-MeOH_ and _StartInjBl_ will group together)\n",
    "    3. If your group name is not between _ 11 and 12 you can adjust those values in the split commands below.  All other (non-controlledvocab) groups will be created from that field.\n",
    "2. STEP 2: Create the groups variable after checking the output from STEP 1\n",
    "3. STEP 3: <br />\n",
    "    Option A: If everything looks fine the group names and short names, Store groups once you know you have files in correct groups by running and checking the output of STEPS 1 and 2.<br />\n",
    "    Option B (optional): If you would like to edit the groups, uncomment the options B-I and B-II. Run Option B-I to export a prefilled tab infosheet. Edit the file and then run Option B-II to import the new groups and save it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#STEP 1: View the groups\n",
    "\n",
    "files = metob.retrieve('lcmsruns',experiment='%ENTERSTRING%',username='*')\n",
    "controlled_vocab = ['QC','InjBl','ISTD'] #add _ to beginning. It will be stripped if at begining\n",
    "version_identifier = 'vs1'\n",
    "exclude_files = []  # Exclude files containing a substring (list) Eg., ['peas']\n",
    "file_dict = {}\n",
    "groups_dict = {}\n",
    "for f in files:\n",
    "    if not any(map(f.name.__contains__, exclude_files)):\n",
    "        k = f.name.split('.')[0]\n",
    "        #     get index if any controlled vocab in filename\n",
    "        indices = [i for i, s in enumerate(controlled_vocab) if s.lower() in k.lower()]\n",
    "        prefix = '_'.join(k.split('_')[:11])\n",
    "        if len(indices)>0:\n",
    "            short_name = controlled_vocab[indices[0]].lstrip('_')\n",
    "            group_name = '%s_%s_%s'%(prefix,version_identifier,short_name)\n",
    "            short_name = k.split('_')[9]+'_'+short_name # Prepending POL to short_name\n",
    "        else:\n",
    "            short_name = k.split('_')[12]\n",
    "            group_name = '%s_%s_%s'%(prefix,version_identifier,short_name)\n",
    "            short_name = k.split('_')[9]+'_'+k.split('_')[12]  # Prepending POL to short_name\n",
    "        file_dict[k] = {'file':f,'group':group_name,'short_name':short_name}\n",
    "        groups_dict[group_name] = {'items':[],'name':group_name,'short_name':short_name}\n",
    "df = pd.DataFrame(file_dict).T\n",
    "df.index.name = 'filename'\n",
    "df.reset_index(inplace=True)#['group'].unique()\n",
    "df.drop(columns=['file'],inplace=True)\n",
    "for ug in groups_dict.keys():\n",
    "    for file_key,file_value in file_dict.items():\n",
    "        if file_value['group'] == ug:\n",
    "            groups_dict[ug]['items'].append(file_value['file'])\n",
    "df.head(100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#STEP 2: create the groups variable, if the above looks OK\n",
    "\n",
    "groups = []\n",
    "for group_key,group_values in groups_dict.items():\n",
    "    g = metob.Group(name=group_key,items=group_values['items'],short_name=group_values['short_name'])\n",
    "    groups.append(g)        \n",
    "    for item in g.items:\n",
    "        print(g.name,g.short_name,item.name)\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# STEP 3 Option A: store the groups variable content in the DB (currently only the long group name is stored)\n",
    "metob.store(groups)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## STEP 3 Option B-I: OPTIONAL: Export groups to csv file for editing (filename, short_name, group, description)\n",
    "#dp.make_prefilled_fileinfo_sheet(groups,os.path.join(output_dir,'prefilled_fileinfo.tab'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## STEP 3 Option B-II: Import groups from csv file after editing the prefilled_fileinfo.tab\n",
    "#groups = dp.make_groups_from_fileinfo_sheet(os.path.join(output_dir,'prefilled_fileinfo.tab'), filetype='tab', store=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# OPTION B: Register LCMS Runs into categorical groups from a file.\n",
    "\n",
    "Typically, you will make one fileinfo sheet with all of your files (pos and neg) for this experiment.  At a minimum, group names MUST contain the first 11 underscore delimited fields (DATE_NORTHENLABINITIALS_COLLABINITIALS_PROJ_EXP_SAMPSET_SYSTEM_COLUMN-method_SERIAL_POL_ACQ) and the 'SAMPLEGROUP' field.\n",
    "\n",
    "Files can be from multiple folders at nersc.\n",
    "\n",
    "1. STEP 1: select files\n",
    "    1. Edit the experiment and name fields to find the files you want.\n",
    "2. STEP 2: create and save a .tab file to your project directory.\n",
    "    1. After running the block, find the .tab file in your project directory.\n",
    "    2. Open in excel or other spreadsheet editor.\n",
    "    3. Fill out the group names as per above in the editor.\n",
    "    4. Save the file as filled_fileinfo.txt\n",
    "3. STEP 3: Create groups from spreadsheet\n",
    "    1. Transfer the .txt. file back to your project directory \n",
    "    2. Run the make groups block using store=False\n",
    "4. STEP 4: CHECK groups\n",
    "    1. Run the next block 'metob.to_dataframe(g) and check that the information looks correct\n",
    "5. If it is correct, rerun the STEP 2 make groups block, using store=True.  If not, fix your file in excel and redo Steps 2&3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#STEP 1: Select files\n",
    "files = dp.get_metatlas_files(experiment = '%ENTERSTRING%',name = '%ENTERSTRING%',most_recent = True)\n",
    "# ^ edit the text string in experiment and name fields"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#STEP 2: Save spreadsheet file\n",
    "dp.make_empty_fileinfo_sheet('%s%s' % (output_dir,'empty_fileinfo.tab'),files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#STEP 3: create groups from file\n",
    "g = dp.make_groups_from_fileinfo_sheet('%s%s' % (output_dir,'filled_fileinfo.txt'),\n",
    "                                       filetype='tab',\n",
    "                                       store=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# STEP 4: check groups\n",
    "metob.to_dataframe(g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make data frame of short filenames and samplenames\n",
    "Uncomment the below 2 blocks to make short file names and smaple names.<br>\n",
    "This creates a dataframe and a csv file which can be edited, exported and imported. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make short_filename and short_samplename \n",
    "files = metob.retrieve('lcmsruns',experiment='%ENTERSTRING%',username='*')\n",
    "short_filename_delim_ids = [0,2,4,5,7,9,14]\n",
    "short_samplename_delim_ids = [9,12,13,14]\n",
    "short_names_df = pd.DataFrame(columns=['sample_treatment','short_filename','short_samplename'])\n",
    "ctr = 0\n",
    "for f in files:\n",
    "    short_filename = []\n",
    "    short_samplename = []\n",
    "    tokens = f.name.split('.')[0].split('_')\n",
    "    for id in short_filename_delim_ids:\n",
    "        short_filename.append(str(tokens[id]))\n",
    "    for id in short_samplename_delim_ids:\n",
    "        short_samplename.append(str(tokens[id]))\n",
    "    short_filename = \"_\".join(short_filename)\n",
    "    short_samplename = \"_\".join(short_samplename)\n",
    "    short_names_df.loc[ctr, 'full_filename'] = f.name.split('.')[0]\n",
    "    short_names_df.loc[ctr, 'sample_treatment'] = str(tokens[12]) # delim 12\n",
    "    short_names_df.loc[ctr, 'short_filename'] = short_filename\n",
    "    short_names_df.loc[ctr, 'short_samplename'] = short_samplename\n",
    "    short_names_df.loc[ctr, 'last_modified'] = pd.to_datetime(f.last_modified,unit='s')\n",
    "    ctr +=1\n",
    "short_names_df.sort_values(by='last_modified', inplace=True)\n",
    "short_names_df.drop(columns=['last_modified'], inplace=True)\n",
    "short_names_df.drop_duplicates(subset=['full_filename'], keep='last', inplace=True)\n",
    "short_names_df.set_index('full_filename', inplace=True)\n",
    "short_names_df.to_csv(os.path.join(output_dir, 'short_names.csv'), sep=',', index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optional import edited short_names.csv \n",
    "short_names_df = pd.read_csv(os.path.join(output_dir, 'short_names.csv'), sep=',', index_col='full_filename')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Select groups of files to operate on\n",
    "\n",
    "Here, you will assign your database groups to a local variable which will be used downstream in the notebook for analyzing your data with an atlas.\n",
    "\n",
    "1. in block below, fill out the fields for name, include_list and exclude_list using text strings from the group names you created in the previous step.  The include/exlcude lists do not need wildcards.  Name is a string unique to all of your groups (ex. fields 0-11 of your filenames)\n",
    "\n",
    "### Typically, you will run one polarity at a time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "polarity = 'POS'  #IMPORTANT: Please make sure you set the correct polarity for the analysis\n",
    "\n",
    "groups = dp.select_groups_for_analysis(name = '%ENTERSEARCHSTRING%',    # <- edit text search string here\n",
    "                                       most_recent = True,\n",
    "                                       remove_empty = True,\n",
    "                                       include_list = [], exclude_list = ['NEG','QC','InjBl'])# ex. ['QC','Blank'])\n",
    "print(\"sorted groups\")\n",
    "groups = sorted(groups, key=operator.attrgetter('name'))\n",
    "for i,a in enumerate(groups):\n",
    "    print(i, a.name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to view metadata about your groups, run the block below\n",
    "metob.to_dataframe(groups)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Create new Atlas entries in the Metatlas database from a csv file\n",
    "\n",
    "## QC, IS, and EMA template atlases are available on the google drive.\n",
    "\n",
    "1. Create your atlas as a csv file, check that it looks correct (has all the correct headers and no blank values in rows; all columns are the correct data type\n",
    "2. Save it with the type of atlas (EMA, QC or IS), your initials, the experiment name, the polarity, and the version or timestamp\n",
    "3. Upload it to your nersc project directory (the one you named above).  (If it doesn't work, double check your file permissions are set to at least rw-rw----).\n",
    "4. Run blocks below to create the DB entries for negative and positive mode atlases\n",
    "5. WARNING: Don't run this block over and over again - it will create multiple new DB entries with the same atlas name\n",
    "\n",
    "Required Atlas file headers:\n",
    "\n",
    "inchi_key,label,rt_min,rt_max,rt_peak,mz,mz_tolerance,adduct,polarity,identification_notes\n",
    "\n",
    "values in rows must be completed for all fields except inchi_key (leaving this blank will not allow you to perform MSMS matching below), and identification notes\n",
    "\n",
    "INFO: store=True will register your atlas in the database.  If you are not sure if your atlas structure is correct, set store=False for the first time your run the block to check if you get an error.  If there is no error, then rerun it with store=True."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NEGATIVE MODE ATLAS UPLOAD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "atlasfilename='%ENTERSTRING%'  # <- enter the exact name of your csv file without the file extension\n",
    "\n",
    "names = dp.make_atlas_from_spreadsheet('%s%s%s' % (pathtoatlas,atlasfilename,'.csv'),  # <- DO NOT EDIT THIS LINE\n",
    "                                       atlasfilename,\n",
    "                                       filetype='csv',\n",
    "                                       sheetname='',\n",
    "                                       polarity = 'negative',\n",
    "                                       store=True,\n",
    "                                       mz_tolerance = 12\n",
    "                                      )   \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# POSITIVE MODE ATLAS UPLOAD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "atlasfilename='%ENTERSTRING%'   # <- enter the exact name of your csv file without the file extension\n",
    "\n",
    "names = dp.make_atlas_from_spreadsheet('%s%s%s' % (pathtoatlas,atlasfilename,'.csv'),  # <- DO NOT EDIT THIS LINE\n",
    "                                       atlasfilename,\n",
    "                                       filetype='csv',\n",
    "                                       sheetname='',\n",
    "                                       polarity = 'positive',\n",
    "                                       store=True,\n",
    "                                       mz_tolerance = 12\n",
    "                                      )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Select Atlas to use\n",
    "\n",
    "1. The first block will retrieve a list of atlases matching the 'name' string that you enter.  Also, you must enter your username.\n",
    "2. The next block will select one from the list, using the index number.  Make sure to enter the index number for the atlas you want to use for your analysis by setting in this line: my_atlas = atlases[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "atlases = metob.retrieve('Atlas',name='%ENTERSTRING%',username='YOUR-NERSC-USERNAME')\n",
    "names = []\n",
    "for i,a in enumerate(atlases):\n",
    "    print(i,a.name,pd.to_datetime(a.last_modified,unit='s'))#len(a.compound_identifications)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_atlas = atlases[-1]\n",
    "atlas_df = ma_data.make_atlas_df(my_atlas)\n",
    "atlas_df['label'] = [cid.name for cid in my_atlas.compound_identifications]\n",
    "print(my_atlas.name)\n",
    "metob.to_dataframe([my_atlas])\n",
    "# the first line of the output will show the dimensions of the atlas dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# OPTIONAL: to view your atlas, run this block\n",
    "print(my_atlas.name)\n",
    "atlas_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. Get EICs and MSMS for all files in your groups, using all compounds in your atlas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### This block builds the metatlas_dataset variable.  This holds your EIC data (mz, rt, intensity values within your mz and rt ranges).\n",
    "\n",
    "The EIC data contains mz, intensity and RT values across your RT range.  There are two parameters that you will need to edit: extra_time and extra_mz.  Extra time will collect mz, intensity and RT values from outside of your atlas defined min and max rt values.  For example if your rt_min = 1.0, and rt_max = 2.0 and you set extra_time to 0.3, then your new rt range will be 0.7 to 2.3.  This is helpful for checking if you have nearby peaks at the same m/z.  Extra_mz should only be used for troubleshooting.  You should keep this at 0 unless you believe you have poor mass accuracy during your run.  Other ways to address this issue is by changing the mz_tolerance values in your atlas.  Before changing this value, you should check in with a metatlas experienced lab member to discuss when/how to use this value.\n",
    "\n",
    "1. Change the value in \"extra_time = 0.0\" to something like 0.5 to 1.0 for the first EMA runthrough on your files.  This will take longer but collect msms outside your retention windows which allows you to check the msms of nearby peaks before adjusting your rt bounds around the correct peak.\n",
    "2. extra_mz should almost always be set to 0.0   If you need to troubleshoot a low mz compound you could potentially use this value to run it back through with a larger mz error window than what was specified in your atlas (ppm tolerance).\n",
    "\n",
    ">On Your final runthrough, set extra_time to 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_files = []\n",
    "for my_group in groups:\n",
    "    for my_file in my_group.items:\n",
    "        extra_time = 0.75          # NOTE: 0.75 for the first run, 0.5 for final \n",
    "        extra_mz = 0.00\n",
    "        all_files.append((my_file,my_group,atlas_df,my_atlas,extra_time,extra_mz))\n",
    "pool = mp.Pool(processes=min(4, len(all_files)))\n",
    "t0 = time.time()\n",
    "metatlas_dataset = pool.map(ma_data.get_data_for_atlas_df_and_file, all_files)\n",
    "pool.close()\n",
    "pool.terminate()\n",
    "print(time.time() - t0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make data sources tables (atlas_metadata.tab, groups_metadata.tab, groups.tab and [atlasname]_originalatlas.tab within data_sources subfolder)\n",
    "ma_data.make_data_sources_tables(groups, my_atlas, output_dir) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6b Optional: Filter atlas for compounds with no or low signals\n",
    "\n",
    "Uncomment the below 3 blocks to filter the atlas.\n",
    "Please ensure that correct polarity is used for the atlases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dp = reload(dp)\n",
    "# num_data_points_passing = 5\n",
    "# peak_height_passing = 4e5\n",
    "# atlas_df_passing = dp.filter_atlas(atlas_df=atlas_df, input_dataset=metatlas_dataset, num_data_points_passing = num_data_points_passing, peak_height_passing = peak_height_passing)\n",
    "# print(\"# Compounds in Atlas: \"+str(len(atlas_df)))\n",
    "# print(\"# Compounds passing filter: \"+str(len(atlas_df_passing)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create new atlas and store in database\n",
    "This block creates a filtered atlas with a new name !!\n",
    "Automatically selects this atlas for processing. \n",
    "Make sure to use this atlas for downstream analyses. (NOTE: If you restart kernel or come back to the analysis, you need to reselect this newly created filtered atlas for processing)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# atlas_passing = my_atlas.name+'_filteredby-datapnts'+str(num_data_points_passing)+'-pkht'+str(peak_height_passing)\n",
    "# myAtlas_passing = dp.make_atlas_from_spreadsheet(atlas_df_passing,\n",
    "#                           atlas_passing,\n",
    "#                           filetype='dataframe',\n",
    "#                           sheetname='',\n",
    "#                           polarity = 'positive',\n",
    "#                           store=True,\n",
    "#                           mz_tolerance = 12)\n",
    "\n",
    "# atlases = dp.get_metatlas_atlas(name=atlas_passing,do_print = True, most_recent=True)\n",
    "\n",
    "# myAtlas = atlases[-1]\n",
    "# atlas_df = ma_data.make_atlas_df(myAtlas)\n",
    "# atlas_df['label'] = [cid.name for cid in myAtlas.compound_identifications]\n",
    "# print(myAtlas.name)\n",
    "# print(myAtlas.username)\n",
    "# metob.to_dataframe([myAtlas])# "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# all_files = []\n",
    "# for my_group in groups:\n",
    "#     for my_file in my_group.items:\n",
    "#         all_files.append((my_file,my_group,atlas_df,myAtlas))\n",
    "        \n",
    "# pool = mp.Pool(processes=min(4, len(all_files)))\n",
    "# t0 = time.time()\n",
    "# metatlas_dataset = pool.map(ma_data.get_data_for_atlas_df_and_file, all_files)\n",
    "# pool.close()\n",
    "# pool.terminate()\n",
    "# #If you're code crashes here, make sure to terminate any processes left open.\n",
    "#(print time.time() - t0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### One of the two blocks below builds the hits variable.  This holds your MSMS spectra (from within your mz, and rt ranges, and within the extra time indicated above).\n",
    "\n",
    "There are two options for generating the hits variable:\n",
    "1. block A: use when your files have msms. It create the hits variable and also saves a binary (pickled) serialized hits file to the output directory.\n",
    "2. block B: only run if your files were collected in MS1 mode\n",
    "3. If you have already run block A and then the kernel dies, you can skip block A and directly unplickle the binary hits file from the output directory. Skip block A, uncomment the Optional block and run it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##BLOCK A\n",
    "import warnings; warnings.simplefilter('ignore')\n",
    "t0 = time.time()\n",
    "\n",
    "hits=dp.get_msms_hits(metatlas_dataset,extra_time=True,keep_nonmatches=True, frag_mz_tolerance=0.01, ref_loc='/global/project/projectdirs/metatlas/projects/spectral_libraries/msms_refs_v3.tab')\n",
    "pickle.dump(hits, open(os.path.join(output_dir,polarity+'_hits.pkl'), \"wb\"))\n",
    "\n",
    "print(time.time() - t0)\n",
    "print('%s%s' % (len(hits),' <- total number of MSMS spectra found in your files'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## BLOCK B (uncomment lines below to run this.  Only use when all data files are MS1)\n",
    "#hits=pd.DataFrame([], columns=['database','id','file_name','msms_scan', u'score', u'num_matches', u'msv_query_aligned', u'msv_ref_aligned', u'name', u'adduct', u'inchi_key', u'precursor_mz', u'measured_precursor_mz'])\n",
    "#hits.set_index(['database','id','file_name','msms_scan'], inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optional: If you already have a pickled hits file and do not need to run get_msms_hits again, uncomment this block\n",
    "# hits = pickle.load(open(os.path.join(output_dir,polarity+'_hits.pkl'), \"rb\")) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. Adjust Retention Times. \n",
    "\n",
    "This block creates an interactive plot.  The top panel displays MSMS from within the two green RT bounds selected below (rt min and max, initially set in atlas).  When the database holds reference spectra, mirror plots are generated with the reference spectra inverted below the sample spectra.  The lower panel displays the EICs overlayed for all of the files in your selected groups.  You can highlight your groups different colors.  It is recommended that you do this, at least, for your extraction blank (or if not available, use a solvent injection blank).  This plot also displays radio buttons that can be interactively selected; the values will be exported in your final identifications table and in your atlas export.  Use these to mark peak/MSMS quality.\n",
    "\n",
    "How to use:\n",
    "1. STEP 1: Set peak flag radio buttons\n",
    "    1. OPTION A (custom flags): fill out the peak flags list (list of strings) \n",
    "        peak_flag_list = ('A','B') some recommendations are below.  \n",
    "    2. OPTION B (default flags): comment out the custom peak_flag_list line.  Uncomment the default peak_flags = \"\". \n",
    "        Flags default to: keep, remove, unresolvable isomers, check.\n",
    "2. STEP 2: Set EIC colors\n",
    "    1. Option A (custom EIC colors): fill out the colorlist in the format of below\n",
    "   > ***\n",
    "   > colorlist = [['color1nameorhexadec','partialgroupstring1'],\n",
    "   >         ['color2nameorhexadec','partialgroupstring2']]\n",
    "   > ***\n",
    "    <ul><li>You can add more comma delimited colors/groups as needed.</li>\n",
    "    <li>These are partial strings that match to you file names (not your group names).</li>\n",
    "    <li>The order they are listed in your list is the order they are displayed in the overlays (first is front, last is back)</li>\n",
    "    <li>Named colors available in matplotlib are here: https://matplotlib.org/3.1.0/gallery/color/named_colors.html\n",
    "          or use hexadecimal values '#000000'</li></ul>\n",
    "    B. Option B (default EIC colors): comment out the custom colorlist lines and uncomment the default colorlist = \"\". \n",
    "        Colors all default to black.\n",
    "3. User the right/left buttons on your keyboard to cycle through compounds in your atlas.\n",
    "4. Use the up/down buttons on your keyboard to cycle through MSMS spectra within the RT bounds of the lower plot.\n",
    "5. Use the horizontal rt min and rt max bars below the plots to adjust the rt bounds around your peak.  If there are multiple peaks, select one at a time and then click up/down to update the msms available in that new RT range.  If necessary evaluate your data in an external program such as mzmine to make sure you are selecting the correct peak.\n",
    "\n",
    "TIPS: use compound_idx = 0 in step 3 to change to a different compound in your atlas using the index number.  If your plot does not fit in your browser window, adjust height and width values.  Use alpha to change the transparency of the lines this is a value 0 (transparent) to 1 (opaque).\n",
    "\n",
    "DO NOT change your RT theoretical peak (the purple line).  It is locked from editing (unless you change a hidden parameter) and only to be changed in special cases.  The measured retention times of your peaks will be calculated and exported in your output files.  These will be compared with the RT theoreticals and used in your evidence of identification table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###STEP 1: Set the peak flag radio buttons using one of the two lines below, for custom flags or default flags\n",
    "import warnings; warnings.simplefilter('ignore')\n",
    "peak_flag_list=('','L1+ - 1 pk, good RT&MSMS','L1+ - known isomer overlap','L1+ - 1 pk, good RT, MSMS ok (coisolated mz/partial match/low int)','L1+ - 1 pk, good RT&MSMS from external library','L1 - 1 pk, correct RT, no MSMS or int too low for matching','L1 - 1 pk, good RT, very low intensity/poor pk shape','L2 put comp','L3 putative class','Remove - background/noise','Remove - bad EMA MSMS','Remove - bad MSMS NIST/MONA/Metlin')\n",
    "msms_flags_list = \"\" #\n",
    "#peak_flag_list =\"\"       # this will default to ('keep','remove','unresolvable isomers','poor peak shape')\n",
    "\n",
    "###STEP 2: Set the EIC line colors using on of the two lines below, for custom colors or default \n",
    "colorlist= [['red','ExCtrl'],                                                     \n",
    "           ['green','TxCtrl'],\n",
    "           ['blue','InjBl']]\n",
    "#colorlist=\"\"   # this will default to black\n",
    "\n",
    "###STEP 3\n",
    "a = dp.adjust_rt_for_selected_compound(metatlas_dataset, msms_hits=hits, peak_flags=peak_flag_list, msms_flags=msms_flags_list, color_me = colorlist, compound_idx=0,alpha=0.5,width=15,height=4.5)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Create filtered atlas excluding compounds marked removed\n",
    "\n",
    "Re-run the following before filtering atlas\n",
    "1. Get Groups (include InjBl)\n",
    "2. Get Atlas\n",
    "3. Get Data\n",
    "4. Get MSMS Hits"
{'fraction_non_alphanumeric': 0.13952673728967108, 'fraction_numerical': 0.008274245956490754, 'mean_word_length': 2.9288154897494305, 'pattern_counts': {'":': 496, '<': 83, '<?xml version=': 0, '>': 83, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 80, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5713489', 'n_tokens_mistral': 27626, 'n_tokens_neox': 26837, 'n_words': 8624}