diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
index 663f575..0bc4b52 100644
--- a/.github/workflows/release.yml
+++ b/.github/workflows/release.yml
@@ -10,13 +10,21 @@ jobs:
   create-container:
     uses: ./.github/workflows/container.yml
 
-  run-all-tests:
+  run-all-tests-sage:
     needs: create-container
+    env:
+      SEARCH: sage
+    uses: ./.github/workflows/tests_all.yml
+
+  run-all-tests-msgf:
+    needs: create-container
+    env:
+      SEARCH: msgf
     uses: ./.github/workflows/tests_all.yml
 
   create-release:  
     runs-on: ubuntu-latest
-    needs: run-all-tests
+    needs: [run-all-tests-sage, run-all-tests-msgf]
     steps:
       - uses: ncipollo/release-action@v1
         with:
diff --git a/.github/workflows/single_test.yml b/.github/workflows/single_test.yml
index 8f9f102..44ce0c2 100644
--- a/.github/workflows/single_test.yml
+++ b/.github/workflows/single_test.yml
@@ -20,4 +20,4 @@ jobs:
             uses: nf-core/setup-nextflow@v2
 
           - name: Run single test
-            run: bash run_tests.sh tmt16_fast
+            run: bash run_tests.sh sage tmt16_fast
diff --git a/.github/workflows/tests_all.yml b/.github/workflows/tests_all.yml
index 8c9b75e..c08cb3e 100644
--- a/.github/workflows/tests_all.yml
+++ b/.github/workflows/tests_all.yml
@@ -15,7 +15,7 @@ jobs:
             uses: nf-core/setup-nextflow@v2
 
           - name: Run tests
-            run: bash run_tests.sh tmt16
+            run: bash run_tests.sh $SEARCH tmt16
 
     test-tmt18:
         runs-on: ubuntu-latest
@@ -27,7 +27,7 @@ jobs:
             uses: nf-core/setup-nextflow@v2
 
           - name: Run tests
-            run: bash run_tests.sh tmt18
+            run: bash run_tests.sh $SEARCH tmt18
 
     test-labelfree:
         runs-on: ubuntu-latest
@@ -39,7 +39,7 @@ jobs:
             uses: nf-core/setup-nextflow@v2
 
           - name: Run tests
-            run: bash run_tests.sh labelfree
+            run: bash run_tests.sh $SEARCH labelfree
 
     test-lf_phos:
         runs-on: ubuntu-latest
@@ -51,7 +51,7 @@ jobs:
             uses: nf-core/setup-nextflow@v2
 
           - name: Run tests
-            run: bash run_tests.sh labelfree_phos
+            run: bash run_tests.sh $SEARCH labelfree_phos
 
     test-tims:
         runs-on: ubuntu-latest
@@ -63,7 +63,7 @@ jobs:
             uses: nf-core/setup-nextflow@v2
 
           - name: Run tests
-            run: bash run_tests.sh tims
+            run: bash run_tests.sh $SEARCH tims
 
     test-mixed-tmt:
         runs-on: ubuntu-latest
@@ -75,4 +75,4 @@ jobs:
             uses: nf-core/setup-nextflow@v2
 
           - name: Run tests
-            run: bash run_tests.sh tmt16_18
+            run: bash run_tests.sh $SEARCH tmt16_18
diff --git a/CHANGELOG.md b/CHANGELOG.md
index c2c5e1b..bfa3e4f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,16 @@
 # lehtiolab/ddamsproteomics: Changelog
 
+## Version 3.1 [2025-02-03]
+- Search engine choice of Sage / MSGF
+- QC plots changes in colors, max PSM use in isobaric step, PCA for sets even when sample groups
+- Tables for PSM results in report
+- Possible to remove (e.g. empty) channel for specific set
+- Fixed volcano plots
+- Error early on incorrect totalproteome / remove channel set names
+- Output ExplainedIonCurrentRatio in MSGF PSM table (msstitch to 3.17)
+- OpenMS to 3.2, MSGFPlus to 2024.03.26, Percolator to 3.6.5
+
+
 ## Version 3.0 [2024-09-19]
 - Rewrote in DSL2
 - QC report with plotly
