if(libgen_method_is_enabled("dda",params)){
    deconvolution_methods += [output: { dda_files_ch } ]
}

if(libgen_method_is_enabled("dda",params)){
    Channel.fromPath(params.ddafiles).tap(dda_files_ch)
}

// fragment tolarance for diam2spep in ppm
// (other tools require it in dalton)
params.diams2pep_fragment_tolerance = null

if(libgen_method_is_enabled("diams2pep",params)){
    deconvolution_methods += [output: { diams2pep_pseudospectra},
                              input: { diams2pep_input_mzml}]
}

if(libgen_method_is_enabled("diams2pep",params)){

process convert_for_DIAMS2PEP {
    input:
    file mzml from diams2pep_input_mzml
    output:
    // there is no good way in nextflow that makes a UUUID that persists across -resume things
    // task.hash is forgotten in resume, as is task.id.
    tuple val("${mzml.baseName}"), path(ofile) into diams2pep_mgf_mzml, diams2pep_window_mzml, diams2pep_for_pseudo_mzml
    script:
    ofile="converted/${mzml.baseName}.mzML"
    """
    mkdir -p converted
    msconvert --mzML --mz64 --zlib --inten64 --simAsSpectra --filter "peakPicking cwt msLevel=1-2" --outdir converted $mzml
    """
}

process convert_mgf_for_DIAMS2PEP {
    input:
    tuple val(hash), path(mzml) from diams2pep_mgf_mzml
    output:
    tuple val(hash), path("${mzml.baseName}.mgf") into diams2pep_mgf

    """
    msconvert --mgf $mzml
    """
}

process DIAMS2PEP_window {
    input:
    tuple val(hash), path(mzml) from  diams2pep_window_mzml
    output:
    tuple val(hash), path("${mzml}.DIA_acquisition_window.txt") into diams2pep_window

    """
    DIA_acquistion_window_generator.pl $mzml
    """
}

if(params.diams2pep_fragment_tolerance == null)
    raise RunTimeError("DIAMSM2PEP enabled but no diams2pep_fragment_tolerance specified.")
process DIAMS2PEP_generate_pseudo {
    input:
    tuple val(hash), path(mzml), path(mgf), path(acq_window) from diams2pep_for_pseudo_mzml.join(diams2pep_mgf).join(diams2pep_window)
    val tolerance from Channel.value(params.diams2pep_fragment_tolerance)
    output:
    file "mgf-output/*.mgf" into diams2pep_pseudospectra_mgf mode flatten
    """
    mkdir -p mgf-output
    DIA_pesudo_MS2_multiforks.pl ${mzml.baseName} mgf-output $tolerance ${task.cpus}
    """
}

// this will use the default container because we need msconvert
process MgfToMzml_DIAMS2PEP {
    input:
    file mgf from diams2pep_pseudospectra_mgf
    output:
    file "*.mzXML" into diams2pep_pseudospectra
    """
    msconvert --mzXML $mgf
    """
}

} // end of diams2pep guard

params.search_engines = ["comet","xtandem"]

// Float or Int; in ppm ; eg. params.precursor_mass_tolerance=10
params.precursor_mass_tolerance=null  
// Float or Int; in Dalton; e.g. parames.fragment_mass_tolerance=0.2
params.fragment_mass_tolerance=null

// Int, if you are using the comet, this can by at maximum 5,
params.max_missed_cleavages=null

max_missed_cleavages = Channel.value(params.max_missed_cleavages)

params.search_engines = ["comet"]
params.max_missed_cleavages= 2

# comet_version 2022.01 rev. 0
# Comet MS/MS search engine parameters file.
# Everything following the '#' symbol is treated as a comment.

database_name = @DDA_DB_FILE@
decoy_search = 0                       # 0=no (default), 1=concatenated search, 2=separate search
peff_format = 0                        # 0=no (normal fasta, default), 1=PEFF PSI-MOD, 2=PEFF Unimod
peff_obo =                             # path to PSI Mod or Unimod OBO file

num_threads = 0                        # 0=poll CPU to set num threads; else specify num threads directly (max 128)

#
# masses
#
peptide_mass_tolerance = @PRECURSOR_MASS_TOLERANCE@
peptide_mass_units = 2                 # 0=amu, 1=mmu, 2=ppm
mass_type_parent = 1                   # 0=average masses, 1=monoisotopic masses
mass_type_fragment = 1                 # 0=average masses, 1=monoisotopic masses
precursor_tolerance_type = 1           # 0=MH+ (default), 1=precursor m/z; only valid for amu/mmu tolerances
isotope_error = 3                      # 0=off, 1=0/1 (C13 error), 2=0/1/2, 3=0/1/2/3, 4=-8/-4/0/4/8 (for +4/+8 labeling)

