<commit_before>from __future__ import print_function
import numpy as np
import sys
import json
import pandas as pd
import re
import glob as glob
import os
from collections import defaultdict
from xml.etree import cElementTree as ET
import multiprocessing as mp

import matplotlib.pyplot as plt

import metatlas.metatlas_objects as metob
from metatlas.helpers import metatlas_get_data_helper_fun as ma_data
from metatlas.helpers import dill2plots as dp

try:
    basestring
except NameError:  # python3
    basestring = str

# setting this very high easily causes out of memory
NUM_THREADS = 4

BATCH_FILE_PATH = '/project/projectdirs/metatlas/projects/mzmine_parameters/batch_files/'
PYTHON_BINARY = '/global/common/software/m2650/python-cori/bin/python'
BINARY_PATH = '/project/projectdirs/metatlas/projects/mzmine_parameters/MZmine'

#we don't need to request haswell on genepool partition of Cori
#remove this line
#SBATCH -C haswell


SLURM_HEADER = """#!/bin/bash
#SBATCH -N 1 -c 64
#SBATCH --exclusive
#SBATCH --error="slurm.err"
#SBATCH --output="slurm.out"
#SBATCH --qos=genepool
#SBATCH -A pkscell
#SBATCH -t 24:00:00
#SBATCH -L project

export MPLBACKEND="agg"
export HDF5_USE_FILE_LOCKING=FALSE
# cori specific tells it to not allocate memory on a per thread basis
export MALLOC_ARENA_MAX=1

#module load java

# echo every command and terminate script if there is an error
set -ex

env | grep -i java | sort

# to see the resources a job used:
# sacct -j <job_id> --format jobidraw,jobname,maxrss,maxvmsize --unit G

"""

def see_if_plots_work():
    print('I am here')
    fig = plt.figure()
    print('I made a figure')
    fig.savefig('fig.pdf')
    print('I saved a figure')


def make_figures_from_filtered_data(params,all_files,my_atlas):

    pool = mp.Pool(processes=min(NUM_THREADS, len(all_files)))
    print("# Made Pool")
    print('# Getting data 2')
    metatlas_dataset = pool.map(ma_data.get_data_for_atlas_df_and_file, all_files)
    print('# Getting data 3')
    pool.close()
    pool.join()
    # pool.terminate()
    print('# Acquired Data for plotting')
    output_dir = os.path.join(params['basedir'],'%s_%s'%(params['basename'],params['polarity']))
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)
    atlas_identifications = dp.export_atlas_to_spreadsheet(my_atlas,os.path.join(output_dir,'atlas_export.csv'))
    # atlas_identifications = dp.export_atlas_to_spreadsheet(myAtlas,'%s/sheets/%s.csv'%(plot_location_label,myAtlas.name))
    peak_height = dp.make_output_dataframe(input_fname = '',input_dataset = metatlas_dataset,include_lcmsruns = [],exclude_lcmsruns = [], fieldname='peak_height' , output_loc=os.path.join(output_dir,'sheets'))
    peak_area = dp.make_output_dataframe(input_fname = '',input_dataset = metatlas_dataset,include_lcmsruns = [],exclude_lcmsruns = [], fieldname='peak_area' , output_loc=os.path.join(output_dir,'sheets'))
    mz_peak = dp.make_output_dataframe(input_fname = '',input_dataset = metatlas_dataset,include_lcmsruns = [],exclude_lcmsruns = [], fieldname='mz_peak' , output_loc=os.path.join(output_dir,'sheets'))
    rt_peak = dp.make_output_dataframe(input_fname = '', input_dataset = metatlas_dataset,include_lcmsruns = [],exclude_lcmsruns = [],fieldname='rt_peak' , output_loc=os.path.join(output_dir,'sheets'))
    mz_centroid = dp.make_output_dataframe(input_fname = '',input_dataset = metatlas_dataset,include_lcmsruns = [],exclude_lcmsruns = [], fieldname='mz_centroid' , output_loc=os.path.join(output_dir,'sheets'))
    rt_centroid = dp.make_output_dataframe(input_fname = '',input_dataset = metatlas_dataset,include_lcmsruns = [],exclude_lcmsruns = [], fieldname='rt_centroid' , output_loc=os.path.join(output_dir,'sheets'))
    print('$$$$ Making ID FIGURES')
    dp.make_identification_figure_v2(input_dataset = metatlas_dataset, input_fname = '', include_lcmsruns = [],exclude_lcmsruns = params['blank_str'].append('QC'), output_loc=os.path.join(output_dir,'identification'))
    print('$$$$ Done Making FIGURES')