diff --git a/README.md b/README.md
index 8df41f4..5bab1dc 100644
--- a/README.md
+++ b/README.md
@@ -18,16 +18,16 @@ The pipeline is built using [Nextflow](https://www.nextflow.io), a workflow tool
 - run pipeline:
 
 ```
-nextflow run lehtiolab/ddamsproteomics --input /path/to/input_definition.txt --tdb /path/to/proteins.fa --mods 'oxidation;carbamidomethylation' -profile standard,docker
+nextflow run lehtiolab/ddamsproteomics --input /path/to/input_definition.txt --sage --tdb /path/to/proteins.fa --mods 'oxidation;carbamidomethylation' -profile standard,docker
 ```
 
 Or for two sample sets of isobaric data you can:
 
 ```
-nextflow run lehtiolab/ddamsproteomics --input /path/to/input_definition.txt --tdb /path/to/proteins.fa --mods 'oxidation;carbamidomethylation --isobaric 'setA:tmt10plex:126 setB:tmt10plex:127N'
+nextflow run lehtiolab/ddamsproteomics --input /path/to/input_definition.txt --sage --tdb /path/to/proteins.fa --mods 'oxidation;carbamidomethylation --isobaric 'setA:tmt10plex:126 setB:tmt10plex:127N'
 ```
 
-For more elaborate examples covering fractionation, PTMs, and more, the lehtiolab/ddamsproteomics pipeline comes with documentation about the pipeline, found in the `docs/` directory:
+For more elaborate examples covering fractionation, MSGF+, PTMs, and more, the lehtiolab/ddamsproteomics pipeline comes with documentation about the pipeline, found in the `docs/` directory:
 
 - [Running the pipeline](docs/usage.md)
 - [Output and how to interpret the results](docs/output.md)
diff --git a/assets/report.html b/assets/report.html
index f4c5de7..86cf1ea 100644
--- a/assets/report.html
+++ b/assets/report.html
@@ -28,9 +28,12 @@ <h4 class="subtitle is-4">
           {% if ptmtables.summary %}
           <li id="ptmqctab"> <a> PTMs </a> </li>
           {% endif %}
-          {% if expplots.pca %}
+          {% if expplots.pcagroup %}
           <li id="deqctab"><a>Expression</a> </li>
           {% endif %}
+          {% if expplots.pcaset %}
+          <li id="pcatab"><a>PCA</a> </li>
+          {% endif %}
           <li id="rundatatab"><a>Run data</a></li>
         </ul>
       </div>
@@ -95,6 +98,52 @@ <h5 class="title is-5">{{ ftitle }}</h5>
         {% endfor %}
       </div>
 
+      <h4 class="title is-4">PSMs</h4>
+      <div class="columns">
+        <div class="column">
+          <h5 class="title is-5">IDs</h5>
+          <table class="table is-striped is-narrow is-hoverable is-size-7">
+            <thead>
+              <th>Sample/plate</th>
+              <th>Scans</th>
+              <th>PSMs</th>
+              <th>% ID</th>
+            </thead>
+            <tbody>
+              {% for plate, scans, psms, pc in psmtables.ids %}
+              <tr>
+                <td> {{ plate }} </td>
+                <td> {{ scans }} </td>
+                <td> {{ psms }} </td>
+                <td> {{ pc }} </td>
+              </tr>
+              {% endfor %}
+            </tbody>
+          </table>
+        </div>
+        <div class="column">
+          <h5 class="title is-5">Missed cleavages</h5>
+          <table class="table is-striped is-narrow is-hoverable is-size-7">
+            <thead>
+              <th>Sample/plate</th>
+              <th># missed cleavages</th>
+              <th># PSMs</th>
+              <th>% PSMs</th>
+            </thead>
+            <tbody>
+              {% for plate, mc, psms, pc in psmtables.miscleav %}
+              <tr>
+                <td> {{ plate }} </td>
+                <td> {{ mc }} </td>
+                <td> {{ psms }} </td>
+                <td> {{ pc }} </td>
+              </tr>
+              {% endfor %}
+            </tbody>
+          </table>
+        </div>
+      </div>
+
       {% if ptmtables.summary %}
       <h4 class="title is-4">PTMs</h4>
       <div class="columns">
@@ -360,7 +409,7 @@ <h5 class="title is-5">{{ ftitle }}</h5>
     </div>
     {% endif %}
 
-    {% if expplots.pca %}
+    {% if expplots.pcagroup %}
     <div id="deqcpage" class="mt-5 container">
       {% for plotname, plottitle in expplotnames %}
       <h4 class="title is-4">{{ plottitle }}</h4>
@@ -397,6 +446,24 @@ <h5 class="title is-5">{{ ftitle }}</h5>
     </div>
     {% endif %}
 
+    {% if expplots.pcaset %}
+    <div id="pcapage" class="mt-5 container">
+      {% for plotname, plottitle in pcaplotnames %}
+      <h4 class="title is-4">{{ plottitle }}</h4>
+      <div class="columns">
+        {% for fname, ftitle in featnames %}
+        {% if expplots[plotname][fname] %}
+        <div class="column">
+          <h5 class="title is-5">{{ ftitle }}</h5>
+          {{ expplots[plotname][fname] }}
+        </div>
+        {% endif %}
+        {% endfor %}
+      </div>
+      {% endfor %}
+    </div>
+    {% endif %}
+
     <script type="text/javascript">
       let tabs = [
         ['results', document.getElementById('resultstab'), document.getElementById('resultspage')],
@@ -406,9 +473,12 @@ <h5 class="title is-5">{{ ftitle }}</h5>
         {% if ptmtables.summary %}
         ['ptmqc', document.getElementById('ptmqctab'), document.getElementById('ptmqcpage')],
         {% endif %}
-        {% if expplots.pca %}
+        {% if expplots.pcagroup %}
         ['deqc', document.getElementById('deqctab'), document.getElementById('deqcpage')],
         {% endif %}
+        {% if expplots.pcaset %}
+        ['pca', document.getElementById('pcatab'), document.getElementById('pcapage')],
+        {% endif %}
       ];
       tabs.forEach(t => {
         t[1].addEventListener('click', e => toggleTab(t[0]));
diff --git a/assets/sage.json b/assets/sage.json
new file mode 100644
index 0000000..40ef5d9
--- /dev/null
+++ b/assets/sage.json
@@ -0,0 +1,46 @@
+{
+  "database": {
+    "bucket_size": 8192,
+    "enzyme": {
+      "missed_cleavages": 2,
+      "min_len": 7,
+      "max_len": 50,
+      "cleave_at": "KR",
+      "restrict": "P",
+      "c_terminal": true
+    },
+    "fragment_min_mz": 200.0,
+    "fragment_max_mz": 2000.0,
+    "peptide_min_mass": 500.0,
+    "peptide_max_mass": 5000.0,
+    "ion_kinds": ["a", "b", "y"],
+    "min_ion_index": 2,
+    "static_mods": {
+      "C": 57.0215
+    },
+    "variable_mods": {
+      "M": [15.9949]
+    },
+    "max_variable_mods": 2,
+    "decoy_tag": "decoy_",
+    "generate_decoys": false
+  },
+  "precursor_tol": {
+    "ppm": [-10, 10]
+  },
+  "fragment_tol": {
+    "ppm": [-20, 20]
+  },
+  "precursor_charge": [2, 6],
+  "isotope_errors": [ -1, 2 ],
+  "deisotope": false,
+  "chimera": false,
+  "wide_window": false,
+  "predict_rt": false,
+  "min_peaks": 7,
+  "max_peaks": 150,
+  "min_matched_peaks": 3,
+  "max_fragment_charge": 2,
+  "report_psms": 1,
+  "output_directory": "./"
+}
diff --git a/bin/create_modfile.py b/bin/create_modfile.py
index 90c89fb..c7f810b 100755
--- a/bin/create_modfile.py
+++ b/bin/create_modfile.py
@@ -1,22 +1,29 @@
 #!/usr/bin/env python3
 
 import sys
+import json
 from mods import Mods
 
 
 def main():
     nummods = sys.argv[1]
     modlibfile = sys.argv[2]
-    passed_mods = sys.argv[3].split(';')
+    search_engine = sys.argv[3]
+    passed_mods = sys.argv[4].split(';')
     fixedmods = {}
     varmods = []
     # Parse modifications passed
-    msgfmods = Mods()
-    msgfmods.parse_msgf_modfile(modlibfile, passed_mods)
-    with open('mods.txt', 'w') as fp:
-        fp.write('NumMods={}'.format(nummods))
-        for modline in msgfmods.get_msgf_modlines():
-            fp.write(f'\n{modline}')
+    parsemods = Mods()
+    parsemods.parse_msgf_modfile(modlibfile, passed_mods)
+    modfn = 'mods.txt'
+    if search_engine == 'msgf':
+        with open(modfn, 'w') as fp:
+            fp.write('NumMods={}'.format(nummods))
+            for modline in parsemods.get_msgf_modlines():
+                fp.write(f'\n{modline}')
+    elif search_engine == 'sage':
+        with open(modfn, 'w') as fp:
+            json.dump(parsemods.get_sage_mods(), fp)
 
 
 if __name__ == '__main__':
diff --git a/bin/luciphor_parse.py b/bin/luciphor_parse.py
index 109dd9c..77dbd85 100755
--- a/bin/luciphor_parse.py
+++ b/bin/luciphor_parse.py
@@ -10,15 +10,19 @@
 
 
 # PTM input/output header fields
-TOPPTM = 'Top luciphor PTM'
-TOPFLR = 'Top PTM FLR'
-TOPSCORE = 'Top PTM score'
-OTHERPTMS = 'High-scoring PTMs'
-SE_PEPTIDE = 'SearchEnginePeptide'
-PEPTIDE = 'Peptide'
-MASTER_PROTEIN = 'Master protein(s)'
-FLANKING_SEQS = 'PTM flanking seq(s)'
-PTMFIELDS = [SE_PEPTIDE, TOPPTM, TOPSCORE, TOPFLR, OTHERPTMS, FLANKING_SEQS]
+def get_headerfields(search_engine):
+    fields = {'TOPPTM': 'Top luciphor PTM', 'TOPFLR': 'Top PTM FLR', 'TOPSCORE': 'Top PTM score',
+            'OTHERPTMS': 'High-scoring PTMs', 'SE_PEPTIDE': 'SearchEnginePeptide',
+            'MASTER_PROTEIN': 'Master protein(s)', 'FLANKING_SEQS': 'PTM flanking seq(s)',
+
+            'PEPTIDE': {'sage': 'peptide', 'msgf': 'Peptide'}[search_engine],
+            'CHARGE': {'sage': 'charge', 'msgf': 'Charge'}[search_engine],
+            'FILENAME': {'sage': 'filename', 'msgf': '#SpecFile'}[search_engine],
+            'SCANNR': {'sage': 'scannr', 'msgf': 'ScanNum'}[search_engine],
+            }
+    fields['PTMFIELDS'] = [fields['SE_PEPTIDE'], fields['TOPPTM'], fields['TOPSCORE'], fields['TOPFLR'], fields['OTHERPTMS'], fields['FLANKING_SEQS']]
+    return fields
+
 
 def main():
     parser = argparse.ArgumentParser(formatter_class=argparse.RawTextHelpFormatter)
@@ -32,18 +36,20 @@ def main():
     parser.add_argument('--stabileptms', nargs='+', default=[])
     parser.add_argument('--mods', nargs='+', default=[])
     parser.add_argument('--fasta')
+    parser.add_argument('--searchengine', default='msgf', type=str)
     args = parser.parse_args(sys.argv[1:])
 
     minscore_high = args.minscore
     labileptms = [x.lower() for x in args.labileptms]
     stabileptms = [x.lower() for x in args.stabileptms]
     mods = [x.lower() for x in args.mods]
+    fields = get_headerfields(args.searchengine)
 
     # First prepare a residue + PTM weight -> PTM name dict for naming mods
     msgfmods = Mods()
     msgfmods.parse_msgf_modfile(args.modfile, [*args.labileptms, *args.stabileptms, *args.mods])
     luci_modmap = msgfmods.lucimass_mod_dict()
-    msgf_mod_map = msgfmods.msgfmass_mod_dict()
+    msgf_mod_map = msgfmods.msgfmass_mod_dict(args.searchengine)
 
     # load sequences
     tdb = SeqIO.index(args.fasta, 'fasta')
@@ -81,30 +87,34 @@ def main():
     # because that needs parsing from PSM table (luciphor throws them out)
     with open(args.psms) as psms, open(args.outfile, 'w') as wfp:
         psmheader = next(psms).strip('\n').split('\t')
-        outheader = psmheader + PTMFIELDS
+        outheader = psmheader + fields['PTMFIELDS']
         wfp.write('\t'.join(outheader))
         for psm in psms:
             psm = psm.strip('\n').split('\t')
             psm = {k: v for k,v in zip(psmheader, psm)}
             msgfpsm = PSM()
-            msgfpsm.parse_msgf_peptide(psm[PEPTIDE], msgf_mod_map,
-                    labileptms, stabileptms)
+            msgfpsm.parse_search_peptide(psm[fields['PEPTIDE']], msgf_mod_map,
+                    labileptms, stabileptms, args.searchengine)
             if not msgfpsm.has_labileptms():
                 continue
-            psmid = '{}.{}.{}.{}'.format(os.path.splitext(psm['#SpecFile'])[0], psm['ScanNum'], psm['ScanNum'], psm['Charge'])
+            if args.searchengine == 'sage':
+                scan = re.sub('.*scan=', '', psm[fields['SCANNR']])
+            else:
+                scan = psm[fields['SCANNR']]
+            psmid = '{}.{}.{}.{}'.format(os.path.splitext(psm[fields['FILENAME']])[0], scan, scan, psm[fields['CHARGE']])
             if psmid in lucptms:
                 luciptm = lucptms[psmid]
                 ptm = {
-                    TOPFLR: luciptm.top_flr,
-                    TOPSCORE: luciptm.top_score,
-                    OTHERPTMS: luciptm.format_alt_ptm_locs(),
+                    fields['TOPFLR']: luciptm.top_flr,
+                    fields['TOPSCORE']: luciptm.top_score,
+                    fields['OTHERPTMS']: luciptm.format_alt_ptm_locs(),
                     }
             else:
                 luciptm = False
                 ptm = {
-                    TOPFLR: 'NA',
-                    TOPSCORE: 'NA',
-                    OTHERPTMS: 'NA',
+                    fields['TOPFLR']: 'NA',
+                    fields['TOPSCORE']: 'NA',
+                    fields['OTHERPTMS']: 'NA',
                     }
             # If there are PTMs which are on fixed mod residues, we need to parse
             # those from the residues from sequences in the PSM table because luciphor
@@ -114,22 +124,22 @@ def main():
                 luciptm = msgfpsm
             elif msgfmods.has_varmods_on_fixmod_residues:
                 luciptm.add_ptms_from_psm(msgfpsm.mods)
-            ptm[TOPPTM] = luciptm.topptm_output()
+            ptm[fields['TOPPTM']] = luciptm.topptm_output()
 
             # Get protein site location of mods
-            if MASTER_PROTEIN in psm:
-                annotate_protein_and_flanks(psm, luciptm, tdb, [*labileptms, *stabileptms])
+            if fields['MASTER_PROTEIN'] in psm:
+                annotate_protein_and_flanks(psm, luciptm, tdb, [*labileptms, *stabileptms], fields)
             outpsm = {k: v for k,v in psm.items()}
             outpsm.update(ptm)
-            outpsm[SE_PEPTIDE] = outpsm.pop(PEPTIDE)
-            outpsm[PEPTIDE] = '{}_{}'.format(luciptm.sequence, ptm[TOPPTM])
+            outpsm[fields['SE_PEPTIDE']] = outpsm.pop(fields['PEPTIDE'])
+            outpsm[fields['PEPTIDE']] = '{}_{}'.format(luciptm.sequence, ptm[fields['TOPPTM']])
             wfp.write('\n{}'.format('\t'.join([outpsm[k] for k in outheader])))
 
 
-def annotate_protein_and_flanks(psm, ptmpsm, tdb, ptmnames):
+def annotate_protein_and_flanks(psm, ptmpsm, tdb, ptmnames, fields):
     flankseqs = set()
     proteins_peploc = {}
-    for p in psm[MASTER_PROTEIN].split(';'):
+    for p in psm[fields['MASTER_PROTEIN']].split(';'):
         proteins_peploc[p] = [x.start() for x in re.finditer(ptmpsm.sequence, str(tdb[p].seq))]
     proteins_loc = {p: [] for p, peplocs in proteins_peploc.items() if len(peplocs)}
     for p, peplocs in proteins_peploc.items():
@@ -145,8 +155,8 @@ def annotate_protein_and_flanks(psm, ptmpsm, tdb, ptmnames):
             flankpos = [(max(x-7, 0) , min(x+8, len(protseq))) for x in site_protlocs]
             flankseqs.update([str(protseq[x[0]:x[1]]) for x in flankpos])
             proteins_loc[p].append('{}:{}'.format(ptm['name'], ','.join(protptms)))
-    psm[MASTER_PROTEIN] = ';'.join(['{}__{}'.format(p, '_'.join(ptmloc)) for p, ptmloc in proteins_loc.items()])
-    psm[FLANKING_SEQS] = ';'.join(flankseqs)
+    psm[fields['MASTER_PROTEIN']] = ';'.join(['{}__{}'.format(p, '_'.join(ptmloc)) for p, ptmloc in proteins_loc.items()])
+    psm[fields['FLANKING_SEQS']] = ';'.join(flankseqs)
 
 
 if __name__ == '__main__':
diff --git a/bin/luciphor_prep.py b/bin/luciphor_prep.py
index 0460d14..d189d7a 100755
--- a/bin/luciphor_prep.py
+++ b/bin/luciphor_prep.py
@@ -43,14 +43,18 @@ def get_mod_dict(self, residue, sitenum, modptm, labileptmnames, stableptmnames)
                 'mass': modptm['mass'], 'name': modptm['name'],
                 'name_lower': modptm['name_lower'], 'adjusted_mass': modptm['adjusted_mass']}
 
-    def parse_msgf_peptide(self, msgfseq, msgf_mods, labileptmnames, stableptmnames):
+    def parse_search_peptide(self, searchseq, msgf_mods, labileptmnames, stableptmnames, search_engine):
         self.mods = []
         barepep = ''
         start = 0
-        for x in re.finditer(r'([A-Z]){0,1}([0-9\.+\-]+)', msgfseq):
+        if search_engine == 'sage':
+            regex = r'([A-Z]){0,1}([\[\]0-9\.+\-]+)'
+        else:
+            regex = r'([A-Z]){0,1}([0-9\.+\-]+)'
+        for x in re.finditer(regex, searchseq):
             if x.group(1) is not None:
                 # mod is on a residue
-                barepep = f'{barepep}{msgfseq[start:x.start()+1]}'
+                barepep = f'{barepep}{searchseq[start:x.start()+1]}'
                 residue = barepep[-1]
                 sitenum = len(barepep) - 1
             else:
@@ -60,10 +64,19 @@ def parse_msgf_peptide(self, msgfseq, msgf_mods, labileptmnames, stableptmnames)
             # TODO cterm = 100, ']'
             start = x.end()
             for mass in re.findall(r'[\+\-][0-9.]+', x.group(2)):
-                mod = msgf_mods[float(mass)][0] # only take first, contains enough info
+                try:
+                    mod = msgf_mods[float(mass)][0] # only take first, contains enough info
+                except KeyError:
+                    if search_engine == 'sage':
+                        # Sage turns 57.021464 into 57.021465 (carbamidomethyl)
+                        # Not sure about others or why, but try to counteract here
+                        alt_sagemass = float(f'{float(mass):.4f}')
+                        mod = msgf_mods[alt_sagemass][0]
+                    else:
+                        raise
                 self.mods.append(self.get_mod_dict(residue, sitenum, mod, labileptmnames,
                     stableptmnames))
-        self.sequence = f'{barepep}{msgfseq[start:]}'
+        self.sequence = f'{barepep}{searchseq[start:]}'
 
     def parse_luciphor_peptide(self, luciline, ptms_map, labileptmnames, stableptmnames):
         '''From a luciphor sequence, create a peptide with PTMs
@@ -144,6 +157,8 @@ def main():
     parser.add_argument('--modfile')
     parser.add_argument('--labileptms', nargs='+', default=[])
     parser.add_argument('--mods', nargs='+', default=[])
+    parser.add_argument('--searchengine', default='msgf', type=str)
+
     args = parser.parse_args(sys.argv[1:])
 
     labileptms = [x.lower() for x in args.labileptms]
@@ -196,22 +211,41 @@ def main():
     # translation table needed...
     # But how to spec in luciphor, it also wants fixed/var/target mods? Does it apply fixed regardless?
 
-    msgf_mod_map = msgfmods.msgfmass_mod_dict()
+    msgf_mod_map = msgfmods.msgfmass_mod_dict(args.searchengine)
+    colummap = {'sage': {
+        'pep': 'peptide',
+        'fn': 'filename',
+        'ch': 'charge',
+        'scan': 'scannr',
+        'value': 'PSM q-value',
+        }, 'msgf': {
+        'pep': 'Peptide',
+        'fn': '#SpecFile',
+        'ch': 'Charge',
+        'scan': 'ScanNum',
+        'value': 'PSM q-value',
+            }
+        }
+            
     with open(args.psmfile) as fp, open(args.lucipsms, 'w') as wfp:
         header = next(fp).strip('\n').split('\t')
-        pepcol = header.index('Peptide')
-        spfile = header.index('#SpecFile')
-        charge = header.index('Charge')
-        scan = header.index('ScanNum')
-        evalue = header.index('PSM q-value')
+        pepcol = header.index(colummap[args.searchengine]['pep'])
+        spfile = header.index(colummap[args.searchengine]['fn'])
+        charge = header.index(colummap[args.searchengine]['ch'])
+        scancol = header.index(colummap[args.searchengine]['scan'])
+        evalue = header.index(colummap[args.searchengine]['value'])
         wfp.write('srcFile\tscanNum\tcharge\tPSMscore\tpeptide\tmodSites')
         for line in fp:
             line = line.strip('\n').split('\t')
             psm = PSM()
-            psm.parse_msgf_peptide(line[pepcol], msgf_mod_map, labileptms, othermods)
+            psm.parse_search_peptide(line[pepcol], msgf_mod_map, labileptms, othermods, args.searchengine)
+            if args.searchengine == 'sage':
+                scan = re.sub('.*scan=', '', line[scancol])
+            else:
+                scan = line[scancol]
             # TODO add C-terminal mods (rare)
             if psm.has_labileptms():
-                wfp.write('\n{}\t{}\t{}\t{}\t{}\t{}'.format(line[spfile], line[scan], line[charge], line[evalue], psm.sequence, psm.luciphor_input_sites()))
+                wfp.write('\n{}\t{}\t{}\t{}\t{}\t{}'.format(line[spfile], scan, line[charge], line[evalue], psm.sequence, psm.luciphor_input_sites()))
 
 
 if __name__ == '__main__':
diff --git a/bin/mods.py b/bin/mods.py
index da3a2cb..c21f17b 100644
--- a/bin/mods.py
+++ b/bin/mods.py
@@ -1,3 +1,5 @@
+import sys
+
 NON_BLOCKING_MODS = {
         'GG': ['TMTpro', 'TMT6plex', 'iTRAQ8plex', 'iTRAQ4plex'],
         }
@@ -111,7 +113,7 @@ def parse_msgf_modfile(self, modfile, mods_passed):
                     adjustment += fmod['mass']
             mod['adjusted_mass'] = round(-(adjustment - mod['mass']), 5)
 
-    def get_msgf_modlines(self):
+    def get_grouped_mods_by_mass(self):
         grouped = {}
         for mod in self.mods:
             mass = self.get_mass_or_adj(mod)
@@ -123,7 +125,10 @@ def get_msgf_modlines(self):
                 if len(check_names) != 1:
                     print('Cannot have two modifications of the same mass but different names')
                     sys.exit(1)
+        return grouped
 
+    def get_msgf_modlines(self):
+        grouped = self.get_grouped_mods_by_mass()
         for mass, mods in grouped.items():
             name = mods[0]['name']
             line_res = {}
@@ -139,6 +144,27 @@ def get_msgf_modlines(self):
                 vf = 'opt' if int(var) else 'fix'
                 yield f'{mass},{"".join(residues)},{vf},{pos},{name}'
 
+    def get_sage_mods(self):
+        grouped = self.get_grouped_mods_by_mass()
+        outmods = {'static_mods': {}, 'variable_mods': {}}
+        posmap = {'N-term': '^',
+                'C-term': '$',
+                }
+        mtypemap = {False: 'static_mods', True: 'variable_mods'}
+        for mass, mods in grouped.items():
+            line_res = {}
+            for mod in mods:
+                mtype = mtypemap[mod['var']]
+                pos = '' if mod['pos'] == 'any' else posmap[mod['pos']]
+                res = '' if mod['residue'] == '*' else mod['residue']
+                sage_residue = f'{pos}{res}'
+                try:
+                    outmods[mtype][sage_residue].append(mass)
+                except KeyError:
+                    outmods[mtype][sage_residue] = [mass]
+        outmods['static_mods'] = {k: sum(v) for k,v in outmods['static_mods'].items()}
+        return outmods
+
     def get_mass_or_adj(self, mod):
         return mod['adjusted_mass'] or mod['mass']
  
@@ -153,14 +179,18 @@ def get_luci_input_mod_line(self, mod):
             residue = mod['residue']
         return f'{residue} {self.get_mass_or_adj(mod)}'
 
-    def msgfmass_mod_dict(self):
-        '''Create MSGF output mass (round(x,3) ) to mod lookup'''
+    def msgfmass_mod_dict(self, search_engine):
+        '''Create MSGF output mass (round(x,3) ), or sage (round(x,5), round(x,4) ) to mod lookup'''
         mod_map = {}
+        nrdigits = {'sage': [4,5], 'msgf': [3]}[search_engine]
         for mod in self.mods:
-            try:
-                mod_map[round(self.get_mass_or_adj(mod), 3)].append(mod)
-            except KeyError:
-                mod_map[round(self.get_mass_or_adj(mod), 3)] = [mod]
+            mass = self.get_mass_or_adj(mod)
+            for digits in nrdigits:
+                key = round(mass, digits)
+                try:
+                    mod_map[key].append(mod)
+                except KeyError:
+                    mod_map[key] = [mod]
         return mod_map
 
     def lucimass_mod_dict(self):
diff --git a/bin/nonlabile_ptm_columns.py b/bin/nonlabile_ptm_columns.py
index 86f003f..3557d82 100755
--- a/bin/nonlabile_ptm_columns.py
+++ b/bin/nonlabile_ptm_columns.py
@@ -18,6 +18,7 @@ def main():
     parser.add_argument('--stabileptms', nargs='+', default=[])
     parser.add_argument('--mods', nargs='+', default=[])
     parser.add_argument('--fasta')
+    parser.add_argument('--searchengine', default='msgf', type=str)
     args = parser.parse_args(sys.argv[1:])
 
     ptms = [x.lower() for x in args.stabileptms]
@@ -26,26 +27,29 @@ def main():
 
     msgfmods = Mods()
     msgfmods.parse_msgf_modfile(args.modfile, [*locptms, *ptms, *mods])
-    msgf_mod_map = msgfmods.msgfmass_mod_dict()
+    msgf_mod_map = msgfmods.msgfmass_mod_dict(args.searchengine)
     tdb = SeqIO.index(args.fasta, 'fasta')
+    fields = lucp.get_headerfields(args.searchengine)
+    print(fields)
     with open(args.psms) as fp, open(args.outfile, 'w') as wfp:
         header = next(fp).strip('\n').split('\t')
-        outheader = header + lucp.PTMFIELDS
+        outheader = header + fields['PTMFIELDS']
         wfp.write('\t'.join(outheader))
 
         for psm in fp:
             psm = psm.strip('\n').split('\t')
             psm = {k: v for k,v in zip(header, psm)}
-            psm.update({x: 'NA' for x in lucp.PTMFIELDS})
+            psm.update({x: 'NA' for x in fields['PTMFIELDS']})
             ptmpsm = PSM()
-            ptmpsm.parse_msgf_peptide(psm[lucp.PEPTIDE], msgf_mod_map, locptms, ptms)
+            ptmpsm.parse_search_peptide(psm[fields['PEPTIDE']], msgf_mod_map, locptms, ptms,
+                    args.searchengine)
             if ptmpsm.has_stableptms() and not ptmpsm.has_labileptms():
-                psm[lucp.TOPPTM] = ptmpsm.topptm_output()
+                psm[fields['TOPPTM']] = ptmpsm.topptm_output()
                 # Add protein location annotation
-                if lucp.MASTER_PROTEIN in psm:
-                    lucp.annotate_protein_and_flanks(psm, ptmpsm, tdb, ptms)
-                psm[lucp.SE_PEPTIDE] = psm.pop(lucp.PEPTIDE)
-                psm[lucp.PEPTIDE] = '{}_{}'.format(ptmpsm.sequence, psm[lucp.TOPPTM])
+                if fields['MASTER_PROTEIN'] in psm:
+                    lucp.annotate_protein_and_flanks(psm, ptmpsm, tdb, ptms, fields)
+                psm[fields['SE_PEPTIDE']] = psm.pop(fields['PEPTIDE'])
+                psm[fields['PEPTIDE']] = '{}_{}'.format(ptmpsm.sequence, psm[fields['TOPPTM']])
                 wfp.write('\n{}'.format('\t'.join([psm[k] for k in outheader])))
 
 
diff --git a/bin/qc_protein.R b/bin/qc_protein.R
index a4d93aa..bf7b383 100755
--- a/bin/qc_protein.R
+++ b/bin/qc_protein.R
@@ -46,6 +46,27 @@ if (is_isobaric) {
   nrpsmscols = colnames(feats)[grep('_Fully.quanted.PSM.count', colnames(feats))]
 }
 
+tmtcolors = c(
+  "#b15928",
+  "#e31a1c",
+  "#fb9a99",
+  "#ff7f00",
+  "#fdbf6f",
+  "gold4",
+  "gold1",
+  "#33a02c",
+  "#b2df8a",
+  "#01665e",
+  "#80cdc1",
+  "#1f78b4",
+  "#a6cee3",
+  "#6a3d9a",
+  "#cab2d6",
+  "#c51b7d",
+  "#f1b6da",
+  "magenta4"
+  )
+
 # nrpsms used for isobaric quant
 if (is_isobaric) {
   nrpsms = pivot_longer(feats, any_of(nrpsmscols), names_to='psmset', values_to='psmcount') 
@@ -202,7 +223,10 @@ if (is_isobaric) {
 	    group_by(Set, channel) %>%
 	    boxplot_stats(value)
   # Plot using geom_crossbar, geom_errorbar
-  ggp = ggplot(summary_stats, aes(x=Set, y=middle, fill=channel)) +
+  allchannels = unique(summary_stats$channel)
+  ggp = ggplot(summary_stats, aes(x=Set, y=middle, fill=channel, width=0.35)) +
+    scale_fill_manual(values=tmtcolors[1:length(allchannels)]) +
+    scale_colour_manual(values=tmtcolors[1:length(allchannels)]) +
     coord_flip() + ylab('Fold change') +
     xlab('Channels') + theme_bw() + 
     theme(axis.title=element_text(size=10), axis.text=element_text(size=10)) + 
@@ -226,6 +250,7 @@ if (is_isobaric) {
     norms$variable = sub('[a-z].*[0-9]*plex_', '', norms$variable)
     ggp = ggplot(norms, aes(Set, value, group=variable)) + 
       geom_col(aes(fill=variable), position=position_dodge(width=1)) +
+      scale_fill_manual(values=tmtcolors[1:length(allchannels)]) +
       geom_text(aes(y=min(value), label=round(value, 4)), position=position_dodge(width=1), colour="gray20", size=3, hjust=0) +
       coord_flip() + ylab('Normalizing factor') + xlab('Channels') + theme_bw() + 
       theme(axis.title=element_text(size=15), axis.text=element_text(size=10)) + 
@@ -367,21 +392,8 @@ if (length(deqpval_cols)) {
   }
 }
 
-
-# PCA
-if (is_isobaric && use_sampletable) {
-  topca = na.omit(feats[,tmtcols])
-  if (nrow(topca)) {
-    pca_ana <- prcomp(t(topca), scale. = TRUE)
-    score.df <- as.data.frame(pca_ana$x)
-    rownames(score.df) = sub('_[a-z0-9]*plex', '', rownames(score.df))
-    if (length(sampletable[sampletable$group != 'NO__GROUP', 'group'])) {
-        colortype = 'group'
-    } else {
-        colortype = 'set'
-    }
-    score.df$type = sampletable[rownames(score.df), colortype]
-  
+run_pca = function(scoredf, colortype) {
+    scoredf$type = sampletable[rownames(scoredf), colortype]
     #Scree plot
     contributions <- data.frame(contrib=round(summary(pca_ana)$importance[2,] * 100, 2)[1:20])
     contributions$pc = sub('PC', '', rownames(contributions))
@@ -391,10 +403,10 @@ if (is_isobaric && use_sampletable) {
       ylab("Contribution (%)") + xlab('PC (ranked by contribution)')
     p = ggplotly(ggp, width=400) %>%
             layout(legend = list(title='', orientation = 'h', x = 0, y = 1.1, xanchor='left', yanchor='bottom'))
-    htmlwidgets::saveWidget(p, glue('scree.html'), selfcontained=F)
+    htmlwidgets::saveWidget(p, glue('scree_{colortype}.html'), selfcontained=F)
 
     # PCA plot
-    ggp = ggplot(data=score.df, aes(x=PC1, y=PC2, label=rownames(score.df), colour=type)) +
+    ggp = ggplot(data=scoredf, aes(x=PC1, y=PC2, label=rownames(scoredf), colour=type)) +
       geom_hline(yintercept = 0, colour = "gray65") +
       geom_vline(xintercept = 0, colour = "gray65") +
       geom_point(size=4) +
@@ -403,6 +415,24 @@ if (is_isobaric && use_sampletable) {
       xlab(sprintf("PC1 (%s%%)", contributions$contrib[1])) + ylab(sprintf("PC2 (%s%%)", contributions$contrib[2]))
     p = ggplotly(ggp, width=400) %>%
             layout(legend = list(title='', orientation = 'h', x = 0, y = 1.1, xanchor='left', yanchor='bottom'))
-    htmlwidgets::saveWidget(p, glue('pca.html'), selfcontained=F)
+    htmlwidgets::saveWidget(p, glue('pca_{colortype}.html'), selfcontained=F)
+}
+
+
+# PCA
+if (is_isobaric && use_sampletable) {
+  topca = na.omit(feats[,tmtcols])
+  if (nrow(topca)) {
+    pca_ana <- prcomp(t(topca), scale. = TRUE)
+    score.df <- as.data.frame(pca_ana$x)
+    rownames(score.df) = sub('_[a-z0-9]*plex', '', rownames(score.df))
+    if (length(sampletable[sampletable$group != 'NO__GROUP', 'group'])) {
+        colortypes = c('group', 'set')
+    } else {
+        colortypes = c('set')
+    }
+    for (colortype in colortypes) {
+      run_pca(score.df, colortype)
+    }
   }
 }
diff --git a/bin/qc_psms.R b/bin/qc_psms.R
index 4342f94..59e2bdc 100755
--- a/bin/qc_psms.R
+++ b/bin/qc_psms.R
@@ -7,19 +7,52 @@ library(stringr)
 
 args = commandArgs(trailingOnly=TRUE)
 has_fractions = args[1] == TRUE
+search_engine = args[2]
 
 # Fraction needs to be a factor and not numeric, since it can contain strings (e.g. A2)
 # That would be fine but when mixing oldmzmls from a rerun in, the join on Fraction op may fail
 # when there is both the factor 02, and the numeric 2 in different tables
 feats = read.table("psms_clean", colClasses=c('Fraction'='factor'), header=T, sep="\t", comment.char = "", quote = "")
 
-scancol = 'SpecID' # MSGF (sage has scannr)
-miscleavcol = 'missed_cleavage'
-rtcol = 'Retention.time.min.' # 'rt'
-ioncol = 'Ion.injection.time.ms.' # ?
-precerrcol = 'PrecursorError.ppm.' # precursor_ppm
-scorecol = 'MSGFScore' # 'sage_discriminant_score'
-filenamecol = 'SpectraFile' # 'filename'
+if (search_engine == 'sage') {
+  scancol = 'scannr'
+  miscleavcol = 'missed_cleavages'
+  rtcol = 'rt'
+  ioncol = 'Ion.injection.time.ms.'
+  precerrcol = 'precursor_ppm'
+  scorecol = 'sage_discriminant_score'
+  filenamecol = 'filename'
+
+} else if (search_engine == 'msgf') {
+  scancol = 'SpecID' # MSGF (sage has scannr)
+  miscleavcol = 'missed_cleavage'
+  rtcol = 'Retention.time.min.' # 'rt'
+  ioncol = 'Ion.injection.time.ms.' # ?
+  precerrcol = 'PrecursorError.ppm.' # precursor_ppm
+  scorecol = 'MSGFScore' # 'sage_discriminant_score'
+  filenamecol = 'SpectraFile' # 'filename'
+}
+
+tmtcolors = c(
+  "#b15928",
+  "#e31a1c",
+  "#fb9a99",
+  "#ff7f00",
+  "#fdbf6f",
+  "gold4",
+  "gold1",
+  "#33a02c",
+  "#b2df8a",
+  "#01665e",
+  "#80cdc1",
+  "#1f78b4",
+  "#a6cee3",
+  "#6a3d9a",
+  "#cab2d6",
+  "#c51b7d",
+  "#f1b6da",
+  "magenta4"
+  )
 
 boxplot_stats = function(data, col) {
   summary_stats = data %>%
@@ -83,8 +116,8 @@ p = ggplotly(ggp, width=400, height=vert_height) %>%
         layout(legend = list(orientation = 'h', x = 0, y = 1.1, xanchor='left', yanchor='bottom'))
 # Work around since plotly does not honor above legend.title=element_blank call
 p$x$layout$legend$title$text = ''
-
 htmlwidgets::saveWidget(p, 'amount_psms.html', selfcontained=F)
+write.table(amount_id, 'psms_ids.txt', row.names=F, quote=F, sep='\t')
 
 # Missing isobaric values
 if (length(grep('plex', names(feats)))) {
@@ -96,8 +129,13 @@ if (length(grep('plex', names(feats)))) {
     psm_empty = aggregate(value~get(xcol)+ name, psm_empty, sum)
     names(psm_empty) = c(xcol, 'channels', 'nr_missing_values')
     psm_empty$channels = sub('.*plex_', '', psm_empty$channels)
+    allchannels = unique(psm_empty$channels)
     ggp = ggplot(psm_empty, aes(x=.data[[ xcol ]], y=nr_missing_values, fill=channels)) + 
-      geom_bar(stat='identity', position="dodge") + ylab('# PSMs without quant') + coord_flip() + theme_bw() + theme(axis.title.x=element_text(size=15), axis.title.y=element_blank(), axis.text=element_text(size=10), axis.text.y=element_text(angle=90), legend.position="top", legend.text=element_text(size=10), legend.title=element_blank())
+      geom_bar(stat='identity', position="dodge") + 
+      scale_fill_manual(values=tmtcolors[1:length(allchannels)]) +
+      ylab('# PSMs without quant') + coord_flip() + theme_bw() +
+      theme(axis.title.x=element_text(size=15), axis.title.y=element_blank(), axis.text=element_text(size=10), axis.text.y=element_text(angle=90), legend.position="top", legend.text=element_text(size=10), legend.title=element_blank()
+    )
     p = ggplotly(ggp, width=400, height=vert_height) %>%
         layout(legend = list(orientation = 'h', x = 0, y = 1.1, xanchor='left', yanchor='bottom'))
   p$x$layout$legend$title$text = ''
@@ -107,16 +145,16 @@ if (length(grep('plex', names(feats)))) {
 
 # Missed cleavages
 mcl = aggregate(get(scancol)~get(xcol)+get(miscleavcol), feats, length)
-colnames(mcl) = c(xcol, 'missed_cleavage', 'nrscan')
+colnames(mcl) = c(xcol, 'missed_cleavage', 'nrpsms')
 mcl_am = subset(merge(mcl, amount_psms, by=xcol), missed_cleavage %in% c(0,1,2))
-mcl_am$percent = mcl_am$nrscan / mcl_am$"PSMs IDed" * 100
+mcl_am$percent = mcl_am$nrpsms / mcl_am$"PSMs IDed" * 100
 mc_text_y = max(mcl_am$percent) * 2/6
 mcl_am$missed_cleavage = as.factor(mcl_am$missed_cleavage)
 
 mcplot = ggplot(mcl_am) +
     geom_bar(aes(x=.data[[ xcol ]], y=percent, fill=missed_cleavage, group=missed_cleavage), position='dodge', stat='identity') +
     # 0.9 is the default dodge (90% of 1, 1 used bc all same value) but when not spec -> no dodge at all?
-    geom_text(position=position_dodge(width=0.9), aes(x=.data[[xcol]], y=mc_text_y, group=missed_cleavage, label=glue('{nrscan} PSMs')), colour="black", size=4, inherit.aes=T) +
+    geom_text(position=position_dodge(width=0.9), aes(x=.data[[xcol]], y=mc_text_y, group=missed_cleavage, label=glue('{nrpsms} PSMs')), colour="black", size=4, inherit.aes=T) +
     ylim(c(0, 100)) + ylab('% of PSMs') +
     theme_bw() +
     theme(axis.title.x=element_text(size=15), axis.title.y=element_blank(), axis.text=element_text(size=10), axis.text.y=element_text(angle=90), legend.position="top", legend.text=element_text(size=10), legend.title=element_blank()) +
@@ -125,6 +163,7 @@ p = ggplotly(mcplot, width=400, height=vert_height) %>%
         layout(legend = list(orientation = 'h', x = 0, y = 1.1, xanchor='left', yanchor='bottom'))
 p$x$layout$legend$title$text = ''
 htmlwidgets::saveWidget(p, 'missed_cleavages.html', selfcontained=F)
+write.table(mcl_am, 'miscleav.txt', row.names=F, quote=F, sep='\t')
 
 
 # Now the per-fraction or per-file stats
diff --git a/bin/qc_ptms.R b/bin/qc_ptms.R
index 05f6a91..4e8dd2a 100755
--- a/bin/qc_ptms.R
+++ b/bin/qc_ptms.R
@@ -9,11 +9,20 @@ library(stringr)
 args = commandArgs(trailingOnly=TRUE)
 psmfn = args[1]
 pepfn = args[2]
+search_engine = args[3]
 
 psms = read.table(psmfn, header=T, sep="\t", comment.char = "", quote = "")
 peptides = read.table(pepfn, header=T, sep="\t", comment.char="", quote="")
 
-scancol = 'SpecID'
+if (search_engine == 'sage') {
+  scancol = 'scannr'
+  pepcol = 'peptide'
+
+} else if (search_engine == 'msgf') {
+  scancol = 'SpecID' # MSGF (sage has scannr)
+  pepcol = 'Peptide'
+}
+
 
 ## Prepare a data frame with all the info for plots
 
@@ -25,7 +34,7 @@ repeater = sapply(multiptms, FUN=length)
 sites = data.frame(sites=unlist(multiptms),
                    specid=rep(psms[[ scancol ]], repeater),
                    bioset=rep(psms$Biological.set, repeater),
-                   peptide=rep(psms$Peptide, repeater)
+                   peptide=rep(psms[[ pepcol ]], repeater)
 )
 if ('Master.protein.s.' %in% names(psms)) {
     sites$protein = rep(psms$Master.protein.s., repeater)
diff --git a/bin/report_tables.py b/bin/report_tables.py
index 6c9222a..b9e3db0 100755
--- a/bin/report_tables.py
+++ b/bin/report_tables.py
@@ -109,11 +109,17 @@ def get_plotly_html(fn):
         featplots[plotname] = False
 
 expplotnames = [
-          ('pca', 'Principal component analysis'),
-          ('scree', 'PCA Scree plot'),
+          ('pcagroup', 'Principal component analysis'),
+          ('screegroup', 'PCA Scree plot'),
           ]
-expplotfns = [('pca', 'pca.html'),
-        ('scree', 'scree.html'),
+pcaplotnames = [
+          ('pcaset', 'Principal component analysis'),
+          ('screeset', 'PCA Scree plot'),
+          ]
+expplotfns = [('pcaset', 'pca_set.html'),
+        ('pcagroup', 'pca_group.html'),
+        ('screeset', 'scree_set.html'),
+        ('screegroup', 'scree_group.html'),
         ]
 expplots = defaultdict(defaultdict)
 for plotname, pfn in expplotfns:
@@ -257,6 +263,27 @@ def get_plotly_html(fn):
     summary_fields = [x for x in summary_field_order if x in fields]
     break
 
+# PSM tables
+psmtables = {'ids': [], 'miscleav': []}
+with open('psmids') as fp:
+    head = next(fp).strip().split()
+    plates = defaultdict(defaultdict)
+    for line in fp:
+        lnmap = {head[ix]: x for ix, x in enumerate(line.strip().split('\t'))}
+        if lnmap['name'] == 'MS2 scans':
+            plates[lnmap['plateID']]['scans'] = lnmap['count']
+        elif lnmap['name'] == 'PSMs IDed':
+            plates[lnmap['plateID']].update({'psms': lnmap['count'], 'pc': lnmap['labeltext']})
+psmtables['ids'] = [[p, nms['scans'], nms['psms'], nms['pc']] for p, nms in plates.items()]
+
+with open('miscleav') as fp:
+    head = next(fp).strip().split()
+    plates = defaultdict()
+    for line in fp:
+        lnmap = {head[ix]: x for ix, x in enumerate(line.strip().split('\t'))}
+        print(lnmap)
+        psmtables['miscleav'].append([lnmap['plateID'], lnmap['missed_cleavage'], lnmap['nrpsms'], lnmap['IDed']])
+
 
 # Overlap
 overlap = defaultdict(dict)
@@ -305,9 +332,11 @@ def get_plotly_html(fn):
         plates=args.plates,
         expplots=expplots,
         expplotnames=expplotnames,
+        pcaplotnames=pcaplotnames,
         deqmsplots=deqmsplots,
         deqmscomps=deqmscomps,
         tabletitles=tabletitles,
+        psmtables=psmtables,
         summary_fields=summary_fields,
         summary_table=summary_table,
         overlap=overlap,
diff --git a/conf/base.config b/conf/base.config
index 5b34df6..ae0463a 100644
--- a/conf/base.config
+++ b/conf/base.config
@@ -16,6 +16,10 @@ process {
     time = { check_max( 6.h * task.attempt, 'time' ) }
     memory = { (db.size() >> 30) < 1 ? 16.GB : "${db.size() * 16}B"  }
   }
+  withName: sage {
+    cpus = 2
+    memory = { (db.size() >> 30) < 1 ? 16.GB : "${db.size() * 16}B"  }
+  }
   withName: luciphorPTMLocalizationScoring {
     cpus = 2
   }
diff --git a/conf/test.config b/conf/test.config
index 07c2cc1..698554d 100644
--- a/conf/test.config
+++ b/conf/test.config
@@ -8,6 +8,11 @@ process {
     time = { check_max( 0.5.h * task.attempt, 'time' ) }
     memory = 4.GB
   }
+  withName: sage {
+    cpus = 2
+    time = { check_max( 0.5.h * task.attempt, 'time' ) }
+    memory = 8.GB
+  }
 }
 
 docker {
diff --git a/containers.json b/containers.json
index f4c28c0..78ae30f 100644
--- a/containers.json
+++ b/containers.json
@@ -1,8 +1,8 @@
 {
   "msstitch": {
-    "version": "3.16",
-    "singularity": "https://depot.galaxyproject.org/singularity/msstitch:3.16--pyhdfd78af_0",
-    "docker": "quay.io/biocontainers/msstitch:3.16--pyhdfd78af_0"
+    "version": "3.17",
+    "singularity": "https://depot.galaxyproject.org/singularity/msstitch:3.17--pyhdfd78af_0",
+    "docker": "quay.io/biocontainers/msstitch:3.17--pyhdfd78af_0"
   },
   "python": {
     "version": "3.12",
@@ -15,9 +15,9 @@
       "singularity": "docker://proteowizard/pwiz-skyline-i-agree-to-the-vendor-licenses:3.0.24172-63d00b1"
   },
   "openms": {
-      "version": "3.1.0",
-      "docker": "quay.io/biocontainers/openms:3.1.0--h191ead1_4",
-      "singularity": "https://depot.galaxyproject.org/singularity/openms:3.1.0--h191ead1_4"
+      "version": "3.3.0",
+      "docker": "quay.io/biocontainers/openms:3.3.0--h0656172_8",
+      "singularity": "https://depot.galaxyproject.org/singularity/openms:3.3.0--h0656172_8"
   },
   "dinosaur": {
       "version": "1.2.0",
@@ -50,14 +50,24 @@
       "singularity": "docker://ghcr.io/lehtiolab/ddamsproteomics:3.0"
   },
   "msgfplus": {
-      "version": "2023.01.1202",
-      "docker": "quay.io/biocontainers/msgf_plus:2023.01.1202--hdfd78af_0",
-      "singularity": "https://depot.galaxyproject.org/singularity/msgf_plus:2023.01.1202--hdfd78af_0"
+      "version": "2024.03.26",
+      "docker": "quay.io/biocontainers/msgf_plus:2024.03.26--hdfd78af_0",
+      "singularity": "https://depot.galaxyproject.org/singularity/msgf_plus:2024.03.26--hdfd78af_0"
+  },
+  "jq": {
+      "version": "1.5",
+      "docker": "quay.io/biocontainers/jq:1.5--4",
+      "singularity": "https://depot.galaxyproject.org/singularity/jq:1.5--4"
+  },
+  "sage": {
+      "version": "0.14.7",
+      "docker": "ghcr.io/lazear/sage:v0.14.7",
+      "singularity": "https://depot.galaxyproject.org/singularity/sage-proteomics:0.14.7--hc1c3326_2"
   },
   "percolator": {
-      "version": "3.5",
-      "docker": "quay.io/biocontainers/percolator:3.5--hfd1433f_1",
-      "singularity": "https://depot.galaxyproject.org/singularity/percolator:3.5--hfd1433f_1"
+      "version": "3.6.5",
+      "docker": "quay.io/biocontainers/percolator:3.6.5--h6351f2a_0",
+      "singularity": "https://depot.galaxyproject.org/singularity/percolator:3.6.5--h6351f2a_0"
   },
   "luciphor2": {
       "version": "2020_04_03",
diff --git a/docs/output.md b/docs/output.md
index 1e310c4..441124e 100644
--- a/docs/output.md
+++ b/docs/output.md
@@ -26,6 +26,7 @@ table is filtered on a PSM and peptide FDR of 0.01.
 * Protein group(s) content (Protein groups have members which can be explained by the representative protein)
 * Total number of matching proteins in group (Those proteins in the group that this PSM matches to)
 * Gene ID, Gene Name, Description (gene identity and protein description information from fasta)
+* Explained ion current ratio (explained ions divided by total current)
 * percolator svm-score
 * PSM q-value and posterior error probability (PEP) (both calculated as TD competition by percolator)
 * peptide q-value and PEP (as PSM q-value but for peptides)
@@ -97,8 +98,13 @@ and processes data using the following steps:
 * [Msstitch](#msstitch) - Post processing, protein inference
 * [DEqMS](#deqms) - Differential expression analysis
 
+
 ## MSGF+
-[MSGF+](https://omics.pnl.gov/software/ms-gf) (aka MSGF+ or MSGFPlus) performs peptide identification by scoring MS/MS spectra against peptides derived from a protein sequence database. [PMID 25358478](https://pubmed.ncbi.nlm.nih.gov/25358478/)
+[MSGF+](https://msgfplus.github.io/)(aka MSGF+ or MSGFPlus) performs peptide identification by scoring MS/MS spectra against peptides derived from a protein sequence database. [PMID 25358478](https://pubmed.ncbi.nlm.nih.gov/25358478/)
+
+
+## Sage
+[Sage](https://sage-docs.vercel.app/) is another search engine to do peptide identification with, and is very fast.
 
 
 ## Percolator
diff --git a/docs/usage.md b/docs/usage.md
index 935970c..d441942 100644
--- a/docs/usage.md
+++ b/docs/usage.md
@@ -41,7 +41,7 @@ NXF_OPTS='-Xms1g -Xmx4g'
 ## Running the pipeline
 The typical command for running the pipeline is as follows:
 ```bash
-nextflow run lehtiolab/ddamsproteomics --input /path/to/input_definition.txt --tdb /path/to/proteins.fa --mods 'oxidation;carbamidomethylation' -profile standard,docker
+nextflow run lehtiolab/ddamsproteomics --input /path/to/input_definition.txt --sage --tdb /path/to/proteins.fa --mods 'oxidation;carbamidomethylation' -profile standard,docker
 ```
 
 This will launch the pipeline with the `docker` configuration profile. See below for more information about profiles.
@@ -95,6 +95,13 @@ Use this parameter to choose a configuration profile. Profiles can give configur
 * `none`
     * No configuration at all. Useful if you want to build your own config from scratch and want to avoid loading in the default `base` config profile (not recommended).
 
+### `--sage`
+Specifies to use the [Sage](https://sage-docs.vercel.app/) search engine in the pipeline
+
+### `--msgf`
+Specifies to use the [MSGF+](https://msgfplus.github.io/)
+) search engine
+
 ### `--input`
 This param passes an mzML definition (txt) file which contains the mzML specifications. This also enables runs with specific fractionation such as HiRIEF or high pH, and the specification of individual instruments per file.
 
@@ -175,6 +182,10 @@ To annotate peptide/protein/gene results after a search, marking hits in another
 one column for each file in the result files, containing the fasta IDs for the record that a peptide (or any peptide from a protein) matched to.
 
 
+### Removing channels
+In case an experiment contains empty isobaric channels, you can remove those from the run by specifying e.g. `--remove_channels setA:127N:129C setB:126`
+
+
 ### Reusing data
 If you have finished a rather large analysis and wish to rerun a part of it or add more fractions, due to e.g. new MS data, you may do so by passing
 
diff --git a/main.nf b/main.nf
index e560050..1407a90 100644
--- a/main.nf
+++ b/main.nf
@@ -1,9 +1,10 @@
 #!/usr/bin/env nextflow
 
-include { paramsSummaryMap } from 'plugin/nf-validation'
+include { paramsSummaryMap } from 'plugin/nf-schema'
 
-include { msgf_info_map; listify; stripchars_infile; get_regex_specialchars; read_header } from './modules.nf' 
+include { msgf_info_map; get_complement_field_nr; listify; stripchars_infile; get_regex_specialchars; read_header } from './modules.nf' 
 include { MSGFPERCO } from './workflows/msgf_perco.nf'
+include { SAGEPERCO } from './workflows/sage_perco.nf'
 include { PTMANALYSIS } from './workflows/ptms.nf'
 include { MATCH_SEQUENCES } from './workflows/match_sequences.nf'
 include { REPORTING } from './workflows/reporting.nf'
@@ -241,7 +242,7 @@ process createNewSpectraLookup {
   container params.__containers[tag][workflow.containerEngine]
 
   input:
-  tuple val(setnames), file(mzmlfiles), val(platenames)
+  tuple val(setnames), path(mzmlfiles), val(platenames)
 
   output:
   path('target_db.sqlite')
@@ -318,7 +319,7 @@ process peptidePiAnnotation {
   container params.__containers[tag][workflow.containerEngine]
 
   input:
-  tuple path('psms'), val(strips), path('hirief_training_pep')
+  tuple path('psms'), val(strips), val(search_engine), path('hirief_training_pep')
 
   output:
   path('target_psmtable.txt')
@@ -326,7 +327,9 @@ process peptidePiAnnotation {
   script:
   """
   echo '${groovy.json.JsonOutput.toJson(strips)}' >> strip.json
-  peptide_pi_annotator.py -i hirief_training_pep -p psms --out target_psmtable.txt --stripcolpattern Strip --pepcolpattern Peptide --fraccolpattern Fraction --stripdef strip.json --ignoremods '*'
+  peptide_pi_annotator.py -i hirief_training_pep -p psms --out target_psmtable.txt \
+    --fraccolpattern Fraction --stripdef strip.json --ignoremods '*' --stripcolpattern Strip \
+    --pepcolpattern ${search_engine == 'sage' ? 'peptide' : 'Peptide'}
   """
 }
 
@@ -337,7 +340,7 @@ process splitPSMs {
   container params.__containers[tag][workflow.containerEngine]
 
   input:
-  tuple val(td), path('psms'), val(setnames)
+  tuple val(td), path('psms'), val(setnames), val(remove_channels)
 
   output:
   tuple val(td), path({listify(setnames).collect { "${it}.tsv" }}) optional true
@@ -345,10 +348,19 @@ process splitPSMs {
   script:
   """
   msstitch split -i psms --splitcol bioset
+  ${td == 'target' ?
+    remove_channels.collect {
+      setn, chs -> chs.collect {
+        ch -> "colnum=${get_complement_field_nr("${setn}.tsv", ch)} && \
+        cut -f \$colnum ${setn}.tsv > tmprm && mv tmprm ${setn}.tsv"
+      }.join(' && ')
+    }.join(' && ')
+  : ''}
   """
 }
 
 
+
 process splitTotalProteomePSMs {
 
   tag 'msstitch'
@@ -358,11 +370,14 @@ process splitTotalProteomePSMs {
   tuple path('tppsms_in'), val(setnames)
 
   output:
-  path({listify(setnames).collect() { "tppsms/${it}.tsv" }}), emit: totalprotpsms_allsets optional true
+  path({listify(setnames).collect() { "tppsms/${it}.tsv" }}) optional true
 
   script:
   """
   mkdir tppsms && msstitch split -i tppsms_in -d tppsms --splitcol bioset
+  ${listify(setnames).collect { set -> "[ -e 'tppsms/${set}.tsv' ] || (\
+    echo Not all sets in the experiment are in the --totalproteomepsms table, \
+    we need: ${listify(setnames).join(', ')}. The total proteome PSM table contains: \$(ls tppsms) && exit 1)" }.join( '&&')}
   """
 }
 
@@ -448,13 +463,20 @@ process sampleTableCheckClean {
   container params.__containers[tag][workflow.containerEngine]
  
   input:
-  tuple path('sampletable'), val(do_deqms)
+  tuple path('sampletable'), val(do_deqms), val(remove_channels)
 
   output:
   tuple path('clean_sampletable'), path('sampletable_no_special_chars')
   
   script:
   """
+  # Remove empty channels
+  ${remove_channels.collect {
+    setch -> setch[1].collect {
+      ch -> "grep -v '^${ch}\t${setch[0]}' sampletable > tmpst && mv tmpst sampletable"
+      }.join(' && ')
+    }.join(' && ')
+  }
   # First add NO__GROUP marker for no-samplegroups clean sampletable from special chars
   awk -v FS="\\t" -v OFS="\\t" \'{if (NF==3) print \$1,\$2,\$3,"NO__GROUP"; else print}\' sampletable > clean_sampletable
   # Check if there are samplegroups at all
@@ -627,6 +649,23 @@ workflow {
   }.collectEntries() {
     x-> [x[0], x[2..-1]]
   } : [:]
+  // Remove channels from specific sets if those are empty: --remove_channels 'setA:126:127 setB:131'
+  rmch = params.remove_channels ? params.remove_channels.tokenize(' ') : false
+  remove_channels_psmtable = rmch ? rmch.collect { y -> y.tokenize(':')
+  }.collect { x -> [x[0], x[1..-1].collect { ch -> "${setisobaric[x[0]]}_${ch}" } ] } : [:]
+  remove_channels_sampletable = rmch ? rmch.collect { y -> y.tokenize(':')
+  }.collect { x -> [x[0], x[1..-1]] } : [:]
+  rm_ch_err = []
+  remove_channels_psmtable.each { sn, chs -> 
+    if (!(sn in setisobaric)) {
+      rm_ch_err.push("Set ${sn} not in --isobaric.")
+    }
+  }
+  if (rm_ch_err) {
+    exit 1, "Errors in --remove_channels: ${rm_ch_err.join(', ')}, please check your isobaric channel input"
+  }
+  
+
   
   do_ms1 = !params.noquant && !params.noms1quant
   do_normalize = (!params.noquant && (params.mediannormalize || params.deqms) && params.isobaric)
@@ -755,34 +794,63 @@ workflow {
   tdb = Channel.fromPath(params.tdb)
   createTargetDecoyFasta(tdb)
 
+  if (params.sage) {
+    search_engine = 'sage'
+  } else if (params.msgf) {
+    search_engine = 'msgf'
+  } else {
+    exit 1, 'Must specify --sage or --msgf'
+  }
+
   if (!is_rerun) {
     search_mods = [params.mods ? params.mods : false,
       params.ptms ?: false,
       params.locptms ?: false,
       ].findAll { it }
       .join(';')
-    MSGFPERCO(mzml_in,
-      createTargetDecoyFasta.out.concatdb,
-      setisobaric,
-      fractionation,
-      mzmls,
-      params.maxvarmods,
-      params.msgfmods,
-      search_mods,
-      params.psmconflvl,
-      params.pepconflvl,
-      params.locptms,
-      params.maxmiscleav,
-      params.enzyme,
-      params.minpeplen,
-      params.maxpeplen,
-      params.mincharge,
-      params.maxcharge,
-    )
-  
-    MSGFPERCO.out.t_tsv
+    if (params.sage) {
+      SAGEPERCO(mzml_in,
+        createTargetDecoyFasta.out.concatdb,
+        setisobaric,
+        fractionation,
+        mzmls,
+        params.maxvarmods,
+        params.msgfmods,
+        search_mods,
+        params.psmconflvl,
+        params.pepconflvl,
+        params.locptms,
+        params.maxmiscleav,
+        params.enzyme,
+        params.minpeplen,
+        params.maxpeplen,
+        params.mincharge,
+        params.maxcharge,
+      ).set { SEARCH }
+    } else {
+      MSGFPERCO(mzml_in,
+        createTargetDecoyFasta.out.concatdb,
+        setisobaric,
+        fractionation,
+        mzmls,
+        params.maxvarmods,
+        params.msgfmods,
+        search_mods,
+        params.psmconflvl,
+        params.pepconflvl,
+        params.locptms,
+        params.maxmiscleav,
+        params.enzyme,
+        params.minpeplen,
+        params.maxpeplen,
+        params.mincharge,
+        params.maxcharge,
+      ).set { SEARCH }
+    }
+
+    SEARCH.t_tsv
     | ifEmpty(['target', nofile])
-    | concat(MSGFPERCO.out.d_tsv)
+    | concat(SEARCH.d_tsv)
     | groupTuple
     | join(specquant_lookups)
     | combine(createTargetDecoyFasta.out.bothdbs)
@@ -806,7 +874,7 @@ workflow {
     hiriefpep = Channel.fromPath(params.hirief)
     psmtables_ch
     | filter { it[0] == 'target' }
-    | map { [it[1], params.strips]  }
+    | map { [it[1], params.strips, search_engine]  }
     | combine(hiriefpep)
     | peptidePiAnnotation
     | map { ['target', it]}
@@ -819,7 +887,7 @@ workflow {
   psmtables_ch
   | filter { it[0] == 'decoy' }
   | concat(target_psmtable)
-  | map { [it[0], it[1], all_setnames] }
+  | map { [it[0], it[1], all_setnames, remove_channels_psmtable] }
   | splitPSMs
   | map{ it -> [it[0], listify(it[1]).collect() { it.baseName.replaceFirst(/\.tsv$/, "") }, it[1]]} // get setname from {setname}.tsv
   | transpose
@@ -854,8 +922,9 @@ workflow {
       setisobaric,
       mzml_in | map { [it.setname, it.mzmlfile] } | groupTuple,
       splitpsms_ch | filter { it[0] == 'target' } | map { it[1..-1] },
+      search_engine,
       tdb,
-      !is_rerun ? MSGFPERCO.out.unfiltered : Channel.empty(),
+      !is_rerun ? SEARCH.unfiltered : Channel.empty(),
       mzml_list,
       params.maxpeplen,
       params.maxcharge,
@@ -917,7 +986,7 @@ workflow {
 
   if (params.sampletable) {
     Channel.fromPath(params.sampletable)
-    | map { [it, params.deqms] }
+    | map { [it, params.deqms, remove_channels_sampletable] }
     | sampleTableCheckClean
     | set { sampletable_ch }
   } else {
@@ -961,6 +1030,7 @@ workflow {
   REPORTING(
     sorted_mzml_in,
     tdspeclookup.map { it[0] },
+    search_engine,
     oldmzmls,
     fractionation,
     target_psmtable.map { it[1] },
@@ -972,7 +1042,7 @@ workflow {
     all_setnames,
     ptm_ch,
     psmwarnings_ch
-      .concat(!is_rerun ? MSGFPERCO.out.warnings: Channel.empty())
+      .concat(!is_rerun ? SEARCH.warnings: Channel.empty())
       .concat(ptmwarn_ch)
       .toList().toList()
       .filter { it[0] }
diff --git a/modules.nf b/modules.nf
index e325cc9..5efd4dd 100644
--- a/modules.nf
+++ b/modules.nf
@@ -11,6 +11,10 @@ def get_field_nr_multi(fn, fieldnames) {
     return "\$(head -n1 ${fn} | tr '\\t' '\\n' | grep -En '(${fieldnames.join('|')})' | cut -f 1 -d':' | tr '\\n' ',' | sed 's/\\,\$//')"
 }
 
+def get_complement_field_nr(fn, fieldname) {
+    /* return field nrs NOT matching input, comma separated like: 1,2,5,9 */
+    return "\$(head -n1 ${fn} | tr '\\t' '\\n' | grep -vwn '^${fieldname}\$' | cut -f 1 -d':' | tr '\\n' ',' | sed 's/\\,\$//')"
+}
 
 def parse_isotype(isobtype) {
   return ['tmt16plex', 'tmt18plex'].any { it == isobtype } ? 'tmtpro' : isobtype
@@ -97,3 +101,20 @@ def msgf_info_map(info_fn) {
   return samples
 }
 
+
+process createMods {
+
+  tag 'python'
+  container params.__containers[tag][workflow.containerEngine]
+
+  input:
+  tuple val(setname), val(isobtype), val(maxvarmods), val(search_engine), path(msgfmods), val(mods)
+
+  output:
+  tuple val(setname), path('mods.txt')
+
+  script:
+  """
+  create_modfile.py $maxvarmods "${msgfmods}" $search_engine "${mods}${isobtype ? ";${parse_isotype(isobtype)}" : ''}"
+  """
+}
diff --git a/nextflow.config b/nextflow.config
index 32dcd18..134a36f 100644
--- a/nextflow.config
+++ b/nextflow.config
@@ -1,7 +1,7 @@
 import groovy.json.JsonSlurper
 
 plugins {
-  id 'nf-validation'
+  id 'nf-schema@2.2.0'
 }
 
 // Dont warn on withName in config without a process in running pipeline
@@ -16,9 +16,11 @@ params {
   external_config_version = 'master'
 
   name = false
-  input = false
+  input = null
   tdb = false
   mods = false
+  sage = false
+  msgf = false
   locptms = false
   ptms = false
   totalproteomepsms = false
@@ -28,6 +30,7 @@ params {
   phospho = false
   maxvarmods = 2
   isobaric = false
+  remove_channels = false
   instrument = 'qe' // Default instrument is Q-Exactive
   prectol = '10.0ppm'
   iso_err = '-1,2'
@@ -149,7 +152,7 @@ manifest {
   description = 'Quantitative DDA MS proteomics pipeline'
   mainScript = 'main.nf'
   nextflowVersion = '>=24.04.4'
-  version = '3.0'
+  version = '3.1'
   doi = '10.5281/zenodo.3548311'
 }
 
diff --git a/nextflow_schema.json b/nextflow_schema.json
index 5bf7100..ef25439 100644
--- a/nextflow_schema.json
+++ b/nextflow_schema.json
@@ -294,6 +294,15 @@
     "overbook_cpus_factor": {
       "type": "integer",
       "default": 1
+    },
+    "sage": {
+      "type": "boolean"
+    },
+    "msgf": {
+      "type": "boolean"
+    },
+    "remove_channels": {
+      "type": "boolean"
     }
   }
 }
diff --git a/run_tests.sh b/run_tests.sh
index cc0d6b5..68a2e52 100644
--- a/run_tests.sh
+++ b/run_tests.sh
@@ -11,6 +11,12 @@ rundir=$(pwd)
 export repodir=$(dirname "$(realpath -s "$0")")
 export testdir="${repodir}/tests/"
 export testdata="${rundir}/static-resources/test-data/ddamsproteomics"
+if [ "$1" != 'sage' ] && [ "$1" != 'msgf' ]
+then
+	echo 'Must run tests with "bash run_tests.sh sage [or msgf]"'
+	exit 1
+fi
+export NXFCMD="nextflow run -resume -profile test ${repodir}/main.nf --$1"
 
 if [ -e "${testdata}" ]
 then
@@ -33,7 +39,8 @@ test_names=(
     tims
 )
 
-if [ -z "$1" ]
+
+if [ -z "$2" ]
 then
 	# Run all tests, do not exit with failures
 	# Manual runs
@@ -54,12 +61,12 @@ then
 else
 	# Single test, exit with error when failing
 	# Used in github actions
-       	bash "${testdir}/$1.sh"
+       	bash "${testdir}/$2.sh"
        	if [[ "$?" == 0 ]]
        	then
-	       	echo ${green}"$1" - SUCCESS ${reset}
+	       	echo ${green}"$2" - SUCCESS ${reset}
        	else
-	       	echo ${red}"$1" - FAIL ${reset}
+	       	echo ${red}"$2" - FAIL ${reset}
 		exit 1
        	fi
 fi
diff --git a/tests/labelfree.sh b/tests/labelfree.sh
index 4def2ae..871bf3e 100644
--- a/tests/labelfree.sh
+++ b/tests/labelfree.sh
@@ -5,7 +5,7 @@ set -eu
 echo Normal labelfree test
 name=lf
 cat ${testdir}/lf_mzmls.txt | envsubst > test_output/mzmldef
-nextflow run -resume -profile test ${repodir}/main.nf --name ${name} \
+$NXFCMD --name ${name} \
     --outdir test_output/${name} \
     --input test_output/mzmldef \
     --genes \
@@ -23,7 +23,7 @@ nextflow run -resume -profile test ${repodir}/main.nf --name ${name} \
 echo Labelfree run with warnings: no decoy in setB, no target in setA
 name=lf_notarget
 cat ${testdir}/lf_mzmls.txt | envsubst > test_output/mzmldef
-nextflow run -resume -profile test ${repodir}/main.nf --name ${name} \
+$NXFCMD --name ${name} \
     --outdir test_output/${name} \
     --input test_output/mzmldef \
     --genes \
@@ -34,7 +34,7 @@ nextflow run -resume -profile test ${repodir}/main.nf --name ${name} \
 echo Labelfree run without fractions, with warnings: no decoy in setB, no target in setA
 name=lf_nofrac_notarget
 cat ${testdir}/lf_mzmls_nofrac.txt | envsubst > test_output/mzmldef
-nextflow run -resume -profile test ${repodir}/main.nf --name ${name} \
+$NXFCMD --name ${name} \
     --outdir test_output/${name} \
     --input test_output/mzmldef \
     --genes \
@@ -46,7 +46,7 @@ nextflow run -resume -profile test ${repodir}/main.nf --name ${name} \
 echo Single file labelfree test
 name=lf_singlefile
 cat <(head -n1 ${testdir}/lf_mzmls.txt) <(grep setA ${testdir}/lf_mzmls.txt) | envsubst > test_output/mzmldef
-nextflow run -resume -profile test ${repodir}/main.nf --name ${name} \
+$NXFCMD --name ${name} \
     --outdir test_output/${name} \
     --input test_output/mzmldef \
     --genes \
diff --git a/tests/labelfree_phos.sh b/tests/labelfree_phos.sh
index ec2870f..1e95415 100644
--- a/tests/labelfree_phos.sh
+++ b/tests/labelfree_phos.sh
@@ -5,7 +5,7 @@ echo Phospho, labelfree, one luciphor setB not enough PSMs
 name=labelfree_phos
 lfphos_dir=test_output/${name}
 cat ${testdir}/lf_mzmls.txt | envsubst > test_output/mzmldef
-nextflow run -resume -profile test ${repodir}/main.nf --name ${name} \
+$NXFCMD --name ${name} \
     --outdir ${lfphos_dir} \
     --input test_output/mzmldef \
     --genes \
@@ -19,7 +19,7 @@ nextflow run -resume -profile test ${repodir}/main.nf --name ${name} \
 # LF phos where we add replace a set with a "new" set
 name=labelfree_phos_addset
 cat <(head -n1 ${testdir}/lf_mzmls.txt) <(grep setA ${testdir}/lf_mzmls.txt) | envsubst > test_output/mzmldef
-nextflow run -resume -profile test ${repodir}/main.nf --name ${name} \
+$NXFCMD --name ${name} \
     --outdir test_output/${name} \
     --input test_output/mzmldef \
     --oldmzmldef ${testdir}/lf_mzmls.txt \
@@ -41,7 +41,7 @@ nextflow run -resume -profile test ${repodir}/main.nf --name ${name} \
 echo LF phos no decoy
 name=labelfree_phos_nodecoy
 cat <(head -n1 ${testdir}/lf_mzmls.txt) <(grep setA ${testdir}/lf_mzmls.txt) | envsubst > test_output/mzmldef
-nextflow run -resume -profile test ${repodir}/main.nf --name ${name} \
+$NXFCMD --name ${name} \
     --outdir test_output/${name} \
     --input test_output/mzmldef \
     --genes \
@@ -55,9 +55,10 @@ set +eu
 echo Run in which no target PSMs are found - crashes
 name=labelfree_phos_notarget
 cat <(head -n1 ${testdir}/lf_mzmls.txt) <(grep setA ${testdir}/lf_mzmls.txt) | envsubst > test_output/mzmldef
-nextflow run -resume -profile test ${repodir}/main.nf --name ${name} \
+$NXFCMD --name ${name} \
     --outdir test_output/${name} \
     --input test_output/mzmldef \
+    --psmconflvl 0.001 \
     --genes \
     --tdb ${testdata}/lf.fa \
     --mods 'carbamidomethyl;oxidation' \
diff --git a/tests/tims.sh b/tests/tims.sh
index 49af58e..11f0530 100644
--- a/tests/tims.sh
+++ b/tests/tims.sh
@@ -7,7 +7,7 @@ echo TIMS TMT16, noms1quant
 name=tims_tmt16
 baseresults=test_output/${name}
 cat "${testdir}/tims_mzmls.txt" | envsubst > test_output/mzmldef
-nextflow run -resume -profile test ${repodir}/main.nf --name ${name} --outdir ${baseresults} \
+$NXFCMD --name ${name} --outdir ${baseresults} \
     --sampletable "${testdir}/tmt16_samples.txt" \
     --noms1quant \
     --input test_output/mzmldef \
diff --git a/tests/tmt16.sh b/tests/tmt16.sh
index 81a5ce2..c17c4e8 100644
--- a/tests/tmt16.sh
+++ b/tests/tmt16.sh
@@ -10,12 +10,14 @@ name=tmt16denomdeq
 resultsdir=test_output/${name}
 mkdir -p $resultsdir
 cat "${testdir}/tmt16_mzmls.txt" | envsubst > ${resultsdir}/mzmldef
-nextflow run -resume -profile test ${repodir}/main.nf --name ${name} --outdir ${resultsdir} \
+$NXFCMD --name ${name} --outdir ${resultsdir} \
     --input ${resultsdir}/mzmldef \
     --sampletable "${testdir}/tmt16_samples.txt" \
     --hardklor --isobaric '0set-A:tmtpro:126:131N' \
     --tdb "${testdata}/tmt16_fa.fa" \
     --mods 'carbamidomethyl;oxidation;43.005814,*,opt,N-term,Unknown' \
+    --maxmiscleav 2 \
+    --pepconflvl 0.05 \
     --deqms --keepnapsmsquant --genes
 
 
@@ -25,7 +27,7 @@ baseresults=${resultsdir}
 resultsdir=test_output/${name}
 mkdir -p $resultsdir
 cat <(head -n1 "${testdir}/tmt16_mzmls.txt") <(grep fr07 ${testdir}/tmt16_mzmls.txt) | envsubst > ${resultsdir}/mzmldef
-nextflow run -resume -profile test ${repodir}/main.nf --name ${name} --outdir test_output/${name} \
+$NXFCMD --name ${name} --outdir test_output/${name} \
     --input ${resultsdir}/mzmldef \
     --sampletable "${testdir}/tmt16_samples.txt" \
     --hardklor --isobaric '0set-A:tmtpro:126:131N' \
@@ -48,7 +50,7 @@ resultsdir=test_output/${name}
 mkdir -p $resultsdir
 cat "${testdir}/tmt16_mzmls.txt" | envsubst > ${resultsdir}/oldmzmls
 cat "${testdir}/tmt16_setB_mzmls.txt" | envsubst > ${resultsdir}/mzmldef
-nextflow run -resume -profile test ${repodir}/main.nf --name ${name} --outdir test_output/${name} \
+$NXFCMD --name ${name} --outdir test_output/${name} \
     --input ${resultsdir}/mzmldef \
     --sampletable "${testdir}/tmt16_setAB_samples.txt" \
     --isobaric '0set-A:tmtpro:126:131N set-B:tmt16plex:126:131N' \
@@ -69,7 +71,7 @@ baseresults=${resultsdir}
 resultsdir=test_output/${name}
 mkdir -p $resultsdir
 cat <(head -n1 "${testdir}/tmt16_mzmls.txt") <(grep fr07 ${testdir}/tmt16_mzmls.txt) <(grep fr08 ${testdir}/tmt16_mzmls.txt | cut -f1-3) | envsubst > ${resultsdir}/mzmldef
-nextflow run -resume -profile test ${repodir}/main.nf --name ${name} --outdir test_output/${name} \
+$NXFCMD --name ${name} --outdir test_output/${name} \
     --input ${resultsdir}/mzmldef \
     --sampletable "${testdir}/tmt16_samples.txt" \
     --hardklor --isobaric '0set-A:tmtpro:126:131N' \
@@ -85,7 +87,7 @@ name=tmt16_singlefile
 resultsdir=test_output/${name}
 mkdir -p $resultsdir
 cat <(head -n1 ${testdir}/tmt16_mzmls.txt) <(grep fr08 ${testdir}/tmt16_mzmls.txt) | envsubst > ${resultsdir}/mzmldef
-nextflow run -resume -profile test ${repodir}/main.nf --name ${name} --outdir test_output/${name} \
+$NXFCMD --name ${name} --outdir test_output/${name} \
     --input ${resultsdir}/mzmldef \
     --sampletable "${testdir}/tmt16_samples.txt" \
     --hardklor --isobaric '0set-A:tmtpro:126:131N' \
diff --git a/tests/tmt16_18.sh b/tests/tmt16_18.sh
index 4c3686e..be4d690 100644
--- a/tests/tmt16_18.sh
+++ b/tests/tmt16_18.sh
@@ -6,7 +6,7 @@ echo TMT 16/18 mix, multi DB run
 name=tmt16_18mix
 baseresults=test_output/${name}
 cat "${testdir}/tmt16_mzmls.txt" <(sed 's/0set-A/20set-A/' "${testdir}/tmt18_mzmls.txt" | tail -n+2) | envsubst > test_output/mzmldef
-nextflow run -resume -profile test ${repodir}/main.nf --name ${name} --outdir ${baseresults} \
+$NXFCMD --name ${name} --outdir ${baseresults} \
     --input test_output/mzmldef \
     --sampletable "${testdir}/tmt18_setAB_samples.txt" \
     --isobaric '0set-A:tmt16plex:126:131N 20set-A:tmt18plex:131' \
diff --git a/tests/tmt16_fast.sh b/tests/tmt16_fast.sh
index 1aa598d..c528a6d 100644
--- a/tests/tmt16_fast.sh
+++ b/tests/tmt16_fast.sh
@@ -1,16 +1,17 @@
 #!/usr/bin/env bash
 
+# Single test to not spend too much time, tests a lot of things but not all.
+
 set -eu
   
 echo TMT16 test with hiRIEF, isobaric and Phospho
-# no hirief, so it can be added in the addsetB test even if not used here
 # Test TMT16, DEqMS w denominator, keepnapsmsquant, implicit normalizing (deqms forces normalize)
 # Warning: no decoys for any set (aka only setA)
 name=tmt16denomdeq
 resultsdir=test_output/${name}
 mkdir -p $resultsdir
 cat "${testdir}/tmt16_mzmls.txt" | envsubst > ${resultsdir}/mzmldef
-nextflow run -resume -profile test ${repodir}/main.nf --name ${name} --outdir ${resultsdir} \
+$NXFCMD --name ${name} --outdir ${resultsdir} \
     --input ${resultsdir}/mzmldef \
     --sampletable "${testdir}/tmt16_samples.txt" \
     --hardklor --isobaric '0set-A:tmtpro:126:131N' \
@@ -19,5 +20,6 @@ nextflow run -resume -profile test ${repodir}/main.nf --name ${name} --outdir ${
     --locptms Phospho \
     --psmconflvl 0.2 --pepconflvl 0.2 \
     --deqms --keepnapsmsquant --genes \
+    --remove_channels '0set-A:127C:128N' \
     --hirief https://github.com/nf-core/test-datasets/raw/6defbf8a92a46b0ac48bb05f9ad96b62716b4a5d/testdata/formatted_known_peptides_ENSUniRefseq_TMT_predpi_20150825.txt
     # FIXME cannot run with carbamyl +43 -> -261 and Phospho, luciprep crash \
diff --git a/tests/tmt18.sh b/tests/tmt18.sh
index 7f90175..29e1453 100644
--- a/tests/tmt18.sh
+++ b/tests/tmt18.sh
@@ -11,7 +11,7 @@ echo TMT18 phos, no MS1 found somehow
 name=tmt18phos
 baseresults=test_output/${name}
 cat "${testdir}/tmt18_mzmls.txt" | envsubst > test_output/mzmldef
-nextflow run -resume -profile test ${repodir}/main.nf --name ${name} --outdir ${baseresults} \
+$NXFCMD --name ${name} --outdir ${baseresults} \
     --input test_output/mzmldef \
     --sampletable "${testdir}/tmt18_samples.txt" \
     --hardklor --isobaric '0set-A:tmt18plex:126:131N' \
@@ -30,7 +30,7 @@ mkdir -p test_output/${name}
 ln -fs "${testdata}/tmt18_fr06_1000.mzML" "${testdata}/linked_tmt18_fr06_1000.mzML"
 cat "${testdir}/tmt18_mzmls.txt" | envsubst > test_output/${name}/oldmzmls
 sed 's/tmt18_fr/linked_tmt18_fr/;s/0set-A/20set-A/' "${testdir}/tmt18_mzmls.txt" | envsubst > test_output/mzmldef
-nextflow run -resume -profile test ${repodir}/main.nf --name ${name} --outdir test_output/${name} \
+$NXFCMD --name ${name} --outdir test_output/${name} \
     --input test_output/mzmldef \
     --sampletable "${testdir}/tmt18_setAB_samples.txt" \
     --hardklor --isobaric '0set-A:tmt18plex:sweep 20set-A:tmt18plex:sweep' \
@@ -52,7 +52,7 @@ echo TMT18 rerun with different settings post PSMs, --noms1quant
 name=tmt18phos_rerun
 mkdir -p test_output/${name}
 cat "${testdir}/tmt18_mzmls.txt" | envsubst > test_output/${name}/oldmzmls
-nextflow run -resume -profile test ${repodir}/main.nf --name ${name} --outdir test_output/${name} \
+$NXFCMD --name ${name} --outdir test_output/${name} \
     --sampletable "${testdir}/tmt18_samples.txt" \
     --isobaric '0set-A:tmt18plex:131' \
     --tdb "${testdata}/tmt18_fa.fa" \
diff --git a/workflows/msgf_perco.nf b/workflows/msgf_perco.nf
index b5ab9fb..947d475 100644
--- a/workflows/msgf_perco.nf
+++ b/workflows/msgf_perco.nf
@@ -1,21 +1,4 @@
-include { listify; stripchars_infile; parse_isotype } from '../modules.nf'
-
-process createMods {
-  tag 'python'
-  container params.__containers[tag][workflow.containerEngine]
-
-  input:
-  tuple val(setname), val(isobtype), val(maxvarmods), path(msgfmods), val(mods)
-
-  output:
-  tuple val(setname), path('mods.txt')
-
-  script:
-  """
-  create_modfile.py $maxvarmods "${msgfmods}" "${mods}${isobtype ? ";${parse_isotype(isobtype)}" : ''}"
-  """
-}
-
+include { createMods; listify; stripchars_infile; parse_isotype } from '../modules.nf'
 
 process msgfPlus {
 
@@ -60,7 +43,7 @@ process percolator {
   container params.__containers[tag][workflow.containerEngine]
 
   input:
-  tuple val(setname), path(mzids), path(tsvs), val(enzyme)
+  tuple val(setname), path(mzids), val(enzyme)
 
   output:
   tuple val(setname), path('perco.xml')
@@ -130,7 +113,7 @@ workflow MSGFPERCO {
   mzml_in
   | map { [it.setname]}
   | groupTuple
-  | map { [it[0], setisobaric && setisobaric[it[0]] ? setisobaric[it[0]] : false, maxvarmods, file(msgfmodfile), mods] }
+  | map { [it[0], setisobaric && setisobaric[it[0]] ? setisobaric[it[0]] : false, maxvarmods, 'msgf', file(msgfmodfile), mods] }
   | createMods
 
   mzml_in
@@ -145,7 +128,7 @@ workflow MSGFPERCO {
   .set { mzid_tsv_2psm }
 
   mzid_tsv_2psm
-  | map {it + [enzyme] }
+  | map {[it[0], it[1], enzyme] }
   | percolator
   | join(mzid_tsv_2psm)
   | map { it + [psmconflvl, pepconflvl, output_unfilt_psms] }
diff --git a/workflows/ptms.nf b/workflows/ptms.nf
index 393a47f..b30312b 100644
--- a/workflows/ptms.nf
+++ b/workflows/ptms.nf
@@ -7,7 +7,7 @@ process luciphorPrep {
   container params.__containers[tag][workflow.containerEngine]
 
   input:
-  tuple val(setname), path(allpsms), val(locptms), val(stab_ptms), val(all_non_ptm_mods), path(msgfmodfile)
+  tuple val(setname), path(allpsms), val(locptms), val(stab_ptms), val(all_non_ptm_mods), path(msgfmodfile), val(search_engine)
 
   output:
   tuple val(setname), path('luciphor_config.txt'), path('lucipsms')
@@ -20,7 +20,7 @@ process luciphorPrep {
   msstitch split -i "${allpsms}" --splitcol TD
   luciphor_prep.py --psmfile target.tsv --template "$baseDir/assets/luciphor2_input_template.txt" \
       --labileptms "${lab_ptms}" --mods ${all_non_ptm_mods} ${stab_ptms} \
-      --modfile "${msgfmodfile}" --lucipsms lucipsms
+      --modfile "${msgfmodfile}" --lucipsms lucipsms --searchengine $search_engine
   """
 }
 
@@ -52,7 +52,8 @@ process luciphorPTMLocalizationScoring {
   export OUTFILE=luciphor.out
   cat "${template}" | envsubst > lucinput.txt
 
-  if luciphor2 -Xmx${task.memory.toGiga()}G lucinput.txt 2>&1 | grep 'not have enough PSMs'
+  luciphor2 -Xmx${task.memory.toGiga()}G lucinput.txt 2> luci_stderr
+  if grep 'not have enough PSMs' luci_stderr
   then
     echo 'Not enough PSMs for luciphor FLR calculation in set ${setname}' > warnings
   fi
@@ -68,7 +69,7 @@ process luciphorParse {
   // is no luciphor data due to errors, it will put NA for luciphor columns
 
   input:
-  tuple val(setname), path(luciphor_out), path('all_scores.debug'), path('psms'), val(ptm_minscore_high), path(msgfmodfile), val(lab_ptms), val(stab_ptms), val(all_non_ptm_mods), path(tdb)
+  tuple val(setname), path(luciphor_out), path('all_scores.debug'), path('psms'), val(ptm_minscore_high), path(msgfmodfile), val(lab_ptms), val(stab_ptms), val(all_non_ptm_mods), val(search_engine), path(tdb)
 
   output:
   tuple val(setname), path('labileptms.txt')
@@ -80,7 +81,7 @@ process luciphorParse {
      --luci_in ${luciphor_file} --luci_scores all_scores.debug --psms psms \
      --modfile "${msgfmodfile}" --labileptms ${lab_ptms.join(' ')} \
      ${stab_ptms ? "--stabileptms ${stab_ptms.join(' ')}": ''} --mods ${all_non_ptm_mods} \
-     --fasta "${tdb}"
+     --fasta "${tdb}" --searchengine $search_engine
   """
 }
 // FIXME msgfmods is false? oxidation so probably never.
@@ -92,7 +93,7 @@ process stabilePTMPrep {
   container params.__containers[tag][workflow.containerEngine]
 
   input:
-  tuple val(setname), val(ptms), path('psms'), path(tdb), val(non_ptm_mods), val(lab_ptms), path(msgfmods)
+  tuple val(setname), val(ptms), path('psms'), path(tdb), val(non_ptm_mods), val(lab_ptms), path(msgfmods), val(search_engine)
   
   output:
   tuple val(setname), path('stabileptms.txt')
@@ -102,7 +103,7 @@ process stabilePTMPrep {
   lab_ptms = listify(lab_ptms).join(' ')
   """
   nonlabile_ptm_columns.py --psms psms -o stabileptms.txt --modfile "${msgfmods}" --fasta "${tdb}" \
-      --stabileptms $stab_ptms \
+      --searchengine $search_engine --stabileptms $stab_ptms \
       ${lab_ptms ? "--labileptms $lab_ptms" : ""} \
       ${non_ptm_mods ? "--mods ${non_ptm_mods}" : ''}
   """
@@ -130,7 +131,7 @@ process createPTMTable {
   # PSM peptide sequences include the PTM site
   cat speclup.sql > ptmlup.sql
   ${has_newptms ? "msstitch concat -i ptms* ${oldptms ?: ''} -o" : "mv ${oldptms}"} concatptmpsms
-  msstitch psmtable -i concatptmpsms --dbfile ptmlup.sql -o ${ptmtable} --spectracol 1
+  msstitch psmtable -i concatptmpsms --dbfile ptmlup.sql -o ${ptmtable}
   msstitch split -i ${ptmtable} --splitcol bioset
   ${setnames.collect() { "test -f '${it}.tsv' || echo 'No PTMs found for set ${it}' >> warnings" }.join(' && ') }
   # No PTMs at all overwrites the per-set messages
@@ -281,6 +282,7 @@ workflow PTMANALYSIS {
   setisobaric
   mzmls
   psms
+  search_engine
   tdb
   unfiltered_psms
   fparams
@@ -304,7 +306,7 @@ workflow PTMANALYSIS {
   nofile = "${baseDir}/assets/NO__FILE"
   
   unfiltered_psms
-  | map { it + [locptms, stab_ptms, get_non_ptms(it[0], setisobaric, othermods), file(msgfmodfile)] }
+  | map { it + [locptms, stab_ptms, get_non_ptms(it[0], setisobaric, othermods), file(msgfmodfile), search_engine] }
   | luciphorPrep
   | join(mzmls)
   | map { it + [fparams.find { x -> x.setname == it[0] }.activation, maxpeplen, maxcharge, minpsms_luci] }
@@ -318,7 +320,7 @@ workflow PTMANALYSIS {
   | join(luciphorPTMLocalizationScoring.out.warnings, remainder: true)
   | map { [it[0], it[1] ?: nofile, it[1] ? it[2] : nofile] }
   | join(psms)
-  | map { it + [ptm_minscore_high, file(msgfmodfile), locptms, stab_ptms, get_non_ptms(it[0], setisobaric, othermods), ]}
+  | map { it + [ptm_minscore_high, file(msgfmodfile), locptms, stab_ptms, get_non_ptms(it[0], setisobaric, othermods), search_engine]}
   | combine(tdb)
   | luciphorParse
 
@@ -326,7 +328,7 @@ workflow PTMANALYSIS {
   | combine(psms)
   | map { [it[1], it[0], it[2]] }
   | combine(tdb)
-  | map { it + [get_non_ptms(it[0], setisobaric, othermods), locptms, file(msgfmodfile)] }
+  | map { it + [get_non_ptms(it[0], setisobaric, othermods), locptms, file(msgfmodfile), search_engine] }
   | stabilePTMPrep
   | concat(luciphorParse.out)
   | toList
diff --git a/workflows/reporting.nf b/workflows/reporting.nf
index 324bc77..08a84ed 100644
--- a/workflows/reporting.nf
+++ b/workflows/reporting.nf
@@ -76,16 +76,16 @@ process PSMQC {
   container params.__containers[tag][workflow.containerEngine]
 
   input:
-  tuple path('psms'), path('filescans'), path('platescans'), val(mzmls), val(fractionation), val(has_newmzmls), val(has_oldmzmls)
+  tuple path('psms'), path('filescans'), path('platescans'), val(mzmls), val(fractionation), val(has_newmzmls), val(has_oldmzmls), val(search_engine)
 
   output:
-  tuple path('platescans'), path('amount_psms_files'), path("psmplothtml"), path('psmtable__summary.txt')
+  tuple path('platescans'), path('amount_psms_files'), path("psmplothtml"), path('psmtable__summary.txt'), path('psms_ids.txt'), path('miscleav.txt')
 
   script:
   """
   ${mzmls.collect { "echo \$'${it.join('\t')}' >> mzmls"  }.join('\n')}
   sed '1s/\\#SpecFile/SpectraFile/' < psms > psms_clean
-  qc_psms.R ${fractionation ? 'TRUE' : 'FALSE'}
+  qc_psms.R ${fractionation ? 'TRUE' : 'FALSE'} $search_engine
   mkdir psmplothtml
   mv *.html psmplothtml/
   """
@@ -117,7 +117,7 @@ process featQC {
 
   mkdir ${htmldir}
   mv *.html ${htmldir}/
-  if [ -e *.png ] ; then mv *.png ${htmldir}/; fi
+  if ls *.png > /dev/null 2>&1 ; then mv *.png ${htmldir}/; fi
   ${parse_normfactors ? "mv allnormfacs ${htmldir}" : ''}
   """
 }
@@ -129,7 +129,7 @@ process PTMQC {
   container params.__containers[tag][workflow.containerEngine]
 
   input:
-  tuple path(ptmpsms), path(peptable)
+  tuple path(ptmpsms), path(peptable), val(search_engine)
 
   output:
   tuple path('ptmplothtml'), path('*.txt')
@@ -137,7 +137,7 @@ process PTMQC {
   script:
   """
   sed '1s/\\#SpecFile/SpectraFile/' < "${ptmpsms}" > psms_clean
-  qc_ptms.R psms_clean "${peptable}"
+  qc_ptms.R psms_clean "${peptable}" $search_engine
   mkdir ptmplothtml
   mv *.html ptmplothtml/
   """
@@ -151,7 +151,7 @@ process summaryReport {
 
 
   input:
-  tuple path(platescans), path(plotlibs), path('psmplots'), path(psm_summary), path(featplots), path(feat_summaries), path(feat_overlaps), path('ptmplots'), path(ptmfiles), path('warnings*')
+  tuple path(platescans), path(plotlibs), path('psmplots'), path(psm_summary), path('psmids'), path('miscleav'), path(featplots), path(feat_summaries), path(feat_overlaps), path('ptmplots'), path(ptmfiles), path('warnings*')
   
   output:
   tuple path('report_groovy_template.html'), path('libs.js')
@@ -171,6 +171,7 @@ workflow REPORTING {
   take:
   mzmls
   speclookup
+  search_engine
   oldmzmls
   fractionation
   plot_psms
@@ -220,7 +221,7 @@ workflow REPORTING {
   plot_psms
   | combine(countMS2sPerPlate.out.counted)
   | combine(all_mzmls_for_psmqc)
-  | map { it + [fractionation, 1,1] }
+  | map { it + [fractionation, 1, 1, search_engine] }
   | PSMQC
 
   pepprotgenes
@@ -238,6 +239,7 @@ workflow REPORTING {
   // Only pick one PTM table for reporting nr of sites etc
   ptms
   | filter { it[1].name.contains('_not_adjusted') }
+  | map { [it, search_engine].flatten() }
   | PTMQC
   | ifEmpty([nofile, nofile])
   | set { ptmqc }
diff --git a/workflows/sage_perco.nf b/workflows/sage_perco.nf
new file mode 100644
index 0000000..3e8d575
--- /dev/null
+++ b/workflows/sage_perco.nf
@@ -0,0 +1,175 @@
+include { createMods; listify; stripchars_infile; parse_isotype } from '../modules.nf'
+
+
+process sagePrepare {
+  tag 'jq'
+  container params.__containers[tag][workflow.containerEngine]
+ 
+  input:
+  tuple val(setname), val(id), val(minlen), val(maxlen), val(mincharge), val(maxcharge), val(maxmiscleav), val(maxvarmods), val(fparams), path('sage.json'), path('mods.json')
+
+  output:
+  tuple val(id), path('config.json')
+  
+  script:
+  prectol_num = fparams.prectol.replaceAll('[^0-9\\.]', '')
+  prectol_unit = fparams.prectol.replaceAll('[0-9\\.]', '')
+  sage_maxmiscleav = maxmiscleav == -1 ? 20 : maxmiscleav
+  """
+  cat sage.json | jq --argjson MODS "\$(cat mods.json)" '.database.enzyme.missed_cleavages=${sage_maxmiscleav} | \
+    .database+=\$MODS | \
+    .database.enzyme.min_len=${minlen} | \
+    .database.enzyme.max_len=${maxlen} | .database.max_variable_mods=${maxvarmods} | \
+    .precursor_tol={$prectol_unit: [-${prectol_num}, ${prectol_num}]} | \
+    .precursor_charge=[${mincharge}, ${maxcharge}]' > config.json
+  """
+    //.fragment_tol.ppm=[-${fparams.fragtol}, ${fparams.fragtol}] | \
+}
+
+
+process sage {
+  tag 'sage'
+  container params.__containers[tag][workflow.containerEngine]
+
+  input:
+  tuple val(setname), path('config.json'), path(specfile), val(fparams), val(instrumenttype), val(fractionation), path(db)
+
+  output:
+  tuple val(setname), path("${specfile.baseName}.pin"), emit: perco
+  tuple val(setname), path("${specfile.baseName}.tsv"), emit: tsv
+
+  script:
+  // Bruker from msconvert has merged= frame= scanStart= scanEnd= as scannr
+  // Bruker from tdf2mzml has index=.. as scannr
+  remove_scan_index_str = instrumenttype == 'bruker'
+  (is_stripped, parsed_infile) = stripchars_infile(specfile)
+  (_, parsed_basename) = stripchars_infile(file(specfile.baseName))
+  """
+  ${is_stripped ? "ln -s ${specfile} '${parsed_infile}'" : ''}
+  export RAYON_NUM_THREADS=${task.cpus}
+  export SAGE_LOG=trace
+  sage --disable-telemetry-i-dont-want-to-improve-sage --write-pin -f $db config.json $parsed_infile
+  ${remove_scan_index_str ? "sed -i 's/index=//' results.sage.pin" : ''}
+  ${remove_scan_index_str ? "sed -E -i 's/merged=([0-9]+) [SEa-z0-9\\ =]+/\\1/' results.sage.pin" : ''}
+
+  # Add set/strip/fraction
+  awk -F \$'\\t' '{OFS=FS ; print \$0, \
+    "Biological set" \
+    ${fractionation ? ', "Strip", "Fraction"' : ''}}' \
+        <( head -n1 results.sage.tsv) > "${specfile.baseName}.tsv"
+  awk -F \$'\\t' '{OFS=FS ; print "${parsed_basename}" \$0, \
+    "$setname" ${fractionation ? ", \"$fparams.plate\", \"$fparams.fraction\"" : ''}}' \
+        <( tail -n+2 results.sage.tsv) >> "${specfile.baseName}.tsv"
+
+  ${remove_scan_index_str ? "sed -i 's/index=//' results.sage.tsv" : ''}
+  ${remove_scan_index_str ? "sed -E -i 's/merged=([0-9]+) [SEa-z0-9\\ =]+/\\1/' results.sage.tsv" : ''}
+
+  head -n1 results.sage.pin > "${specfile.baseName}.pin"
+  awk -F \$'\\t' '{OFS=FS ; print "${parsed_basename}" \$0}' <( tail -n+2 results.sage.pin) >> "${specfile.baseName}.pin"
+  """
+}
+
+
+process percolator {
+  tag 'percolator'
+  container params.__containers[tag][workflow.containerEngine]
+
+  input:
+  tuple val(setname), path(pins)
+
+  output:
+  tuple val(setname), path('perco.xml')
+
+  script:
+  """
+  head -n1 ${listify(pins)[0]} > percoin.tsv
+  cat ${listify(pins).collect { "<(tail -n+2 $it)" }.join(' ')} >> percoin.tsv
+  percolator -j percoin.tsv -X perco.xml -N 500000 --decoy-xml-output -Y --num-threads ${task.cpus}
+  """
+}
+
+
+process percolatorToPsms {
+
+  tag 'msstitch'
+  container params.__containers[tag][workflow.containerEngine]
+
+  input:
+  tuple val(setname), path(perco), path(tsvs), val(psmconf), val(pepconf), val(output_unfiltered)
+
+  output:
+  tuple val('target'), path("${setname}_target.tsv"), emit: tmzidtsv_perco optional true
+  tuple val('decoy'), path("${setname}_decoy.tsv"), emit: dmzidtsv_perco optional true
+  tuple val(setname), path('allpsms'), emit: unfiltered_psms optional true
+  path('warnings'), emit: percowarnings optional true
+
+  script:
+    //--mzids ${listify(mzids).collect(){"${it}"}.join(' ')} \
+  """
+  mkdir outtables
+  msstitch perco2psm --perco perco.xml -d outtables \
+    -i ${listify(tsvs).collect(){"${it}"}.join(' ')} \
+    ${!output_unfiltered ? "--filtpsm ${psmconf} --filtpep ${pepconf}" : ''}
+  msstitch concat -i outtables/* -o allpsms
+  ${output_unfiltered ?
+    "msstitch conffilt -i allpsms -o filtpsm --confcolpattern 'PSM q-value' --confidence-lvl ${psmconf} --confidence-better lower && \
+    msstitch conffilt -i filtpsm -o psms --confcolpattern 'peptide q-value' --confidence-lvl ${pepconf} --confidence-better lower" : 'cp allpsms psms'}
+  msstitch split -i psms --splitcol TD
+  ${['target', 'decoy'].collect() { 
+    "test -f '${it}.tsv' && mv '${it}.tsv' '${setname}_${it}.tsv' || echo 'No ${it} PSMs found for set ${setname} at PSM FDR ${psmconf} and peptide FDR ${pepconf} ${it == 'decoy' ? '(not possible to output protein/gene FDR)' : ''}' >> warnings" }.join(' && ') }
+  """
+}
+
+
+workflow SAGEPERCO {
+  take:
+  mzml_in
+  concatdb
+  setisobaric
+  fractionation
+  fparam_map
+  maxvarmods
+  msgfmodfile
+  mods
+  psmconflvl
+  pepconflvl
+  output_unfilt_psms
+  maxmiscleav
+  enzyme
+  minpeplen
+  maxpeplen
+  mincharge
+  maxcharge
+
+  main: 
+
+  mzml_in
+  | map { [it.setname]}
+  | groupTuple
+  | map { [it[0], setisobaric && setisobaric[it[0]] ? setisobaric[it[0]] : false, maxvarmods, 'sage', file(msgfmodfile), mods] }
+  | createMods
+
+  mzml_in
+  | map { [it.setname, it.id, minpeplen, maxpeplen, mincharge, maxcharge, maxmiscleav, maxvarmods, fparam_map[it.id], file("$baseDir/assets/sage.json")] }
+  | combine(createMods.out, by: 0)
+  | sagePrepare
+  | join(mzml_in.map { [it.id, it.setname, it.mzmlfile, fparam_map[it.id], it.instrument == 'timstof' ? 'bruker' : 'thermo', fractionation] })
+  | map { [it[2], it[1], it[3..-1]].flatten() } // Set, config, rest
+  | combine(concatdb)
+  | sage
+
+  sage.out.perco
+  | groupTuple
+  | percolator
+  | join(sage.out.tsv | groupTuple)
+  | map { it + [psmconflvl, pepconflvl, output_unfilt_psms] }
+  | percolatorToPsms
+
+//  tuple val(id), path('config.json'), path(specfile), val(instrumenttype), path(mods), path(db)
+
+  emit:
+  t_tsv = percolatorToPsms.out.tmzidtsv_perco
+  d_tsv = percolatorToPsms.out.dmzidtsv_perco
+  unfiltered = percolatorToPsms.out.unfiltered_psms
+  warnings = percolatorToPsms.out.percowarnings
+}