#
# search enzyme
#
search_enzyme_number = 1               # choose from list at end of this params file
search_enzyme2_number = 0              # second enzyme; set to 0 if no second enzyme
num_enzyme_termini = 2                 # 1 (semi-digested), 2 (fully digested, default), 8 C-term unspecific , 9 N-term unspecific
allowed_missed_cleavage = @MAX_MISSED_CLEAVAGES@            # maximum value is 5; for enzyme search

#
# Up to 9 variable modifications are supported
# format:  <mass> <residues> <0=variable/else binary> <max_mods_per_peptide> <term_distance> <n/c-term> <required> <neutral_loss>
#     e.g. 79.966331 STY 0 3 -1 0 0 97.976896
#
variable_mod01 = 15.9949 M 0 3 -1 0 0 0.0
variable_mod02 = 0.0 X 0 3 -1 0 0 0.0
variable_mod03 = 0.0 X 0 3 -1 0 0 0.0
variable_mod04 = 0.0 X 0 3 -1 0 0 0.0
variable_mod05 = 0.0 X 0 3 -1 0 0 0.0
variable_mod06 = 0.0 X 0 3 -1 0 0 0.0
variable_mod07 = 0.0 X 0 3 -1 0 0 0.0
variable_mod08 = 0.0 X 0 3 -1 0 0 0.0
variable_mod09 = 0.0 X 0 3 -1 0 0 0.0
max_variable_mods_in_peptide = 5
require_variable_mod = 0

#
# fragment ions
#
# ion trap ms/ms:  1.0005 tolerance, 0.4 offset (mono masses), theoretical_fragment_ions = 1
# high res ms/ms:    0.02 tolerance, 0.0 offset (mono masses), theoretical_fragment_ions = 0, spectrum_batch_size = 15000
#
fragment_bin_tol = @FRAGMENT_MASS_TOLERANCE@              # binning to use on fragment ions
fragment_bin_offset = 0.0              # offset position to start the binning (0.0 to 1.0)
theoretical_fragment_ions = 1          # 0=use flanking peaks, 1=M peak only
use_A_ions = 0
use_B_ions = 1
use_C_ions = 0
use_X_ions = 0
use_Y_ions = 1
use_Z_ions = 0
use_Z1_ions = 0
use_NL_ions = 0                        # 0=no, 1=yes to consider NH3/H2O neutral loss peaks

#
# output
#
output_sqtfile = 0                     # 0=no, 1=yes  write sqt file
output_txtfile = 0                     # 0=no, 1=yes  write tab-delimited txt file
output_pepxmlfile = 1                  # 0=no, 1=yes  write pepXML file
output_mzidentmlfile = 0               # 0=no, 1=yes  write mzIdentML file
output_percolatorfile = 1              # 0=no, 1=yes  write Percolator pin file
print_expect_score = 1                 # 0=no, 1=yes to replace Sp with expect in out & sqt
num_output_lines = 5                   # num peptide results to show

sample_enzyme_number = 1               # Sample enzyme which is possibly different than the one applied to the search.
                                       # Used to calculate NTT & NMC in pepXML output (default=1 for trypsin).

#
# mzXML parameters
#
scan_range = 0 0                       # start and end scan range to search; either entry can be set independently
precursor_charge = 0 0                 # precursor charge range to analyze; does not override any existing charge; 0 as 1st entry ignores parameter
override_charge = 0                    # 0=no, 1=override precursor charge states, 2=ignore precursor charges outside precursor_charge range, 3=see online
ms_level = 2                           # MS level to analyze, valid are levels 2 (default) or 3
activation_method = HCD                # activation method; used if activation method set; allowed ALL, CID, ECD, ETD, ETD+SA, PQD, HCD, IRMPD, SID