def clean_and_filter_mzmine_output(json_filename=None,n_peaks=1000):
    """
    This is a command line function that runs after the first feature finding is complete.
    After mzmine runs, make an atlas, get data, and filter
    The principle is to cast the mzmine peaks into three different categories:
        blank
        bad peak shape
        good no msms
        good with msms
    First, Remove in blank:
        This can be done from mzmine output
        exclude peaks where neither peak height nor peak area is appreciably greater in sample vs blank
    Second, Remove Bad Peak shape:
        seems like rt_peak being near rt@max is a good thing
        seems like looking at intensity max
        you don't want to exclude peaks like leu-ileu.
        many peaks will have multiple-humps
        A simplified logic is: is max intensity near the rt_peak? if yes, is it decaying by 10x +/- 0.5 minutes. Obviously, peaks that are +/- 0.5 minutes appart will get removed. Thus, find local minimum +/- 0.5 minutes from rt_peak. Set those as rt_min and rt_max.
    Third remove no MS/MS
        exclude if no feature in non-blank samples has msms
        
    Last: only keep n_peaks by max intensity in any sample.

    """


    if not json_filename:
        json_filename = sys.argv[1]
    with open(json_filename) as data_file:    
        params = json.load(data_file)
    file_to_convert = os.path.join(params['basedir'],'intermediate_results','%s_%s.csv'%(params['basename'],params['polarity']))
    print('# Working on %s %s'%(params['basename'],params['polarity']))
    if os.path.isfile(file_to_convert):
        # take the comprehensive csv from mzmine and make a peak-height only version of it
        df,original_mzmine = metatlas_formatted_atlas_from_mzmine_output(file_to_convert,params['polarity'],
                                                                                 make_atlas=False,min_rt=0.55,
                                                                                remove_fragments=False,
                                                                                remove_adducts=False,
                                                                                remove_clusters=False,)
        df.to_csv(file_to_convert.replace('.csv','') + '_formatted.csv',index=True) #save a simplified mzmine-like csv as a backup of all features found.
        
        #Filter features found in the blank
        df = clean_up_mzmine_dataframe(df)
        df_blank_compare = df.transpose().groupby(['b' if any([s in g.lower() for s in params['blank_str']]) else 's' for g in df.columns]).max().transpose()
        if 'b' in df_blank_compare.columns:
            df_features_not_in_blank = df_blank_compare[df_blank_compare['s'] > (params['sample_to_blank'] * df_blank_compare['b'])]
            print('# %d features total'%(df.shape[0]))
            print('# %d features removed by blank\n'%(df_blank_compare.shape[0] - df_features_not_in_blank.shape[0]))
        else:
            print('# No files have "blank" or "injbl" in their names.')
            df_features_not_in_blank = df_blank_compare
        df_features_not_in_blank.reset_index(inplace=True)
        print('#There are now %d features not in blank'%df_features_not_in_blank.shape[0])<commit_msg>added testing filter<commit_after>from __future__ import print_function
import numpy as np
import sys
import json
import pandas as pd
import re
import glob as glob
import os
from collections import defaultdict
from xml.etree import cElementTree as ET
import multiprocessing as mp

import matplotlib.pyplot as plt

import metatlas.metatlas_objects as metob
from metatlas.helpers import metatlas_get_data_helper_fun as ma_data
from metatlas.helpers import dill2plots as dp

try:
    basestring
except NameError:  # python3
    basestring = str

# setting this very high easily causes out of memory
NUM_THREADS = 4

BATCH_FILE_PATH = '/project/projectdirs/metatlas/projects/mzmine_parameters/batch_files/'
PYTHON_BINARY = '/global/common/software/m2650/python-cori/bin/python'
BINARY_PATH = '/project/projectdirs/metatlas/projects/mzmine_parameters/MZmine'

#we don't need to request haswell on genepool partition of Cori
#remove this line
#SBATCH -C haswell


SLURM_HEADER = """#!/bin/bash
#SBATCH -N 1 -c 64
#SBATCH --exclusive
#SBATCH --error="slurm.err"
#SBATCH --output="slurm.out"
#SBATCH --qos=genepool
#SBATCH -A pkscell
#SBATCH -t 24:00:00
#SBATCH -L project

export MPLBACKEND="agg"
export HDF5_USE_FILE_LOCKING=FALSE
# cori specific tells it to not allocate memory on a per thread basis
export MALLOC_ARENA_MAX=1

#module load java

# echo every command and terminate script if there is an error
set -ex

env | grep -i java | sort

# to see the resources a job used:
# sacct -j <job_id> --format jobidraw,jobname,maxrss,maxvmsize --unit G

"""

def see_if_plots_work():
    print('I am here')
    fig = plt.figure()
    print('I made a figure')
    fig.savefig('fig.pdf')
    print('I saved a figure')


