{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "567f7581",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---------- Preparations ----------\n",
    "# Load Libraries\n",
    "library(Spectra)\n",
    "library(MsBackendMgf)\n",
    "library(MsBackendHmdb)\n",
    "library(MsCoreUtils)\n",
    "library(MsBackendMsp)\n",
    "library(readr)\n",
    "library(dplyr)\n",
    "# 3 dependencies for latest MassBank version\n",
    "library(rvest)\n",
    "library(stringr)\n",
    "library(xml2)\n",
    "options(warn=-1)\n",
    "library(\"mzR\")\n",
    "#' Download the CompDb database using curl\n",
    "library(curl)\n",
    "#library(CompoundDb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9fee3781",
   "metadata": {},
   "outputs": [],
   "source": [
    "download_specDB_new <- function(input_dir, db = \"all\"){\n",
    "\n",
    "    if (dir.exists(input_dir)){\n",
    "        # Track Time\n",
    "        start_time <- Sys.time()\n",
    "\n",
    "        # only input available as of now\n",
    "        databases <- 'gnps, hmdb, mbank, all'\n",
    "        # creat a summary file, open and store timings of download and version if possible\n",
    "        if (!(file.exists(paste(input_dir, \"/summaryFile.txt\", sep = \"\")))){\n",
    "            summaryFile <- paste(input_dir, \"/summaryFile.txt\", sep = \"\")\n",
    "            file.create(summaryFile, recursive = TRUE)\n",
    "        }\n",
    "        else{\n",
    "            summaryFile <- paste(input_dir, \"/summaryFile.txt\", sep = \"\")\n",
    "        }\n",
    "        file.conn <- file(summaryFile)\n",
    "        open(file.conn, open = \"at\")\n",
    "        # gnps\n",
    "        if (db == \"all\" || db ==\"gnps\"){\n",
    "            # Download file\n",
    "            system(paste(\"wget -P\", \n",
    "                         input_dir, \n",
    "                         \"https://gnps-external.ucsd.edu/gnpslibrary/ALL_GNPS.msp\", \n",
    "                         sep =  \" \"))\n",
    "\n",
    "            # load the spectra into MsBackendMgf\n",
    "            gnpsdb <- Spectra(paste(input_dir, \"/ALL_GNPS.msp\", sep = ''), source = MsBackendMsp())\n",
    "            save(gnpsdb, file = paste(input_dir,\"/gnps.rda\", sep = \"\"))\n",
    "\n",
    "            # delete the database in its format to free up space\n",
    "            system(paste(\"rm\", (paste(input_dir, \"/ALL_GNPS.mgf\", sep = '')), sep = \" \"))\n",
    "\n",
    "            writeLines(paste(\"GNPS saved at\", Sys.time(), sep=\" \"),con=file.conn)\n",
    "        }\n",
    "        #mbank\n",
    "        if (db == \"all\" || db ==\"mbank\"){\n",
    "\n",
    "            #print(\"MassBank WORKS\")\n",
    "\n",
    "            page <- read_html(\"https://github.com/MassBank/MassBank-data/releases\")\n",
    "            page %>%\n",
    "                html_nodes(\"a\") %>%       # find all links\n",
    "                html_attr(\"href\") %>%     # get the url\n",
    "                str_subset(\"MassBank_NIST.msp\") -> tmp # find those that have the name MassBank_NIST.msp\n",
    "\n",
    "            #download file\n",
    "            system(paste(\"wget \",\n",
    "                         \"https://github.com\", tmp[1],\n",
    "                         sep =  \"\"))\n",
    "\n",
    "            mbank <- Spectra(paste(input_dir, \"/MassBank_NIST.msp\", sep = ''), source = MsBackendMsp())\n",
    "            save(mbank, file = paste(input_dir,\"/mbankNIST.rda\", sep = \"\"))\n",
    "\n",
    "            # delete the database in its format to free up space\n",
    "            system(paste(\"rm\", (paste(input_dir, \"/MassBank_NIST.msp\", sep = '')), sep = \" \"))\n",
    "\n",
    "            # obtain the month and year for the database release to add to summary\n",
    "            res <- str_match(tmp[1], \"download/\\\\s*(.*?)\\\\s*/MassBank_NIST\")\n",
    "\n",
    "            writeLines(paste(\"MassBank saved at\", Sys.time(), \"with release version\", res[,2], sep=\" \"),con=file.conn)\n",
    "        }\n",
    "\n",
    "        #mbank\n",
    "        if (db == \"all\" || db ==\"hmdb\"){\n",
    "            # extract HMDB Current version\\n\",\n",
    "            html <- read_html(\"https://hmdb.ca/downloads\")\n",
    "            strings <- html%>% html_elements(\"a\") %>% html_text2()\n",
    "            ls <- unique(strings)\n",
    "            hmdb_curr_ver <- c()\n",
    "            for (i in ls){\n",
    "                if (grepl(\"Current\", i)){\n",
    "                hmdb_curr_ver<- c(i, hmdb_curr_ver)\n",
    "                }\n",
    "            }\n",
    "            dbname <- \"CompDb.Hsapiens.HMDB.5.0.sqlite\"\n",
    "            db_file <- file.path(tempdir(), dbname)\n",
    "            curl_download(\n",
    "                paste0(\"https://github.com/jorainer/MetaboAnnotationTutorials/\",\n",
    "               \"releases/download/2021-11-02/\", dbname),\n",
    "                destfile = db_file)\n",
    "            #' Load a CompDb database with compound annotation from HMDB\n",
    "            cdb <- CompDb(db_file)\n",
    "            hmdb <- Spectra(cdb)\n",
    "            hmdb$collisionEnergy <- as.numeric(hmdb$collisionEnergy)\n",
    "            hmdb <- setBackend(hmdb, backend = MsBackendDataFrame())\n",
    "            save(hmdb, file = paste(input_dir,\"/hmdb.rda\", sep = \"\"))\n",
    "            writeLines(paste(\"HMDB saved at\", Sys.time(), \"with release version\", hmdb_curr_ver, sep=\" \"),con=file.conn)\n",
    "        }\n",
    "\n",
    "        #wrong input error message\n",
    "        else if (!grepl(db, databases, fixed = TRUE)){\n",
    "            stop(\"Wrong db input. Following inputs apply: gnps, hmdb, mbank or all\")\n",
    "        }\n",
    "        close(file.conn)\n",
    "        end_time <- Sys.time()\n",
    "        print(end_time - start_time)\n",
    "    }\n",
    "    else{\n",
    "        stop(\"Your input_dir is incorrect. Please provide the directory where all your input files are stored.\")\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3707280d",
   "metadata": {},
   "outputs": [],
   "source": [
    "##-----------------------------------------------------------------\n",
    "## filter intensity \n",
    "##-----------------------------------------------------------------\n",
    "\n",
    "#' Define a filtering function and remove peaks less than 0.05 of intensity\n",
    "low_int <- function(c, ...) {\n",
    "    c > max(c, na.rm = TRUE) * 0.05\n",
    "}\n",
    "\n",
    "# Usage:\n",
    "# filterIntensity(spectra_object, intensity = low_int)\n",
    "\n",
    "\n",
    "##-----------------------------------------------------------------\n",
    "## normalize intensity \n",
    "##-----------------------------------------------------------------\n",
    "\n",
    "#' Define a function to *normalize* the intensities\n",
    "norm_int <- function(y, ...) {\n",
    "    maxint <- max(y[, \"intensity\"], na.rm = TRUE)\n",
    "    y[, \"intensity\"] <- 100 * y[, \"intensity\"] / maxint\n",
    "    y\n",
    "}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d263f9f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Specifying a function for creating result directories for each input mzml\n",
    "# input for the function:\n",
    "# input directory\n",
    "ms2_rfilename<- function(input_dir){\n",
    "    if (dir.exists(input_dir)){\n",
    "        #list_ms2_files <- intersect(list.files(input_dir, pattern = \"_PRM_\"), list.files(input_dir, pattern = \".mzML\"))\n",
    "        list_ms2_files <- list.files(input_dir, pattern = \".mzML\")\n",
    "        mzml_file <- paste(input_dir, \"/\", list_ms2_files, sep = \"\")\n",
    "\n",
    "        #store the result file names to return to this function as output\n",
    "        mzml_files <- c()\n",
    "        ResultFileNames <- c()\n",
    "        File_id <- c()\n",
    "        nx <- 0\n",
    "        # x is mzML files\n",
    "        for (i in 1:length(mzml_file)){\n",
    "            nx <- nx+1\n",
    "            # remove .mzML to extract just the names\n",
    "            mzml_filex <- str_replace(mzml_file[i], input_dir, \".\")\n",
    "            name_mzmls <- str_remove(as.character(mzml_filex), \".mzML\")\n",
    "            name_mzmlsd <- str_remove(mzml_file[i], \".mzML\")\n",
    "            #name_mzml <- str_replace(name_mzmls, input_dir, \"./\")\n",
    "            #' for each file a subdirectory is created to store all results in that, add working directory\n",
    "            if (!file.exists(name_mzmlsd)){\n",
    "                dir.create(name_mzmlsd) ##create folder\n",
    "            }\n",
    "            ResultFileNames<- c(ResultFileNames, name_mzmls)\n",
    "            mzml_files <- c(mzml_files, mzml_filex)\n",
    "            File_id <- c(File_id, paste(\"file_\", nx, sep = \"\"))\n",
    "        }\n",
    "        input_table <- cbind(mzml_files, ResultFileNames, File_id)\n",
    "\n",
    "        write.csv(input_table, paste(input_dir, \"/input_table.csv\", sep = \"\"))\n",
    "        return(data.frame(input_table))\n",
    "    }\n",
    "    else{\n",
    "        stop(\"Your input_dir is incorrect. Please provide the directory where all your input files are stored. : ) Good Luck\")\n",
    "    }\n",
    "}\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "26952fc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#' All spectra in mzML files preprocessing, return two outputs, pre-processed MS2 spectra and all precursor masses\n",
    "# x is one mzML file\n",
    "spec_Processing <- function(input_dir, x, result_dir){\n",
    "\n",
    "    x <- paste(input_dir, str_remove(x, \".\"), sep = \"\")\n",
    "\n",
    "    result_dir <- paste(input_dir, str_remove(result_dir, \".\"), sep = \"\")\n",
    "    if (file.exists(x) && substring(x, nchar(x)) == \"L\"){\n",
    "        if (dir.exists(result_dir)){\n",
    "            # read the spectra\n",
    "            sps_all <- Spectra(x, backend = MsBackendMzR())\n",
    "            #' Change backend to a MsBackendDataFrame: load data into memory\n",
    "            #sps_all <- setBackend(sps_all, MsBackendDataFrame())\n",
    "            #' Filter Empty Spectra\n",
    "            sps_all <- filterEmptySpectra(sps_all)\n",
    "            #' Extract Precursor m/z(s) in each file\n",
    "            pre_mz <- unique(precursorMz(sps_all))\n",
    "            #' Remove any NAs\n",
    "            pre_mz <- na.omit(pre_mz)\n",
    "            if (!file.exists(paste(result_dir, \"/processedSpectra.mzML\", sep = \"\"))){\n",
    "                export(sps_all, backend = MsBackendMzR(), file = paste(result_dir, \"/processedSpectra.mzML\", sep = \"\"))\n",
    "            }\n",
    "            if (!file.exists(paste(result_dir, \"/premz_list.txt\", sep = \"\"))){\n",
    "                write.table(pre_mz, file = paste(result_dir, \"/premz_list.txt\", sep = \"\"), sep = \"/t\",row.names = FALSE, col.names = FALSE)\n",
    "            }\n",
    "\n",
    "            spsall_pmz <- list(sps_all, pre_mz)\n",
    "            return(spsall_pmz)\n",
    "        }\n",
    "        else{\n",
    "            stop(\"Seems like it is not the result directory of the input .mzML file which is provided as x. Please use the function ms2_rfilename to generate a result directory or create one yourself with the same name as the .mzML input file.\")\n",
    "        }\n",
    "    }\n",
    "    else{\n",
    "        stop(\"Are you sure x is an mzML input file?\")\n",
    "    }\n",
    "\n",
    "}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2f11afa3",
   "metadata": {},
   "outputs": [],
   "source": [
    "spec2_Processing <- function(z, obj, spec = \"spec_all\", ppmx = 15){\n",
    "    if (spec == \"spec_all\"){\n",
    "        #' Subset the dataset to MS2 spectra matching the m/z\n",
    "        sps <- filterPrecursorMzValues(obj, mz = z + ppm(c(-z, z), 10))\n",
    "    } else if (spec == \"gnps\"){\n",
    "        #gnps spectra that contains precursor mass\n",
    "        has_mz <- containsMz(obj, mz = z, ppm = ppmx)\n",
    "        #' Subset the GNPS Spectra\n",
    "        sps <- obj[has_mz]\n",
    "    } else if (spec == \"hmdb\"){\n",
    "        #hmdb spectra that contains precursor mass\n",
    "        has_mz <- containsMz(obj, mz = z, ppm = ppmx)\n",
    "        #' Subset the HMDB Spectra\n",
    "        sps <- obj[has_mz]\n",
    "    } else if (spec == \"mbank\"){\n",
    "        has_mz <- containsMz(obj, mz = z, ppm = ppmx)\n",
    "        #' Subset the MB Spectra\n",
    "        sps <- obj[has_mz]\n",
    "    }\n",
    "\n",
    "    #wrong input error message\n",
    "    else if (!grepl(db, databases, fixed = TRUE)){\n",
    "        stop(\"Wrong db input. Following inputs apply: gnps, hmdb, mbank or all\")\n",
    "    }\n",
    "\n",
    "    if (length(sps)>0){\n",
    "        #' Apply the function to filter the spectra\n",
    "        sps <- filterIntensity(sps, intensity = low_int)\n",
    "        #' *Apply* the function to the data\n",
    "        sps <- addProcessing(sps, norm_int)\n",
    "        # cleaning peaks that are heavier or equal to the precursor mass\n",
    "        pkd <- peaksData(sps)@listData\n",
    "\n",
    "        #obtain the list of peaks that are higher or equal to precursor mass\n",
    "        # y is peaksData from spectra\n",
    "        removePrecursorPeaks <- function(m){\n",
    "            m <- m[m[, \"mz\"] <= z, ]\n",
    "        }\n",
    "        # use lapply to apply the function to the list of peaksData\n",
    "        pkd <- lapply(pkd, removePrecursorPeaks)\n",
    "        #store the indices of spectra with 0 peaks\n",
    "        store_i <- c()\n",
    "        for (i in 1:length(pkd)){\n",
    "            if (is.null(nrow(pkd[[i]]))){\n",
    "                #convert the object of one peak into a matrix\n",
    "                mz <- pkd[[i]][[1]]\n",
    "                intensity <- pkd[[i]][[2]]\n",
    "                mat <- cbind(mz, intensity)\n",
    "                pkd[[i]] <- mat\n",
    "            }else if(nrow(pkd[[i]])==0){\n",
    "                #store indices with 0 peaks\n",
    "                store_i <- c(store_i, i)\n",
    "            }\n",
    "        }\n",
    "        # if 0 peaks, remove the relevant spectra from gnps or hmdb or mbank\n",
    "        if (!(is.null(store_i))){\n",
    "            pkd <- pkd[-(store_i)]\n",
    "            sps <- sps[-(store_i)]\n",
    "            peaksData(sps@backend)<- pkd\n",
    "        }else {\n",
    "            peaksData(sps@backend)<- pkd\n",
    "        }\n",
    "        return(sps)\n",
    "    }\n",
    "    else {\n",
    "        sps <- NULL\n",
    "        return(sps)\n",
    "    }\n",
    "}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "eb1e74e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "##-----------------------------------------------------------------\n",
    "## Extract peaksdata in a dataframe\n",
    "##-----------------------------------------------------------------\n",
    "\n",
    "#' obtain peaksData for each spectral matching between query and database spectra\n",
    "#inputs a is best match from Database, b is best match from query spectra\n",
    "peakdf <- function(a, b, ppmx){\n",
    "\n",
    "    #' obtain peaklists for both query spectra and best matched spectra from MassBank\n",
    "    z<- peaksData(a)[[1]] #from GNPS/HMDB\n",
    "    y <- peaksData(b)[[1]] #from query\n",
    "    if (!(nrow(z)==0)){\n",
    "    #' Since we used 15 ppm, so to find the range, calculate the mass error range\n",
    "    range <-  (ppmx/1000000)*y[,\"mz\"]\n",
    "    y <- cbind(y, range)\n",
    "    low_range <- y[,\"mz\"]-y[,\"range\"] # low range of m/z\n",
    "    high_range <- y[,\"mz\"]+y[,\"range\"] # high range of m/z\n",
    "    y <- cbind(y, low_range, high_range)\n",
    "    #from GNPS/HMDB/MassBank spectra\n",
    "    mz.z <- c()\n",
    "    intensity.z <- c()\n",
    "    #from query spectra\n",
    "    mz.y <- c()\n",
    "    intensity.y <- c()\n",
    "    #difference between their intensity\n",
    "    diff <- c()\n",
    "\n",
    "    #' for all rows of y\n",
    "    for (m in 1:nrow(y)){\n",
    "        #' for all rows of z\n",
    "\n",
    "        for(j in 1:nrow(z)){\n",
    "\n",
    "            ###################################################################\n",
    "\n",
    "            ## IFELSE Statement no.2 -- LOOP 1.1.1.1\n",
    "\n",
    "            #' if the m/z of MB Spectra is within the 20 ppm range, save difference between intensities\n",
    "            if (y[m,\"low_range\"] <= z[j, \"mz\"] && z[j, \"mz\"] <= y[m,\"high_range\"]){\n",
    "\n",
    "                #GNPS/HMDB\n",
    "                mz_z <- as.numeric(z[j, \"mz\"])\n",
    "                mz.z <- c(mz.z, mz_z)\n",
    "\n",
    "                intensity_z <- as.numeric(z[j, \"intensity\"])\n",
    "                intensity.z <- c(intensity.z, intensity_z)\n",
    "\n",
    "                #QUERY\n",
    "                mz_y <- as.numeric(y[m, \"mz\"])\n",
    "                mz.y <- c(mz.y, mz_y)\n",
    "\n",
    "                intensity_y <- as.numeric(y[m, \"intensity\"])\n",
    "                intensity.y <- c(intensity.y, intensity_y)\n",
    "\n",
    "                #Difference between intensities\n",
    "                difference <- as.numeric(abs(z[j, \"intensity\"]-y[m, \"intensity\"]))\n",
    "                diff <- c(diff, difference)\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "    df_peaklists <- cbind(mz.y, intensity.y, mz.z, intensity.z, diff)\n",
    "    return(df_peaklists)\n",
    "    }\n",
    "    else{\n",
    "        df_peaklists <- NULL\n",
    "        return(df_peaklists)\n",
    "    }\n",
    "    #output is a dataframe with mz and intensity from db spectra and query spectra and their difference\n",
    "}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "5be702aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "##-----------------------------------------------------------------\n",
    "## Plotting Mirror Spectra\n",
    "##-----------------------------------------------------------------\n",
    "\n",
    "#' Specifying a function to draw peak labels\n",
    "#label_fun <- function(x) {\n",
    "    #ints <- unlist(intensity(x))\n",
    "    #mzs <- format(unlist(mz(x)), digits = 4)\n",
    "    #mzs[ints < 5] <- \"\"\n",
    "    #mzs\n",
    "#}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "355411dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "spec_dereplication_file <- function(mzml_file, pre_tbl, proc_mzml, db, result_dir, file_id, input_dir, no_of_candidates = 30, ppmx, error = TRUE){\n",
    "    # if the database selected is HMDB or all\n",
    "    # if the database selected is GNPS or all\n",
    "    if (db == \"all\" || db ==\"gnps\"){\n",
    "        if (file.exists(paste(input_dir,\"/gnps.rda\", sep = \"\"))){\n",
    "            # load the gnps spectral database\n",
    "            load(file = paste(input_dir,\"/gnps.rda\", sep = \"\"))\n",
    "        }\n",
    "        else if (file.exists(\"gnps.rda\")){\n",
    "            # load the gnps spectral database\n",
    "            load(\"gnps.rda\")\n",
    "        }\n",
    "        \n",
    "    }\n",
    "    # if the database selected is HMDB or all\n",
    "    if (db == \"all\" || db ==\"hmdb\"){\n",
    "        \n",
    "        if (file.exists(paste(input_dir,\"/hmdb.rda\", sep = \"\"))){\n",
    "            # load the hmdb spectral database\n",
    "            load(file = paste(input_dir,\"/hmdb.rda\", sep = \"\"))\n",
    "        }\n",
    "        else if (file.exists(\"hmdb.rda\")){\n",
    "            # load the gnps spectral database\n",
    "            load(\"hmdb.rda\")\n",
    "        }\n",
    "        \n",
    "    }\n",
    "    # if the database selected is HMDB or all\n",
    "    if (db == \"all\" || db == \"mbank\"){\n",
    "        if (file.exists(paste(input_dir,\"/mbankNIST.rda\", sep = \"\"))){\n",
    "            # load the mbank spectral database\n",
    "            load(file = paste(input_dir,\"/mbankNIST.rda\", sep = \"\"))\n",
    "        }\n",
    "        else if (file.exists(\"mbankNIST.rda\")){\n",
    "            # load the gnps spectral database\n",
    "            load(\"mbankNIST.rda\")\n",
    "        }\n",
    "        \n",
    "    }\n",
    "\n",
    "    # read spectra object\n",
    "    sps_all <- Spectra(proc_mzml, source = MsBackendMzR())\n",
    "\n",
    "    # extract precursor m/z\n",
    "    tbl <- read.table(pre_tbl)\n",
    "    pre_mz <- tbl[[1]]\n",
    "\n",
    "    # common feature information\n",
    "\n",
    "    id_X <- c() # id\n",
    "    premz <- c() # precursor mz\n",
    "    rtmin <- c() # stores rtmin\n",
    "    rtmax <- c() # stores rtmax\n",
    "    rtmed <- c() # stores calculated median of rtmin and rtmax\n",
    "    rtmean <- c() # stores calculated mean of rtmin and rtmax\n",
    "    col_eng <- c() # stores collision energy\n",
    "    pol <- c() # stores polarity\n",
    "    int <- c() # store intensity\n",
    "    source_file <- c() # source file\n",
    "    nx <- 0 # numbering the ids\n",
    "\n",
    "    pre_mzs <- listenv() # list for holding pre_mz futures\n",
    "\n",
    "    # for each pre mass\n",
    "    for (x in pre_mz){\n",
    "        print(x)\n",
    "        # to name the file\n",
    "        nx <- nx+1\n",
    "        # filter spectra based on precusror m/z\n",
    "        # this is done to extract all common information for id_X\n",
    "        spsrt <- filterPrecursorMzRange(sps_all, x)\n",
    "        # id based on file id,\n",
    "        id_Xx <- paste(file_id,  \"M\",  as.character(round(x, digits = 0)),\n",
    "                        \"R\", as.character(round(median(spsrt$rtime, na.rm = TRUE), digits = 0)),\n",
    "                        \"ID\", as.character(nx), sep = '')\n",
    "        id_X <- c(id_X, id_Xx)\n",
    "        # pre_mas\n",
    "        pre <- x\n",
    "        premz <- c(premz, pre)\n",
    "\n",
    "        # rt min\n",
    "        rti <- min(spsrt$rtime)\n",
    "        rtmin <- c(rtmin, rti)\n",
    "        #rt max\n",
    "        rtx <- max(spsrt$rtime)\n",
    "        rtmax <- c(rtmax, rtx)\n",
    "        #rt median\n",
    "        rtmd <- median(spsrt$rtime, na.rm = TRUE)\n",
    "        rtmed <- c(rtmed, rtmd)\n",
    "        #rt mean\n",
    "        rtmn <- mean(spsrt$rtime, na.rm = TRUE)\n",
    "        rtmean <- c(rtmean, rtmn)\n",
    "        #collision energy\n",
    "        ce <- max(spsrt$collisionEnergy)\n",
    "        col_eng <- c(col_eng, ce)\n",
    "        #polarity\n",
    "        pl <- max(spsrt$polarity)\n",
    "        if (pl == 1){\n",
    "            px <- 'pos'\n",
    "            pol <- c(pol, px)\n",
    "        }\n",
    "        else {\n",
    "            px <- 'neg'\n",
    "            pol <- c(pol, px)\n",
    "        }\n",
    "        #int\n",
    "        ints <- max(spsrt$precursorIntensity)\n",
    "        int <- c(int, ints)\n",
    "        #mzmlfile\n",
    "        source_file <- c(source_file, mzml_file)\n",
    "        # after all the common infromation is stored,\n",
    "        # move to extracting matching candidates with input spectra\n",
    "        pre_mzs[[x]] <- future({\n",
    "            sps <- spec2_Processing(x, sps_all, spec = \"spec_all\")\n",
    "            ####-------------------------------------------------------------\n",
    "            #### Dereplication with all or GNPS ----\n",
    "            ####-------------------------------------------------------------\n",
    "            # define variables for result dataframe\n",
    "\n",
    "            # if the database selected is GNPS or all\n",
    "            f_gnps <- future(\n",
    "            if (db == \"all\" || db ==\"gnps\"){\n",
    "                GNPSmax_similarity <- c() # dot product score\n",
    "                GNPSmzScore <- c() # similar m/z score\n",
    "                GNPSintScore <- c() # similar int score\n",
    "                GQMatchingPeaks <- c() # matching peaks between gnps candidate and input spectra\n",
    "                GNPSTotalPeaks <- c() # total peaks in gnps candidate\n",
    "                gQueryTotalPeaks<- c() # total peaks in input spectra\n",
    "                GNPSSMILES <- c() # smiles of gnps candidate\n",
    "                #GNPSspectrumID <- c() # spectrum id of gnps candidate\n",
    "                GNPScompound_name <- c() # compound name of gnps candidate\n",
    "                #GNPSmirrorSpec <- c() # path for mirror spectra between gnps candidate and input of gnps candidate\n",
    "                Source <- c() # GNPS as source of result\n",
    "                #### GNPS spec with pre_mz\n",
    "                gnps_with_mz <- spec2_Processing(x, gnpsdb, spec = \"gnps\", ppmx) # change here later\n",
    "\n",
    "                # define the directoyr name to store all GNPS results\n",
    "                dir_name <- paste(input_dir, str_remove(paste(result_dir, \"/spectral_dereplication/GNPS/\", sep = \"\"), \".\"), sep =\"\")\n",
    "                if (!file.exists(dir_name)){\n",
    "                    dir.create(dir_name, recursive = TRUE)\n",
    "                }\n",
    "                if (length(sps) != 0 && length(gnps_with_mz) !=0){\n",
    "                     #' Compare experimental spectra against GNPS\n",
    "                    res <- compareSpectra(sps, gnps_with_mz, ppm = 15, FUN = MsCoreUtils::gnps, MAPFUN = joinPeaksGnps)\n",
    "\n",
    "                    # first condition for GNPS\n",
    "                    # if more input spectra and more candidates have been extracted from GNPS\n",
    "                    if (length(sps) > 1 && length(gnps_with_mz) >1){\n",
    "                        # given threshold of 0.85 for GNPS, extract top candidates\n",
    "                        res_top <- which(res > res[res>0.85], arr.ind = TRUE)\n",
    "\n",
    "                        # if there are some compounds from GNPS detected\n",
    "                        if (length(res_top) > 0){\n",
    "                            res_topdf <- data.frame(res_top)\n",
    "                            # to store the scores to add to res_topdf\n",
    "                            gnps_scores <- c()\n",
    "\n",
    "                            # for all rows and columns in res_topdf\n",
    "                            for (i in 1:nrow(res_topdf)){\n",
    "\n",
    "                                # store the scores\n",
    "                                gnps_scores <- c(gnps_scores, res[(res_topdf[i, \"row\"]), (res_topdf[i, \"col\"])])\n",
    "                            }\n",
    "                            if (length(gnps_scores)>0){\n",
    "                                # add the score column to res_top\n",
    "                                gnps_res <- cbind(res_top, gnps_scores)\n",
    "                                gnps_res <- data.frame(gnps_res)\n",
    "\n",
    "                                # sort in descending order\n",
    "                                ordered_gnps_res <- gnps_res[order(-gnps_res[,\"gnps_scores\"]),]\n",
    "                                df_ord_gnps_res <- data.frame(ordered_gnps_res)\n",
    "                                if (nrow(df_ord_gnps_res)>no_of_candidates){\n",
    "                                    df_ord_gnps_res <- df_ord_gnps_res[1:no_of_candidates,]\n",
    "                                }\n",
    "                                #for each candidate from GNPS\n",
    "                                for (k in 1:nrow(df_ord_gnps_res)){\n",
    "                                    # take each component from df_ord_gnps_res\n",
    "                                    idv <- df_ord_gnps_res[k,]\n",
    "                                    df_peaklists <- peakdf(gnps_with_mz[idv[[2]]], sps[idv[[1]]], ppmx)\n",
    "\n",
    "                                    if (!(is.null(df_peaklists))){\n",
    "                                        GNPSscore <- idv[1, \"gnps_scores\"]\n",
    "                                        GNPSmax_similarity <- c(GNPSmax_similarity, GNPSscore)\n",
    "\n",
    "                                        GNPSmz <- (nrow(df_peaklists)*2)/(nrow(peaksData(gnps_with_mz[idv[[2]]])[[1]])+nrow(peaksData(sps[idv[[1]]])[[1]]))\n",
    "                                        GNPSmzScore <- c(GNPSmzScore, GNPSmz)\n",
    "\n",
    "                                        GNPSint <- mean(1-(df_peaklists[,\"diff\"]/100))\n",
    "                                        GNPSintScore <- c(GNPSintScore, GNPSint)\n",
    "\n",
    "                                        GQMatPeaks <- nrow(df_peaklists)\n",
    "                                        GQMatchingPeaks <- c(GQMatchingPeaks, GQMatPeaks)\n",
    "\n",
    "                                        GNPSTPeaks <- nrow(peaksData(gnps_with_mz[idv[[2]]])[[1]])\n",
    "                                        GNPSTotalPeaks <- c(GNPSTotalPeaks, GNPSTPeaks)\n",
    "\n",
    "                                        gQTPeaks<- nrow(peaksData(sps[idv[[1]]])[[1]])\n",
    "                                        gQueryTotalPeaks <- c(gQueryTotalPeaks, gQTPeaks)\n",
    "\n",
    "                                        GNPS_SMILES <- gnps_with_mz[idv[[2]]]$smiles\n",
    "                                        GNPSSMILES <- c(GNPSSMILES, GNPS_SMILES)\n",
    "\n",
    "                                        #GNPSID <- gnps_with_mz[idv[[2]]]$SPECTRUMID\n",
    "                                        #GNPSspectrumID <- c(GNPSspectrumID, GNPSID)\n",
    "\n",
    "                                        GNPSname <- gnps_with_mz[idv[[2]]]$name\n",
    "                                        GNPScompound_name <- c(GNPScompound_name, GNPSname)\n",
    "                                        Src <- \"GNPS\"\n",
    "                                        Source <- c(Source, Src)\n",
    "\n",
    "                                    }# if df_peaklists isnt empty\n",
    "                                }# for each candidate\n",
    "\n",
    "                            }# gnps_score exists\n",
    "                        }# if res_top has some good candidates\n",
    "                    }# first condition\n",
    "                    # if only one sepctrum from input and more candidates from GNPS\n",
    "                    else if (length(sps) == 1 && length(gnps_with_mz) >1){\n",
    "                        # given threshold of 0.85 for GNPS, extract top candidates\n",
    "                        res_top <- which(res > res[res>0.85], arr.ind = TRUE)\n",
    "                        # if there are candidates with good score\n",
    "                        if (length(res_top) > 0){\n",
    "                            res_topdf <- data.frame(res_top)\n",
    "\n",
    "                            # top store the scores to add to res_topdf\n",
    "                            gnps_scores <- c()\n",
    "                            # for all rows and columns in res_topdf\n",
    "\n",
    "                            for (i in 1:nrow(res_topdf)){\n",
    "                                # store the scores\n",
    "                                gnps_scores <- c(gnps_scores, res[(res_topdf[i, \"res_topdf\"])])\n",
    "                            }\n",
    "                            if (length(gnps_scores)>0){\n",
    "                                # add the score column to res_top\n",
    "                                gnps_res <- cbind(res_top, gnps_scores)\n",
    "                                gnps_res <- data.frame(gnps_res)\n",
    "\n",
    "                                # sort in descending order\n",
    "                                ordered_gnps_res <- gnps_res[order(-gnps_res[,\"gnps_scores\"]),]\n",
    "\n",
    "                                df_ord_gnps_res <- data.frame(ordered_gnps_res)\n",
    "                                if (nrow(df_ord_gnps_res)>no_of_candidates){\n",
    "                                    df_ord_gnps_res <- df_ord_gnps_res[1:no_of_candidates,]\n",
    "                                }\n",
    "                                # for each candidate\n",
    "                                for (k in 1:nrow(df_ord_gnps_res)){\n",
    "                                    # take each candidate\n",
    "                                    idv <- df_ord_gnps_res[k,]\n",
    "                                    df_peaklists <- peakdf(gnps_with_mz[idv[[1]]], sps, ppmx)\n",
    "\n",
    "                                    # if there are matchingpeaks\n",
    "                                    if (!(is.null(df_peaklists))){\n",
    "                                        GNPSscore <- idv[1, \"gnps_scores\"]\n",
    "                                        GNPSmax_similarity <- c(GNPSmax_similarity, GNPSscore)\n",
    "\n",
    "                                        GNPSmz <- (nrow(df_peaklists)*2)/(nrow(peaksData(gnps_with_mz[idv[[1]]])[[1]])+nrow(peaksData(sps)[[1]]))\n",
    "                                        GNPSmzScore <- c(GNPSmzScore, GNPSmz)\n",
    "\n",
    "                                        GNPSint <- mean(1-(df_peaklists[,\"diff\"]/100))\n",
    "                                        GNPSintScore <- c(GNPSintScore, GNPSint)\n",
    "\n",
    "                                        GQMatPeaks <- nrow(df_peaklists)\n",
    "                                        GQMatchingPeaks <- c(GQMatchingPeaks, GQMatPeaks)\n",
    "\n",
    "                                        GNPSTPeaks <- nrow(peaksData(gnps_with_mz[idv[[1]]])[[1]])\n",
    "                                        GNPSTotalPeaks <- c(GNPSTotalPeaks, GNPSTPeaks)\n",
    "\n",
    "                                        gQTPeaks<- nrow(peaksData(sps)[[1]])\n",
    "                                        gQueryTotalPeaks <- c(gQueryTotalPeaks, gQTPeaks)\n",
    "\n",
    "                                        GNPS_SMILES <- gnps_with_mz[idv[[1]]]$smiles\n",
    "                                        GNPSSMILES <- c(GNPSSMILES, GNPS_SMILES)\n",
    "\n",
    "                                        #GNPSID <- gnps_with_mz[idv[[1]]]$SPECTRUMID\n",
    "                                        #GNPSspectrumID <- c(GNPSspectrumID, GNPSID)\n",
    "\n",
    "                                        GNPSname <- gnps_with_mz[idv[[1]]]$name\n",
    "                                        GNPScompound_name <- c(GNPScompound_name, GNPSname)\n",
    "\n",
    "                                        Src <- \"GNPS\"\n",
    "                                        Source <- c(Source, Src)\n",
    "                                    }# if df_peaklists isnt empty\n",
    "                                }# for each candidate\n",
    "                            }# gnps_score exists\n",
    "                        }# if res_top has some good candidates\n",
    "                    }# second condition\n",
    "                    # if there are more input spectra and one candidate from GNPS\n",
    "                    else if (length(sps) > 1 && length(gnps_with_mz) == 1){\n",
    "                        # given threshold of 0.85 for GNPS, extract top candidates\n",
    "                        res_top <- which(res > res[res>0.85], arr.ind = TRUE)\n",
    "\n",
    "                        # if there are good matching candidates\n",
    "                        if (length(res_top) > 0){\n",
    "                            res_topdf <- data.frame(res_top)\n",
    "\n",
    "                            # top store the scores to add to res_topdf\n",
    "                            gnps_scores <- c()\n",
    "                            # for all rows and columns in res_topdf\n",
    "\n",
    "                            # for all candidates\n",
    "                            for (i in 1:nrow(res_topdf)){\n",
    "                                # store the scores\n",
    "                                gnps_scores <- c(gnps_scores, res[(res_topdf[i, \"res_topdf\"])])\n",
    "                            }\n",
    "                            if (length(gnps_scores)>0){\n",
    "                                # add the score column to res_top\n",
    "                                gnps_res <- cbind(res_top, gnps_scores)\n",
    "                                gnps_res <- data.frame(gnps_res)\n",
    "\n",
    "                                # sort in descending order\n",
    "                                ordered_gnps_res <- gnps_res[order(-gnps_res[,\"gnps_scores\"]),]\n",
    "\n",
    "                                df_ord_gnps_res <- data.frame(ordered_gnps_res)\n",
    "                                if (nrow(df_ord_gnps_res)>no_of_candidates){\n",
    "                                    df_ord_gnps_res <- df_ord_gnps_res[1:no_of_candidates,]\n",
    "                                }\n",
    "                                # for each candidate match\n",
    "                                for (k in 1:nrow(df_ord_gnps_res)){\n",
    "                                    # take each candidate\n",
    "                                    idv <- df_ord_gnps_res[k,]\n",
    "\n",
    "                                    df_peaklists <- peakdf(gnps_with_mz, sps[idv[[1]]], ppmx)\n",
    "\n",
    "                                    # if there are matching peaks\n",
    "                                    if (!(is.null(df_peaklists))){\n",
    "                                        GNPSscore <- idv[1, \"gnps_scores\"]\n",
    "                                        GNPSmax_similarity <- c(GNPSmax_similarity, GNPSscore)\n",
    "\n",
    "                                        GNPSmz <- (nrow(df_peaklists)*2)/(nrow(peaksData(gnps_with_mz)[[1]])+nrow(peaksData(sps[idv[[1]]])[[1]]))\n",
    "                                        GNPSmzScore <- c(GNPSmzScore, GNPSmz)\n",
    "\n",
    "                                        GNPSint <- mean(1-(df_peaklists[,\"diff\"]/100))\n",
    "                                        GNPSintScore <- c(GNPSintScore, GNPSint)\n",
    "\n",
    "                                        GQMatPeaks <- nrow(df_peaklists)\n",
    "                                        GQMatchingPeaks <- c(GQMatchingPeaks, GQMatPeaks)\n",
    "\n",
    "                                        GNPSTPeaks <- nrow(peaksData(gnps_with_mz)[[1]])\n",
    "                                        GNPSTotalPeaks <- c(GNPSTotalPeaks, GNPSTPeaks)\n",
    "\n",
    "                                        gQTPeaks<- nrow(peaksData(sps[idv[[1]]])[[1]])\n",
    "                                        gQueryTotalPeaks <- c(gQueryTotalPeaks, gQTPeaks)\n",
    "\n",
    "                                        GNPS_SMILES <- gnps_with_mz$smiles\n",
    "                                        GNPSSMILES <- c(GNPSSMILES, GNPS_SMILES)\n",
    "\n",
    "                                        #GNPSID <- gnps_with_mz$SPECTRUMID\n",
    "                                        #GNPSspectrumID <- c(GNPSspectrumID, GNPSID)\n",
    "\n",
    "                                        GNPSname <- gnps_with_mz$name\n",
    "                                        GNPScompound_name <- c(GNPScompound_name, GNPSname)\n",
    "\n",
    "                                        Src <- \"GNPS\"\n",
    "                                        Source <- c(Source, Src)\n",
    "                                    }# if df_peaklists isnt empty\n",
    "                                }# for each candidate\n",
    "                            }# gnps_score exists \n",
    "                        }# if res_top has some good candidates\n",
    "                    }# third condition\n",
    "                    else if (length(sps) == 1 && length(gnps_with_mz) == 1){\n",
    "                        if (res>= 0.85){\n",
    "                            #take that one candidate\n",
    "                            gnps_best_match <- gnps_with_mz\n",
    "\n",
    "                            df_peaklists <- peakdf(gnps_best_match, sps, ppmx)\n",
    "\n",
    "                            # if there are matching peaks\n",
    "                            if (!(is.null(df_peaklists))){\n",
    "\n",
    "                                GNPSscore <- max(res)\n",
    "                                GNPSmax_similarity <- c(GNPSmax_similarity, GNPSscore)\n",
    "\n",
    "\n",
    "                                GNPSmz <- (nrow(df_peaklists)*2)/(nrow(peaksData(gnps_best_match)[[1]])+nrow(peaksData(sps)[[1]]))\n",
    "                                GNPSmzScore <- c(GNPSmzScore, GNPSmz)\n",
    "\n",
    "                                GNPSint <- mean(1-(df_peaklists[,\"diff\"]/100))\n",
    "                                GNPSintScore <- c(GNPSintScore, GNPSint)\n",
    "\n",
    "                                GQMatPeaks <- NA\n",
    "                                GQMatchingPeaks <- c(GQMatchingPeaks, GQMatPeaks)\n",
    "\n",
    "\n",
    "                                GNPSTPeaks <- nrow(peaksData(gnps_best_match)[[1]])\n",
    "                                GNPSTotalPeaks <- c(GNPSTotalPeaks, GNPSTPeaks)\n",
    "\n",
    "                                gQTPeaks<- nrow(peaksData(sps)[[1]])\n",
    "                                gQueryTotalPeaks <- c(gQueryTotalPeaks, gQTPeaks)\n",
    "\n",
    "\n",
    "                                GNPS_SMILES <- gnps_best_match$smiles\n",
    "                                GNPSSMILES <- c(GNPSSMILES, GNPS_SMILES)\n",
    "\n",
    "                                GNPSname <- gnps_best_match$name\n",
    "                                GNPScompound_name <- c(GNPScompound_name, GNPSname)\n",
    "\n",
    "\n",
    "                                #GNPSID <- gnps_best_match$SPECTRUMID\n",
    "                                #GNPSspectrumID <- c(GNPSspectrumID, GNPSID)\n",
    "\n",
    "\n",
    "                                Src <- \"GNPS\"\n",
    "                                Source <- c(Source, Src)\n",
    "                            }# if candidate exixts\n",
    "                        }#  if res_top has some good candidates\n",
    "\n",
    "                    }# fourth condition\n",
    "                }# if sps and gnps has some candidate\n",
    "                gnps_x <- data.frame(cbind(GNPSmax_similarity, GNPSmzScore,\n",
    "                                       GNPSintScore, GQMatchingPeaks,\n",
    "                                       GNPSTotalPeaks, gQueryTotalPeaks,\n",
    "                                       GNPSSMILES, GNPScompound_name, Source))\n",
    "                write.csv(gnps_x, file = paste(dir_name, \"/gnps_results_for_\", id_Xx, \".csv\", sep = \"\"))\n",
    "            })# gnps ends\n",
    "\n",
    "\n",
    "            ####-------------------------------------------------------------\n",
    "            #### Dereplication with all or HMDB ----\n",
    "            ####-------------------------------------------------------------\n",
    "             # if the database selected is HMDB or all\n",
    "            f_hmdb <- future(\n",
    "            if (db == \"all\" || db ==\"hmdb\"){\n",
    "                # hmdb\n",
    "                HMDBmax_similarity <- c()\n",
    "                HMDBmzScore <- c()\n",
    "                HMDBintScore <- c()\n",
    "                HQMatchingPeaks <- c()\n",
    "                HMDBTotalPeaks <- c()\n",
    "                hQueryTotalPeaks<- c()\n",
    "                HMDBcompoundID <- c()\n",
    "                Source <- c()\n",
    "\n",
    "                #### HMDB spec with pre_mz\n",
    "                hmdb_with_mz <- spec2_Processing(x, hmdb, spec = \"hmdb\", ppmx) # change here later\n",
    "\n",
    "                # directory name for HMDB results\n",
    "                dir_name <- paste(input_dir, str_remove(paste(result_dir, \"/spectral_dereplication/HMDB/\", sep = \"\"), \".\"), sep =\"\")\n",
    "                if (!file.exists(dir_name)){\n",
    "                    dir.create(dir_name, recursive = TRUE)\n",
    "                }\n",
    "                if (length(sps) != 0 && length(hmdb_with_mz) !=0){\n",
    "                    #' Compare experimental spectra against HMDB\n",
    "                    res <- compareSpectra(sps, hmdb_with_mz, ppm = 15)\n",
    "\n",
    "                    # if there are more input spectra and more candidates from GNPS\n",
    "                    if (length(sps) > 1 && length(hmdb_with_mz) >1){\n",
    "                        # given threshold of 0.70 for HMDB, extract top candidates\n",
    "                        res_top <- which(res > res[res>0.70], arr.ind = TRUE)\n",
    "                        if (length(res_top) > 0){\n",
    "                            res_topdf <- data.frame(res_top)\n",
    "                            # to store the scores to add to res_topdf\n",
    "                            hmdb_scores <- c()\n",
    "                            # for all rows and columns in res_topdf\n",
    "                            for (i in 1:nrow(res_topdf)){\n",
    "\n",
    "                                # store the scores\n",
    "                                hmdb_scores <- c(hmdb_scores, res[(res_topdf[i, \"row\"]), (res_topdf[i, \"col\"])])\n",
    "                            }\n",
    "                            if (length(hmdb_scores)>0){\n",
    "                                # add the score column to res_top\n",
    "                                hmdb_res <- cbind(res_top, hmdb_scores)\n",
    "                                hmdb_res <- data.frame(hmdb_res)\n",
    "                                # sort in descending order\n",
    "                                ordered_hmdb_res <- hmdb_res[order(-hmdb_res[,\"hmdb_scores\"]),]\n",
    "                                df_ord_hmdb_res <- data.frame(ordered_hmdb_res)\n",
    "                                if (nrow(df_ord_hmdb_res)>no_of_candidates){\n",
    "                                    df_ord_hmdb_res <- df_ord_hmdb_res[1:no_of_candidates,]\n",
    "                                }\n",
    "                                for (k in 1:nrow(df_ord_hmdb_res)){\n",
    "                                    idv <- df_ord_hmdb_res[k,]\n",
    "                                    df_peaklists <- peakdf(hmdb_with_mz[idv[[2]]], sps[idv[[1]]], ppmx)\n",
    "                                    if (!(is.null(df_peaklists))){\n",
    "                                        HMDBscore <- idv[1, \"hmdb_scores\"]\n",
    "                                        HMDBmax_similarity <- c(HMDBmax_similarity, HMDBscore)\n",
    "\n",
    "                                        HMDBmz <- (nrow(df_peaklists)*2)/(nrow(peaksData(hmdb_with_mz[idv[[2]]])[[1]])+nrow(peaksData(sps[idv[[1]]])[[1]]))\n",
    "                                        HMDBmzScore <- c(HMDBmzScore, HMDBmz)\n",
    "\n",
    "                                        HMDBint <- mean(1-(df_peaklists[,\"diff\"]/100))\n",
    "                                        HMDBintScore <- c(HMDBintScore, HMDBint)\n",
    "\n",
    "\n",
    "                                        HQMatPeaks <- nrow(df_peaklists)\n",
    "                                        HQMatchingPeaks <- c(HQMatchingPeaks, HQMatPeaks)\n",
    "\n",
    "\n",
    "                                        HMDBTPeaks <- nrow(peaksData(hmdb_with_mz[idv[[2]]])[[1]])\n",
    "                                        HMDBTotalPeaks <- c(HMDBTotalPeaks, HMDBTPeaks)\n",
    "\n",
    "                                        hQTPeaks<- nrow(peaksData(sps[idv[[1]]])[[1]])\n",
    "                                        hQueryTotalPeaks<- c(hQueryTotalPeaks, hQTPeaks)\n",
    "\n",
    "\n",
    "                                        HMDBID <- hmdb_with_mz[idv[[2]]]$compound_id\n",
    "                                        HMDBcompoundID <- c(HMDBcompoundID, HMDBID)\n",
    "\n",
    "                                        Src <- \"HMDB\"\n",
    "                                        Source <- c(Source, Src)\n",
    "                                    }# if df_peaklists is not empty\n",
    "                                }# for each candidate\n",
    "                            }# if hmdb_scores exist\n",
    "                        }# if there are top candidadates with good scores\n",
    "                    }#first condition ends\n",
    "                    else if (length(sps) == 1 && length(hmdb_with_mz) >1){\n",
    "                        # given threshold of 0.70 for HMDB, extract top candidates\n",
    "                        res_top <- which(res > res[res>0.70], arr.ind = TRUE)\n",
    "                        if (length(res_top) > 0){\n",
    "                            res_topdf <- data.frame(res_top)\n",
    "\n",
    "                            # top store the scores to add to res_topdf\n",
    "                            hmdb_scores <- c()\n",
    "                            # for all rows and columns in res_topdf\n",
    "                            for (i in 1:nrow(res_topdf)){\n",
    "                                # store the scores\n",
    "                                hmdb_scores <- c(hmdb_scores, res[(res_topdf[i, \"res_topdf\"])])\n",
    "                            }\n",
    "                            if (length(hmdb_scores > 0)){\n",
    "                                # add the score column to res_top\n",
    "                                hmdb_res <- cbind(res_top, hmdb_scores)\n",
    "                                hmdb_res <- data.frame(hmdb_res)\n",
    "                                # sort in descending order\n",
    "                                ordered_hmdb_res <- hmdb_res[order(-hmdb_res[,\"hmdb_scores\"]),]\n",
    "\n",
    "                                df_ord_hmdb_res <- data.frame(ordered_hmdb_res)\n",
    "                                if (nrow(df_ord_hmdb_res)>no_of_candidates){\n",
    "                                    df_ord_hmdb_res <- df_ord_hmdb_res[1:no_of_candidates,]\n",
    "                                }\n",
    "                                for (k in 1:nrow(df_ord_hmdb_res)){\n",
    "                                    idv <- df_ord_hmdb_res[k,]\n",
    "                                    df_peaklists <- peakdf(hmdb_with_mz[idv[[1]]], sps, ppmx)\n",
    "                                    if (!(is.null(df_peaklists))){\n",
    "                                        HMDBscore <- idv[1, \"hmdb_scores\"]\n",
    "                                        HMDBmax_similarity <- c(HMDBmax_similarity, HMDBscore)\n",
    "\n",
    "                                        HMDBmz <- (nrow(df_peaklists)*2)/(nrow(peaksData(hmdb_with_mz[idv[[1]]])[[1]])+nrow(peaksData(sps)[[1]]))\n",
    "                                        HMDBmzScore <- c(HMDBmzScore, HMDBmz)\n",
    "\n",
    "                                        HMDBint <- mean(1-(df_peaklists[,\"diff\"]/100))\n",
    "                                        HMDBintScore <- c(HMDBintScore, HMDBint)\n",
    "\n",
    "\n",
    "                                        HQMatPeaks <- nrow(df_peaklists)\n",
    "                                        HQMatchingPeaks <- c(HQMatchingPeaks, HQMatPeaks)\n",
    "\n",
    "\n",
    "                                        HMDBTPeaks <- nrow(peaksData(hmdb_with_mz[idv[[1]]])[[1]])\n",
    "                                        HMDBTotalPeaks <- c(HMDBTotalPeaks, HMDBTPeaks)\n",
    "\n",
    "                                        hQTPeaks<- nrow(peaksData(sps)[[1]])\n",
    "                                        hQueryTotalPeaks<- c(hQueryTotalPeaks, hQTPeaks)\n",
    "\n",
    "\n",
    "                                        HMDBID <- hmdb_with_mz[idv[[1]]]$compound_id\n",
    "                                        HMDBcompoundID <- c(HMDBcompoundID, HMDBID)\n",
    "\n",
    "                                        Src <- \"HMDB\"\n",
    "                                        Source <- c(Source, Src)\n",
    "                                    }#df_peaklists is not null\n",
    "                                }# for each candidate\n",
    "                            }# hmdb_scores exists\n",
    "                        }\n",
    "                    }# second condition ends\n",
    "                    else if (length(sps) > 1 && length(hmdb_with_mz) == 1){\n",
    "                        # given threshold of 0.70 for HMDB, extract top candidates\n",
    "                        res_top <- which(res > res[res>0.70], arr.ind = TRUE)\n",
    "                        if (length(res_top) > 0){\n",
    "                            res_topdf <- data.frame(res_top)\n",
    "                            # top store the scores to add to res_topdf\n",
    "                            hmdb_scores <- c()\n",
    "                            # for all rows and columns in res_topdf\n",
    "\n",
    "                            for (i in 1:nrow(res_topdf)){\n",
    "                                # store the scores\n",
    "                                hmdb_scores <- c(hmdb_scores, res[(res_topdf[i, \"res_topdf\"])])\n",
    "                            }\n",
    "                            if (length(hmdb_scores)>0){\n",
    "                                # add the score column to res_top\n",
    "                                hmdb_res <- cbind(res_top, hmdb_scores)\n",
    "                                hmdb_res <- data.frame(hmdb_res)\n",
    "\n",
    "                                # sort in descending order\n",
    "                                ordered_hmdb_res <- hmdb_res[order(-hmdb_res[,\"hmdb_scores\"]),]\n",
    "\n",
    "                                df_ord_hmdb_res <- data.frame(ordered_hmdb_res)\n",
    "                                if (nrow(df_ord_hmdb_res)>no_of_candidates){\n",
    "                                    df_ord_hmdb_res <- df_ord_hmdb_res[1:no_of_candidates,]\n",
    "                                }\n",
    "                                for (k in 1:nrow(df_ord_hmdb_res)){\n",
    "                                    idv <- df_ord_hmdb_res[k,]\n",
    "                                    df_peaklists <- peakdf(hmdb_with_mz, sps[idv[[1]]], ppmx)\n",
    "                                    if (!(is.null(df_peaklists))){\n",
    "                                        HMDBscore <- idv[1, \"hmdb_scores\"]\n",
    "                                        HMDBmax_similarity <- c(HMDBmax_similarity, HMDBscore)\n",
    "\n",
    "                                        HMDBmz <- (nrow(df_peaklists)*2)/(nrow(peaksData(hmdb_with_mz)[[1]])+nrow(peaksData(sps[idv[[1]]])[[1]]))\n",
    "                                        HMDBmzScore <- c(HMDBmzScore, HMDBmz)\n",
    "\n",
    "                                        HMDBint <- mean(1-(df_peaklists[,\"diff\"]/100))\n",
    "                                        HMDBintScore <- c(HMDBintScore, HMDBint)\n",
    "\n",
    "\n",
    "                                        HQMatPeaks <- nrow(df_peaklists)\n",
    "                                        HQMatchingPeaks <- c(HQMatchingPeaks, HQMatPeaks)\n",
    "\n",
    "\n",
    "                                        HMDBTPeaks <- nrow(peaksData(hmdb_with_mz)[[1]])\n",
    "                                        HMDBTotalPeaks <- c(HMDBTotalPeaks, HMDBTPeaks)\n",
    "\n",
    "                                        hQTPeaks<- nrow(peaksData(sps[idv[[1]]])[[1]])\n",
    "                                        hQueryTotalPeaks<- c(hQueryTotalPeaks, hQTPeaks)\n",
    "\n",
    "\n",
    "                                        HMDBID <- hmdb_with_mz$compound_id\n",
    "                                        HMDBcompoundID <- c(HMDBcompoundID, HMDBID)\n",
    "\n",
    "                                        Src <- \"HMDB\"\n",
    "                                        Source <- c(Source, Src)\n",
    "                                    }# if df_peaklists isn't null\n",
    "                                }# for each candidate\n",
    "                            }#hmdb_score exists\n",
    "                        }\n",
    "                    }#third condition ends\n",
    "                    else if (length(sps) == 1 && length(hmdb_with_mz) == 1){\n",
    "                        if (res>=0.70){\n",
    "                            hmdb_best_match <- hmdb_with_mz\n",
    "                            df_peaklists <- peakdf(hmdb_best_match, sps, ppmx)\n",
    "                            if (!(is.null(df_peaklists))){\n",
    "                                HMDBscore <- max(res)\n",
    "                                HMDBmax_similarity <- c(HMDBmax_similarity, HMDBscore)\n",
    "\n",
    "                                HMDBmz <- (nrow(df_peaklists)*2)/(nrow(peaksData(hmdb_best_match)[[1]])+nrow(peaksData(sps)[[1]]))\n",
    "                                HMDBmzScore <- c(HMDBmzScore, HMDBmz)\n",
    "\n",
    "                                HMDBint <- mean(1-(df_peaklists[,\"diff\"]/100))\n",
    "                                HMDBintScore <- c(HMDBintScore, HMDBint)\n",
    "\n",
    "\n",
    "                                HQMatPeaks <- nrow(df_peaklists)\n",
    "                                HQMatchingPeaks <- c(HQMatchingPeaks, HQMatPeaks)\n",
    "\n",
    "\n",
    "                                HMDBTPeaks <- nrow(peaksData(hmdb_best_match)[[1]])\n",
    "                                HMDBTotalPeaks <- c(HMDBTotalPeaks, HMDBTPeaks)\n",
    "\n",
    "                                hQTPeaks<- nrow(peaksData(sps)[[1]])\n",
    "                                hQueryTotalPeaks<- c(hQueryTotalPeaks, hQTPeaks)\n",
    "\n",
    "\n",
    "                                HMDBID <- hmdb_best_match$compound_id\n",
    "                                HMDBcompoundID <- c(HMDBcompoundID, HMDBID)\n",
    "                                 Src <- \"HMDB\"\n",
    "                                Source <- c(Source, Src)\n",
    "                            }\n",
    "                        }\n",
    "                    }# fourth condition\n",
    "                }# if both sps and hmdb has some matching candidates\n",
    "                hmdb_x <- data.frame(cbind(HMDBmax_similarity, HMDBmzScore,\n",
    "                                       HMDBintScore, HQMatchingPeaks,\n",
    "                                       HMDBTotalPeaks, hQueryTotalPeaks,\n",
    "                                       HMDBcompoundID, Source))\n",
    "                write.csv(hmdb_x, file = paste(dir_name, \"/hmdb_results_for_\", id_Xx, \".csv\", sep = \"\"))\n",
    "            })#hmdb ends here\n",
    "\n",
    "\n",
    "\n",
    "            ####-------------------------------------------------------------\n",
    "            #### Dereplication with all or MassBank ----\n",
    "            ####-------------------------------------------------------------\n",
    "            # define variables for result dataframe\n",
    "            f_mbank <- future(\n",
    "                # if the database selected is MassBank or all\n",
    "                if (db == \"all\" || db ==\"mbank\"){\n",
    "                    # mbank\n",
    "                    MBmax_similarity <- c()\n",
    "                    MBmzScore <- c()\n",
    "                    MBintScore <- c()\n",
    "                    MQMatchingPeaks <- c()\n",
    "                    MBTotalPeaks <- c()\n",
    "                    mQueryTotalPeaks<- c()\n",
    "                    MBformula <- c()\n",
    "                    MBSMILES <- c()\n",
    "                    MBspectrumID <- c()\n",
    "                    MBcompound_name <- c()\n",
    "                    Source <- c()\n",
    "\n",
    "                    mbank_with_mz <- spec2_Processing(x, mbank, spec = \"mbank\", ppmx) # change here later\n",
    "\n",
    "                    dir_name <- paste(input_dir, str_remove(paste(result_dir, \"/spectral_dereplication/MassBank/\", sep = \"\"), \".\"), sep =\"\")\n",
    "                    if (!file.exists(dir_name)){\n",
    "                        dir.create(dir_name, recursive = TRUE)\n",
    "                    }\n",
    "                    if (length(sps) != 0 && length(mbank_with_mz) !=0){\n",
    "                        #' Compare experimental spectra against MassBank\n",
    "                        res <- compareSpectra(sps, mbank_with_mz, ppm = 15)\n",
    "                        if (length(sps) > 1 && length(mbank_with_mz) >1){\n",
    "                            # given threshold of 0.70 for MassBank, extract top candidates\n",
    "                            res_top <- which(res > res[res>0.70], arr.ind = TRUE)\n",
    "                            if (length(res_top) > 0){\n",
    "                                res_topdf <- data.frame(res_top)\n",
    "                                # to store the scores to add to res_topdf\n",
    "                                mbank_scores <- c()\n",
    "                                # for all rows and columns in res_topdf\n",
    "                                for (i in 1:nrow(res_topdf)){\n",
    "                                    # store the scores\n",
    "                                    mbank_scores <- c(mbank_scores, res[(res_topdf[i, \"row\"]), (res_topdf[i, \"col\"])])\n",
    "                                }\n",
    "                                if (length(mbank_scores)>0){\n",
    "                                    # add the score column to res_top\n",
    "                                    mbank_res <- cbind(res_top, mbank_scores)\n",
    "                                    mbank_res <- data.frame(mbank_res)\n",
    "                                    # sort in descending order\n",
    "                                    ordered_mbank_res <- mbank_res[order(-mbank_res[,\"mbank_scores\"]),]\n",
    "                                    df_ord_mbank_res <- data.frame(ordered_mbank_res)\n",
    "                                    if (nrow(df_ord_mbank_res)>no_of_candidates){\n",
    "                                        df_ord_mbank_res <- df_ord_mbank_res[1:no_of_candidates,]\n",
    "                                    }\n",
    "                                    for (k in 1:nrow(df_ord_mbank_res)){\n",
    "                                        idv <- df_ord_mbank_res[k,]\n",
    "                                        df_peaklists <- peakdf(mbank_with_mz[idv[[2]]], sps[idv[[1]]], ppmx)\n",
    "                                        if (!(is.null(df_peaklists))){\n",
    "                                            mbscore <- idv[1, \"mbank_scores\"]\n",
    "                                            MBmax_similarity<- c(MBmax_similarity, mbscore)\n",
    "\n",
    "                                            MBmz <- (nrow(df_peaklists)*2)/(nrow(peaksData(mbank_with_mz[idv[[2]]])[[1]])+nrow(peaksData(sps[idv[[1]]])[[1]]))\n",
    "                                            MBmzScore <- c(MBmzScore, MBmz)\n",
    "\n",
    "                                            MBint <- mean(1-(df_peaklists[,\"diff\"]/100))\n",
    "                                            MBintScore <- c(MBintScore, MBint)\n",
    "\n",
    "                                            MQMatPeaks <- nrow(df_peaklists)\n",
    "                                            MQMatchingPeaks <- c(MQMatchingPeaks, MQMatPeaks)\n",
    "\n",
    "                                            MBTPeaks <- nrow(peaksData(mbank_with_mz[idv[[2]]])[[1]])\n",
    "                                            MBTotalPeaks<- c(MBTotalPeaks, MBTPeaks)\n",
    "\n",
    "                                            mQTPeaks<- nrow(peaksData(sps[idv[[1]]])[[1]])\n",
    "                                            mQueryTotalPeaks<- c(mQueryTotalPeaks, mQTPeaks)\n",
    "\n",
    "\n",
    "                                            MBfor <- mbank_with_mz[idv[[2]]]$Formula\n",
    "                                            MBformula<- c(MBformula, MBfor)\n",
    "\n",
    "                                            MBS <- mbank_with_mz[idv[[2]]]$smiles\n",
    "                                            MBSMILES <- c(MBSMILES, MBS)\n",
    "\n",
    "                                            MBID <- mbank_with_mz[idv[[2]]]$accession\n",
    "                                            MBspectrumID<- c(MBspectrumID, MBID)\n",
    "\n",
    "                                            MBname <- mbank_with_mz[idv[[2]]]$Name\n",
    "                                            MBcompound_name <- c(MBcompound_name, MBname)\n",
    "\n",
    "                                            Src <- \"MassBank\"\n",
    "                                            Source <- c(Source, Src)\n",
    "                                        }# if df_peaklists is not empty\n",
    "                                    }# for each candidate\n",
    "                                }# if mbank_score exists\n",
    "                            }# if there are in candidates in res_top\n",
    "                        }# first condition\n",
    "                        else if (length(sps) == 1 && length(mbank_with_mz) >1){\n",
    "                            # given threshold of 0.70 for MassBank, extract top candidates\n",
    "                            res_top <- which(res > res[res>0.70], arr.ind = TRUE)\n",
    "                            if (length(res_top) > 0){\n",
    "                                res_topdf <- data.frame(res_top)\n",
    "\n",
    "                                # top store the scores to add to res_topdf\n",
    "                                mbank_scores <- c()\n",
    "                                # for all rows and columns in res_topdf\n",
    "\n",
    "                                for (i in 1:nrow(res_topdf)){\n",
    "                                    # store the scores\n",
    "                                    mbank_scores <- c(mbank_scores, res[(res_topdf[i, \"res_top\"])])\n",
    "                                }\n",
    "                                if (length(mbank_scores)>0){\n",
    "                                    # add the score column to res_top\n",
    "                                    mbank_res <- cbind(res_top, mbank_scores)\n",
    "                                    mbank_res <- data.frame(mbank_res)\n",
    "\n",
    "                                    # sort in descending order\n",
    "                                    ordered_mbank_res <- mbank_res[order(-mbank_res[,\"mbank_scores\"]),]\n",
    "\n",
    "                                    df_ord_mbank_res <- data.frame(ordered_mbank_res)\n",
    "                                    if (nrow(df_ord_mbank_res)>no_of_candidates){\n",
    "                                        df_ord_mbank_res <- df_ord_mbank_res[1:no_of_candidates,]\n",
    "                                    }\n",
    "\n",
    "                                    for (k in 1:nrow(df_ord_mbank_res)){\n",
    "                                        idv <- df_ord_mbank_res[k,]\n",
    "                                        df_peaklists <- peakdf(mbank_with_mz[idv[[1]]], sps, ppmx)\n",
    "                                        if (!(is.null(df_peaklists))){\n",
    "                                            mbscore <- idv[1, \"mbank_scores\"]\n",
    "                                            MBmax_similarity<- c(MBmax_similarity, mbscore)\n",
    "\n",
    "                                            MBmz <- (nrow(df_peaklists)*2)/(nrow(peaksData(mbank_with_mz[idv[[1]]])[[1]])+nrow(peaksData(sps)[[1]]))\n",
    "                                            MBmzScore <- c(MBmzScore, MBmz)\n",
    "\n",
    "                                            MBint <- mean(1-(df_peaklists[,\"diff\"]/100))\n",
    "                                            MBintScore <- c(MBintScore, MBint)\n",
    "\n",
    "                                            MQMatPeaks <- nrow(df_peaklists)\n",
    "                                            MQMatchingPeaks <- c(MQMatchingPeaks, MQMatPeaks)\n",
    "\n",
    "                                            MBTPeaks <- nrow(peaksData(mbank_with_mz[idv[[1]]])[[1]])\n",
    "                                            MBTotalPeaks<- c(MBTotalPeaks, MBTPeaks)\n",
    "\n",
    "                                            mQTPeaks<- nrow(peaksData(sps)[[1]])\n",
    "                                            mQueryTotalPeaks<- c(mQueryTotalPeaks, mQTPeaks)\n",
    "\n",
    "\n",
    "                                            MBfor <- mbank_with_mz[idv[[1]]]$Formula\n",
    "                                            MBformula<- c(MBformula, MBfor)\n",
    "\n",
    "                                            MBS <-  mbank_with_mz[idv[[1]]]$smiles\n",
    "                                            MBSMILES <- c(MBSMILES, MBS)\n",
    "\n",
    "                                            MBID <- mbank_with_mz[idv[[1]]]$accession\n",
    "                                            MBspectrumID<- c(MBspectrumID, MBID)\n",
    "\n",
    "                                            MBname <- mbank_with_mz[idv[[1]]]$Name\n",
    "                                            MBcompound_name <- c(MBcompound_name, MBname)\n",
    "\n",
    "                                            Src <- \"MassBank\"\n",
    "                                            Source <- c(Source, Src)\n",
    "                                        }# if df_peaklists is not empty\n",
    "                                    }# for each candidate\n",
    "                                }# if mbank_score exists\n",
    "\n",
    "                            }# res_top has some candidates\n",
    "                        }#second condition ends\n",
    "                        else if (length(sps) > 1 && length(mbank_with_mz) == 1){\n",
    "                            # given threshold of 0.70 for MassBank, extract top candidates\n",
    "                            res_top <- which(res > res[res>0.70], arr.ind = TRUE)\n",
    "                            if (length(res_top) > 0){\n",
    "                                res_topdf <- data.frame(res_top)\n",
    "\n",
    "                                # top store the scores to add to res_topdf\n",
    "                                mbank_scores <- c()\n",
    "                                # for all rows and columns in res_topdf\n",
    "\n",
    "                                for (i in 1:nrow(res_topdf)){\n",
    "                                    # store the scores\n",
    "                                    mbank_scores <- c(mbank_scores, res[(res_topdf[i, \"res_top\"])])\n",
    "                                }\n",
    "                                if (length(mbank_scores)>0){\n",
    "                                    # add the score column to res_top\n",
    "                                    mbank_res <- cbind(res_top, mbank_scores)\n",
    "                                    mbank_res <- data.frame(mbank_res)\n",
    "                                    # sort in descending order\n",
    "                                    ordered_mbank_res <- mbank_res[order(-mbank_res[,\"mbank_scores\"]),]\n",
    "\n",
    "                                    df_ord_mbank_res <- data.frame(ordered_mbank_res)\n",
    "                                    if (nrow(df_ord_mbank_res)>no_of_candidates){\n",
    "                                        df_ord_mbank_res <- df_ord_mbank_res[1:no_of_candidates,]\n",
    "                                    }\n",
    "                                    for (k in 1:nrow(df_ord_mbank_res)){\n",
    "                                        idv <- df_ord_mbank_res[k,]\n",
    "                                        df_peaklists <- peakdf(mbank_with_mz, sps[idv[[1]]], ppmx)\n",
    "                                        if (!(is.null(df_peaklists))){\n",
    "                                            mbscore <- idv[1, \"mbank_scores\"]\n",
    "                                            MBmax_similarity<- c(MBmax_similarity, mbscore)\n",
    "\n",
    "                                            MBmz <- (nrow(df_peaklists)*2)/(nrow(peaksData(mbank_with_mz)[[1]])+nrow(peaksData(sps[idv[[1]]])[[1]]))\n",
    "                                            MBmzScore <- c(MBmzScore, MBmz)\n",
    "\n",
    "                                            MBint <- mean(1-(df_peaklists[,\"diff\"]/100))\n",
    "                                            MBintScore <- c(MBintScore, MBint)\n",
    "\n",
    "                                            MQMatPeaks <- nrow(df_peaklists)\n",
    "                                            MQMatchingPeaks <- c(MQMatchingPeaks, MQMatPeaks)\n",
    "\n",
    "                                            MBTPeaks <- nrow(peaksData(mbank_with_mz)[[1]])\n",
    "                                            MBTotalPeaks<- c(MBTotalPeaks, MBTPeaks)\n",
    "\n",
    "                                            mQTPeaks<- nrow(peaksData(sps[idv[[1]]])[[1]])\n",
    "                                            mQueryTotalPeaks<- c(mQueryTotalPeaks, mQTPeaks)\n",
    "\n",
    "\n",
    "                                            MBfor <- mbank_with_mz$Formula\n",
    "                                            MBformula<- c(MBformula, MBfor)\n",
    "\n",
    "                                            MBS <- mbank_with_mz$smiles\n",
    "                                            MBSMILES <- c(MBSMILES, MBS)\n",
    "                                            MBID <- mbank_with_mz$accession\n",
    "                                            MBspectrumID<- c(MBspectrumID, MBID)\n",
    "\n",
    "                                            MBname <- mbank_with_mz$Name\n",
    "                                            MBcompound_name <- c(MBcompound_name, MBname)\n",
    "\n",
    "                                            Src <- \"MassBank\"\n",
    "                                            Source <- c(Source, Src)\n",
    "                                        }# if df_peaklists is not empty\n",
    "                                    }# for each candidate\n",
    "                                }# if mbank_score exists\n",
    "                            }# res_top has some candidates\n",
    "                        }#third condition ends\n",
    "                        else if (length(sps) == 1 && length(mbank_with_mz) == 1){\n",
    "                            if (res>0.70){\n",
    "                                mbank_best_match <- mbank_with_mz\n",
    "                                df_peaklists <- peakdf(mbank_best_match, sps, ppmx)\n",
    "                                if (!(is.null(df_peaklists))){\n",
    "                                    mbscore <- max(res)\n",
    "                                    MBmax_similarity<- c(MBmax_similarity, mbscore)\n",
    "\n",
    "                                    MBmz <- (nrow(df_peaklists)*2)/(nrow(peaksData(mbank_with_mz)[[1]])+nrow(peaksData(sps)[[1]]))\n",
    "                                    MBmzScore <- c(MBmzScore, MBmz)\n",
    "\n",
    "                                    MBint <- mean(1-(df_peaklists[,\"diff\"]/100))\n",
    "                                    MBintScore <- c(MBintScore, MBint)\n",
    "\n",
    "                                    MQMatPeaks <- nrow(df_peaklists)\n",
    "                                    MQMatchingPeaks <- c(MQMatchingPeaks, MQMatPeaks)\n",
    "\n",
    "                                    MBTPeaks <- nrow(peaksData(mbank_with_mz)[[1]])\n",
    "                                    MBTotalPeaks<- c(MBTotalPeaks, MBTPeaks)\n",
    "\n",
    "                                    mQTPeaks<- nrow(peaksData(sps)[[1]])\n",
    "                                    mQueryTotalPeaks<- c(mQueryTotalPeaks, mQTPeaks)\n",
    "\n",
    "\n",
    "                                    MBfor <- mbank_with_mz$Formula\n",
    "                                    MBformula<- c(MBformula, MBfor)\n",
    "\n",
    "                                    MBS <- mbank_with_mz$smiles\n",
    "                                    MBSMILES <- c(MBSMILES, MBS)\n",
    "\n",
    "                                    MBID <- mbank_with_mz$accession\n",
    "                                    MBspectrumID<- c(MBspectrumID, MBID)\n",
    "\n",
    "                                    MBname <- mbank_with_mz$Name\n",
    "                                    MBcompound_name <- c(MBcompound_name, MBname)\n",
    "\n",
    "                                    Src <- \"MassBank\"\n",
    "                                    Source <- c(Source, Src)\n",
    "\n",
    "                                }# if df_peaklists is not empty\n",
    "                            }# res_top has some candidates\n",
    "                        }#fourth condition ends\n",
    "                    }# both sps and mbank have candidates\n",
    "                    mbank_x <- data.frame(cbind(MBmax_similarity, MBmzScore,\n",
    "                                            MBintScore, MQMatchingPeaks,\n",
    "                                            MBTotalPeaks, mQueryTotalPeaks,\n",
    "                                            MBformula, MBSMILES, MBspectrumID,\n",
    "                                            MBcompound_name, Source))\n",
    "                    write.csv(mbank_x, file = paste(dir_name, \"/mbank_results_for_\", id_Xx, \".csv\", sep = \"\"))\n",
    "\n",
    "                }\n",
    "            )\n",
    "            v <- c(future::value(f_gnps), future::value(f_hmdb), future::value(f_mbank))\n",
    "\n",
    "        })\n",
    "    }\n",
    "    pre_mzs <- as.list(pre_mzs)\n",
    "    v_pre_mzs <- future::value(pre_mzs)\n",
    "    result_dir_spectra <- paste(input_dir, str_remove(paste(result_dir, \"/spectral_dereplication\", sep = \"\"), \".\"), sep = \"\")\n",
    "    spectra_input <- data.frame(cbind(id_X, premz, rtmin,\n",
    "                                      rtmax, rtmed, rtmean,\n",
    "                                      col_eng, pol, int, source_file))\n",
    "    write.csv(spectra_input, file = paste(result_dir_spectra, \"/spectral_results_for_\", file_id, \".csv\", sep = \"\"))\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "ab55dd7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#' Extract MS2 Fragment peaks\n",
    "# This functon returns a dataframe and stores a csv file\n",
    "    # the directory for csv file is input_dir + /insilico/MS2DATA.csv\n",
    "# input is from spec_Processing and result directory for each mzML input file\n",
    "\n",
    "ms2_peaks <- function(pre_tbl, proc_mzml, input_dir, result_dir, file_id){\n",
    "\n",
    "\n",
    "    sps_all <- Spectra(proc_mzml, backend = MsBackendMzR())\n",
    "\n",
    "    tbl <- read.table(pre_tbl)\n",
    "    pre_mz <- tbl[[1]]\n",
    "\n",
    "\n",
    "    ## Define variables\n",
    "    premz <- c() # stores mz\n",
    "    rtmin <- c() # stores rtmin\n",
    "    rtmax <- c() # stores rtmax\n",
    "    rtmed <- c() # stores calculated median of rtmin and rtmax\n",
    "    rtmean <- c() # stores calculated mean of rtmin and rtmax\n",
    "    col_eng <- c() # stores collision energy\n",
    "    pol <- c() # stores polarity\n",
    "    ms2Peaks <- c() # stores the peak list file directory\n",
    "    id_X <- c() # creates a unique ID based on mz, rt and also the index\n",
    "        #(since the mz and rt can be similar in some cases)\n",
    "    no_of_ms2_peaks <- c()\n",
    "    int <- c() # stores intensity of the MS1 feature\n",
    "    nx <- 0 # stores number for the ID\n",
    "    indeX <- 0 # stores number to name the peaklist files\n",
    "\n",
    "    # pre_mz is a list of precursor m/z\n",
    "    for (i in pre_mz){\n",
    "\n",
    "\n",
    "        #filter based on pre mz; sps_all is preprocessed spectra\n",
    "        sps <- filterPrecursorMzRange(sps_all, i)\n",
    "        #sps <- filterIntensity(sps, intensity = low_int)\n",
    "\n",
    "        if (length(sps)>0){\n",
    "\n",
    "            #ids\n",
    "            nx <- nx+1\n",
    "            id_Xx <- paste(file_id,  \"M\",  as.character(round(i, digits = 0)),\n",
    "                              \"R\", as.character(round(median(sps$rtime, na.rm = TRUE), digits = 0)),\n",
    "                              \"ID\", as.character(nx), sep = '')\n",
    "            id_X <- c(id_X, id_Xx)\n",
    "\n",
    "\n",
    "            #mz\n",
    "            premz <- c(premz, i)\n",
    "\n",
    "            #rtmin\n",
    "            rn <- min(sps$rtime)\n",
    "            rtmin <- c(rtmin, rn)\n",
    "\n",
    "            #rtmax\n",
    "            rx <- max(sps$rtime)\n",
    "            rtmax <- c(rtmax, rx)\n",
    "\n",
    "            #rtmedian\n",
    "            rtm <- median(sps$rtime, na.rm = TRUE)\n",
    "            rtmed <- c(rtmed, rtm)\n",
    "\n",
    "            #rtmean\n",
    "            rtme <- mean(sps$rtime, na.rm = TRUE)\n",
    "            rtmean <- c(rtmean, rtme)\n",
    "\n",
    "\n",
    "            #collision energy\n",
    "            ce <- max(sps$collisionEnergy)\n",
    "            col_eng <- c(col_eng, ce)\n",
    "\n",
    "            #polarity\n",
    "            pl <- max(sps$polarity)\n",
    "            if (pl == 1){\n",
    "                px <- 'pos'\n",
    "                pol <- c(pol, px)\n",
    "            }\n",
    "            else {\n",
    "                px <- 'neg'\n",
    "                pol <- c(pol, px)\n",
    "            }\n",
    "\n",
    "            #int\n",
    "            ints <- max(sps$precursorIntensity)\n",
    "            int <- c(int, ints)\n",
    "\n",
    "\n",
    "            #peak lists\n",
    "            # variable for name\n",
    "            names <- c()\n",
    "\n",
    "            # create a new directory to store all the peak list txt files\n",
    "            dir_name <- paste(input_dir, str_remove(paste(result_dir, \"/insilico/peakfiles_ms2\", sep =\"\"), \".\"), sep = \"\")\n",
    "            if (!file.exists(dir_name)){\n",
    "                dir.create(dir_name, recursive = TRUE)\n",
    "            }\n",
    "\n",
    "            for (j in 1:length(sps)){\n",
    "                nam <- paste('pk', j, sep = '') ## name of variable\n",
    "                assign(nam, cbind(mz = unlist(mz(sps[j])),intensity = unlist(intensity(sps[j])))) ## assign name to peaklist\n",
    "                names <- c(names, nam) ## save names in another variable\n",
    "\n",
    "                ## at the end of each list, extract the peak list via combinePeaks function\n",
    "                if (j == length(sps)){\n",
    "                    n <- paste(names, collapse = ', ') #paste names at the end\n",
    "                    func <- eval(parse(text = paste('combinePeaks(list(',n,'))', sep = ''))) #write the function and then run it\n",
    "                    indeX <- indeX+1\n",
    "                    Y <- as.character(indeX)# numbering for naming peak lists\n",
    "                    #create separate folder for peaklists files\n",
    "                    fileN <- paste(dir_name, '/Peaks_0', Y, '.txt', sep = '')\n",
    "                    write.table(func, fileN, row.names = FALSE, col.names = FALSE)\n",
    "                    fileN1 <- str_replace(fileN, input_dir, \".\")\n",
    "                    ms2Peaks <- c(ms2Peaks, fileN1)\n",
    "                }\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "\n",
    "\n",
    "\n",
    "    first_list <- data.frame(cbind(id_X, premz, rtmed, rtmean, int ,col_eng, pol, ms2Peaks))\n",
    "    write.csv(first_list, file = paste(input_dir, str_remove(paste(result_dir,'/insilico/MS2DATA.csv', sep = \"\"), \".\"), sep =\"\"))\n",
    "    return(first_list)\n",
    "}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "07e308df",
   "metadata": {},
   "outputs": [],
   "source": [
    "#cam_funcMode <- function(path, pattern = \".mzML\"){\n",
    "    #library(\"CAMERA\")\n",
    "    # List all files present in QC folder\n",
    "    #files_QC_N <- list.files(path, pattern = pattern ,full.names=TRUE)\n",
    "\n",
    "    #for (i in 1:length(files_QC_N)){\n",
    "\n",
    "        # read each file using Spectra\n",
    "        #sps_all <- Spectra(files_QC_N[i], backend = MsBackendMzR())\n",
    "\n",
    "\n",
    "        #if (length(unique(sps_all$polarity)) == 1){\n",
    "            #if(unique(sps_all$polarity) == 1){\n",
    "                #Read the same file with MS1 information; note CAMERA reads xcmsSet object\n",
    "                #xs <- xcmsSet(file = as.character(files_QC_N[i]),\n",
    "                              #profmethod = \"bin\", profparam = list(), lockMassFreq=FALSE,\n",
    "                              #mslevel= 1, progressCallback=NULL, polarity=\"positive\",\n",
    "                              #scanrange = NULL, BPPARAM = bpparam(),\n",
    "                              #stopOnError = TRUE)\n",
    "                # Create an xsAnnotate object\n",
    "                #an <- xsAnnotate(xs)\n",
    "                # Group based on RT\n",
    "                #anF <- groupFWHM(an, perfwhm = 0.6)\n",
    "                # Annotate isotopes\n",
    "                #anI <- findIsotopes(anF, mzabs = 0.01)\n",
    "                # Verify grouping\n",
    "                #anIC <- groupCorr(anI, cor_eic_th = 0.75)\n",
    "                #Annotate adducts\n",
    "                #anFA <- findAdducts(anIC, polarity=\"positive\")\n",
    "                #get a feature list\n",
    "                #peaklist <- getPeaklist(anFA)\n",
    "                # add file_origin information\n",
    "                #peaklist$file_origin <- as.character(files_QC_N[i])\n",
    "                #file_name <- paste(path,\"/posCAMERA_Results_\", i, \".csv\", sep = \"\")\n",
    "                # write individual QC files which are in pos mode (later code will combine them)\n",
    "                #write.csv(peaklist, file = file_name)\n",
    "            #}else if (unique(sps_all$polarity) == 0){\n",
    "                #Read the same file with MS1 information; note CAMERA reads xcmsSet object\n",
    "                #xs <- xcmsSet(file = as.character(files_QC_N[i]),\n",
    "                              #profmethod = \"bin\", profparam = list(), lockMassFreq=FALSE,\n",
    "                              #mslevel= 1, progressCallback=NULL, polarity=\"negative\",\n",
    "                              #scanrange = NULL, BPPARAM = bpparam(),\n",
    "                              #stopOnError = TRUE)\n",
    "                # Create an xsAnnotate object\n",
    "                #an <- xsAnnotate(xs)\n",
    "                # Group based on RT\n",
    "                #anF <- groupFWHM(an, perfwhm = 0.6)\n",
    "                # Annotate isotopes\n",
    "                #anI <- findIsotopes(anF, mzabs = 0.01)\n",
    "                # Verify grouping\n",
    "                #anIC <- groupCorr(anI, cor_eic_th = 0.75)\n",
    "                #Annotate adducts\n",
    "                #anFA <- findAdducts(anIC, polarity=\"negative\")\n",
    "                #get a feature list\n",
    "                #peaklist <- getPeaklist(anFA)\n",
    "                # add file_origin information\n",
    "                #peaklist$file_origin <- as.character(files_QC_N[i])\n",
    "                #file_name <- paste(path, \"/negCAMERA_Results_\", i, \".csv\", sep = \"\")\n",
    "                # write individual QC files which are in pos mode (later code will combine them)\n",
    "                #write.csv(peaklist, file = file_name)\n",
    "            #}\n",
    "        #}\n",
    "    #else{\n",
    "        #pos <- sps_all[sps_all$polarity == 1]\n",
    "        #neg <- sps_all[sps_all$polarity == 0]\n",
    "        #file_p <- paste(path, \"/QC_280k_pos\", i, \".mzML\", sep = \"\")\n",
    "        #file_n <- paste(path, \"/QC_280k_neg\", i, \".mzML\", sep = \"\")\n",
    "\n",
    "        # create new mzML QC files for pos and neg modes from each common QC file\n",
    "        #export(pos, backend = MsBackendMzR(), file = file_p)\n",
    "        #export(neg, backend = MsBackendMzR(), file = file_n)\n",
    "        #Read the same file with MS1 information; note CAMERA reads xcmsSet object\n",
    "        #xs <- xcmsSet(file = as.character(file_p),\n",
    "                        #profmethod = \"bin\", profparam = list(), lockMassFreq=FALSE,\n",
    "                        #mslevel= 1, progressCallback=NULL, polarity=\"positive\",\n",
    "                        #scanrange = NULL, BPPARAM = bpparam(),\n",
    "                        #stopOnError = TRUE)\n",
    "        # Create an xsAnnotate object\n",
    "        #an <- xsAnnotate(xs)\n",
    "        # Group based on RT\n",
    "        #anF <- groupFWHM(an, perfwhm = 0.6)\n",
    "        # Annotate isotopes\n",
    "        #anI <- findIsotopes(anF, mzabs = 0.01)\n",
    "        # Verify grouping\n",
    "        #anIC <- groupCorr(anI, cor_eic_th = 0.75)\n",
    "        #Annotate adducts\n",
    "        #anFA <- findAdducts(anIC, polarity=\"positive\")\n",
    "        #get a feature list\n",
    "        #peaklist <- getPeaklist(anFA)\n",
    "        # add file_origin information\n",
    "        #peaklist$file_origin <- as.character(files_QC_N[i])\n",
    "        #file_name <- paste(path, \"/posCAMERA_Results_\", i, \".csv\", sep = \"\")\n",
    "        # write individual QC files which are in pos mode (later code will combine them)\n",
    "        #write.csv(peaklist, file = file_name)\n",
    "\n",
    "        #Read the same file with MS1 information; note CAMERA reads xcmsSet object\n",
    "        #xs <- xcmsSet(file = as.character(file_n),\n",
    "                        #profmethod = \"bin\", profparam = list(), lockMassFreq=FALSE,\n",
    "                        #mslevel= 1, progressCallback=NULL, polarity=\"negative\",\n",
    "                        #scanrange = NULL, BPPARAM = bpparam(),\n",
    "                        #stopOnError = TRUE)\n",
    "        # Create an xsAnnotate object\n",
    "        #an <- xsAnnotate(xs)\n",
    "        # Group based on RT\n",
    "        #anF <- groupFWHM(an, perfwhm = 0.6)\n",
    "        # Annotate isotopes\n",
    "        #anI <- findIsotopes(anF, mzabs = 0.01)\n",
    "        # Verify grouping\n",
    "        #anIC <- groupCorr(anI, cor_eic_th = 0.75)\n",
    "        #Annotate adducts\n",
    "        #anFA <- findAdducts(anIC, polarity=\"negative\")\n",
    "        #get a feature list\n",
    "        #peaklist <- getPeaklist(anFA)\n",
    "        # add file_origin information\n",
    "        #peaklist$file_origin <- as.character(files_QC_N[i])\n",
    "        #file_name <- paste(path, \"/negCAMERA_Results_\", i, \".csv\", sep = \"\")\n",
    "        # write individual QC files which are in pos mode (later code will combine them)\n",
    "        #write.csv(peaklist, file = file_name)\n",
    "    #}\n",
    "\n",
    "    #}\n",
    "\n",
    "    #detach(\"package:CAMERA\", unload=TRUE)\n",
    "#}\n",
    "\n",
    "# ---------- merge_qc ----------\n",
    "\n",
    "#merge_qc<- function(path){\n",
    "    # combine all QC which are in positive mode\n",
    "    #df_pos <- list.files(path, pattern = \"posCAMERA_Results_\", full.names = TRUE) %>%\n",
    "        #lapply(read_csv) %>%\n",
    "        #bind_rows\n",
    "    # remove any duplicated rows\n",
    "    #df_pos <- as.data.frame(df_pos[!duplicated(df_pos), ])\n",
    "\n",
    "    #extract isotope column numbers, the numbers represent the group of isotope\n",
    "    #nm_p <- regmatches(df_pos[, \"isotopes\"],gregexpr(\"[[:digit:]]+\\\\.*[[:digit:]]*\",df_pos[, \"isotopes\"]))\n",
    "\n",
    "    # for all the numbers, extract only first number, since it is the group number,\n",
    "    # second number can be charge\n",
    "    #for (i in 1:length(nm_p)){\n",
    "       # y <- as.numeric(unlist(nm_p[i]))\n",
    "        #df_pos[i,'istops'] = y[1]\n",
    "    #}\n",
    "\n",
    "    # write csv for the combined_camera_pos results\n",
    "    #write.csv(df_pos, paste(path, \"/Combined_Camera_pos.csv\", sep = \"\"))\n",
    "\n",
    "    # combine all QC which are in negative mode\n",
    "    #df_neg <- list.files(path, pattern = \"negCAMERA_Results_\", full.names = TRUE) %>%\n",
    "        #lapply(read_csv) %>%\n",
    "        #bind_rows\n",
    "    # remove any duplicated rows based on mz\n",
    "    #df_neg <- as.data.frame(df_neg[!duplicated(df_neg), ])\n",
    "\n",
    "    #extract isotope column numbers, the numbers represent the group of isotope\n",
    "    #nm_n <- regmatches(df_neg[, \"isotopes\"],gregexpr(\"[[:digit:]]+\\\\.*[[:digit:]]*\",df_neg[, \"isotopes\"]))\n",
    "\n",
    "    # for all the numbers, extract only first number, since it is the group number,\n",
    "    # second number can be charge\n",
    "    #for (i in 1:length(nm_n)){\n",
    "        #y <- as.numeric(unlist(nm_n[i]))\n",
    "        #df_neg[i,'istops'] = y[1]\n",
    "    #}\n",
    "    # write csv for the combined_camera_neg results\n",
    "    #write.csv(df_neg, paste(path, \"/Combined_Camera_neg.csv\", sep = \"\"))\n",
    "#}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e9f9b0ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "cam_func <- function(input_dir, f, ms2features){\n",
    "    modes_file <- read_csv(ms2features)\n",
    "    mode = unique(modes_file[\"pol\"])\n",
    "    fl <- paste(input_dir, str_remove(f, \".\"), sep =\"\")\n",
    "    if(mode == \"pos\"){\n",
    "        library(\"CAMERA\")\n",
    "        xs <- xcmsSet(file = fl,profmethod = \"bin\",\n",
    "              profparam = list(), lockMassFreq=FALSE,\n",
    "              mslevel= 1, progressCallback=NULL, polarity=\"positive\",\n",
    "              scanrange = NULL, BPPARAM = bpparam(),stopOnError = TRUE)\n",
    "        # Create an xsAnnotate object\n",
    "        an <- xsAnnotate(xs)\n",
    "        # Group based on RT\n",
    "        anF <- groupFWHM(an, perfwhm = 0.6)\n",
    "        # Annotate isotopes\n",
    "        anI <- findIsotopes(anF, mzabs = 0.01)\n",
    "        # Verify grouping\n",
    "        anIC <- groupCorr(anI, cor_eic_th = 0.75)\n",
    "        #Annotate adducts\n",
    "        anFA <- findAdducts(anIC, polarity=\"positive\")\n",
    "        peaklist <- getPeaklist(anFA)\n",
    "        peaklist$file_origin <- fl\n",
    "        #extract isotope column numbers, the numbers represent the group of isotope\n",
    "        nm_po <- regmatches(peaklist[, \"isotopes\"],gregexpr(\"[[:digit:]]+\\\\.*[[:digit:]]*\",peaklist[, \"isotopes\"]))\n",
    "        # for all the numbers in v, extract only first number, since it is the group number,\n",
    "        # second number can be charge\n",
    "        for (i in 1:length(nm_po)){\n",
    "            y <- as.numeric(unlist(nm_po[i]))\n",
    "            peaklist[i,'istops'] = y[1]\n",
    "        }\n",
    "        name <- str_remove(f, \".mzML\")\n",
    "        write.csv(peaklist, file = paste(input_dir, str_remove(name, \".\"), \"/CAMERAResults\",\".csv\", sep = \"\"))\n",
    "        unloadNamespace(\"CAMERA\")\n",
    "        unloadNamespace(\"xcms\")\n",
    "        unloadNamespace(\"MSnBase\")\n",
    "\n",
    "\n",
    "    }else{\n",
    "        library(\"CAMERA\")\n",
    "        xs <- xcmsSet(file = fl,profmethod = \"bin\",\n",
    "              profparam = list(), lockMassFreq=FALSE,\n",
    "              mslevel= 1, progressCallback=NULL, polarity=\"negative\",\n",
    "              scanrange = NULL, BPPARAM = bpparam(),stopOnError = TRUE)\n",
    "        # Create an xsAnnotate object\n",
    "        an <- xsAnnotate(xs)\n",
    "        # Group based on RT\n",
    "        anF <- groupFWHM(an, perfwhm = 0.6)\n",
    "        # Annotate isotopes\n",
    "        anI <- findIsotopes(anF, mzabs = 0.01)\n",
    "        # Verify grouping\n",
    "        anIC <- groupCorr(anI, cor_eic_th = 0.75)\n",
    "        #Annotate adducts\n",
    "        anFA <- findAdducts(anIC, polarity=\"negative\")\n",
    "        peaklist <- getPeaklist(anFA)\n",
    "        peaklist$file_origin <- fl\n",
    "        #extract isotope column numbers, the numbers represent the group of isotope\n",
    "        nm_ne <- regmatches(peaklist[, \"isotopes\"],gregexpr(\"[[:digit:]]+\\\\.*[[:digit:]]*\",peaklist[, \"isotopes\"]))\n",
    "        # for all the numbers in v, extract only first number, since it is the group number,\n",
    "        # second number can be charge\n",
    "        for (i in 1:length(nm_ne)){\n",
    "            y <- as.numeric(unlist(nm_ne[i]))\n",
    "            peaklist[i,'istops'] = y[1]\n",
    "        }\n",
    "        name <- str_remove(f, \".mzML\")\n",
    "        write.csv(peaklist, file = paste(input_dir, str_remove(name, \".\"), \"/CAMERAResults\", \".csv\", sep = \"\"))\n",
    "        unloadNamespace(\"CAMERA\")\n",
    "        unloadNamespace(\"xcms\")\n",
    "        unloadNamespace(\"MSnBase\")\n",
    "    }\n",
    "    return(peaklist)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "6228bb1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract isotopic peaks for each pre_mz\n",
    "# The input is x = first_list (from ms2peaks function) and y = camera results\n",
    "\n",
    "ms1_peaks <- function(x, y, result_dir, input_dir, QCfile){\n",
    "    # store the ms1_peak list path here\n",
    "    ms1Peaks <- c()\n",
    "    x = read.csv(x)\n",
    "\n",
    "    if (QCfile){\n",
    "\n",
    "        dir_name <- paste(input_dir, str_remove(paste(result_dir, \"/insilico/peakfiles_ms1\", sep =\"\"), \".\"), sep = \"\")\n",
    "        # create a new directory to store all the peak list txt files\n",
    "        if (!file.exists(dir_name)){\n",
    "            dir.create(dir_name, recursive = TRUE)\n",
    "        }\n",
    "\n",
    "        # read the CAMERA results\n",
    "        y = read.csv(y)\n",
    "\n",
    "\n",
    "        # for all indices in the ms2 features table\n",
    "        for (i in 1:nrow(x)){\n",
    "            #store the indices of CAMERA that have same mz and rt as the ms2 features table\n",
    "            store_c <- c()\n",
    "            # for all indices in CAMERA Results\n",
    "            for (j in 1:nrow(y)){\n",
    "                # if mz and rt from ms2 features are within the range of CAMERA mz and rt\n",
    "                if (x[i, 'premz'] <= y[j, \"mzmax\"]   && y[j, \"mzmin\"] <= x[i, 'premz'] && x[i, 'rtmed'] <= y[j, \"rtmax\"] && y[j, \"rtmin\"] <= x[i, 'rtmed']){\n",
    "                    store_c <- c(store_c, j)\n",
    "                }\n",
    "            }\n",
    "            # indices with same pre m/z and same rt\n",
    "            df_y <- y[store_c, ]\n",
    "            df_y <- as.data.frame(df_y)\n",
    "\n",
    "            #if there was only one index\n",
    "            if (nrow(df_y) == 1){\n",
    "\n",
    "                # if there was no isotope annotation for that one index\n",
    "                if (is.na(df_y[1, \"istops\"])){\n",
    "\n",
    "                    mz <- df_y[1, \"mz\"] # save mz\n",
    "                    int <- df_y[1, \"into\"] # save intensity\n",
    "                    no_isotop <- cbind(mz, int) # save as table\n",
    "                    name_file <- paste(dir_name, \"/ms1_peaks_\", x[i, 'premz'], \"_no_isotopes.txt\", sep = \"\") # save name of the peaklist\n",
    "                    write.table(no_isotop, name_file, row.names = FALSE, col.names = FALSE) # save peak list\n",
    "                    name_file1 <- str_replace(name_file, input_dir, \".\")\n",
    "                    ms1Peaks <- c(ms1Peaks, name_file1) # add the path of the peak list to a list\n",
    "                }\n",
    "\n",
    "                # if there was an isotope annotation\n",
    "                else{\n",
    "\n",
    "                    df_x <- y[which(y[, \"file_origin\"] ==df_y[1, \"file_origin\"]), ] # extract camera results from one file origin\n",
    "                    df_x <- df_x[which(df_x[, 'istops'] == df_y[1, 'istops']), ] # extract only certain isotope annotation group\n",
    "                    mz <- df_x[, \"mz\"] # save mz\n",
    "                    int <- df_x[, \"into\"] # save intensity\n",
    "                    no_isotop <- cbind(mz, int) # save as table\n",
    "                    name_file <- paste(dir_name, \"/ms1_peaksISOTOPE_\", x[i, 'premz'], \"_isotopeNum_\", df_x[1, \"istops\"], \".txt\", sep = \"\")\n",
    "                    write.table(no_isotop, name_file, row.names = FALSE, col.names = FALSE)\n",
    "                    name_file1 <- str_replace(name_file, input_dir, \".\")\n",
    "                    ms1Peaks <- c(ms1Peaks, name_file1)\n",
    "                }\n",
    "            }\n",
    "            # if there are more indices for df_y\n",
    "            else if(nrow(df_y) > 1){\n",
    "                # if all enteries have no isotope annotation\n",
    "                if(all(is.na(df_y[, 'istops']))){\n",
    "\n",
    "                    df_z <- df_y[which(df_y[,\"into\"] == max(df_y[,\"into\"])), ] # extract the ms1 peak with highest intensity\n",
    "                    mz <- df_z[1, \"mz\"] # save mz\n",
    "                    int <- df_z[1, \"into\"] # save intensity\n",
    "                    no_isotop <- cbind(mz, int) # save as table\n",
    "                    name_file <- paste(dir_name, \"/ms1_peaks_\", x[i, 'premz'], \"_no_isotopes.txt\", sep = \"\") # save name of the peaklist\n",
    "                    write.table(no_isotop, name_file, row.names = FALSE, col.names = FALSE) # save peak list\n",
    "                    name_file1 <- str_replace(name_file, input_dir, \".\")\n",
    "                    ms1Peaks <- c(ms1Peaks, name_file1) # add the path of the peak list to a list\n",
    "                }\n",
    "                # if not all isotope annotations are NA\n",
    "                else if (!(all(is.na(df_y[, 'istops'])))){\n",
    "\n",
    "                    df_y <- df_y[!is.na(df_y$'istops'),] # Remove the NA isotope annotations\n",
    "                    df_z <- df_y[which(df_y[,\"into\"] == max(df_y[,\"into\"])), ] # Select the MS1 peak with highest intensity\n",
    "                    df_z1 <- y[which(y[, \"file_origin\"] == df_z[1, \"file_origin\"]), ]  # extract camera results from one file origin\n",
    "                    df_z1 <- df_z1[which(df_z1[, 'istops'] == df_z[1, 'istops']), ] # extract only certain isotope annotation group\n",
    "                    mz <- df_z1[, \"mz\"] # save mz\n",
    "                    int <- df_z1[, \"into\"] # save intensity\n",
    "                    no_isotop <- cbind(mz, int) # save as table\n",
    "                    name_file <- paste(dir_name, \"/ms1_peaksISOTOPE_\", x[i, 'premz'], \"_isotopeNum_\", df_z1[1, 'istops'],\".txt\", sep = \"\") # save name of the peaklist\n",
    "                    write.table(no_isotop, name_file, row.names = FALSE, col.names = FALSE) # save peak list\n",
    "                    name_file1 <- str_replace(name_file, input_dir, \".\")\n",
    "                    ms1Peaks <- c(ms1Peaks, name_file1) # add the path of the peak list to a list\n",
    "                }\n",
    "            }\n",
    "            else if (nrow(df_y)==0){\n",
    "                ms1Peaks <- c(ms1Peaks, 'no ms1 peaks in QC')\n",
    "            }\n",
    "        }\n",
    "        second_list <- data.frame(cbind(x, ms1Peaks))\n",
    "        write.csv(second_list, file = paste(input_dir, str_remove(paste(result_dir,'/insilico/MS1DATA.csv', sep = \"\"), \".\"), sep =\"\"))\n",
    "        return(second_list)\n",
    "    }\n",
    "    else{\n",
    "\n",
    "        ms1Peaks <- c(ms1Peaks, 'no ms1 peaks in QC')\n",
    "        second_list <- data.frame(cbind(x, ms1Peaks))\n",
    "        write.csv(second_list, file = paste(input_dir, str_remove(paste(result_dir,'/insilico/MS1DATA.csv', sep = \"\"), \".\"), sep =\"\"))\n",
    "        return(second_list)\n",
    "    }\n",
    "\n",
    "}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "5ede9d69",
   "metadata": {},
   "outputs": [],
   "source": [
    "# sirius_postprocess <- function(x, SL = TRUE){\n",
    "#     feat_scale <- function(p) {\n",
    "#     (p - min(p)) / (max(p) - min(p))\n",
    "#     }\n",
    "#     #the result directory name for each file\n",
    "#     dir_name <- paste(x, \"/insilico/SIRIUS\", sep = '')\n",
    "#     ##result json files in each result directory\n",
    "#     #for suspect list\n",
    "#     SL_Param <- list.files(dir_name, pattern = 'SList.json', full.names = TRUE)\n",
    "#     #for all DB\n",
    "#     Param <- list.files(dir_name, pattern = 'param.json', full.names = TRUE)\n",
    "#     #DATA FRAME of both directories for each feature\n",
    "#     parameter_json <- cbind(SL_Param, Param)\n",
    "#     #read the msdata csv file that contains all features and their data\n",
    "#     msdata <- as.data.frame(read.csv(paste(x, \"/insilico/MS1DATA.csv\", sep = '')))\n",
    "#     #add new empty columns to store information from SIRIUS results\n",
    "#     msdata$Adducts <- NA           #adducts\n",
    "#     msdata$name <- NA              #name of compound\n",
    "#     msdata$PubChemIDs <- NA        #all pubchem ids\n",
    "#     msdata$SMILES <- NA            #smiles of compound\n",
    "#     msdata$Formula <- NA           #formula of compound\n",
    "#     msdata$FormulaRank <- NA       #formula rank\n",
    "#     msdata$SIRIUSscore <- NA       #Formula score\n",
    "#     msdata$CSIFingerIDscore <- NA  #Structure score\n",
    "#     msdata$SMILESforMCSS <- NA   #SMILES of top scoring candidates; to calculate their tanimoto later\n",
    "#     msdata$exp_int <- NA           #explained intensity of formula\n",
    "#     msdata$dir <- NA               #directory for results either the strcuture or formula candidates\n",
    "#     msdata$Result <- NA             #the type of candidate, either from Suspect list, strcuture command or formula command\n",
    "#     # for all entries in the json result files for each each feature\n",
    "#     for (i in 1:nrow(parameter_json)){\n",
    "\n",
    "#         # find the number of the row corresponding to the row in msdata that has the same precursor m/z as json result files\n",
    "#         rowMS <- msdata[grepl(str_match(as.character(parameter_json[i, 'Param']), \"MS1p_\\\\s*(.*?)\\\\s*_SIRIUS\")[2] ,msdata$premz), ]\n",
    "\n",
    "#         if (SL){\n",
    "\n",
    "#             # file path for strcuture candidate from Suspect List json folder\n",
    "#             str_canS <- paste(list.dirs(parameter_json[i,'SL_Param'])[2], '/structure_candidates.tsv', sep = '')\n",
    "#             # file path for formula candidate from Param json folder\n",
    "#             for_canS <- paste(list.dirs(parameter_json[i,'SL_Param'])[2], '/formula_candidates.tsv', sep = '')\n",
    "#             # file path for strcuture candidate from Suspect List json folder\n",
    "#             str_can <- paste(list.dirs(parameter_json[i,'Param'])[2], '/structure_candidates.tsv', sep = '')\n",
    "#             # file path for formula candidate from Param json folder\n",
    "#             for_can <- paste(list.dirs(parameter_json[i,'Param'])[2], '/formula_candidates.tsv', sep = '')\n",
    "\n",
    "#             # if the strcuture candidate file exists\n",
    "#             if (file.exists(str_canS) && file.exists(str_can)){\n",
    "\n",
    "#                 # read the corresponding structure and formula candidate files\n",
    "#                 str_canSL <- as.data.frame(read_tsv(str_canS))\n",
    "#                 for_canSL <- as.data.frame(read_tsv(for_canS))\n",
    "\n",
    "#                 # read the corresponding structure and formula candidate files\n",
    "#                 str_canP <- as.data.frame(read_tsv(str_can))\n",
    "#                 for_canP <- as.data.frame(read_tsv(for_can))\n",
    "\n",
    "#                 # if the strcuture candidate file contains 1 or more rows, it has detected a candidate from suspect list, add relevant info\n",
    "\n",
    "#                 if (nrow(str_canSL) >= 1){\n",
    "\n",
    "#                     if (str_canSL[1, 'CSI:FingerIDScore'] > str_canP[1, 'CSI:FingerIDScore']){\n",
    "#                         # information from structure candidate file\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- str_canSL[1, 'adduct']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'name'] <- str_canSL[1, 'name']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'PubChemIDs'] <- str_canSL[1, 'pubchemids']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'SMILES'] <- str_canSL[1, 'smiles']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'Formula'] <- str_canSL[1, 'molecularFormula']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'FormulaRank'] <- str_canSL[1, 'formulaRank']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'CSIFingerIDscore'] <- str_canSL[1, 'CSI:FingerIDScore']\n",
    "\n",
    "#                         # information from formula candidate file\n",
    "#                         formulaRow <- which(for_canSL[,'rank'] == str_canSL[1, 'formulaRank'])\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'SIRIUSscore'] <- for_canSL[formulaRow, 'SiriusScore']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'exp_int'] <- for_canSL[formulaRow, 'explainedIntensity']\n",
    "\n",
    "#                         # other info\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'Result'] <- 'SIRIUS_SL'\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'dir'] <- str_canS\n",
    "#                     }\n",
    "#                     else{\n",
    "\n",
    "#                         # if the structure candidate file contains 1 row, it has detected a candidate from all DBs, add relevant info\n",
    "#                         if (nrow(str_canP) == 1){\n",
    "\n",
    "#                             # information from structure candidate file\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- str_canP[1, 'adduct']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'name'] <- str_canP[1, 'name']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'PubChemIDs'] <- str_canP[1, 'pubchemids']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'SMILES'] <- str_canP[1, 'smiles']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'Formula'] <- str_canP[1, 'molecularFormula']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'FormulaRank'] <- str_canP[1, 'formulaRank']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'CSIFingerIDscore'] <- str_canP[1, 'CSI:FingerIDScore']\n",
    "\n",
    "#                             # information from formula candidate file\n",
    "#                             formulaRow1 <- which(for_canP[,'rank'] == str_canP[1, 'formulaRank'])\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'SIRIUSscore'] <- for_canP[formulaRow1, 'SiriusScore']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'exp_int'] <- for_canP[formulaRow1, 'explainedIntensity']\n",
    "\n",
    "#                             # other info\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'Result'] <- 'SIRIUS_STR'\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'dir'] <- str_can\n",
    "#                         }\n",
    "#                         # if the structure candidate file contains more rows, extract SMILES of top candidates and check their similarity later\n",
    "#                         else if (nrow(str_canP) > 1){\n",
    "\n",
    "#                             # information from structure candidate file\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- str_canP[1, 'adduct']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'name'] <- str_canP[1, 'name']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'PubChemIDs'] <- str_canP[1, 'pubchemids']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'SMILES'] <- str_canP[1, 'smiles']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'Formula'] <- str_canP[1, 'molecularFormula']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'FormulaRank'] <- str_canP[1, 'formulaRank']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'CSIFingerIDscore'] <- str_canP[1, 'CSI:FingerIDScore']\n",
    "\n",
    "#                             # information from formula candidate file, take info from the formula rank that corresponds to the formula rank with the top strcuture candidate\n",
    "#                             formulaRow2 <- which(for_canP[,'rank'] == str_canP[1, 'formulaRank'])\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'SIRIUSscore'] <- for_canP[formulaRow2, 'SiriusScore']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'exp_int'] <- for_canP[formulaRow2, 'explainedIntensity']\n",
    "\n",
    "#                             # other info\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'Result'] <- 'SIRIUS_STR'\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'dir'] <- str_can\n",
    "\n",
    "#                             # normalize the CSI:FingerIDScores\n",
    "#                             norm_score <- feat_scale(str_canP[,\"CSI:FingerIDScore\"])\n",
    "#                             # store the upper quartile\n",
    "#                             upper_quartile <- str_canP[which(norm_score > as.numeric(quantile(norm_score)[4])), \"smiles\"]\n",
    "\n",
    "#                             # if the upper quartile has more than 5 candidates, then just take the top 5 candidates\n",
    "#                             if (length(upper_quartile) > 5){\n",
    "#                                 upper_quartile <- upper_quartile[1:5]\n",
    "#                             }\n",
    "#                             # save the top candidates SMILES, to check similarity later with rdkit in Python\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'SMILESforMCSS'] <- paste(upper_quartile, collapse = '|')\n",
    "\n",
    "#                         }\n",
    "#                         # if the structure candidate file is empty, take information from just the formula candidate file\n",
    "#                         else if (nrow(str_canP) == 0){\n",
    "\n",
    "#                             # if formula candidate file is not empty\n",
    "#                             if (nrow(for_canP) >= 1){\n",
    "\n",
    "#                                 # information from formula candidate file\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- for_canP[1, 'adduct']\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'Formula'] <- for_canP[1, 'molecularFormula']\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'FormulaRank'] <- for_canP[1, 'rank']\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'exp_int'] <- for_canP[1, 'explainedIntensity']\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'SIRIUSscore'] <- for_canP[1, 'SiriusScore']\n",
    "#                                 # other info\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'Result'] <- 'SIRIUS_FOR'\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'dir'] <- for_can\n",
    "\n",
    "#                             }\n",
    "#                         }\n",
    "\n",
    "#                         # if the structure candidate from all DBs does not exist\n",
    "#                         else{\n",
    "#                             # check if the formula candidate file exists\n",
    "#                             if (file.exists(for_can)){\n",
    "#                                 for_canF1 <- as.data.frame(read_tsv(for_can))\n",
    "\n",
    "#                                 # if formula candidate file is not empty\n",
    "#                                 if (nrow(for_canF1)>= 1){\n",
    "\n",
    "#                                     # information from formula candidate file\n",
    "#                                     msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- for_canF1[1, 'adduct']\n",
    "#                                     msdata[as.numeric(rownames(rowMS)), 'Formula'] <- for_canF1[1, 'molecularFormula']\n",
    "#                                     msdata[as.numeric(rownames(rowMS)), 'FormulaRank'] <- for_canF1[1, 'rank']\n",
    "#                                     msdata[as.numeric(rownames(rowMS)), 'exp_int'] <- for_canF1[1, 'explainedIntensity']\n",
    "#                                     msdata[as.numeric(rownames(rowMS)), 'SIRIUSscore'] <- for_canF1[1, 'SiriusScore']\n",
    "#                                     # other info\n",
    "#                                     msdata[as.numeric(rownames(rowMS)), 'Result'] <- 'SIRIUS_FOR'\n",
    "#                                     msdata[as.numeric(rownames(rowMS)), 'dir'] <- for_can\n",
    "#                                 }\n",
    "#                             }\n",
    "#                         }\n",
    "#                     }\n",
    "#                 }\n",
    "#                 # if it's empty, move onto the All DB result folder called PARAM here\n",
    "#                 else{\n",
    "\n",
    "#                     # if the structure candidate file contains 1 row, it has detected a candidate from all DBs, add relevant info\n",
    "#                     if (nrow(str_canP) == 1){\n",
    "\n",
    "#                         # information from structure candidate file\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- str_canP[1, 'adduct']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'name'] <- str_canP[1, 'name']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'PubChemIDs'] <- str_canP[1, 'pubchemids']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'SMILES'] <- str_canP[1, 'smiles']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'Formula'] <- str_canP[1, 'molecularFormula']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'FormulaRank'] <- str_canP[1, 'formulaRank']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'CSIFingerIDscore'] <- str_canP[1, 'CSI:FingerIDScore']\n",
    "\n",
    "#                         # information from formula candidate file\n",
    "#                         formulaRow1 <- which(for_canP[,'rank'] == str_canP[1, 'formulaRank'])\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'SIRIUSscore'] <- for_canP[formulaRow1, 'SiriusScore']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'exp_int'] <- for_canP[formulaRow1, 'explainedIntensity']\n",
    "\n",
    "#                         # other info\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'Result'] <- 'SIRIUS_STR'\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'dir'] <- str_can\n",
    "#                     }\n",
    "#                     # if the structure candidate file contains more rows, extract SMILES of top candidates and check their similarity later\n",
    "#                     else if (nrow(str_canP) > 1){\n",
    "\n",
    "#                         # information from structure candidate file\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- str_canP[1, 'adduct']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'name'] <- str_canP[1, 'name']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'PubChemIDs'] <- str_canP[1, 'pubchemids']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'SMILES'] <- str_canP[1, 'smiles']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'Formula'] <- str_canP[1, 'molecularFormula']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'FormulaRank'] <- str_canP[1, 'formulaRank']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'CSIFingerIDscore'] <- str_canP[1, 'CSI:FingerIDScore']\n",
    "\n",
    "#                         # information from formula candidate file, take info from the formula rank that corresponds to the formula rank with the top strcuture candidate\n",
    "#                         formulaRow2 <- which(for_canP[,'rank'] == str_canP[1, 'formulaRank'])\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'SIRIUSscore'] <- for_canP[formulaRow2, 'SiriusScore']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'exp_int'] <- for_canP[formulaRow2, 'explainedIntensity']\n",
    "\n",
    "#                         # other info\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'Result'] <- 'SIRIUS_STR'\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'dir'] <- str_can\n",
    "\n",
    "#                         # normalize the CSI:FingerIDScores\n",
    "#                         norm_score <- feat_scale(str_canP[,\"CSI:FingerIDScore\"])\n",
    "#                         # store the upper quartile\n",
    "#                         upper_quartile <- str_canP[which(norm_score > as.numeric(quantile(norm_score)[4])), \"smiles\"]\n",
    "\n",
    "#                         # if the upper quartile has more than 5 candidates, then just take the top 5 candidates\n",
    "#                         if (length(upper_quartile) > 5){\n",
    "#                             upper_quartile <- upper_quartile[1:5]\n",
    "#                         }\n",
    "#                         # save the top candidates SMILES, to check similarity later with rdkit in Python\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'SMILESforMCSS'] <- paste(upper_quartile, collapse = '|')\n",
    "\n",
    "#                     }\n",
    "#                     # if the structure candidate file is empty, take information from just the formula candidate file\n",
    "#                     else if (nrow(str_canP) == 0){\n",
    "\n",
    "#                         # if formula candidate file is not empty\n",
    "#                         if (nrow(for_canP) >= 1){\n",
    "\n",
    "#                             # information from formula candidate file\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- for_canP[1, 'adduct']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'Formula'] <- for_canP[1, 'molecularFormula']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'FormulaRank'] <- for_canP[1, 'rank']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'exp_int'] <- for_canP[1, 'explainedIntensity']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'SIRIUSscore'] <- for_canP[1, 'SiriusScore']\n",
    "#                             # other info\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'Result'] <- 'SIRIUS_FOR'\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'dir'] <- for_can\n",
    "\n",
    "#                         }\n",
    "#                     }\n",
    "\n",
    "#                     # if the structure candidate from all DBs does not exist\n",
    "#                     else{\n",
    "#                         # check if the formula candidate file exists\n",
    "#                         if (file.exists(for_can)){\n",
    "#                             for_canF1 <- as.data.frame(read_tsv(for_can))\n",
    "\n",
    "#                             # if formula candidate file is not empty\n",
    "#                             if (nrow(for_canF1)>= 1){\n",
    "\n",
    "#                                 # information from formula candidate file\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- for_canF1[1, 'adduct']\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'Formula'] <- for_canF1[1, 'molecularFormula']\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'FormulaRank'] <- for_canF1[1, 'rank']\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'exp_int'] <- for_canF1[1, 'explainedIntensity']\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'SIRIUSscore'] <- for_canF1[1, 'SiriusScore']\n",
    "#                                 # other info\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'Result'] <- 'SIRIUS_FOR'\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'dir'] <- for_can\n",
    "#                             }\n",
    "#                         }\n",
    "#                     }\n",
    "#                 }\n",
    "#             }\n",
    "#             else{\n",
    "#                 # directory for structure candidate file\n",
    "#                 str_can1 <- paste(list.dirs(parameter_json[i,'Param'])[2], '/structure_candidates.tsv', sep = '')\n",
    "#                 # directory for formula candidate file\n",
    "#                 for_can1 <- paste(list.dirs(parameter_json[i,'Param'])[2], '/formula_candidates.tsv', sep = '')\n",
    "\n",
    "#                 # if the structure candidate file exists\n",
    "#                 if (file.exists(str_can1)){\n",
    "\n",
    "#                     # read the structure file from All Dbs (PARAM) json file\n",
    "#                     str_canP1 <- as.data.frame(read_tsv(str_can1))\n",
    "#                     # read the formula file from All Dbs (PARAM) json file\n",
    "#                     for_canP1 <- as.data.frame(read_tsv(for_can1))\n",
    "\n",
    "#                     #if the structure candidate file has one candidate, it has detected a candidate from all DBs, add relevant info\n",
    "#                     if (nrow(str_canP1) == 1){\n",
    "\n",
    "#                         # information from structure candidate file\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- str_canP1[1, 'adduct']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- str_canP1[1, 'adduct']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'name'] <- str_canP1[1, 'name']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'PubChemIDs'] <- str_canP1[1, 'pubchemids']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'SMILES'] <- str_canP1[1, 'smiles']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'Formula'] <- str_canP1[1, 'molecularFormula']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'FormulaRank'] <- str_canP1[1, 'formulaRank']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'CSIFingerIDscore'] <- str_canP1[1, 'CSI:FingerIDScore']\n",
    "#                         # information from formula candidate file\n",
    "#                         formulaRow3 <- which(for_canP1[,'rank'] == str_canP1[1, 'formulaRank'])\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'SIRIUSscore'] <- for_canP1[formulaRow3, 'SiriusScore']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'exp_int'] <- for_canP1[formulaRow3, 'explainedIntensity']\n",
    "#                         # other info\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'Result'] <- 'SIRIUS_STR'\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'dir'] <- str_can1\n",
    "\n",
    "#                     }\n",
    "#                     # if the strcuture cabdidate file has more than 1 candidates, it has detected candidates from all DBs, add relevant info\n",
    "#                     else if (nrow(str_canP1) > 1){\n",
    "\n",
    "#                         # information from structure candidate file\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- str_canP1[1, 'adduct']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'name'] <- str_canP1[1, 'name']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'PubChemIDs'] <- str_canP1[1, 'pubchemids']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'SMILES'] <- str_canP1[1, 'smiles']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'Formula'] <- str_canP1[1, 'molecularFormula']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'FormulaRank'] <- str_canP1[1, 'formulaRank']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'CSIFingerIDscore'] <- str_canP1[1, 'CSI:FingerIDScore']\n",
    "#                         # information from formula candidate file\n",
    "#                         formulaRow4 <- which(for_canP1[,'rank'] == str_canP1[1, 'formulaRank'])\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'SIRIUSscore'] <- for_canP1[formulaRow4, 'SiriusScore']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'exp_int'] <- for_canP1[formulaRow4, 'explainedIntensity']\n",
    "#                         # other info\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'Result'] <- 'SIRIUS_STR'\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'dir'] <- str_can1\n",
    "\n",
    "#                         # normalize the CSI:FingerIDScores\n",
    "#                         norm_score1 <- feat_scale(str_canP1[,\"CSI:FingerIDScore\"])\n",
    "#                         # store the upper quartile\n",
    "#                         upper_quartile1 <- str_canP1[which(norm_score1 > as.numeric(quantile(norm_score1)[4])), \"smiles\"]\n",
    "#                         # if the upper quartile has more than 5 candidates, then just take the top 5 candidates\n",
    "#                         if (length(upper_quartile1) > 5){\n",
    "#                             upper_quartile1 <- upper_quartile1[1:5]\n",
    "#                         }\n",
    "#                         # save the top candidates SMILES, to check similarity later with rdkit in Python\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'SMILESforMCSS'] <- paste(upper_quartile1, collapse = '|')\n",
    "#                     }\n",
    "#                     #\n",
    "#                     else if (nrow(str_canP1) == 0){\n",
    "#                         if (file.exists(for_can1)){\n",
    "#                             for_canP2 <- as.data.frame(read_tsv(for_can1))\n",
    "#                             if (nrow(for_canP2)>= 1){\n",
    "\n",
    "#                                 # information from formula candidate file\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- for_canP2[1, 'adduct']\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'Formula'] <- for_canP2[1, 'molecularFormula']\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'FormulaRank'] <- for_canP2[1, 'rank']\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'exp_int'] <- for_canP2[1, 'explainedIntensity']\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'SIRIUSscore'] <- for_canP2[1, 'SiriusScore']\n",
    "#                                 # other info\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'Result'] <- 'SIRIUS_FOR'\n",
    "#                                 msdata[as.numeric(rownames(rowMS)), 'dir'] <- for_can1\n",
    "#                             }\n",
    "#                         }\n",
    "#                     }\n",
    "#                 }\n",
    "\n",
    "#                 # if the structure candidate file doesn't exists (and no str_sl exists)\n",
    "#                 else{\n",
    "#                     # if formula candidate file exists\n",
    "#                     if (file.exists(for_can1)){\n",
    "#                         for_canP3 <- as.data.frame(read_tsv(for_can1))\n",
    "#                         if (nrow(for_canP3) >= 1){\n",
    "#                             # information from formula candidate file\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- for_canP3[1, 'adduct']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'Formula'] <- for_canP3[1, 'molecularFormula']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'FormulaRank'] <- for_canP3[1, 'rank']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'exp_int'] <- for_canP3[1, 'explainedIntensity']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'SIRIUSscore'] <- for_canP3[1, 'SiriusScore']\n",
    "#                             # other info\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'Result'] <- 'SIRIUS_FOR'\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'dir'] <- for_can1\n",
    "#                         }\n",
    "#                     }\n",
    "#                 }\n",
    "#             }\n",
    "#         }# SL IS false\n",
    "#         else{\n",
    "#             # directory for structure candidate file\n",
    "#             str_can1 <- paste(list.dirs(parameter_json[i,'Param'])[2], '/structure_candidates.tsv', sep = '')\n",
    "#             # directory for formula candidate file\n",
    "#             for_can1 <- paste(list.dirs(parameter_json[i,'Param'])[2], '/formula_candidates.tsv', sep = '')\n",
    "\n",
    "#             # if the structure candidate file exists\n",
    "#             if (file.exists(str_can1)){\n",
    "\n",
    "#                 # read the structure file from All Dbs (PARAM) json file\n",
    "#                 str_canP1 <- as.data.frame(read_tsv(str_can1))\n",
    "#                 # read the formula file from All Dbs (PARAM) json file\n",
    "#                 for_canP1 <- as.data.frame(read_tsv(for_can1))\n",
    "\n",
    "#                 #if the structure candidate file has one candidate, it has detected a candidate from all DBs, add relevant info\n",
    "#                 if (nrow(str_canP1) == 1){\n",
    "\n",
    "#                     # information from structure candidate file\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- str_canP1[1, 'adduct']\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- str_canP1[1, 'adduct']\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'name'] <- str_canP1[1, 'name']\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'PubChemIDs'] <- str_canP1[1, 'pubchemids']\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'SMILES'] <- str_canP1[1, 'smiles']\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'Formula'] <- str_canP1[1, 'molecularFormula']\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'FormulaRank'] <- str_canP1[1, 'formulaRank']\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'CSIFingerIDscore'] <- str_canP1[1, 'CSI:FingerIDScore']\n",
    "#                     # information from formula candidate file\n",
    "#                     formulaRow3 <- which(for_canP1[,'rank'] == str_canP1[1, 'formulaRank'])\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'SIRIUSscore'] <- for_canP1[formulaRow3, 'SiriusScore']\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'exp_int'] <- for_canP1[formulaRow3, 'explainedIntensity']\n",
    "#                     # other info\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'Result'] <- 'SIRIUS_STR'\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'dir'] <- str_can1\n",
    "\n",
    "#                 }\n",
    "#                 # if the strcuture cabdidate file has more than 1 candidates, it has detected candidates from all DBs, add relevant info\n",
    "#                 else if (nrow(str_canP1) > 1){\n",
    "\n",
    "#                     # information from structure candidate file\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- str_canP1[1, 'adduct']\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'name'] <- str_canP1[1, 'name']\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'PubChemIDs'] <- str_canP1[1, 'pubchemids']\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'SMILES'] <- str_canP1[1, 'smiles']\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'Formula'] <- str_canP1[1, 'molecularFormula']\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'FormulaRank'] <- str_canP1[1, 'formulaRank']\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'CSIFingerIDscore'] <- str_canP1[1, 'CSI:FingerIDScore']\n",
    "#                     # information from formula candidate file\n",
    "#                     formulaRow4 <- which(for_canP1[,'rank'] == str_canP1[1, 'formulaRank'])\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'SIRIUSscore'] <- for_canP1[formulaRow4, 'SiriusScore']\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'exp_int'] <- for_canP1[formulaRow4, 'explainedIntensity']\n",
    "#                     # other info\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'Result'] <- 'SIRIUS_STR'\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'dir'] <- str_can1\n",
    "\n",
    "#                     # normalize the CSI:FingerIDScores\n",
    "#                     norm_score1 <- feat_scale(str_canP1[,\"CSI:FingerIDScore\"])\n",
    "#                     # store the upper quartile\n",
    "#                     upper_quartile1 <- str_canP1[which(norm_score1 > as.numeric(quantile(norm_score1)[4])), \"smiles\"]\n",
    "#                     # if the upper quartile has more than 5 candidates, then just take the top 5 candidates\n",
    "#                     if (length(upper_quartile1) > 5){\n",
    "#                         upper_quartile1 <- upper_quartile1[1:5]\n",
    "#                     }\n",
    "#                     # save the top candidates SMILES, to check similarity later with rdkit in Python\n",
    "#                     msdata[as.numeric(rownames(rowMS)), 'SMILESforMCSS'] <- paste(upper_quartile1, collapse = '|')\n",
    "#                 }\n",
    "#                 #\n",
    "#                 else if (nrow(str_canP1) == 0){\n",
    "#                     if (file.exists(for_can1)){\n",
    "#                         for_canP2 <- as.data.frame(read_tsv(for_can1))\n",
    "#                         if (nrow(for_canP2)>= 1){\n",
    "\n",
    "#                             # information from formula candidate file\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- for_canP2[1, 'adduct']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'Formula'] <- for_canP2[1, 'molecularFormula']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'FormulaRank'] <- for_canP2[1, 'rank']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'exp_int'] <- for_canP2[1, 'explainedIntensity']\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'SIRIUSscore'] <- for_canP2[1, 'SiriusScore']\n",
    "#                             # other info\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'Result'] <- 'SIRIUS_FOR'\n",
    "#                             msdata[as.numeric(rownames(rowMS)), 'dir'] <- for_can1\n",
    "#                         }\n",
    "#                     }\n",
    "#                 }\n",
    "#             }\n",
    "\n",
    "#             # if the structure candidate file doesn't exists (and no str_sl exists)\n",
    "#             else{\n",
    "#                 # if formula candidate file exists\n",
    "#                 if (file.exists(for_can1)){\n",
    "#                     for_canP3 <- as.data.frame(read_tsv(for_can1))\n",
    "#                     if (nrow(for_canP3) >= 1){\n",
    "#                         # information from formula candidate file\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- for_canP3[1, 'adduct']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'Formula'] <- for_canP3[1, 'molecularFormula']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'FormulaRank'] <- for_canP3[1, 'rank']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'exp_int'] <- for_canP3[1, 'explainedIntensity']\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'SIRIUSscore'] <- for_canP3[1, 'SiriusScore']\n",
    "#                         # other info\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'Result'] <- 'SIRIUS_FOR'\n",
    "#                         msdata[as.numeric(rownames(rowMS)), 'dir'] <- for_can1\n",
    "#                     }\n",
    "#                 }\n",
    "#             }\n",
    "#         }\n",
    "#     }\n",
    "#     write.csv(msdata, paste(x, \"/insilico/MS1DATAsirius.csv\", sep = ''))\n",
    "#     return(msdata)\n",
    "# }\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "297cdddb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# input x is result from either ms1_peaks\n",
    "# SL is if a suspect list is present\n",
    "\n",
    "sirius_param <- function(x, result_dir, input_dir, SL = FALSE, collision_info = FALSE){\n",
    "\n",
    "    dir_name <- paste(input_dir, str_remove(paste(result_dir, \"/insilico/SIRIUS\", sep =\"\"), \".\"), sep = \"\")\n",
    "    if (!file.exists(dir_name)){\n",
    "        dir.create(dir_name, recursive = TRUE) ##create folder\n",
    "    }\n",
    "    isotopes <- c() #NA or isotope group number\n",
    "    sirius_param_file <- c() #input for SIRIUS\n",
    "    outputNames <- c() #output for SIRIUS with all db\n",
    "    outputNamesSL <- c() #output for SIRIUS with suspect list as db\n",
    "\n",
    "    parameter_file <- c()\n",
    "    par <- 0\n",
    "\n",
    "    x <- read.csv(x)\n",
    "\n",
    "    a <-0 # counting\n",
    "    y <- 0 # counting\n",
    "    z <- 0 # counting\n",
    "    for (i in 1:nrow(x)){\n",
    "\n",
    "        par <- par+1\n",
    "        para <- as.character(par) # for numbering\n",
    "\n",
    "        #no MS1 PEAKS and no ISOTOPES\n",
    "\n",
    "        if (x[i, \"ms1Peaks\"] == 'no ms1 peaks in QC'){\n",
    "\n",
    "            #INPUT FILE NAME\n",
    "            fileR <- paste(dir_name, \"/\" ,para, \"_NA_iso_NA_MS1p_\", x[i, \"premz\"], \"_SIRIUS_param.ms\", sep = \"\")\n",
    "\n",
    "            sirius_param_file <- c(sirius_param_file, fileR)\n",
    "            #ISOTOPE Information\n",
    "            isotopes <- c(isotopes, NA)\n",
    "            #OUTPUT\n",
    "            fileSR <- paste(str_sub(fileR, end=-4),'.json', sep = '')\n",
    "            fileSRS <- paste(str_sub(fileR, end=-4), 'SList.json', sep = '')\n",
    "            outputNames <- c(outputNames, fileSR)\n",
    "            outputNamesSL <- c(outputNamesSL, fileSRS)\n",
    "            file.create(fileR, recursive = TRUE)\n",
    "            file.conn <- file(fileR)\n",
    "            open(file.conn, open = \"at\")\n",
    "\n",
    "            #compound\n",
    "            writeLines(paste(\">compound\", x[i,\"id_X\"], sep=\" \"),con=file.conn)\n",
    "            #parentmass\n",
    "            writeLines(paste(\">parentmass\", x[i,\"premz\"], sep=\" \"),con=file.conn)\n",
    "            ##charge\n",
    "            if (x[i,\"pol\"] == \"pos\"){\n",
    "                writeLines(paste(\">charge\", \"+1\" ,sep=\" \"),con=file.conn)\n",
    "            }\n",
    "            else{\n",
    "                writeLines(paste(\">charge\", \"-1\" ,sep=\" \"),con=file.conn)\n",
    "            }\n",
    "            #rt\n",
    "            writeLines(paste(\">rt\", paste(x[i,\"rtmed\"], \"s\", sep =''), sep=\" \"),con=file.conn)\n",
    "\n",
    "            #ms1\n",
    "#             writeLines(\">ms1\",con=file.conn)\n",
    "#             writeLines(paste(x[i,\"premz\"], x[i,\"int\"] ,sep=\" \"),con=file.conn)\n",
    "\n",
    "            #ms2\n",
    "            if (collision_info){\n",
    "                writeLines(paste(\">collision\", paste(x[i,\"col_eng\"],\"eV\", sep =''),sep=\" \"),con=file.conn)\n",
    "            }else{\n",
    "                writeLines(\">ms2\" ,con=file.conn)\n",
    "            }\n",
    "            \n",
    "            ms2pk <- paste(input_dir, str_remove(x[i,\"ms2Peaks\"], \".\"), sep =\"\")\n",
    "\n",
    "            peak<- read.table(ms2pk)\n",
    "            for (k in 1:length(peak[,1])){\n",
    "                writeLines(paste(as.character(peak[k,1]),as.character(peak[k,2]), sep =\" \"), con=file.conn)\n",
    "            }\n",
    "            close(file.conn)\n",
    "            parameter_file <- c(parameter_file,file.conn)\n",
    "        }\n",
    "\n",
    "        # MS1 PEAKS and no ISOTOPES\n",
    "\n",
    "        else if (grepl(\"_no_isotopes.txt\", x[i, \"ms1Peaks\"], fixed=TRUE)){\n",
    "\n",
    "            #INPUT FILE NAME\n",
    "            fileR <- paste(dir_name, \"/\", para, \"_NA_iso_MS1p_\", x[i, \"premz\"], \"_SIRIUS_param.ms\", sep = \"\")\n",
    "            sirius_param_file <- c(sirius_param_file, fileR)\n",
    "            #ISOTOPE Information\n",
    "            isotopes <- c(isotopes, NA)\n",
    "            #OUTPUT\n",
    "            fileSR <- paste(str_sub(fileR, end=-4),'.json', sep = '')\n",
    "            fileSRS <- paste(str_sub(fileR, end=-4), 'SList.json', sep = '')\n",
    "            outputNames <- c(outputNames, fileSR)\n",
    "            outputNamesSL <- c(outputNamesSL, fileSRS)\n",
    "            file.create(fileR, recursive = TRUE)\n",
    "            file.conn <- file(fileR)\n",
    "            open(file.conn, open = \"at\")\n",
    "\n",
    "            #compound\n",
    "            writeLines(paste(\">compound\", x[i,\"id_X\"], sep=\" \"),con=file.conn)\n",
    "            #parentmass\n",
    "            writeLines(paste(\">parentmass\", x[i,\"premz\"], sep=\" \"),con=file.conn)\n",
    "            ##charge\n",
    "            if (x[i,\"pol\"] == \"pos\"){\n",
    "                writeLines(paste(\">charge\", \"+1\" ,sep=\" \"),con=file.conn)\n",
    "            }\n",
    "            else{\n",
    "                writeLines(paste(\">charge\", \"-1\" ,sep=\" \"),con=file.conn)\n",
    "            }\n",
    "            #rt\n",
    "            writeLines(paste(\">rt\", paste(x[i,\"rtmed\"], \"s\", sep =''), sep=\" \"),con=file.conn)\n",
    "\n",
    "            #ms1\n",
    "            writeLines(\">ms1\",con=file.conn)\n",
    "\n",
    "            ms1pk <- paste(input_dir, str_remove(x[i,\"ms1Peaks\"], \".\"), sep =\"\")\n",
    "            peakms1<- read.table(ms1pk)\n",
    "\n",
    "            for (l in 1:length(peakms1[,1])){\n",
    "                writeLines(paste(as.character(peakms1[l,1]),as.character(peakms1[l,2]), sep =\" \"), con=file.conn)\n",
    "            }\n",
    "\n",
    "            #ms2\n",
    "            if (collision_info){\n",
    "                writeLines(paste(\">collision\", paste(x[i,\"col_eng\"],\"eV\", sep =''),sep=\" \"),con=file.conn)\n",
    "            }else{\n",
    "                writeLines(\">ms2\" ,con=file.conn)\n",
    "            }\n",
    "\n",
    "            ms2pk <- paste(input_dir, str_remove(x[i,\"ms2Peaks\"], \".\"), sep =\"\")\n",
    "\n",
    "            peakms2<- read.table(ms2pk)\n",
    "\n",
    "            for (k in 1:length(peakms2[,1])){\n",
    "                writeLines(paste(as.character(peakms2[k,1]),as.character(peakms2[k,2]), sep =\" \"), con=file.conn)\n",
    "            }\n",
    "\n",
    "            close(file.conn)\n",
    "            parameter_file <- c(parameter_file,file.conn)\n",
    "        }\n",
    "\n",
    "        # MS1 PEAKS and ISOTOPES\n",
    "\n",
    "        else if (grepl(\"_isotopeNum_\", x[i, \"ms1Peaks\"], fixed=TRUE)){\n",
    "\n",
    "            #INPUT FILE NAME\n",
    "            fileR <- paste(dir_name, \"/\", para, \"_iso_MS1p_\", as.character(x[i, \"premz\"]), \"_SIRIUS_param.ms\", sep = \"\")\n",
    "            sirius_param_file <- c(sirius_param_file, fileR)\n",
    "            #ISOTOPE Information\n",
    "            isotopes <- c(isotopes, \"present\")\n",
    "            #OUTPUT\n",
    "            fileSR <- paste(str_sub(fileR, end=-4),'.json', sep = '')\n",
    "            fileSRS <- paste(str_sub(fileR, end=-4), 'SList.json', sep = '')\n",
    "            outputNames <- c(outputNames, fileSR)\n",
    "            outputNamesSL <- c(outputNamesSL, fileSRS)\n",
    "            file.create(fileR, recursive = TRUE)\n",
    "            file.conn <- file(fileR)\n",
    "            open(file.conn, open = \"at\")\n",
    "\n",
    "             #compound\n",
    "            writeLines(paste(\">compound\", x[i,\"id_X\"], sep=\" \"),con=file.conn)\n",
    "            #parentmass\n",
    "            writeLines(paste(\">parentmass\", x[i,\"premz\"], sep=\" \"),con=file.conn)\n",
    "            ##charge\n",
    "            if (x[i,\"pol\"] == \"pos\"){\n",
    "                writeLines(paste(\">charge\", \"+1\" ,sep=\" \"),con=file.conn)\n",
    "            }\n",
    "            else{\n",
    "                writeLines(paste(\">charge\", \"-1\" ,sep=\" \"),con=file.conn)\n",
    "            }\n",
    "            #rt\n",
    "            writeLines(paste(\">rt\", paste(x[i,\"rtmed\"], \"s\", sep =''), sep=\" \"),con=file.conn)\n",
    "\n",
    "\n",
    "            #ms1\n",
    "            writeLines(\">ms1\",con=file.conn)\n",
    "\n",
    "            ms1pk <- paste(input_dir, str_remove(x[i,\"ms1Peaks\"], \".\"), sep =\"\")\n",
    "            peakms1<- read.table(ms1pk)\n",
    "\n",
    "            for (l in 1:length(peakms1[,1])){\n",
    "                writeLines(paste(as.character(peakms1[l,1]),as.character(peakms1[l,2]), sep =\" \"), con=file.conn)\n",
    "            }\n",
    "\n",
    "            #ms2\n",
    "            if (collision_info){\n",
    "                writeLines(paste(\">collision\", paste(x[i,\"col_eng\"],\"eV\", sep =''),sep=\" \"),con=file.conn)\n",
    "            }else{\n",
    "                writeLines(\">ms2\" ,con=file.conn)\n",
    "            }\n",
    "\n",
    "            ms2pk <- paste(input_dir, str_remove(x[i,\"ms2Peaks\"], \".\"), sep =\"\")\n",
    "\n",
    "            peakms2<- read.table(ms2pk)\n",
    "\n",
    "            for (k in 1:length(peakms2[,1])){\n",
    "                writeLines(paste(as.character(peakms2[k,1]),as.character(peakms2[k,2]), sep =\" \"), con=file.conn)\n",
    "            }\n",
    "\n",
    "            close(file.conn)\n",
    "            parameter_file <- c(parameter_file,file.conn)\n",
    "\n",
    "        }\n",
    "    }\n",
    "    if (SL){\n",
    "\n",
    "        in_out_file <- data.frame(cbind(sirius_param_file, outputNames, outputNamesSL, isotopes))\n",
    "\n",
    "        write.table(in_out_file, paste(input_dir, str_remove(paste(result_dir,'/insilico/MS1DATA_SiriusPandSL.tsv', sep = \"\"), \".\"), sep =\"\"), sep = \"\\t\")\n",
    "        return(in_out_file)\n",
    "\n",
    "    }\n",
    "    else{\n",
    "        in_out_file <- data.frame(cbind(sirius_param_file, outputNames, isotopes))\n",
    "\n",
    "        write.table(in_out_file, paste(input_dir, str_remove(paste(result_dir,'/insilico/MS1DATA_SiriusP.tsv', sep = \"\"), \".\"), sep =\"\"), sep = \"\\t\")\n",
    "        return(in_out_file)\n",
    "\n",
    "    }\n",
    "\n",
    "}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "1ec877be",
   "metadata": {},
   "outputs": [],
   "source": [
    "run_sirius <- function(files, ppm_max = 5, ppm_max_ms2 = 15, QC = TRUE, SL = TRUE, SL_path, candidates = 30, profile, db){\n",
    "    files <- read.csv(files, sep = \"\\t\")\n",
    "\n",
    "    for (b in 1:nrow(files)){\n",
    "\n",
    "        if (QC){\n",
    "            if (is.na(files[b, \"isotopes\"])){\n",
    "                system(paste(\"sirius --input\", files[b, \"sirius_param_file\"], \"--output\", files[b, \"outputNames\"],\n",
    "                             \"formula --profile\", profile, \"--no-isotope-filter --no-isotope-score --candidates\", candidates, \"--ppm-max\", ppm_max, \"--ppm-max-ms2\",  ppm_max_ms2,\"structure --database\", db ,\"canopus\",\n",
    "                             sep = \" \"))\n",
    "                if(!file.exists(files[b, \"outputNames\"])){\n",
    "                    system(paste(\"sirius --input\", files[b, \"sirius_param_file\"], \"--output\", files[b, \"outputNames\"],\n",
    "                                 \"formula --profile\", profile, \"--no-isotope-filter --no-isotope-score --candidates\", candidates, \"--ppm-max\", ppm_max, \"--ppm-max-ms2\",  ppm_max_ms2,\"structure --database\", db ,\"canopus\",\n",
    "                                 sep = \" \"))\n",
    "                }\n",
    "                if(SL){\n",
    "                    system(paste(\"sirius --input\", files[b, \"sirius_param_file\"], \"--output\", files[b, \"outputNamesSL\"],\n",
    "                             \"formula --profile\", profile, \"--no-isotope-filter --no-isotope-score --candidates\", candidates, \"--ppm-max\", ppm_max, \"--ppm-max-ms2\",  ppm_max_ms2,\"structure --database\",  SL_path ,\"canopus\",\n",
    "                             sep = \" \"))\n",
    "                }\n",
    "\n",
    "            }\n",
    "            else if(files[b, \"isotopes\"] == \"present\"){\n",
    "                system(paste(\"sirius --input\", files[b, \"sirius_param_file\"], \"--output\", files[b, \"outputNames\"],\n",
    "                             \"formula --profile\", profile, \"--candidates\", candidates, \"--ppm-max\", ppm_max,\"--ppm-max-ms2\", ppm_max_ms2,\"structure --database\", db ,\"canopus\",\n",
    "                             sep = \" \"))\n",
    "                if(!file.exists(files[b, \"outputNames\"])){\n",
    "                     system(paste(\"sirius --input\", files[b, \"sirius_param_file\"], \"--output\", files[b, \"outputNames\"],\n",
    "                                 \"formula --profile\", profile, \"--candidates\", candidates, \"--ppm-max\", ppm_max,\"--ppm-max-ms2\", ppm_max_ms2,\"structure --database\", db ,\"canopus\",\n",
    "                                 sep = \" \"))\n",
    "                }\n",
    "                if(SL){\n",
    "                    system(paste(\"sirius --input\", files[b, \"sirius_param_file\"], \"--output\", files[b, \"outputNamesSL\"],\n",
    "                             \"formula --profile\", profile, \"--candidates\", candidates, \"--ppm-max\", ppm_max, \"--ppm-max-ms2\",  ppm_max_ms2,\"structure --database\",  SL_path ,\"canopus\",\n",
    "                             sep = \" \"))\n",
    "                }\n",
    "            }\n",
    "        }\n",
    "        else{\n",
    "            system(paste(\"sirius --input\", files[b, \"sirius_param_file\"], \"--output\", files[b, \"outputNames\"],\n",
    "                        \"formula --profile\", profile, \"--no-isotope-filter --no-isotope-score --candidates\", candidates, \"--ppm-max\", ppm_max, \"--ppm-max-ms2\",  ppm_max_ms2,\"structure --database\", db ,\"canopus\",\n",
    "                        sep = \" \"))\n",
    "            if(!file.exists(files[b, \"outputNames\"])){\n",
    "                system(paste(\"sirius --input\", files[b, \"sirius_param_file\"], \"--output\", files[b, \"outputNames\"],\n",
    "                             \"formula --profile\", profile, \"--no-isotope-filter --no-isotope-score --candidates\", candidates, \"--ppm-max\", ppm_max, \"--ppm-max-ms2\",  ppm_max_ms2,\"structure --database\", db ,\"canopus\",\n",
    "                             sep = \" \"))\n",
    "            }\n",
    "            if(SL){\n",
    "                system(paste(\"sirius --input\", files[b, \"sirius_param_file\"], \"--output\", files[b, \"outputNamesSL\"],\n",
    "                            \"formula --profile\", profile, \"--no-isotope-filter --no-isotope-score --candidates\", candidates, \"--ppm-max\", ppm_max, \"--ppm-max-ms2\",  ppm_max_ms2,\"structure --database\",  SL_path ,\"canopus\",\n",
    "                            sep = \" \"))\n",
    "            }\n",
    "        }\n",
    "        Sys.sleep(5)\n",
    "\n",
    "    }\n",
    "\n",
    "}\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b82776b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "sirius_adduct <- function(input_dir, x, SL = TRUE){\n",
    "\n",
    "    feat_scale <- function(p) {\n",
    "    (p - min(p)) / (max(p) - min(p))\n",
    "    }\n",
    "\n",
    "    #the result directory name for each file\n",
    "    dir_name <- paste(input_dir, str_remove(x, \".\"),'/insilico/SIRIUS', sep =\"\")\n",
    "\n",
    "    #for all DB\n",
    "    Param <- list.files(dir_name, pattern = 'param.json', full.names = TRUE)\n",
    "    #read the msdata csv file that contains all features and their data\n",
    "    msdatacsv <- read.csv(paste(input_dir, str_remove(x, \".\"), \"/insilico/MS1DATA.csv\", sep = ''))\n",
    "    msdata <- as.data.frame(msdatacsv)\n",
    "    #add new empty columns to store information from SIRIUS results\n",
    "    msdata$Adducts <- NA           #adducts\n",
    "\n",
    "    # for all entries in the json result files for each each feature\n",
    "    for (i in Param){\n",
    "\n",
    "        # find the number of the row corresponding to the row in msdata that has the same precursor m/z as json result files\n",
    "        rowMS <- msdata[grepl(str_match(as.character(i), \"MS1p_\\\\s*(.*?)\\\\s*_SIRIUS\")[2] ,msdata$premz), ]\n",
    "\n",
    "        # directory for formula candidate file\n",
    "        for_can1 <- paste(list.dirs(i)[2], '/formula_candidates.tsv', sep = '')\n",
    "        if (file.exists(for_can1)){\n",
    "            # read the formula file from All Dbs (PARAM) json file\n",
    "            for_canP1 <- as.data.frame(read_tsv(for_can1))\n",
    "\n",
    "            if (nrow(for_canP1)>= 1){\n",
    "\n",
    "                # information from formula candidate file\n",
    "                msdata[as.numeric(rownames(rowMS)), 'Adducts'] <- for_canP1[1, 'adduct']\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "    write.csv(msdata, paste(input_dir, str_remove(x, \".\"), \"/insilico/MS1DATAsirius.csv\", sep = ''))\n",
    "    return(msdata)\n",
    "}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4402ab15",
   "metadata": {},
   "outputs": [],
   "source": [
    "metfrag_param <- function(x, result_dir, input_dir, sl_mtfrag, SL = TRUE, ppm_max = 5, ppm_max_ms2= 15){\n",
    "\n",
    "\n",
    "\n",
    "    if (file.exists(paste(input_dir, \"/MetFrag_AdductTypes.csv\", sep = \"\"))){\n",
    "        adducts <- paste(input_dir, \"/MetFrag_AdductTypes.csv\", sep = \"\")\n",
    "    }\n",
    "    else{\n",
    "        system(paste(\"wget -P\",\n",
    "                    input_dir,\n",
    "                    \"https://github.com/schymane/ReSOLUTION/raw/master/inst/extdata/MetFrag_AdductTypes.csv\",\n",
    "                    sep = \" \"))\n",
    "        adducts <- paste(input_dir, \"/MetFrag_AdductTypes.csv\", sep = \"\")\n",
    "    }\n",
    "\n",
    "    x <- read.csv(x)\n",
    "\n",
    "    dir_name <- paste(input_dir, str_remove(paste(result_dir, \"/insilico/MetFrag\", sep =\"\"), \".\"), sep = \"\")\n",
    "\n",
    "    if (!file.exists(dir_name)){\n",
    "        dir.create(dir_name, recursive = TRUE) ##create folder\n",
    "    }\n",
    "\n",
    "    db <- \"PubChem\"\n",
    "\n",
    "    parameter_file <- c()\n",
    "    par <- 0\n",
    "    metfrag_param_file <- c()\n",
    "\n",
    "    AdductsMF <- data.frame(read.csv(adducts))\n",
    "\n",
    "    for (j in 1:nrow(x)){\n",
    "        if (!(is.na(x[j, 'Adducts']))){\n",
    "            for (k in db){\n",
    "                par <- par+1\n",
    "                para <- as.character(par)\n",
    "                fileR <- paste(dir_name, \"/\", para, \"_id_\", x[j, 'id_X'], \"_mz_\", x[j, 'premz'], \"_rt_\", x[j, 'rtmed'], \"_db_\", k, \".txt\", sep = '')\n",
    "                metfrag_param_file <- c(metfrag_param_file, fileR)\n",
    "\n",
    "                file.create(fileR, recursive = TRUE)\n",
    "                file.conn <- file(fileR)\n",
    "                open(file.conn, open = \"at\")\n",
    "\n",
    "\n",
    "                peakspath <- str_replace(x[j, \"ms2Peaks\"], \".\", input_dir)\n",
    "\n",
    "\n",
    "                #writeLines(paste(\"PeakListPath = \",as.character(peakspath),sep=\"\"),con=file.conn)\n",
    "                writeLines(paste(\"PeakListPath = \",peakspath, sep=\"\"),con=file.conn)\n",
    "                writeLines(paste(\"IonizedPrecursorMass = \", x[j, \"premz\"]),con = file.conn)\n",
    "\n",
    "                # write code here\n",
    "\n",
    "\n",
    "                PrecursorIonMode <- AdductsMF[which(AdductsMF[, \"PrecursorIonType\"] == gsub(\"[[:blank:]]\", \"\", x[j, 'Adducts'])), \"PrecursorIonMode\"]\n",
    "                IsPositiveIonMode <- AdductsMF[which(AdductsMF[, \"PrecursorIonType\"] == gsub(\"[[:blank:]]\", \"\", x[j, 'Adducts'])), \"IsPositiveIonMode\"]\n",
    "\n",
    "\n",
    "                writeLines(paste(\"PrecursorIonMode = \", PrecursorIonMode, sep = ''), con = file.conn)\n",
    "                writeLines(paste(\"IsPositiveIonMode = \", IsPositiveIonMode, sep = ''), con = file.conn)\n",
    "\n",
    "                writeLines(paste(\"MetFragDatabaseType = \", k),con = file.conn)\n",
    "\n",
    "                writeLines(paste(\"DatabaseSearchRelativeMassDeviation = \", ppm_max, sep = ''),con=file.conn)\n",
    "                writeLines(\"FragmentPeakMatchAbsoluteMassDeviation = 0.001\",con=file.conn)\n",
    "                writeLines(paste(\"FragmentPeakMatchRelativeMassDeviation = \", ppm_max_ms2, sep = ''),con=file.conn)\n",
    "\n",
    "                if (SL){\n",
    "                    writeLines(paste(\"ScoreSuspectLists = \", sl_mtfrag),con=file.conn)\n",
    "\n",
    "                    writeLines(\"MetFragScoreTypes = FragmenterScore, SuspectListScore\",con=file.conn)\n",
    "                    writeLines(\"MetFragScoreWeights = 1.0, 1.0\", con=file.conn)\n",
    "                }\n",
    "                else{\n",
    "                    writeLines(\"MetFragScoreTypes = FragmenterScore\", con=file.conn)\n",
    "                    writeLines(\"MetFragScoreWeights = 1.0\", con=file.conn)\n",
    "                }\n",
    "\n",
    "\n",
    "                writeLines(\"MetFragCandidateWriter = CSV\",con=file.conn)\n",
    "\n",
    "                writeLines(paste(\"SampleName = \", para, \"_id_\", x[j, 'id_X'], \"_mz_\", x[j, 'premz'], \"_rt_\", x[j, 'rtmed'], \"_db_\", k, sep = ''),con=file.conn)\n",
    "                resultspath <- str_replace(result_dir, \".\", input_dir)\n",
    "                writeLines(paste(\"ResultsPath = \", resultspath, \"/insilico/MetFrag/\", sep = ''),con=file.conn)\n",
    "\n",
    "                writeLines(\"MetFragPreProcessingCandidateFilter = UnconnectedCompoundFilter\",con=file.conn)\n",
    "                writeLines(\"MetFragPostProcessingCandidateFilter = InChIKeyFilter\",con=file.conn)\n",
    "                writeLines(\"#MaximumTreeDepth = 2\",con=file.conn)\n",
    "                writeLines(\"#NumberThreads = 1\",con=file.conn)\n",
    "\n",
    "                close(file.conn)\n",
    "                parameter_file <- c(parameter_file,file.conn)\n",
    "\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "\n",
    "    write.table(metfrag_param_file, file = paste(input_dir, str_remove(paste(result_dir, \"/insilico/metparam_list.txt\", sep =\"\"), \".\"), sep = \"\"), sep = \"/t\", row.names = FALSE, col.names = FALSE)\n",
    "    return(metfrag_param_file)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c4f64c50",
   "metadata": {},
   "outputs": [],
   "source": [
    "run_metfrag <- function(met_param){\n",
    "    #line <- system.time({\n",
    "\n",
    "        if (file.exists(paste(input_dir, \"/MetFragCommandLine-2.5.0.jar\", sep = \"\"))){\n",
    "            MetFragjarFile <- paste(input_dir, \"/MetFragCommandLine-2.5.0.jar\", sep = \"\")\n",
    "        }\n",
    "        else{\n",
    "            system(paste(\"wget -P\",\n",
    "                        input_dir,\n",
    "                        \"https://github.com/ipb-halle/MetFragRelaunched/releases/download/v.2.5.0/MetFragCommandLine-2.5.0.jar\",\n",
    "                        sep = \" \"))\n",
    "            MetFragjarFile <- paste(input_dir, \"/MetFragCommandLine-2.5.0.jar\", sep = \"\")\n",
    "        }\n",
    "\n",
    "        start_time <- Sys.time()\n",
    "\n",
    "        filesmet_param <- read.table(met_param)\n",
    "\n",
    "        for (files in filesmet_param[[1]]){\n",
    "            system(paste(\"java -jar\",  MetFragjarFile , files))\n",
    "            Sys.sleep(5)\n",
    "        }\n",
    "    #})\n",
    "\n",
    "#     if (file.exists(paste(input_dir, \"/summaryFile.txt\", sep = \"\"))){\n",
    "\n",
    "#         line1 = paste(\"Compound Database Dereplication with MetFrag took\", line[1], \"of seconds.\", sep = \" \")\n",
    "#         write(line1 ,file=paste(input_dir, \"/summaryFile.txt\", sep = \"\"),append=TRUE)\n",
    "\n",
    "#     }\n",
    "#     else{\n",
    "#         fileConn<-file(paste(input_dir, \"/summaryFile.txt\", sep = \"\"))\n",
    "\n",
    "#         writeLines(paste(\"Compound Database Dereplication with Metfrag took\", line[1], \"of seconds.\", sep = \" \"), fileConn)\n",
    "\n",
    "#         close(fileConn)\n",
    "#     }\n",
    "\n",
    "}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22f73f7f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