#
# misc parameters
#
digest_mass_range = 600.0 5000.0       # MH+ peptide mass range to analyze
peptide_length_range = 5 63            # minimum and maximum peptide length to analyze (default 1 63; max length 63)
num_results = 100                      # number of search hits to store internally
max_duplicate_proteins = 20            # maximum number of additional duplicate protein names to report for each peptide ID; -1 reports all duplicates
max_fragment_charge = 3                # set maximum fragment charge state to analyze (allowed max 5)
max_precursor_charge = 6               # set maximum precursor charge state to analyze (allowed max 9)
nucleotide_reading_frame = 0           # 0=proteinDB, 1-6, 7=forward three, 8=reverse three, 9=all six
clip_nterm_methionine = 0              # 0=leave protein sequences as-is; 1=also consider sequence w/o N-term methionine
spectrum_batch_size = 15000            # max. # of spectra to search at a time; 0 to search the entire scan range in one loop
decoy_prefix = DECOY_                  # decoy entries are denoted by this string which is pre-pended to each protein accession
equal_I_and_L = 1                      # 0=treat I and L as different; 1=treat I and L as same
output_suffix =                        # add a suffix to output base names i.e. suffix "-C" generates base-C.pep.xml from base.mzXML input
mass_offsets =                         # one or more mass offsets to search (values substracted from deconvoluted precursor mass)
precursor_NL_ions =                    # one or more precursor neutral loss masses, will be added to xcorr analysis

#
# spectral processing
#
minimum_peaks = 10                     # required minimum number of peaks in spectrum to search (default 10)
minimum_intensity = 0                  # minimum intensity value to read in
remove_precursor_peak = 0              # 0=no, 1=yes, 2=all charge reduced precursor peaks (for ETD), 3=phosphate neutral loss peaks
remove_precursor_tolerance = 1.5       # +- Da tolerance for precursor removal
clear_mz_range = 0.0 0.0               # for iTRAQ/TMT type data; will clear out all peaks in the specified m/z range

#
# additional modifications
#

add_Cterm_peptide = 0.0
add_Nterm_peptide = 0.0
add_Cterm_protein = 0.0
add_Nterm_protein = 0.0

add_G_glycine = 0.0000                 # added to G - avg.  57.0513, mono.  57.02146
add_A_alanine = 0.0000                 # added to A - avg.  71.0779, mono.  71.03711
add_S_serine = 0.0000                  # added to S - avg.  87.0773, mono.  87.03203
add_P_proline = 0.0000                 # added to P - avg.  97.1152, mono.  97.05276
add_V_valine = 0.0000                  # added to V - avg.  99.1311, mono.  99.06841
add_T_threonine = 0.0000               # added to T - avg. 101.1038, mono. 101.04768
add_C_cysteine = 57.021464             # added to C - avg. 103.1429, mono. 103.00918
add_L_leucine = 0.0000                 # added to L - avg. 113.1576, mono. 113.08406
add_I_isoleucine = 0.0000              # added to I - avg. 113.1576, mono. 113.08406
add_N_asparagine = 0.0000              # added to N - avg. 114.1026, mono. 114.04293
add_D_aspartic_acid = 0.0000           # added to D - avg. 115.0874, mono. 115.02694
add_Q_glutamine = 0.0000               # added to Q - avg. 128.1292, mono. 128.05858
add_K_lysine = 0.0000                  # added to K - avg. 128.1723, mono. 128.09496
add_E_glutamic_acid = 0.0000           # added to E - avg. 129.1140, mono. 129.04259
add_M_methionine = 0.0000              # added to M - avg. 131.1961, mono. 131.04048
add_H_histidine = 0.0000               # added to H - avg. 137.1393, mono. 137.05891
add_F_phenylalanine = 0.0000           # added to F - avg. 147.1739, mono. 147.06841
add_U_selenocysteine = 0.0000          # added to U - avg. 150.0379, mono. 150.95363
add_R_arginine = 0.0000                # added to R - avg. 156.1857, mono. 156.10111
add_Y_tyrosine = 0.0000                # added to Y - avg. 163.0633, mono. 163.06333
add_W_tryptophan = 0.0000              # added to W - avg. 186.0793, mono. 186.07931
add_O_pyrrolysine = 0.0000             # added to O - avg. 237.2982, mono  237.14773
add_B_user_amino_acid = 0.0000         # added to B - avg.   0.0000, mono.   0.00000
add_J_user_amino_acid = 0.0000         # added to J - avg.   0.0000, mono.   0.00000
add_X_user_amino_acid = 0.0000         # added to X - avg.   0.0000, mono.   0.00000
add_Z_user_amino_acid = 0.0000         # added to Z - avg.   0.0000, mono.   0.00000