def make_figures_from_filtered_data(params,all_files,my_atlas):

    pool = mp.Pool(processes=min(NUM_THREADS, len(all_files)))
    print("# Made Pool")
    print('# Getting data 2')
    metatlas_dataset = pool.map(ma_data.get_data_for_atlas_df_and_file, all_files)
    print('# Getting data 3')
    pool.close()
    pool.join()
    # pool.terminate()
    print('# Acquired Data for plotting')
    output_dir = os.path.join(params['basedir'],'%s_%s'%(params['basename'],params['polarity']))
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)
    atlas_identifications = dp.export_atlas_to_spreadsheet(my_atlas,os.path.join(output_dir,'atlas_export.csv'))
    # atlas_identifications = dp.export_atlas_to_spreadsheet(myAtlas,'%s/sheets/%s.csv'%(plot_location_label,myAtlas.name))
    peak_height = dp.make_output_dataframe(input_fname = '',input_dataset = metatlas_dataset,include_lcmsruns = [],exclude_lcmsruns = [], fieldname='peak_height' , output_loc=os.path.join(output_dir,'sheets'))
    peak_area = dp.make_output_dataframe(input_fname = '',input_dataset = metatlas_dataset,include_lcmsruns = [],exclude_lcmsruns = [], fieldname='peak_area' , output_loc=os.path.join(output_dir,'sheets'))
    mz_peak = dp.make_output_dataframe(input_fname = '',input_dataset = metatlas_dataset,include_lcmsruns = [],exclude_lcmsruns = [], fieldname='mz_peak' , output_loc=os.path.join(output_dir,'sheets'))
    rt_peak = dp.make_output_dataframe(input_fname = '', input_dataset = metatlas_dataset,include_lcmsruns = [],exclude_lcmsruns = [],fieldname='rt_peak' , output_loc=os.path.join(output_dir,'sheets'))
    mz_centroid = dp.make_output_dataframe(input_fname = '',input_dataset = metatlas_dataset,include_lcmsruns = [],exclude_lcmsruns = [], fieldname='mz_centroid' , output_loc=os.path.join(output_dir,'sheets'))
    rt_centroid = dp.make_output_dataframe(input_fname = '',input_dataset = metatlas_dataset,include_lcmsruns = [],exclude_lcmsruns = [], fieldname='rt_centroid' , output_loc=os.path.join(output_dir,'sheets'))
    print('$$$$ Making ID FIGURES')
    dp.make_identification_figure_v2(input_dataset = metatlas_dataset, input_fname = '', include_lcmsruns = [],exclude_lcmsruns = params['blank_str'].append('QC'), output_loc=os.path.join(output_dir,'identification'))
    print('$$$$ Done Making FIGURES')

def clean_and_filter_mzmine_output(json_filename=None,n_peaks=1000,do_test=False):
    """
    This is a command line function that runs after the first feature finding is complete.
    After mzmine runs, make an atlas, get data, and filter
    The principle is to cast the mzmine peaks into three different categories:
        blank
        bad peak shape
        good no msms
        good with msms
    First, Remove in blank:
        This can be done from mzmine output
        exclude peaks where neither peak height nor peak area is appreciably greater in sample vs blank
    Second, Remove Bad Peak shape:
        seems like rt_peak being near rt@max is a good thing
        seems like looking at intensity max
        you don't want to exclude peaks like leu-ileu.
        many peaks will have multiple-humps
        A simplified logic is: is max intensity near the rt_peak? if yes, is it decaying by 10x +/- 0.5 minutes. Obviously, peaks that are +/- 0.5 minutes appart will get removed. Thus, find local minimum +/- 0.5 minutes from rt_peak. Set those as rt_min and rt_max.
    Third remove no MS/MS
        exclude if no feature in non-blank samples has msms
        
    Last: only keep n_peaks by max intensity in any sample.

    """


    if not json_filename:
        json_filename = sys.argv[1]
    with open(json_filename) as data_file:    
        params = json.load(data_file)
    file_to_convert = os.path.join(params['basedir'],'intermediate_results','%s_%s.csv'%(params['basename'],params['polarity']))
    print('# Working on %s %s'%(params['basename'],params['polarity']))
    if os.path.isfile(file_to_convert):
        # take the comprehensive csv from mzmine and make a peak-height only version of it
        df,original_mzmine = metatlas_formatted_atlas_from_mzmine_output(file_to_convert,params['polarity'],
                                                                                 make_atlas=False,min_rt=0.55,
                                                                                remove_fragments=False,
                                                                                remove_adducts=False,
                                                                                remove_clusters=False,)
        if do_test == True:
            #if only a test, just do a couple of features
            df = df.head(10)
            
        df.to_csv(file_to_convert.replace('.csv','') + '_formatted.csv',index=True) #save a simplified mzmine-like csv as a backup of all features found.
        
        #Filter features found in the blank
        df = clean_up_mzmine_dataframe(df)
        df_blank_compare = df.transpose().groupby(['b' if any([s in g.lower() for s in params['blank_str']]) else 's' for g in df.columns]).max().transpose()
        if 'b' in df_blank_compare.columns:
            df_features_not_in_blank = df_blank_compare[df_blank_compare['s'] > (params['sample_to_blank'] * df_blank_compare['b'])]
            print('# %d features total'%(df.shape[0]))
            print('# %d features removed by blank\n'%(df_blank_compare.shape[0] - df_features_not_in_blank.shape[0]))
        else:
            print('# No files have "blank" or "injbl" in their names.')
            df_features_not_in_blank = df_blank_compare
        df_features_not_in_blank.reset_index(inplace=True)
        print('#There are now %d features not in blank'%df_features_not_in_blank.shape[0])
{'fraction_non_alphanumeric': 0.1033830778580593, 'fraction_numerical': 0.0056722263488419205, 'mean_word_length': 3.4987849331713243, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 2}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1538124', 'n_tokens_mistral': 4775, 'n_tokens_neox': 4619, 'n_words': 1180}