#
# COMET_ENZYME_INFO _must_ be at the end of this parameters file
#
[COMET_ENZYME_INFO]
0.  Cut_everywhere         0      -           -
1.  Trypsin                1      KR          P
2.  Trypsin/P              1      KR          -
3.  Lys_C                  1      K           P
4.  Lys_N                  0      K           -
5.  Arg_C                  1      R           P
6.  Asp_N                  0      D           -
7.  CNBr                   1      M           -
8.  Glu_C                  1      DE          P
9.  PepsinA                1      FL          P
10. Chymotrypsin           1      FWYL        P
11. No_cut                 1      @           @

params.comet_template="comet_template.txt"

process MakeCometConfig {
    // should we instead return a tuple here of fastadb and config
    // because the config.txt refers to it?
    input:
    val max_missed_cleavages
    file fastadb_with_decoy from joined_fasta_with_decoys_ch.first()
    path template from Channel.fromPath(params.comet_template)
    output:
    file "comet_config.txt" into comet_config_ch
    """
    sed 's/@DDA_DB_FILE@/$fastadb_with_decoy/g;s/@FRAGMENT_MASS_TOLERANCE@/$params.fragment_mass_tolerance/g;s/@PRECURSOR_MASS_TOLERANCE@/$params.precursor_mass_tolerance/g;s/@MAX_MISSED_CLEAVAGES@/$max_missed_cleavages/g' $template > comet_config.txt 
    """

}

process Comet {
    // we probably also want to publish thees
    memory { 5.GB * 2 *  task.attempt }
    errorStrategy { task.exitStatus in 137..137 ? 'retry' : 'terminate' }
    maxRetries 2
    input:
    file comet_config from comet_config_ch.first()
    // future dev: we can .mix with DDA here?
    // though we might need to tag for DDA / Pseudo
    // so that xinteract 
    file mzxml from maybespectra_ch.comet
    file fastadb_with_decoy from joined_fasta_with_decoys_ch.first()
    output:
    file("${mzxml.baseName}.pep.xml") into comet_pepxml_ch
    file mzxml into xinteract_comet_mzxml_ch
    when:
    params.search_engines.contains("comet")

    """
    if command -v command-ms;
    then
      comet-ms -P$comet_config $mzxml
    else
      comet -P$comet_config $mzxml
    fi
    """
}

process XinteractComet {
    memory '16 GB'
    time '5h'
    // memory usage scales with the number of input files
    // find the correct usage per input file or size
    // also for xinteractxtandem
    // usage there seems to be a lot smaller
    // as input files seems to be smaller
    input:
    file pepxmls from comet_pepxml_ch.toSortedList()
    // the filename of needed fastdadb was defined in cometcfg
    // and stored in pepxml in the comet-ms step
    // -a suppplies the absulute path to the data directory where the mzxmls
    // rather than reading wherer the mfrom the xmls
    // where the mzxml are, because its not very
    // nextflow to look outside the cwd.
    file fastadb_with_decoy from joined_fasta_with_decoys_ch.first()
    file mzxmls from  xinteract_comet_mzxml_ch.toSortedList()
    output: 
    file "interact_comet.pep.xml" into comet_search_results_ch
    when:
    pepxmls.size() > 0
    """
    xinteract -a\$PWD -OARPd -dDECOY_ -Ninteract_comet.pep.xml $pepxmls
    """
}

<?xml version="1.0"?>
<bioml label="x! taxon-to-file matching list">
  <taxon label="DB">
    <file format="peptide" URL="%s" />
  </taxon>
</bioml>

<?xml version="1.0"?>
<?xml-stylesheet type="text/xsl" href="tandem-input-style.xsl"?>
<bioml>
<note>list path parameters</note>

<note>spectrum parameters</note>
        <note type="input" label="spectrum, fragment monoisotopic mass error">@FRAGMENT_MASS_TOLERANCE@</note>
        <note type="input" label="spectrum, parent monoisotopic mass error plus">@PRECURSOR_MASS_TOLERANCE@</note>
        <note type="input" label="spectrum, parent monoisotopic mass error minus">@PRECURSOR_MASS_TOLERANCE@</note>
        <note type="input" label="spectrum, parent monoisotopic mass isotope error">yes</note>
        <note type="input" label="spectrum, fragment monoisotopic mass error units">Daltons</note>
        <note>The value for this parameter may be 'Daltons' or 'ppm': all other values are ignored</note>
        <note type="input" label="spectrum, parent monoisotopic mass error units">ppm</note>
                <note>The value for this parameter may be 'Daltons' or 'ppm': all other values are ignored</note>
        <note type="input" label="spectrum, fragment mass type">monoisotopic</note>
                <note>values are monoisotopic|average </note>

<note>spectrum conditioning parameters</note>
        <note type="input" label="spectrum, dynamic range">100.0</note>
                <note>The peaks read in are normalized so that the most intense peak
                is set to the dynamic range value. All peaks with values of less that
                1, using this normalization, are not used. This normalization has the
                overall effect of setting a threshold value for peak intensities.</note>
        <note type="input" label="spectrum, total peaks">50</note> 
                <note>If this value is 0, it is ignored. If it is greater than zero (lets say 50),
                then the number of peaks in the spectrum with be limited to the 50 most intense
                peaks in the spectrum. X! tandem does not do any peak finding: it only
                limits the peaks used by this parameter, and the dynamic range parameter.</note>
        <note type="input" label="spectrum, maximum parent charge">4</note>
        <note type="input" label="spectrum, use noise suppression">yes</note>
        <note type="input" label="spectrum, minimum parent m+h">500.0</note>
        <note type="input" label="spectrum, minimum fragment mz">150.0</note>
        <note type="input" label="spectrum, minimum peaks">15</note> 
        <note type="input" label="spectrum, threads">40</note>
        <note type="input" label="spectrum, sequence batch size">1000</note>

<note>residue modification parameters</note>
        <note type="input" label="residue, modification mass">57.022@C</note>
                <note>The format of this parameter is m@X, where m is the modfication
                mass in Daltons and X is the appropriate residue to modify. Lists of
                modifications are separated by commas. For example, to modify M and C
                with the addition of 16.0 Daltons, the parameter line would be
                +16.0@M,+16.0@C
                Positive and negative values are allowed.
                </note>
        <note type="input" label="residue, potential modification mass">16@M</note>
                <note>The format of this parameter is the same as the format
                for residue, modification mass (see above).</note>
        <note type="input" label="residue, potential modification motif"></note>
                <note>The format of this parameter is similar to residue, modification mass,
                with the addition of a modified PROSITE notation sequence motif specification.
                For example, a value of 80@[ST!]PX[KR] indicates a modification
                of either S or T when followed by P, and residue and the a K or an R.
                A value of 204@N!{P}[ST]{P} indicates a modification of N by 204, if it
                is NOT followed by a P, then either an S or a T, NOT followed by a P.
                Positive and negative values are allowed.
                </note>

<note>protein parameters</note>
        <note type="input" label="protein, taxon">other mammals</note>
                <note>This value is interpreted using the information in taxonomy.xml.</note>
        <note type="input" label="protein, cleavage site">[RK]|{P}</note>
                <note>this setting corresponds to the enzyme trypsin. The first characters
                in brackets represent residues N-terminal to the bond - the '|' pipe -
                and the second set of characters represent residues C-terminal to the
                bond. The characters must be in square brackets (denoting that only
                these residues are allowed for a cleavage) or french brackets (denoting
                that these residues cannot be in that position). Use UPPERCASE characters.
                To denote cleavage at any residue, use [X]|[X] and reset the 
                scoring, maximum missed cleavage site parameter (see below) to something like 50.
                </note>
        <note type="input" label="protein, modified residue mass file"></note>
        <note type="input" label="protein, cleavage C-terminal mass change">+17.002735</note>
        <note type="input" label="protein, cleavage N-terminal mass change">+1.007825</note>
        <note type="input" label="protein, N-terminal residue modification mass">0.0</note>
        <note type="input" label="protein, C-terminal residue modification mass">0.0</note>
        <note type="input" label="protein, homolog management">no</note>
        <note>if yes, an upper limit is set on the number of homologues kept for a particular spectrum</note>
        <note type="input" label="protein, quick acetyl">no</note>
        <note type="input" label="protein, quick pyrolidone">no</note>

<note>model refinement parameters</note>
        <note type="input" label="refine">yes</note>
        <note type="input" label="refine, modification mass"></note>
        <note type="input" label="refine, sequence path"></note>
        <note type="input" label="refine, tic percent">20</note>
        <note type="input" label="refine, spectrum synthesis">yes</note>
        <note type="input" label="refine, maximum valid expectation value">0.1</note>
        <note type="input" label="refine, potential N-terminus modifications">+42.010565@[</note>


        <note type="input" label="refine, potential C-terminus modifications"></note>
        <note type="input" label="refine, unanticipated cleavage">yes</note>
        <note type="input" label="refine, potential modification mass"></note>
        <note type="input" label="refine, point mutations">no</note>
        <note type="input" label="refine, use potential modifications for full refinement">no</note>
        <note type="input" label="refine, point mutations">no</note>
        <note type="input" label="refine, potential modification motif"></note>
        <note>The format of this parameter is similar to residue, modification mass,
                with the addition of a modified PROSITE notation sequence motif specification.
                For example, a value of 80@[ST!]PX[KR] indicates a modification
                of either S or T when followed by P, and residue and the a K or an R.
                A value of 204@N!{P}[ST]{P} indicates a modification of N by 204, if it
                is NOT followed by a P, then either an S or a T, NOT followed by a P.
                Positive and negative values are allowed.
                </note>

<note>scoring parameters</note>
        <note type="input" label="scoring, minimum ion count">4</note>
        <note type="input" label="scoring, maximum missed cleavage sites">@MAX_MISSED_CLEAVAGES@</note>
        <note type="input" label="scoring, x ions">no</note>
        <note type="input" label="scoring, y ions">yes</note>
        <note type="input" label="scoring, z ions">no</note>
        <note type="input" label="scoring, a ions">no</note>
        <note type="input" label="scoring, b ions">yes</note>
        <note type="input" label="scoring, c ions">no</note>
        <note type="input" label="scoring, cyclic permutation">no</note>
                <note>if yes, cyclic peptide sequence permutation is used to pad the scoring histograms</note>
        <note type="input" label="scoring, include reverse">no</note>
                <note>if yes, then reversed sequences are searched at the same time as forward sequences</note>
        <note type="input" label="scoring, cyclic permutation">no</note>
        <note type="input" label="scoring, include reverse">no</note>

<note>output parameters</note>
        <note type="input" label="output, log path"></note>
        <note type="input" label="output, message">testing 1 2 3</note>
        <note type="input" label="output, one sequence copy">no</note>
        <note type="input" label="output, sequence path"></note>
        <note type="input" label="output, path">output.xml</note>
        <note type="input" label="output, sort results by">protein</note>
                <note>values = protein|spectrum (spectrum is the default)</note>
        <note type="input" label="output, path hashing">no</note>
                <note>values = yes|no</note>
        <note type="input" label="output, xsl path">tandem-style.xsl</note>
        <note type="input" label="output, parameters">yes</note>
                <note>values = yes|no</note>
        <note type="input" label="output, performance">yes</note>
                <note>values = yes|no</note>
        <note type="input" label="output, spectra">yes</note>
                <note>values = yes|no</note>
        <note type="input" label="output, histograms">yes</note>
                <note>values = yes|no</note>
        <note type="input" label="output, proteins">yes</note>
                <note>values = yes|no</note>
        <note type="input" label="output, sequences">yes</note>
                <note>values = yes|no</note>
        <note type="input" label="output, one sequence copy">no</note>
                <note>values = yes|no, set to yes to produce only one copy of each protein sequence in the output xml</note>
        <note type="input" label="output, results">valid</note>
                <note>values = all|valid|stochastic</note>
        <note type="input" label="output, maximum valid expectation value">0.1</note>
                <note>value is used in the valid|stochastic setting of output, results</note>
        <note type="input" label="output, histogram column width">30</note>
                <note>values any integer greater than 0. Setting this to '1' makes cutting and pasting histograms
                into spread sheet programs easier.</note>
<note type="description">ADDITIONAL EXPLANATIONS</note>
        <note type="description">Each one of the parameters for X! tandem is entered as a labeled note
                        node. In the current version of X!, keep those note nodes
                        on a single line.
        </note>
        <note type="description">The presence of the type 'input' is necessary if a note is to be considered
                        an input parameter.
        </note>
        <note type="description">Any of the parameters that are paths to files may require alteration for a 
                        particular installation. Full path names usually cause the least trouble,
                        but there is no reason not to use relative path names, if that is the
                        most convenient.
        </note>
        <note type="description">Any parameter values set in the 'list path, default parameters' file are
                        reset by entries in the normal input file, if they are present. Otherwise,
                        the default set is used.
        </note>
        <note type="description">The 'list path, taxonomy information' file must exist.
                </note>
        <note type="description">The directory containing the 'output, path' file must exist: it will not be created.
                </note>
        <note type="description">The 'output, xsl path' is optional: it is only of use if a good XSLT style sheet exists.
                </note> 
</bioml>

<?xml version="1.0"?>
<bioml>
        <note>
        Each one of the parameters for x! tandem is entered as a labeled note node. 
        Any of the entries in the default_input.xml file can be over-ridden by
        adding a corresponding entry to this file. This file represents a minimum
        input file, with only entries for the default settings, the output file
        and the input spectra file name. 
        See the taxonomy.xml file for a description of how FASTA sequence list 
        files are linked to a taxon name.
        </note>

        <note type="input" label="list path, default parameters">%s</note>
        <note type="input" label="list path, taxonomy information">%s</note>

        <note type="input" label="protein, taxon">DB</note>

        <note type="input" label="spectrum, path">%s</note>

        <note type="input" label="output, path">%s</note>
</bioml>

process MakeXtandemConfig {
    input:
    file template from Channel.fromPath(params.xtandem_template)
    file fastadb_with_decoy from joined_fasta_with_decoys_ch.first()
    val max_missed_cleavages
    output:
    file "xtandem_config.xml" into xtandem_config_ch
    """
    sed 's/@DDA_DB_FILE@/$fastadb_with_decoy/g;s/@FRAGMENT_MASS_TOLERANCE@/$params.fragment_mass_tolerance/g;s/@PRECURSOR_MASS_TOLERANCE@/$params.precursor_mass_tolerance/g;s/@MAX_MISSED_CLEAVAGES@/$max_missed_cleavages/g' $template > xtandem_config.xml
    """
}


process XTandem {
    when:
    params.search_engines.contains("xtandem")

    input:
    file mzxml from maybespectra_ch.xtandem
    file tandem_config from xtandem_config_ch.first()
    file fastadb_with_decoy from joined_fasta_with_decoys_ch.first()
    output:
    file("${mzxml.baseName}.tandem.pep.xml") into xtandem_pepxml_ch
    file mzxml into xinteract_xtandem_mzxml_ch
    """
    printf '
    <<taxonomy-template>>'  $fastadb_with_decoy | tail -n+2 > xtandem_taxonomy.xml

    printf '
    <<xtandem-input-template>>' $tandem_config xtandem_taxonomy.xml $mzxml ${mzxml.baseName}.TANDEM.OUTPUT.xml | tail -n+2 > input.xml
    tandem input.xml
    Tandem2XML ${mzxml.baseName}.TANDEM.OUTPUT.xml ${mzxml.baseName}.tandem.pep.xml 
    """
}

process XinteractXTandem {
    memory '16 GB'
    input:
    file pepxmls from xtandem_pepxml_ch.toSortedList()
    // the filename of needed fastdadb was defined in cometcfg
    // and stored in pepxml in the comet-ms step
    // -a suppplies the absulute path to the data directory where the mzxmls
    // rather than reading wherer the mfrom the xmls
    // where the mzxml are, because its not very
    // nextflow to look outside the cwd.
    file fastadb_with_decoy from joined_fasta_with_decoys_ch.first()
    file mzxmls from  xinteract_xtandem_mzxml_ch.toSortedList()
    output: 
    file "interact_xtandem.pep.xml" into xtandem_search_results_ch
    when:
    pepxmls.size() > 0 
    """
    xinteract -a\$PWD -OARPd -dDECOY_ -Ninteract_xtandem.pep.xml $pepxmls
    """
}

params.xtandem_template="xtandem-template.xml"

// we handle the one or two engines case
// DSL2 incompat
// would be in workflow body

if (params.search_engines.size() > 1) {  
    process CombineSearchResults {
        publishDir "${params.outdir}/speclib"
        when:

        input:
        file xtandem_search_results from xtandem_search_results_ch
        file comet_search_results from comet_search_results_ch
        output:
        file "lib_iprophet.peps.xml" into combined_search_results_ch
        """
        InterProphetParser DECOY=DECOY_ THREADS=${task.cpus} $xtandem_search_results $comet_search_results lib_iprophet.peps.xml
        """
    }
} else if (params.search_engines.contains("comet")) {
    combined_search_results_ch = comet_search_results_ch
} else if (params.search_engines.contains("xtandem")) {
    combined_search_results_ch =xtandem_search_results_ch
} else {
    combined_search_results_ch = Channel.create()
}

if (params.pyprophet_use_legacy)
    process pyprophet_legacy {
    publishDir "${params.outdir}/pyprophet/", pattern: "*.csv"
    publishDir "${params.outdir}/reports/pyprophet/", pattern: "*.pdf"
    input:
    file dia_tsv from openswath_tsv_ch
    output:
    file dscore_csv into pyprophet_legacy_ch
    // just for publishing
    file "${dia_tsv.baseName}_report.pdf" 
    script:
    seed="928418756933579397"

    dscore_csv="${dia_tsv.baseName}_with_dscore.csv"
    """
    pyprophet --random_seed=${seed} --delim=tab --export.mayu ${dia_tsv} --ignore.invalid_score_columns
    """
}

// The ratio (0,1] to subsample by in pyprophet.
// leave to null to use 1 / nr_samples
params.pyprophet_subsample_ratio = null 

if(params.pyprophet_subsample_ratio == null){
    Channel.value(1 / Math.max(Channel.fromPath(params.diafiles).toSortedList().size().getVal(), 1))
        .set{subsample_ratio}
} else {
    Channel.value(params.pyprophet_subsample_ratio).set{subsample_ratio}
}

params.pyprophet_fixed_seed=true
params.pyprophet_use_legacy=false

if (!params.pyprophet_use_legacy)
    {

process pyprophet_subsample {
    input:
    file dia_osw_file from openswath_osw_ch.pyprophet_subsample
    val subsample_ratio
    output:
    file subsampled_osw
    script:
    subsampled_osw="${dia_osw_file.baseName}.osws"
    pyprophet_seed_flag=(params.pyprophet_fixed_seed ? "--test" : "--no-test")
    """
    pyprophet subsample $pyprophet_seed_flag --in=$dia_osw_file --out=$subsampled_osw --subsample_ratio=$subsample_ratio
    """
}

process pyprophet_learn_classifier {
    input:
    file subsampled_osws from subsampled_osw.toSortedList()
    file openswath_transitions from openswath_transitions_ch_for_pyprophet.first()
    output:
    file osw_model
    script:
    pyprophet_seed_flag=(params.pyprophet_fixed_seed ? "--test" : "--no-test")
    osw_model="model.osw"
    """
    pyprophet merge --template=$openswath_transitions --out=$osw_model $subsampled_osws
    pyprophet score  $pyprophet_seed_flag --in=$osw_model --level=ms1ms2
    """
}

scored_osw_indirect_ch =Channel.create()
scored_osw_indirect_ch.multiMap{it ->
    reduce: backpropagate: it}.set{scored_osw_ch}

process pyprophet_apply_classifier {
    input:
    file osw_model from osw_model.first()
    file osw from openswath_osw_ch.pyprophet_score
    output:
    file scored_osw into scored_osw_indirect_ch
    script:
    pyprophet_seed_flag=(params.pyprophet_fixed_seed ? "--test" : "--no-test")
    scored_osw="${osw.baseName}.scored.${osw.Extension}"
    """
    pyprophet score $pyprophet_seed_flag --in=$osw --out=$scored_osw --apply_weights=$osw_model --level=ms1ms2
    """
}

process pyprophet_reduce {
    input:
    file scored_osw from scored_osw_ch.reduce
    output:
    file reduced_scored_osw 
    script:
    reduced_scored_osw="${file(scored_osw.baseName).baseName}.${scored_osw.Extension}r"
    """
    pyprophet reduce --in=$scored_osw --out=$reduced_scored_osw
    """
}

process pyprophet_control_error {
    input:
    file reduced_scored_osws from reduced_scored_osw.toSortedList()
    file osw_model from osw_model.first()
    output:
    file osw_global_model
    script:
    osw_global_model="model_global.osw"
    """
    pyprophet merge --template=$osw_model --out=$osw_global_model $reduced_scored_osws
    pyprophet peptide --context=global --in=$osw_global_model
    pyprophet protein --context=global --in=$osw_global_model
    """
}

process pyprophet_backpropagate {
    input:
    file osw_scored from scored_osw_ch.backpropagate
    file osw_global_model from osw_global_model.first()
    output:
    file dscore_tsv into pyprophet_nonlegacy_ch
    script:
    base="${file(osw_scored.baseName).baseName}"
    backproposw="${base}.backprop.osw"
    dscore_tsv="${base}.tsv"
    """
    pyprophet backpropagate --in="$osw_scored" --apply_scores="$osw_global_model" --out=$backproposw
    # we supply --format=legacy_merged so that pyprophet export respect the --out parameter
    pyprophet export --in=$backproposw --format=legacy_merged --out=$dscore_tsv
    """
}

}

if (params.pyprophet_use_legacy)
    pyprophet_legacy_ch.set{pyprophet_ch}
else
    pyprophet_nonlegacy_ch.set{pyprophet_ch}

("swath2stats"
 "r-minimal"
 "r-swath2stats")

#!/usr/bin/env Rscript

.libPaths("/opt/Rlibs/")

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(version = "3.15",ask=FALSE,update=FALSE)
BiocManager::install("SWATH2stats", ask=FALSE,update=FALSE)